Compressibility correction for SST model

Common/include/option_structure.hpp

@@ -990,6 +990,8 @@ enum class SST_OPTIONS {
   V,           /*!< \brief Menter k-w SST model with vorticity production terms. */
   KL,          /*!< \brief Menter k-w SST model with Kato-Launder production terms. */
   UQ,          /*!< \brief Menter k-w SST model with uncertainty quantification modifications. */
+  COMP_Wilcox, /*!< \brief Menter k-w SST model with Compressibility correction of Wilcox. */
+  COMP_Sarkar, /*!< \brief Menter k-w SST model with Compressibility correction of Sarkar. */
 };
 static const MapType<std::string, SST_OPTIONS> SST_Options_Map = {
   MakePair("NONE", SST_OPTIONS::NONE)
@@ -1002,6 +1004,8 @@ static const MapType<std::string, SST_OPTIONS> SST_Options_Map = {
   MakePair("VORTICITY", SST_OPTIONS::V)
   MakePair("KATO-LAUNDER", SST_OPTIONS::KL)
   MakePair("UQ", SST_OPTIONS::UQ)
+  MakePair("COMPRESSIBILITY-WILCOX", SST_OPTIONS::COMP_Wilcox)
+  MakePair("COMPRESSIBILITY-SARKAR", SST_OPTIONS::COMP_Sarkar)
 };
 
 /*!
@@ -1013,6 +1017,8 @@ struct SST_ParsedOptions {
   bool sust = false;                          /*!< \brief Bool for SST model with sustaining terms. */
   bool uq = false;                            /*!< \brief Bool for using uncertainty quantification. */
   bool modified = false;                      /*!< \brief Bool for modified (m) SST model. */
+  bool compWilcox = false;                    /*!< \brief Bool for compressibility correction of Wilcox. */
+  bool compSarkar = false;                    /*!< \brief Bool for compressibility correction of Sarkar. */
 };
 
 /*!
@@ -1048,6 +1054,8 @@ inline SST_ParsedOptions ParseSSTOptions(const SST_OPTIONS *SST_Options, unsigne
   const bool sst_v = IsPresent(SST_OPTIONS::V);
   const bool sst_kl = IsPresent(SST_OPTIONS::KL);
   const bool sst_uq = IsPresent(SST_OPTIONS::UQ);
+  const bool sst_compWilcox = IsPresent(SST_OPTIONS::COMP_Wilcox);
+  const bool sst_compSarkar = IsPresent(SST_OPTIONS::COMP_Sarkar);
 
   if (sst_1994 && sst_2003) {
     SU2_MPI::Error("Two versions (1994 and 2003) selected for SST_OPTIONS. Please choose only one.", CURRENT_FUNCTION);
@@ -1066,11 +1074,17 @@ inline SST_ParsedOptions ParseSSTOptions(const SST_OPTIONS *SST_Options, unsigne
     SSTParsedOptions.production = SST_OPTIONS::KL;
   } else if (sst_uq) {
     SSTParsedOptions.production = SST_OPTIONS::UQ;
+  } else if (sst_compWilcox) {
+    SSTParsedOptions.production = SST_OPTIONS::COMP_Wilcox;
+  } else if (sst_compSarkar) {
+    SSTParsedOptions.production = SST_OPTIONS::COMP_Sarkar;
   }
 
   SSTParsedOptions.sust = sst_sust;
   SSTParsedOptions.modified = sst_m;
   SSTParsedOptions.uq = sst_uq;
+  SSTParsedOptions.compWilcox = sst_compWilcox;
+  SSTParsedOptions.compSarkar = sst_compSarkar;
   return SSTParsedOptions;
 }
 

Common/src/CConfig.cpp

@@ -6170,6 +6170,12 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) {
                 cout << "\nperturbing the Reynold's Stress Matrix towards " << eig_val_comp << " component turbulence";
                 if (uq_permute) cout << " (permuting eigenvectors)";
                 break;
+              case SST_OPTIONS::COMP_Wilcox:
+                cout << " with compressibility correction of Wilcox";
+                break;
+              case SST_OPTIONS::COMP_Sarkar:
+                cout << " with compressibility correction of Sarkar";
+                break;
               default:
                 cout << " with no production modification";
                 break;

SU2_CFD/include/numerics/turbulent/turb_sources.hpp

@@ -740,6 +740,7 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
     AD::SetPreaccIn(Vorticity_i, 3);
     AD::SetPreaccIn(V_i[idx.Density()], V_i[idx.LaminarViscosity()], V_i[idx.EddyViscosity()]);
     AD::SetPreaccIn(V_i[idx.Velocity() + 1]);
+    AD::SetPreaccIn(V_i[idx.SoundSpeed()]);
 
     Density_i = V_i[idx.Density()];
     Laminar_Viscosity_i = V_i[idx.LaminarViscosity()];
@@ -778,6 +779,8 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
 
       const su2double VorticityMag = GeometryToolbox::Norm(3, Vorticity_i);
       su2double P_Base = 0;
+      su2double zetaFMt = 0.0, Dilatation_Sarkar = 0.0;
+      const su2double Mt = sqrt(2.0 * ScalarVar_i[0]) / V_i[idx.SoundSpeed()];
 
       /*--- Apply production term modifications ---*/
       switch (sstParsedOptions.production) {
@@ -795,6 +798,20 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
           P_Base = sqrt(StrainMag_i*VorticityMag);
           break;
 
+        case SST_OPTIONS::COMP_Wilcox:
+          P_Base = StrainMag_i;
+          if (Mt >= 0.25) {
+            zetaFMt = 2.0 * (Mt * Mt - 0.25 * 0.25);
+          }
+          break;
+
+        case SST_OPTIONS::COMP_Sarkar:
+          P_Base = StrainMag_i;
+          if( Mt >= 0.25){
+            zetaFMt = 0.5 * (Mt * Mt);
+          }
+          break;
+
         default:
           /*--- Base production term for SST-1994 and SST-2003 ---*/
           P_Base = StrainMag_i;
@@ -831,10 +848,15 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
         pw = max(pw, sust_w);
       }
 
