Interface between the moving and stationary zones generates wrong results as the moving zone rotates

a3adam commented 7 months ago

Summary

The fluid should freely pass a fluid interface between RANS zones. However, in some unsteady 2D simulations, I encountered the interface generating unexpected flow solutions which is observed through the increased eddy viscosity and unphysical velocity as the inner zone rotates with respect to the stationary outer zone.

I have tried using

a finer mesh with more cells on the zone interface,
the incompressible and compressible RANS solvers,
smaller time steps and more pseudo-time iterations,
different interface interpolation methods (NEAREST_NEIGHBOR and WEIGHTED_AVERAGE),
different convective numerical methods (upwind/central and with and without MUSCL_FLOW)

however the problem persisted.

Setup

Mesh

File: vawt_ami.su2

Configuration files

su2Config.cfg:

SOLVER= INC_RANS
KIND_TURB_MODEL= SST
MATH_PROBLEM= DIRECT
CONFIG_LIST= (zone_sta.cfg, zone_rot.cfg)

INNER_ITER= 1
OUTER_ITER= 20
TIME_ITER= 3160

TIME_DOMAIN= YES
TIME_MARCHING= DUAL_TIME_STEPPING-2ND_ORDER
TIME_STEP = 0.000593412
MAX_TIME= 50.0

INC_DENSITY_INIT= 1.225
INC_VELOCITY_INIT= ( 10.0, 0.0, 0.0 )
INC_TEMPERATURE_INIT= 288.15
INC_INLET_TYPE= VELOCITY_INLET
INC_OUTLET_TYPE= PRESSURE_OUTLET

REF_ORIGIN_MOMENT_X = 0.00
REF_ORIGIN_MOMENT_Y = 0.00
REF_ORIGIN_MOMENT_Z = 0.00
REF_LENGTH= 0.85
REF_AREA= 1.7

VISCOSITY_MODEL= CONSTANT_VISCOSITY
MU_CONSTANT= 1.8375E-5

MARKER_ZONE_INTERFACE= ( AMI1, AMI2 )
MARKER_FLUID_INTERFACE= ( AMI1, AMI2 )
KIND_INTERPOLATION= NEAREST_NEIGHBOR

NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES
CFL_NUMBER= 20.0
MAX_DELTA_TIME= 1E6

SLOPE_LIMITER_FLOW=  VENKATAKRISHNAN
MUSCL_TURB= NO
SLOPE_LIMITER_TURB= VENKATAKRISHNAN
VENKAT_LIMITER_COEFF= 0.05
REF_SHARP_EDGES = 3.0
SENS_REMOVE_SHARP = NO
LIMITER_ITER= 999999

LINEAR_SOLVER= FGMRES
LINEAR_SOLVER_PREC= ILU
LINEAR_SOLVER_ERROR= 1E-6
LINEAR_SOLVER_ITER= 10

CONV_NUM_METHOD_FLOW= JST
TIME_DISCRE_FLOW= EULER_IMPLICIT

CONV_NUM_METHOD_TURB= SCALAR_UPWIND
TIME_DISCRE_TURB= EULER_IMPLICIT

MULTIZONE= YES
MULTIZONE_MESH= YES

MESH_FILENAME= vawt_ami.su2
MESH_FORMAT= SU2

_zonesta.cfg:

GRID_MOVEMENT= NONE

MARKER_FAR= ( upper, lower )
MARKER_INLET= ( inlet, 288.15, 10.0, 1.0, 0.0, 0.0 )
MARKER_OUTLET= ( outlet, 0.0 )
MARKER_FLUID_INTERFACE = (AMI1)

_zonerot.cfg:

GRID_MOVEMENT= RIGID_MOTION
MACH_MOTION= 0.0294
MOTION_ORIGIN= 0.0 0.0 0.0
ROTATION_RATE = 0.0 0.0 29.41176471

MARKER_HEATFLUX= ( blade1, 0.0, blade2, 0.0, blade3, 0.0  )
MARKER_FLUID_INTERFACE = (AMI2)
MARKER_PLOTTING = ( blade1, blade2, blade3 )
MARKER_MONITORING = ( blade1, blade2, blade3 )

Simulation results

Velocity

Eddy viscosity

Desktop

OS: Linux (CentOS 7.9)
C++ compiler and version: g++ (GCC) 12.1
MPI implementation and version: OpenMPI 4.1.4
SU2 Version: v8.0.0 "Harrier"

joshkellyjak commented 7 months ago

Are you able to reproduce this behaviour with any other cases, say a rotating cylinder perhaps? To me it looks like the strange behaviour develops from the area in your mesh where the mesh density along the interface is visibily different (north and south of the image shown). This could be a discretization issue rather than a bug with the solver?

