Potential non-physical solution on Caradonna-Tung test case

[jphkun]

Describe the bug
Reference case and expectations
After many, many tries, the Caradonna-Tung test case seems to "converge" towards a non-physical solution. The flow creates a vortex around the blade instead of generating a downward jet as one can expect (https://www.nas.nasa.gov/SC08/v22_2.html). The case simulates a propeller with a compressible Euler solver.

Main issue
The outlet, pressure outlet, or farfield boundary conditions seems to destabilize after some iterations due to the high speed gradient on the outlet face. To get into more details, there is a vortex forming ring forming near of the outlet creating a inflow at the center. The BC seems to build on this inflow up to a point where the flow destabilizes the solution and diverges or pushes the flow to make a ring around the propeller.

Case setup and context
To put some context the current case tries to simulate the flow-field generated by a blade/propeller in a hover state in order to compute the thrust.

The blade is surrounded by a cylinder. The markers for each surfaces are:
-"inlet" for the top of the cylinder
-"farfield" for the sides of the cylinder
-"outlet" for the bottom of the cylinder.

The blade has been divided into 4 surfaces which have the following self explanatory marker names:

• "propellerUpperSurface"
• "propellerLowerSurface"
• "propellerRightTip"
• "propellerLeftTip"

Many combinations of the possible boundary conditions setup were tested for the "inlet", "farfield" and "outlet" markers. Furthermore the Inlet BC were tried with many different vertical speeds (vx=0/0.001/0.01/0.1/1/10). Simulations with the inlet, farfield and outlet as TOTAL_INLET, FARFIELD, OUTLET and EULER and some (not all) combination between these BC were tried.

The mesh was generated with Salome-Meca. It is a tetrahedral mesh. I have a script that permits to change the domain size with ease. So from the mesh perspective it is fairly easy to test different mesh sizes and domain sizes.
Generally the outer domain radius was between 1.5propRadius to 10propRadus. Domain height between 3propRadius to 15propRadius. The prop was at 2/3 of the total height (1/2 and 3/4 were also tested).

Some simulations approximated from a far experimental results (thrust +/-10% error, torque +/-15%). Some completely diverged some converge towards something way out of the 10% error range.

To Reproduce

You can find the configuration file below.
I cannot provide the mesh at the moment.
If needed I will change the geometry for "reproduction".

`%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% SU2 configuration file %
% Case description: Caradonna & Tung rotor in hover (still air) %
% Author: Thomas D. Economon %
% Institution: Stanford University %
% Date: 2020.05.24 %
% File Version 7.3.1 "Blackbird" %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------%
%
% Physical governing equations (EULER, NAVIER_STOKES,
% WAVE_EQUATION, HEAT_EQUATION, FEM_ELASTICITY,
% POISSON_EQUATION)
SOLVER= EULER
%
% Mathematical problem (DIRECT, CONTINUOUS_ADJOINT)
MATH_PROBLEM= DIRECT
%
% Restart solution (NO, YES)
RESTART_SOL= NO

% ----------- COMPRESSIBLE AND INCOMPRESSIBLE FREE-STREAM DEFINITION ----------%
%
% Mach number (non-dimensional, based on the free-stream values)
MACH_NUMBER= 0.0
%
% Angle of attack (degrees)
AOA= 0.0
%
% Side-slip angle (degrees)
SIDESLIP_ANGLE= 0.0
%
% Free-stream pressure (101325.0 N/m^2 by default, only Euler flows)
FREESTREAM_PRESSURE= 101325.0
%
% Free-stream temperature (288.15 K by default)
FREESTREAM_TEMPERATURE= 288.15

% ---------------------- REFERENCE VALUE DEFINITION ---------------------------%
%
% Reference origin for moment computation
REF_ORIGIN_MOMENT_X = 0.25
REF_ORIGIN_MOMENT_Y = 0.00
REF_ORIGIN_MOMENT_Z = 0.00
%
% Reference length for pitching, rolling, and yawing non-dimensional moment
REF_LENGTH= 1.0
%
% Reference area for force coefficients (0 implies automatic calculation)
REF_AREA= 14.098904

% ----------------------- DYNAMIC MESH DEFINITION -----------------------------%
%
% Type of dynamic mesh (NONE, RIGID_MOTION, DEFORMING, ROTATING_FRAME,
% MOVING_WALL, AEROELASTIC, EXTERNAL)
GRID_MOVEMENT= ROTATING_FRAME
%
% Motion mach number (non-dimensional). Used for intitializing a viscous flow
% with the Reynolds number and for computing force coeffs. with dynamic meshes.
MACH_MOTION= 0.001
%
% Coordinates of the motion origin
MOTION_ORIGIN= 0.0 0.0 0.0
%
% Angular velocity vector (rad/s) about the motion origin
ROTATION_RATE = 335.103216 0.0 0.0
% 2000 RPM -> 209.43950999999998
% 2200 RPM -> 230.38346099999998
% 2400 RPM -> 251.32741199999998
% 2600 RPM -> 272.271363
% 2800 RPM -> 293.215314
% 3000 RPM -> 314.159265
% 3000 RPM -> 335.103216

% -------------------- BOUNDARY CONDITION DEFINITION --------------------------%
%
% Marker of the Euler boundary (NONE = no marker)
MARKER_EULER= ( innerCylinder, propellerUpperSurface, propellerLowerSurface, propellerRightTip, propellerLeftTip )
%
% Marker of the far field (NONE = no marker)
%
% Read inlet profile from a file (YES, NO) default: NO
%SPECIFIED_INLET_PROFILE= YES
%
% File specifying inlet profile
INLET_TYPE= TOTAL_CONDITIONS
MARKER_INLET = (inlet, 288.15, 101325, 0.0, 0.0, 0.0 )
MARKER_FAR= (farfield)
MARKER_OUTLET = ( outlet, 101325 )
%
% Marker of the surface which is going to be plotted or designed
MARKER_PLOTTING= ( propellerUpperSurface, propellerLowerSurface, propellerRightTip, propellerLeftTip )
%
% Marker of the surface where the functional (Cd, Cl, etc.) will be evaluated
MARKER_MONITORING= ( propellerUpperSurface, propellerLowerSurface, propellerRightTip, propellerLeftTip )

% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------%
%
% Numerical method for spatial gradients (GREEN_GAUSS, LEAST_SQUARES,
% WEIGHTED_LEAST_SQUARES)
NUM_METHOD_GRAD= GREEN_GAUSS
%
% Courant-Friedrichs-Lewy condition of the finest grid
CFL_NUMBER= 1
%
% Adaptive CFL number (NO, YES)
CFL_ADAPT= NO
%
% Parameters of the adaptive CFL number (factor down, factor up, CFL min value,
% CFL max value )
CFL_ADAPT_PARAM= ( 0.1, 2.0, 15, 1e6 )
%
% Runge-Kutta alpha coefficients
RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 )
%
% Number of total iterations
ITER= 99999

% ------------------------ LINEAR SOLVER DEFINITION ---------------------------%
%
% Linear solver for implicit formulations (BCGSTAB, FGMRES)
LINEAR_SOLVER= FGMRES
%
% Preconditioner of the Krylov linear solver (JACOBI, LINELET, LU_SGS)
LINEAR_SOLVER_PREC= ILU
%
% Minimum error of the linear solver for implicit formulations
LINEAR_SOLVER_ERROR= 1E-10
%
% Max number of iterations of the linear solver for the implicit formulation
LINEAR_SOLVER_ITER= 10

% -------------------------- MULTIGRID PARAMETERS -----------------------------%
%
% Multi-Grid Levels (0 = no multi-grid)
MGLEVEL= 3
%
% Multi-grid cycle (V_CYCLE, W_CYCLE, FULLMG_CYCLE)
MGCYCLE= W_CYCLE
%
% Multi-grid pre-smoothing level
MG_PRE_SMOOTH= ( 1, 2, 3, 3 )
%
% Multi-grid post-smoothing level
MG_POST_SMOOTH= ( 0, 0, 0, 0 )
%
% Jacobi implicit smoothing of the correction
MG_CORRECTION_SMOOTH= ( 0, 0, 0, 0 )
%
% Damping factor for the residual restriction
MG_DAMP_RESTRICTION= 0.9
%
% Damping factor for the correction prolongation
MG_DAMP_PROLONGATION= 0.9

% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------%
%
% Convective numerical method (JST, LAX-FRIEDRICH, CUSP, ROE, AUSM, HLLC,
% TURKEL_PREC, MSW)
CONV_NUM_METHOD_FLOW= JST
%
% Monotonic Upwind Scheme for Conservation Laws (TVD) in the flow equations.
% Required for 2nd order upwind schemes (NO, YES)
MUSCL_FLOW= YES
%
% Slope limiter (VENKATAKRISHNAN, MINMOD)
SLOPE_LIMITER_FLOW= VENKATAKRISHNAN
%
% 2nd and 4th order artificial dissipation coefficients
JST_SENSOR_COEFF= ( 0.5, 0.02 )
%
% Time discretization (RUNGE-KUTTA_EXPLICIT, EULER_IMPLICIT, EULER_EXPLICIT)
TIME_DISCRE_FLOW= EULER_IMPLICIT

% --------------------------- CONVERGENCE PARAMETERS --------------------------%
%
% Min value of the residual (log10 of the residual)
CONV_RESIDUAL_MINVAL= -12
%
% Start convergence criteria at iteration number
CONV_STARTITER= 10
%
% Number of elements to apply the criteria
CONV_CAUCHY_ELEMS= 100
%
% Epsilon to control the series convergence
CONV_CAUCHY_EPS= 1E-10
%

% ------------------------- INPUT/OUTPUT INFORMATION --------------------------%
%
% Mesh input file
MESH_FILENAME= NLU_Baseline.su2
%
% Mesh input file format (CGNS, NETCDF_ASCII)
MESH_FORMAT= SU2
%
% Mesh output file
MESH_OUT_FILENAME= mesh_caradonna_tung_periodic.su2
%
% Restart flow input file
SOLUTION_FILENAME= solution_flow.dat
%
% Restart adjoint input file
SOLUTION_ADJ_FILENAME= solution_adj.dat
%
TABULAR_FORMAT= CSV
%
% Output file convergence history (w/o extension)
CONV_FILENAME= history
%
% Output file restart flow
RESTART_FILENAME= restart_flow.dat
%
% Output file restart adjoint
RESTART_ADJ_FILENAME= restart_adj.dat
%
% Output file flow (w/o extension) variables
VOLUME_FILENAME= flow
%
% Output file adjoint (w/o extension) variables
VOLUME_ADJ_FILENAME= adjoint
%
% Output Objective function gradient (using continuous adjoint)
GRAD_OBJFUNC_FILENAME= of_grad.dat
%
% Output file surface flow coefficient (w/o extension)
SURFACE_FILENAME= surface_flow
%
% Output file surface adjoint coefficient (w/o extension)
SURFACE_ADJ_FILENAME= surface_adjoint
%
% Writing solution file frequency
OUTPUT_WRT_FREQ= 10
%

SCREEN_OUTPUT= (INNER_ITER, WALL_TIME, RMS_RES, FORCE_X, MOMENT_X)
HISTORY_OUTPUT=(INNER_ITER, WALL_TIME, RMS_RES, AERO_COEFF)`

Desktop (please complete the following information):

• OS: [e.g., Linux (Redhat 6.6)]
• C++ compiler and version: [e.g., g++ (GCC) 4.8.5]
• MPI implementation and version: [e.g., OpenMPI 3.0.0]
• SU2 Version: [e.g., v6.2.0]
  3 OS were used:
  Ubuntu 22.04 on AWS
  WSL with ubuntu 20.04 LST
  Windows

Latest MPICH package was used on the ubuntu versions
(basically the one that come with a 'sudo apt-get install mpich')
Latest mpiexec version is used on windows.

The code was not compiled, the binaries form the website were used (SU2 v7.3.1).

[pcarruscag]

Can you post some pictures?

[jphkun]

Hi Pedro, hi everyone,

Sorry for the closing/reopening. I am still new to Github and did a wrong manipulation. It is my fist issue ever...

Here are some images of the solution every 100 iteration. The python code used to copy the flow file was not working properly and missed some of the iterations. Do not worry if some iterations are missing.
images_100.zip

Here is the config file:
3200_config_Euler.txt

Here is the solution output. Mach number, surface, and chord were changed in or
history.csv

Here is the convergence plot

Do not hesitate to ask for more information. When I have some time I will change the blade to post the mesh in order to replicate the issue. Modified blade mesh should come by the end of the week.

[pcarruscag]

I would use a more recent version of su2.
And you shouldn't mix inlets outlets and far field, you can also try turning off the multi grid, it does not behave very well sometimes.

[clarkpede]

you shouldn't mix inlets outlets and far field,

Forgive my ignorance here-- can you explain why? Isn't that the setup used on the Turbulence Modeling Resource for a flat plate boundary layer?

[pcarruscag]

If one can simply use far field for everything why make it more complicated? Introducing the following risks in so doing:

• Some inconsistency in the values that cause problems (static and total specification).
• Some inconsistency on our end that causes problems even if all inputs are ok.
• I don't think our outlets deal with reversed flow very well.

[clarkpede]

Ah ok. Thanks for the explanation.

[jphkun]

Thank you guys for the explanations! I will re-try an all farfield case this time with the CFL at 1 and without the multi-grid.
I also witnessed the "inflow" problem with the pressure outlet BC. Good to know that it was not just me.

Regarding the SU2 version, I use the pre-compiled version on the website hence the kind of "old" SU2 version. I had some trouble compiling it form source and abandoned that route weeks ago.

I will keep you updated of my endeavors. Mesh is on the way

[jphkun]

Hi again everyone,

After a new day of simulation with the all farfield BCs it looks like the results are better. Sadly the farfield BC is not working as expected. Ideally the flow should be going out with a sort of Poiseuille profile in the center and not recycle back in as it is the case here. You can see the start of the vortex on the image below.

Basically I have the following question:

1. Is there a way to set up a boundary condition that allows the jet to escape the domain?
2. Keep in mind that the farfield BC was set with an infinite speed of Mach 0.0. Is this preventing something from happening?

To put some context I expect the flow to pass thought the bottom of the cylinder as in the figure below. Probably that this simulation has a different setup, the image might not reflect the whole simulated domain, computing real time, modelling turbulence, etc. But nevertheless the fact that the flow goes out of the domain is of primary importance for me, otherwise the simulation risks two things:

1. Measuring ground effect.
2. Having a non 0 inflow which would reduce the lift significantly.

source: https://www.nas.nasa.gov/SC08/v22_2.html

[pcarruscag]

If you are running with cfl of 1 it will take a loonnnng time to wash out that reflexion/vortex, especially because the farfield is too close.
Try to raise the cfl now that it seems more stable.

[jphkun]

Next time i will increase the CFL and see what happens.

When you say that the farfield is too close, do you mean:

1. Sides are too close
2. Top and/or bottom of the domain is too close
3. Both, sides and top/bottom too close from the blade

[pcarruscag]

In that picture everything is way too close, but clearly the bottom is more important, maybe put that at 50x diameter and the rest 25x if you want to save on mesh size.
Depending on whether you are interested in force/moment coefficients or pretty pictures of the wake, you may also coarsen the mesh aggressively away from the rotor, maybe from 5-10x (I'm no expert here, so I recommend testing the sensitivity of the results to this type of choices) perhaps there are good guidelines out there from workshops specific to this type of application.

[jphkun]

Yes, it all makes sense. Thanks a lot!