+      if(sstParsedOptions.production == SST_OPTIONS::COMP_Sarkar) {
+        Dilatation_Sarkar = -0.15 * pk * Mt + 0.2 * beta_star * (1.0 +zetaFMt) * Density_i * ScalarVar_i[1] * ScalarVar_i[0] * Mt * Mt;
+        pk += Dilatation_Sarkar;
+      }
+
       /*--- Dissipation ---*/
 
-      su2double dk = beta_star * Density_i * ScalarVar_i[1] * ScalarVar_i[0];
-      su2double dw = beta_blended * Density_i * ScalarVar_i[1] * ScalarVar_i[1];
+      su2double dk = beta_star * Density_i * ScalarVar_i[1] * ScalarVar_i[0] * (1.0 + zetaFMt);
+      su2double dw = beta_blended * Density_i * ScalarVar_i[1] * ScalarVar_i[1] * (1.0 - 0.09/beta_blended * zetaFMt);
 
       /*--- LM model coupling with production and dissipation term for k transport equation---*/
       if (config->GetKind_Trans_Model() == TURB_TRANS_MODEL::LM) {
@@ -862,10 +884,10 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
 
       /*--- Implicit part ---*/
 
-      Jacobian_i[0][0] = -beta_star * ScalarVar_i[1] * Volume;
-      Jacobian_i[0][1] = -beta_star * ScalarVar_i[0] * Volume;
+      Jacobian_i[0][0] = -beta_star * ScalarVar_i[1] * Volume * (1.0 + zetaFMt);
+      Jacobian_i[0][1] = -beta_star * ScalarVar_i[0] * Volume * (1.0 + zetaFMt);
       Jacobian_i[1][0] = 0.0;
-      Jacobian_i[1][1] = -2.0 * beta_blended * ScalarVar_i[1] * Volume;
+      Jacobian_i[1][1] = -2.0 * beta_blended * ScalarVar_i[1] * Volume * (1.0 - 0.09/beta_blended * zetaFMt);
     }
 
     AD::SetPreaccOut(Residual, nVar);

TestCases/parallel_regression.py

@@ -392,6 +392,22 @@ def main():
     turb_flatplate_species.test_vals = [-4.147727, -0.634899, -1.770894, 1.334987, -3.250340, 9.000000, -6.700853, 5.000000, -6.991055, 10.000000, -6.033829, 0.996033, 0.996033]
     test_list.append(turb_flatplate_species)
 