jblueh commented 7 months ago

I just found a heap buffer overflow in the interface code, mentioning it here in case it is related https://github.com/su2code/SU2/pull/2246#issuecomment-2034404720

a3adam commented 7 months ago

To me it looks like the strange behaviour develops from the area in your mesh where the mesh density along the interface is visibily different (north and south of the image shown).

I have encountered the same issue with a finer mesh where the mesh density along the interface is more or less uniform:

Are you able to reproduce this behaviour with any other cases, say a rotating cylinder perhaps?

I have tried a rotating 2D cylinder case with a similar setting. The only differences are that the time step is reduced to a tenth to prevent any numerical oscillations, and the mesh is generated by using Gmsh (for the VAWT case, it was Pointwise).

The simulation diverged at the end. The discontinuity around the interface is quite apparent, especially on the upstream part.

Mesh:

Axial velocity:

Eddy viscosity:

a3adam commented 6 months ago

I have tried two additional simulations with the same cylinder mesh:

1. Stationary cylinder

The settings are almost identical, except:

su2Config.cfg:

...
TIME_STEP = 0.00593412
...

zone_rot.cfg:

GRID_MOVEMENT= RIGID_MOTION
MOTION_ORIGIN= 0.0 0.0 0.0
ROTATION_RATE = 0.0 0.0 0.0
...

In other words, this is a transient stationary cylinder simulation with two zones. As the axial velocity and eddy viscosity evolutions show below, the simulation runs smoothly and no discontinuities or unphysical flow result is observed around the zone interface.

2. Slowly rotating cylinder

I also tried a very slowly rotating cylinder to see if the error will occur even with a slow rotation. The rotational speed of the inner zone is around 0.03 rad/s, which corresponds to a surface angular speed (0.025 m/s) equal to the 0.25% of the freestream velocity (Uinf=10 m/s).

zone_rot.cfg:

GRID_MOVEMENT= RIGID_MOTION
MOTION_ORIGIN= 0.0 0.0 0.0
ROTATION_RATE = 0.0 0.0 0.02941176471
...

The time step is set as the one used for the stationary cylinder. Each time step corresponds to a 0.01 deg rotation of the inner zone. The below animation shows the movement of the cells on the interface:

The below animations show the axial velocity and eddy viscosity evolution. The results should not deviate much from the stationary cylinder results as the angular speed of the cylinder is quite low compared to the freestream. However, it is not the case for this simulation. Also, the discontinuity on the downstream part of the interface is quite apparent even at the earlier time steps. As the simulation progresses, unphysical flow results (non-zero eddy viscosity at the top part of the interface) appear, and eventually, the simulation diverges.

jblueh commented 6 months ago

I just found a heap buffer overflow in the interface code, mentioning it here in case it is related #2246 (comment)

I ran your test case with the fixes from #2246 and the address sanitizer, it still shows the behaviour that you describe and the address sanitizer did not detect anything else. So it's probably not related to the aforementioned overflow.

bigfooted commented 6 months ago

Does it also occur without turbulence, so for INC_NAVIER_STOKES and Re=50-180?

bigfooted commented 6 months ago

There are 2 cases in the testcases repository, both are using the Euler solver. Maybe also check if your case is running correctly with the settings from those testcases.

a3adam commented 6 months ago

There are 2 cases in the testcases repository, both are using the Euler solver. Maybe also check if your case is running correctly with the settings from those testcases.

Earlier, I was defining the mesh movement for the static zone as:

GRID_MOVEMENT= NONE

since it was stationary. However, after checking the test cases, I found that the stator zones should also be defined as moving as

GRID_MOVEMENT= RIGID_MOTION
MOTION_ORIGIN= 0.0 0.0 0.0
ROTATION_RATE = 0.0 0.0 0.0

even though the mesh stays stationary. In that case, the unphysical flow results I presented earlier vanished. Earlier, I tried the cylinder test case with the inner zone having

GRID_MOVEMENT= RIGID_MOTION
MOTION_ORIGIN= 0.0 0.0 0.0
ROTATION_RATE = 0.0 0.0 0.0

and the outer zone having GRID_MOVEMENT= NONE

where the simulation gave physically viable results. I guess the zone interface may not accurately work if the GRID_MOVEMENT options of the neighbouring zones do not match while one of the zones has a non-zero movement. However, I have not simulated any other test cases.

Also, the sliding_interface/rotating_cylinders test case diverges with the given .cfg files. At around the 775th time step, I received the following error message:

SU2 has diverged (Residual > 10^20 detected).

However, it works just fine if the simulations are done in subsonic flow (e.g. MACH_NUMBER= 0.1) while keeping the rest of the configuration the same.

su2code / SU2