High Lift Prediction Workshop 5

[rois1995]

Hi everyone!

Myself, @tbellosta and @bigfooted are involved in the HLPW-5 (https://hiliftpw.larc.nasa.gov/). Personally, I am within both the TFG group for Fixed and Adaptive Grid RANS, although, at the moment, I plan to be just an observer for the latter.

Currently, we are running the computations for the first test case, which is in conjunction with the High Fidelity CFD Verification Workshop 2024. The geometry is quite simple as it is a half-model of the CRM without flaps, slats or nacelle. The test case consists of a grid convergence analysis performed with the SA-neg-QCR2000-R turbulence model, implemented as in https://doi.org/10.2514/6.2023-1244 (freely available at https://www.google.com/url?sa=t&rct=j&q=&esrc=s&source=web&cd=&cad=rja&uact=8&ved=2ahUKEwiwm_D467P-AhUYQ_EDHV0LCV4QFnoECAoQAQ&url=https%3A%2F%2Fhighfidelitycfdverificationworkshop.github.io%2Fpapers%2Frans.pdf&usg=AOvVaw35mtEWtqFl4ACXiaQ4BUrQ).

I made a few changes to the code to update the current version of SA to match the one on the paper. In the following days, I'll prepare a pull request for a new branch to share such changes.

Any suggestion and/or help is welcome! Moreover, if someone is interested in participating in the Adaptive Grids RANS TFG, please let me know, such that I can update you on the state of the workshop.

[pcarruscag]

Feel free to pull me in, I can help with the settings for convergence if things get challenging. I have some experience running SU2 on the workshop meshes.

[rois1995]

Sure, thank you! Is there a place where I can upload the configs that we are using such that we all are on the same boat (maybe as a testcase)? Or should I create a repository on my github and then share it here?

[pcarruscag]

We could share them here, or create a slack channel for this topic

[bigfooted]

I tried the smallest mesh (1v) on my local machine last weekend, it was just small enough that I could run it with 16Gb of memory. That case converged very well with the existing SA options (in less than 1000 iterations I believe). I tried the second mesh (3v) today on the cluster while I was in the train and I also had some difficulty with that case. I managed to get a converged solution in the end. For me the first thing to check is that the setup is correct and such that we can extract the things we need for the validation.

[bigfooted]

by the way, was there a CL(alpha) curve for exactly this geometry? The computational geometry has changed since v4 and I couldn't find the experimental data for this model and I don't know if we can use the computational results of the last workshop as a reference. Knowing what the relevant CL and CD are will be very helpful.

[tbellosta]

Not sure of any experimental data, but there should be FUN3D mesh convergence results in the paper @rois1995 linked above. Test case 3 in the paper should be the same condition of the HLPW5 case

[bigfooted]

Here is my script below. It converges well for the 3v case. CL and CD are reported as 2.2 and 0.5, which is higher than in the report above, so I guess the scaling has to be done differently.

% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------%
%
SOLVER= RANS
KIND_TURB_MODEL= SA
SA_OPTIONS= NEGATIVE, QCR2000, ROTATION, EXPERIMENTAL
% option 'NONE' converges better, and the CFL can be increased
%SA_OPTIONS= NONE 
MATH_PROBLEM= DIRECT
RESTART_SOL= NO 
% mesh is in meters
SYSTEM_MEASUREMENTS= SI
%
% ---------------- COMPRESSIBLE FLOW CONDITION DEFINITION -------------------%
%
MACH_NUMBER= 0.2
AOA= 11.0
REYNOLDS_NUMBER= 5.6e6
INIT_OPTION= REYNOLDS
FREESTREAM_OPTION= TEMPERATURE_FS
% temperature is 521 R (17 Celsius = 289.44)
FREESTREAM_TEMPERATURE= 289.44
FREESTREAM_PRESSURE= 101325
FREESTREAM_VISCOSITY= 1.8e-5
% Reynolds length is MAC = 275.8 inch = 7.00532 m
REYNOLDS_LENGTH= 7.00532
%
% actual speed of sound is c = 341.057 m/s
% actual velocity is U = 68.21 m/s
%
% ---------------------- REFERENCE VALUE DEFINITION ---------------------------%
%
% Reference origin for moment computation
% moment reference center (MRC)
% 1325.9 inch = 33.67786 m
% 0 inch = 0 m
% 177.95 inch = 4.51993 m
REF_ORIGIN_MOMENT_X = 33.67786
REF_ORIGIN_MOMENT_Y = 0.00
REF_ORIGIN_MOMENT_Z = 4.51993
%
% Reference length for pitching, rolling, and yawing non-dimensional moment
% mean aerodynamic chord (in inches)
% 275.8 inch -> 7.00532 m
REF_LENGTH= 7.00532
%
% Reference area for force coefficients (0 implies automatic calculation)
% This is the reference area of the semi-span model (Sref/2) in inches^2
% Sref = 297360
% S = Sref/2 = 148680
% = 95.922 m^2
REF_AREA= 95.922
%
% Compressible flow non-dimensionalization (DIMENSIONAL, FREESTREAM_PRESS_EQ_ONE,
%                              FREESTREAM_VEL_EQ_MACH, FREESTREAM_VEL_EQ_ONE)
% reference pressure equals free stream pressure: p_ref=p_freestream
REF_DIMENSIONALIZATION= FREESTREAM_PRESS_EQ_ONE

%
% ---- IDEAL GAS, POLYTROPIC, VAN DER WAALS AND PENG ROBINSON CONSTANTS -------%
%
FLUID_MODEL= STANDARD_AIR
GAMMA_VALUE= 1.4

% --------------------------- VISCOSITY MODEL ---------------------------------%
%
VISCOSITY_MODEL= CONSTANT_VISCOSITY
MU_CONSTANT= 1.716E-5
%
% --------------------------- THERMAL CONDUCTIVITY MODEL ----------------------%
%
CONDUCTIVITY_MODEL= CONSTANT_PRANDTL
PRANDTL_LAM= 0.72

% -------------------- BOUNDARY CONDITION DEFINITION --------------------------%
%
%MARKER_HEATFLUX= ( TRI_bdy1, 0.0, TRI_bdy2, 0.0, TRI_bdy3, 0.0, TRI_bdy4, 0.0, TRI_bdy5, 0.0, \
%                   TRI_bdy6, 0.0, TRI_bdy7, 0.0, TRI_bdy8, 0.0, TRI_bdy9, 0.0, TRI_bdy10, 0.0 \
%                   QUAD_bdy3, 0.0, TRI_bdy11, 0.0, TRI_bdy12, 0.0, TRI_bdy13, 0.0, TRI_bdy14, 0.0 \
%                   TRI_bdy15, 0.0 \  
%                 )
% QUAD_bdy3 -> symmetry_quad
% TRI_bdy3 -> symmetry_tri
% TRI_bdy11 -> outlet
% TRI_bdy12 -> inlet
% TRI_bdy13 -> far_back
% TRI_bdy14 -> far_top
% TRI_bdy15 -> far_bottom

MARKER_HEATFLUX= ( TRI_bdy1, 0.0, TRI_bdy2, 0.0, TRI_bdy4, 0.0, TRI_bdy5, 0.0, \
                   TRI_bdy6, 0.0, TRI_bdy7, 0.0, TRI_bdy8, 0.0, TRI_bdy9, 0.0, TRI_bdy10, 0.0 \
                 )
MARKER_FAR= (inlet, far_top, far_bottom, far_back, outlet)
MARKER_SYM= (symmetry_quad, symmetry_tri)
%
MARKER_MONITORING= ( TRI_bdy1, TRI_bdy2, TRI_bdy4, TRI_bdy5, TRI_bdy6, TRI_bdy7, TRI_bdy8, TRI_bdy9, TRI_bdy10)
MARKER_PLOTTING= ( TRI_bdy1 )
MARKER_ANALYZE= ( inlet,outlet )
MARKER_ANALYZE_AVERAGE= MASSFLUX
%
% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------%
%
NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES
%NUM_METHOD_GRAD= GREEN_GAUSS
%ENTROPY_FIX_COEFF= 0.001
%SLOPE_LIMITER_FLOW= VAN_ALBADA_EDGE
 
%
CFL_NUMBER=50.0
%CFL_NUMBER=100.0
CFL_ADAPT= NO 
CFL_ADAPT_PARAM= ( 0.9, 1.1, 50.0, 250.0 ) 
ITER= 10000
%
% ------------------------ LINEAR SOLVER DEFINITION ---------------------------%
%
LINEAR_SOLVER= FGMRES
LINEAR_SOLVER_PREC= ILU
% low memory consumption (not compatible with discrete adjoint):
%LINEAR_SOLVER_PREC= LU_SGS
LINEAR_SOLVER_ERROR= 1E-05
%LINEAR_SOLVER_ITER= 15
% low memory consumption:
LINEAR_SOLVER_ITER= 5
%
% -------------------------- MULTIGRID PARAMETERS -----------------------------%
%

MGLEVEL= 1
%mglevel=2 leads to noise in the momentum residuals for case 3v
%MGLEVEL= 2
%
% Multi-grid cycle (V_CYCLE, W_CYCLE, FULLMG_CYCLE)
MGCYCLE= V_CYCLE
MG_PRE_SMOOTH= ( 1, 2, 3, 3 )
MG_POST_SMOOTH= ( 1, 1, 1, 1 )
MG_CORRECTION_SMOOTH= ( 1, 1, 1, 1 )
MG_DAMP_RESTRICTION= 0.5
MG_DAMP_PROLONGATION= 0.5
%
% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------%
%
CONV_NUM_METHOD_FLOW= ROE
JST_SENSOR_COEFF= (0.5, 0.02)
MUSCL_FLOW= NO 
TIME_DISCRE_FLOW= EULER_IMPLICIT
%
% -------------------- TURBULENT NUMERICAL METHOD DEFINITION ------------------%
%
CONV_NUM_METHOD_TURB= SCALAR_UPWIND
MUSCL_TURB= NO 
%SLOPE_LIMITER_TURB= VENKATAKRISHNAN
SLOPE_LIMITER_TURB= NONE 
TIME_DISCRE_TURB= EULER_IMPLICIT
%
% --------------------------- CONVERGENCE PARAMETERS --------------------------%
%
CONV_FIELD= RMS_DENSITY
CONV_RESIDUAL_MINVAL= -10
CONV_STARTITER= 10
CONV_CAUCHY_ELEMS= 100
CONV_CAUCHY_EPS= 1E-6
%
% ------------------------- INPUT/OUTPUT INFORMATION --------------------------%
%
MESH_FORMAT= SU2
%MESH_FILENAME= mesh_out.su2 
MESH_FILENAME= CRMHL_Wingbody_3v_SI.su2
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TRANSFORM TO SI                      % 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%MESH_FORMAT= CGNS
%MESH_FILENAME= CRMHL_Wingbody_3v.cgns 
%DV_KIND= SCALE_GRID
%OPT_RELAX_FACTOR= 1.0
%DV_PARAM= 1.0
%DV_VALUE= (0.0254)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
SOLUTION_FILENAME= solution
CONV_FILENAME= history
RESTART_FILENAME= restart
WRT_PERFORMANCE= YES
VOLUME_FILENAME= flow
OUTPUT_WRT_FREQ= 50
OUTPUT_FILES= PARAVIEW_MULTIBLOCK RESTART SURFACE_CSV
SCREEN_OUTPUT= INNER_ITER WALL_TIME RMS_DENSITY RMS_MOMENTUM-X RMS_MOMENTUM-Y RMS_MOMENTUM-Z RMS_ENERGY RMS_NU_TILDE LINSOL_ITER LINSOL_RESIDUAL LIFT DRAG
VOLUME_OUTPUT= SOLUTION PRIMITIVE RESIDUAL
HISTORY_OUTPUT= INNER_ITER RMS_RES FLOW_COEFF AERO_COEFF

[tbellosta]

I'm not sure about the /2 factor in the reference area. Maybe that's why results are so different. Also, the config is using a 1st order convection scheme, that may be another reason. We'll upload our config shortly

[rois1995]

The config that we are using is this one:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                              %
% SU2 configuration file                                                       %
% Case description: Turbulent flow, HLPW-5 TC1, Newton-Krylov solver             %
% File Version 7.5.1 "Blackbird"                                               %
%                                                                              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

SOLVER= RANS
KIND_TURB_MODEL= SA
SA_OPTIONS= (NEGATIVE, WITHFT2, QCR2000, ROTATION, EXPERIMENTAL)
MATH_PROBLEM= DIRECT

RESTART_SOL= NO

SYSTEM_MEASUREMENTS= SI

% ------------------------- NEWTON-KRYLOV PARAMETERS --------------------------%
%
% --- Things that matter ---
NEWTON_KRYLOV= YES

% Iterations and tolerance for the Krylov part, it is important not to
% "over solve", tolerance should be as high as possible.
LINEAR_SOLVER_ITER= 10
LINEAR_SOLVER_ERROR= 0.025

% For "n0" iterations or "r0" residual reduction, the normal quasi-Newton iterations
% are used. Then, they become the preconditioner for the NK iterations with "np" linear
% iterations or "tp" tolerance, with "np"=0 the linear preconditioner (e.g. ILU) is
% used directly (this may be enough for unsteady).
% The tolerance for NK iterations is initially relaxed by factor "ft", and reaches
% LINEAR_SOLVER_ERROR after "rf" residual reduction (additional to "r0").
% The Jacobian-free products are based on finite differences with step "e".
NEWTON_KRYLOV_IPARAM= (0, 5, 2) % n0, np, ft
NEWTON_KRYLOV_DPARAM= (-1.0, 0.05, -6.0, 1e-5) % r0, tp, rf, e

CFL_ADAPT= YES % it's needed
CFL_NUMBER= 10
CFL_ADAPT_PARAM= ( 0.1, 1.1, 0.001, 100000) % no point using NK with low CFL values

% It is important (more than usual) to have similar magnitude variables
REF_DIMENSIONALIZATION= FREESTREAM_VEL_EQ_MACH

TIME_DISCRE_FLOW= EULER_IMPLICIT
LINEAR_SOLVER_PREC= ILU

% --- Things that don't matter ---
MGLEVEL= 0 % NK replaces MG
LINEAR_SOLVER= FGMRES

% -------------------- COMPRESSIBLE FREE-STREAM DEFINITION --------------------%
%
MACH_NUMBER= 0.20
AOA= 11.0
SIDESLIP_ANGLE= 0.0

FREESTREAM_TEMPERATURE= 289.44
REYNOLDS_NUMBER= 5.6E6
REYNOLDS_LENGTH= 7.00532

% ---- IDEAL GAS, POLYTROPIC, VAN DER WAALS AND PENG ROBINSON CONSTANTS -------%
%
FLUID_MODEL= STANDARD_AIR
GAMMA_VALUE= 1.4
%GAS_CONSTANT= 287.058
%ACENTRIC_FACTOR= 0.035

% --------------------------- VISCOSITY MODEL ---------------------------------%
%
VISCOSITY_MODEL= SUTHERLAND
%MU_CONSTANT= 1.716E-5
%MU_REF= 1.716E-5
%MU_T_REF= 273.15
%SUTHERLAND_CONSTANT= 110.4

% ---------------------- REFERENCE VALUE DEFINITION ---------------------------%
%
REF_ORIGIN_MOMENT_X = 33.67786
REF_ORIGIN_MOMENT_Y = 0.00
REF_ORIGIN_MOMENT_Z = 4.51993
REF_LENGTH= 7.00532
REF_AREA= 191.844778

% -------------------- BOUNDARY CONDITION DEFINITION --------------------------%
%
MARKER_HEATFLUX= ( wall, 0.0 )
MARKER_FAR= ( farfield )
MARKER_SYM= ( tri_symmetry, quad_symmetry )
MARKER_PLOTTING= ( wall )
MARKER_MONITORING= ( wall )

% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------%
%
NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES
NUM_METHOD_GRAD_RECON= WEIGHTED_LEAST_SQUARES

CONV_NUM_METHOD_FLOW= ROE
MUSCL_FLOW= YES
SLOPE_LIMITER_FLOW= VAN_ALBADA_EDGE

% -------------------- TURBULENT NUMERICAL METHOD DEFINITION ------------------%
%
CONV_NUM_METHOD_TURB= SCALAR_UPWIND
MUSCL_TURB= NO
TIME_DISCRE_TURB= EULER_IMPLICIT

% --------------------------- CONVERGENCE PARAMETERS --------------------------%
%
CONV_RESIDUAL_MINVAL= -10
CONV_STARTITER= 10
ITER= 200000

% ------------------------- INPUT/OUTPUT INFORMATION --------------------------%
%
MESH_FILENAME= mesh.cgns
MESH_FORMAT= CGNS
TABULAR_FORMAT= CSV
SOLUTION_FILENAME= restart.dat
RESTART_FILENAME= restart.dat
SOLUTION_ADJ_FILENAME= solution_adj.dat
RESTART_ADJ_FILENAME= restart_adj.dat
VOLUME_FILENAME= flow
VOLUME_ADJ_FILENAME= adjoint
SURFACE_FILENAME= surface_flow
SURFACE_ADJ_FILENAME= surface_adjoint
OUTPUT_FILES=(RESTART)
OUTPUT_WRT_FREQ= (1000, 500)
SCREEN_OUTPUT = (INNER_ITER, WALL_TIME, RMS_DENSITY, RMS_ENERGY, RMS_NU_TILDE, LIFT, DRAG, LINSOL_ITER, LINSOL_RESIDUAL, AVG_CFL)
CONV_FILENAME= history
HISTORY_OUTPUT= (RMS_RES, WALL_TIME, ITER, LINSOL, CFL_NUMBER, AERO_COEFF)

UPDATED added WITHFT2 to config

Since we scaled the meshes within Pointwise, we also changed the boundary names to have only three boundaries.

The results until now are:

Mesh	CL	CD
v1	0.870151	0.137867
v3	0.907166	0.120253
v5	/	/
v7	1.112564	0.134821

Although CL and CD values are still changing a bit.

[rois1995]

The config that I am using is wrong. There should be the WITHFT2 within the SA_OPTIONS.

UPDATE Even if it is the correct way of using SA-NEG, the residuals shoot up and the simulations diverge.

[pcarruscag]

For the HLCRM this is what I was using https://github.com/pcarruscag/SU2/tree/ENG_NOTE_ADJOINTS/Paper_Data/HLCRM
It took quite a few iterations for the finer grids, up to 40k (without Newton Krylov and around 10k with) on the extra fine I think.
Remember that you can restart from coarse solutions for this, it helps a lot.
I did not run the higher AoA with SU2.

At the time I also got some solutions with the Roe scheme but it took longer to converge and for high lift conditions I don't think it matters much.
I was still using Newton-Krylov but with lower CFL around 25 and NEWTON_KRYLOV_IPARAM= (0, 0, 1) and aiming for ~0.2 linear residual which from memory needed about 10-15 iterations.
The rest of the settings were similar, Green-Gauss, either van Albada or Venkat-Wang, and I did not find any benefit in increasing the entropy correction.

[rois1995]

I tried restarting the solution for the v3 from the v1 but it diverged after the first few iterations. I'll try to change some settings in the config and see if I can make it work.

We also tried 0.25 as linear residual but we had divergence problems. Now we tried lowering it to 0.025 and we achieved somewhat better results.

[aeroamit]

Give a try to SLAU2 scheme. Also do not turn on USE ACCURATE FLUX JACOBIAN option (keep it default to NO with SLAU2), in my experience you can probably go with higher CFL with these (I can't explain this peculiar behaviour at lower Mach numbers).
Also SLAU2 is less dissipative than Roe scheme for low Mach numbers (better here for M=0.2).
Just a suggestion... if you would like to give above a try.

[aeroamit]

Also Multigrid off. Newton Krylov may be interesting, I could do ONERA M6 (3 million grid) case in 2k iterations with that. Low speed case may take more iterations

[rois1995]

I am planning on running some meshes with both JST and SLAU2 to check for convergence differences. Thank you for the tip!

[rois1995]

Tecplot is providing a free license valid until the end of the workshop. It will be useful since we will be given Tecplot scripts to automate the post-processing of the results. If you are interested, please let me know so that I can send your names to the TFG leader.

[rois1995]

The branch that I am using is feature_SA_HLPW5. There are some changes like the one for computing Shat and some corrections on the R variant. Should I create a pull request?

[pcarruscag]

Yes if you found issues in the code we should correct them

[rois1995]

Just a quick update: I cannot make anything past the v5 (included) converge. I tried the following approaches:

1. as numerical schemes ROE, JST, SLAU2.
2. I tried restarting from a converged v3.
3. I tried reducing the CFL for the turbulence to 0.1

It is always the turbulence equation that diverges first. I may try further reducing the CFL, but this is a strange phenomenon for sure.

Anyway, today I cannot follow the Fixed Grid meeting, is anyone able to do it?

[bigfooted]

I now have this:

This is for the Roe scheme without MUSCL, so the Cl,Cd are overpredicted by a lot. These simulations all converge very well. For MUSCL, or JST, I really have to play with the settings to make it converge and I haven't been able to get the entire Cl,Cd convergence for all points yet for another numerical scheme.

[rois1995]

Well that seems nice at least for the first order. Where is the function addclippedsolution? I could not find it anywhere. Also, I will try a relaxation on the linear solver too, maybe that can help.

[bigfooted]

Here it is, I played with multiplying the relaxation by a value between 0.1..0.5 or using fmin(0.5,underrelaxationfactor) :

$ ~/Codes/su2_github/SU2$ grep -niIr 'addclippedsolution' *
SU2_CFD/include/solvers/CScalarSolver.inl:504:          nodes->AddClippedSolution(iPoint, iVar, nodes->GetUnderRelaxation(iPoint) * LinSysSol(iPoint, iVar),
SU2_CFD/include/solvers/CScalarSolver.inl:513:          nodes->AddClippedSolution(iPoint, iVar, nodes->GetUnderRelaxation(iPoint) * LinSysSol(iPoint, iVar),
SU2_CFD/include/variables/CVariable.hpp:452:  inline void AddClippedSolution(unsigned long iPoint, unsigned long iVar, su2double solution,

[bigfooted]

Oh there's also addSolution() in CFVMFlowSolverBase.hpp, where I tried the same, but I think this one is less important.

[bigfooted]

Also got some results for the JST scheme now, which is much closer to the reported values for CL and CD (the red dot in the figure)