+    # Flat plate SST compressibility  correction Wilcox
+    turb_flatplate_CC_Wilcox = TestCase('turb_flatplate_CC_Wilcox')
+    turb_flatplate_CC_Wilcox.cfg_dir   = "rans/flatplate"
+    turb_flatplate_CC_Wilcox.cfg_file  = "turb_SST_flatplate_compressibility_Wilcox.cfg"
+    turb_flatplate_CC_Wilcox.test_iter = 20
+    turb_flatplate_CC_Wilcox.test_vals = [-1.280847, 1.974242, 1.440510, 5.038429, -4.051126, 8.520857]
+    test_list.append(turb_flatplate_CC_Wilcox)
+
+    # Flat plate SST compressibility  correction Sarkar
+    turb_flatplate_CC_Sarkar = TestCase('turb_flatplate_CC_Sarkar')
+    turb_flatplate_CC_Sarkar.cfg_dir   = "rans/flatplate"
+    turb_flatplate_CC_Sarkar.cfg_file  = "turb_SST_flatplate_compressibility_Sarkar.cfg"
+    turb_flatplate_CC_Sarkar.test_iter = 20
+    turb_flatplate_CC_Sarkar.test_vals = [-1.280847, 1.974242, 1.440510, 5.038429, -4.051129, 8.520857]
+    test_list.append(turb_flatplate_CC_Sarkar)
+
     # ONERA M6 Wing
     turb_oneram6           = TestCase('turb_oneram6')
     turb_oneram6.cfg_dir   = "rans/oneram6"

TestCases/rans/flatplate/turb_SST_flatplate_compressibility_Sarkar.cfg

@@ -0,0 +1,150 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%                                                                              %
+% SU2 configuration file                                                       %
+% Case description: Turbulent flow over flat plate(2DZPH)					   %
+% Author: Sunoh. Kang														   %
+% Institution: Pusan National University                                       %
+% Date: 2024.04.30                                                             %
+% File Version 8.0.1 "Harrier"                                                 %
+%                                                                              %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------%
+%
+SOLVER= RANS
+KIND_TURB_MODEL= SST
+SST_OPTIONS= V2003m, COMPRESSIBILITY-SARKAR
+MATH_PROBLEM= DIRECT
+RESTART_SOL= NO
+
+% ------------------------ COMPRESSIBLE FREE-STREAM DEFINITION ----------------%
+%
+MACH_NUMBER= 5.0
+AOA= 0.000
+INIT_OPTION= TD_CONDITIONS
+FREESTREAM_OPTION= TEMPERATURE_FS
+FREESTREAM_TEMPERATURE= 300.3333
+FREESTREAM_PRESSURE= 13753.90558
+FREESTREAM_TURBULENCEINTENSITY = 0.00002
+FREESTREAM_TURB2LAMVISCRATIO = 0.009
+REF_DIMENSIONALIZATION= DIMENSIONAL
+REYNOLDS_NUMBER= 15E6
+
+% ---- IDEAL GAS, POLYTPOPIC, VAN DER WALLS AND PENG ROBINSOPN CNSTANTS ------ %
+%
+FLUID_MODEL= IDEAL_GAS
+GAMMA_VALUE= 1.4
+GAS_CONSTANT= 287.058
+ACENTRIC_FACTOR= 0.035
+SPECIFIC_HEAT_CP= 1004.703
+THERMAL_EXPANSION_COEFF= 0.00347
+MOLECULAR_WEIGHT= 28.96
+
+% --------------------------- VISCOSITY MODEL ---------------------------------%
+%
+
+VISCOSITY_MODEL= SUTHERLAND
+MU_REF= 1.716E-5
+MU_T_REF= 273.15
+SUTHERLAND_CONSTANT= 110.4
+
+% ------------------------ THERMAL CONDUCTIVITY MODEL -------------------------%
+%
+
+CONDUCTIVITY_MODEL= CONSTANT_PRANDTL
+PRANDTL_LAM= 0.72
+TURBULENT_CONDUCTIVITY_MODEL= CONSTANT_PRANDTL_TURB
+PRANDTL_TURB= 0.9
+
+% ---------------------- REFERENCE VALUE DEFINITION ---------------------------%
+%
+
+REF_ORIGIN_MOMENT_X = 0.00
+REF_ORIGIN_MOMENT_Y = 0.00
+REF_ORIGIN_MOMENT_Z = 0.00
+REF_LENGTH= 1.0
+REF_AREA= 1.0
+
+% -------------------- BOUNDARY CONDITION DEFINITION --------------------------%
+%
+
+MARKER_PLOTTING = ( wall )
+MARKER_SUPERSONIC_INLET= ( inlet, 300.3333 , 13753.90558 , 1737.082226 , 0.0 , 0.0 )
+MARKER_FAR= ( farfield )
+MARKER_SYM= ( symmetry )
+MARKER_SUPERSONIC_OUTLET= (outlet, 34384.76396)
+MARKER_ISOTHERMAL= ( wall, 327.36297)
+
+% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------%
+%
+
+NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES
+CFL_NUMBER= 0.5
+CFL_ADAPT= YES
+CFL_ADAPT_PARAM= ( 0.1, 1.5, 0.5, 50.0 )
+
+% --------- SLOPE LIMITER DEFINITION, DISSIPATION SENSOR DEFINITION -----------------%
+%
+
+MUSCL_FLOW= YES 
+SLOPE_LIMITER_FLOW= VENKATAKRISHNAN_WANG
+SLOPE_LIMITER_TURB= NONE
+VENKAT_LIMITER_COEFF= 0.05
+JST_SENSOR_COEFF= (0.5, 0.02)
+
+% ------------------------ LINEAR SOLVER DEFINITION ---------------------------%
+%
+
+LINEAR_SOLVER= FGMRES
+LINEAR_SOLVER_PREC= LU_SGS
+LINEAR_SOLVER_ERROR= 1E-6
+LINEAR_SOLVER_ITER= 5
+
+% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------%
+%
+
+CONV_NUM_METHOD_FLOW= AUSM
+TIME_DISCRE_FLOW= EULER_IMPLICIT
+
+% -------------------- TURBULENT NUMERICAL METHOD DEFINITION ------------------%
+%
+
+CONV_NUM_METHOD_TURB= SCALAR_UPWIND
+TIME_DISCRE_TURB= EULER_IMPLICIT
+
+% ------------------------- SCREEN/HISTORY VOLUME OUTPUT --------------------------%
+%
+
+SCREEN_OUTPUT= (INNER_ITER, RMS_DENSITY, RMS_MOMENTUM-X, RMS_MOMENTUM-Y, RMS_ENERGY, RMS_TKE, RMS_DISSIPATION)
+SCREEN_WRT_FREQ_INNER= 1
+SCREEN_WRT_FREQ_OUTER= 1 
+SCREEN_WRT_FREQ_TIME= 1
+HISTORY_WRT_FREQ_INNER= 1
+HISTORY_WRT_FREQ_OUTER= 1 
+HISTORY_WRT_FREQ_TIME= 1
+OUTPUT_WRT_FREQ= 1000
+
+% ------------------------- INPUT/OUTPUT FILE INFORMATION --------------------------%
+%
+
+MESH_FILENAME= mesh_flatplate_turb_137x97.su2
+MESH_FORMAT= SU2
+MESH_OUT_FILENAME= mesh_out.su2
+SOLUTION_FILENAME= restart_flow.dat
+SOLUTION_ADJ_FILENAME= solution_adj.dat
+OUTPUT_FILES= (RESTART, PARAVIEW, SURFACE_PARAVIEW, TECPLOT, SURFACE_TECPLOT, CSV, SURFACE_CSV )
+CONV_FILENAME= history
+BREAKDOWN_FILENAME= forces_breakdown.dat
+RESTART_FILENAME= restart_flow.dat
+RESTART_ADJ_FILENAME= restart_adj.dat
+
+% ------------------------------- SOLVER CONTROL ------------------------------%
+%
+
+ITER= 100
+OUTER_ITER= 1
+TIME_ITER= 100
+CONV_RESIDUAL_MINVAL= -14
+CONV_STARTITER= 10
+CONV_CAUCHY_ELEMS= 100
+CONV_CAUCHY_EPS= 1E-10