Adding new roughness boundary conditions for SST (#2573)

* Boundary conditions for SST for rough walls

* Fix formatting

* Fix issues about rough SST BCs formulations

* Add SST Limiter Output

* Revert SST Limiter Output changes

* Switch to regression branch

* Fix rough BCs options issues

* Added rough flatplate test cases

* Fixed values regression test rough flat plate

* Fix regression test roughness

* Cleaned regression tests for rough flat plate and moved to parallel regression

* [skip ci]

* Apply suggestions from code review

* Apply suggestions from code review

* [skip ci]

---------

Co-authored-by: Pedro Gomes <38071223+pcarruscag@users.noreply.github.com>

Author: pcarruscag

Common/include/CConfig.hpp

@@ -763,6 +763,7 @@ class CConfig {
   SST_OPTIONS *SST_Options;           /*!< \brief List of modifications/corrections/versions of SST turbulence model.*/
   SA_OPTIONS *SA_Options;             /*!< \brief List of modifications/corrections/versions of SA turbulence model.*/
   LM_OPTIONS *LM_Options;             /*!< \brief List of modifications/corrections/versions of SA turbulence model.*/
+  ROUGHSST_MODEL Kind_RoughSST_Model;         /*!< \brief List of modifications/corrections/versions of rough-wall boundary conditions for SST turbulence model.*/
   unsigned short nSST_Options;        /*!< \brief Number of SST options specified. */
   unsigned short nSA_Options;         /*!< \brief Number of SA options specified. */
   unsigned short nLM_Options;         /*!< \brief Number of SA options specified. */
@@ -1226,6 +1227,7 @@ class CConfig {
   SST_ParsedOptions sstParsedOptions; /*!< \brief Additional parameters for the SST turbulence model. */
   SA_ParsedOptions saParsedOptions;   /*!< \brief Additional parameters for the SA turbulence model. */
   LM_ParsedOptions lmParsedOptions;   /*!< \brief Additional parameters for the LM transition model. */
+  ROUGH_SST_ParsedOptions roughsstParsedOptions; /*!< \brief Additional parameters for the boundary conditions for rough walls for the SST turbulence model. */
   su2double uq_delta_b;         /*!< \brief Parameter used to perturb eigenvalues of Reynolds Stress Matrix */
   unsigned short eig_val_comp;  /*!< \brief Parameter used to determine type of eigenvalue perturbation */
   su2double uq_urlx;            /*!< \brief Under-relaxation factor */
@@ -10173,6 +10175,11 @@ class CConfig {
    */
   LM_ParsedOptions GetLMParsedOptions() const { return lmParsedOptions; }
 
+  /*!
+   * \brief Get parsed rough-wall boundary conditions for SST option data structure.
+   * \return Rough-wall SST option data structure.
+   */
+  ROUGH_SST_ParsedOptions GetROUGHSSTParsedOptions() const { return roughsstParsedOptions; }
 
   /*!
    * \brief Get parsed option data structure for data-driven fluid model.

Common/include/option_structure.hpp

@@ -1101,6 +1101,58 @@ inline SST_ParsedOptions ParseSSTOptions(const SST_OPTIONS *SST_Options, unsigne
   return SSTParsedOptions;
 }
 
+/*!
+ * \brief SST rough-wall boundary conditions Options
+ */
+enum class ROUGHSST_MODEL {
+  NONE,                 /*!< \brief No option / default no surface roughness applied. */
+  WILCOX1998,           /*!< \brief Wilcox 1998 boundary conditions for rough walls. */
+  WILCOX2006,           /*!< \brief Wilcox 2006 boundary conditions for rough walls / default version if roughness is applied. */
+  LIMITER_KNOPP,        /*!< \brief Knopp eddy viscosity limiter. */
+  LIMITER_AUPOIX,       /*!< \brief Aupoix eddy viscosity limiter. */
+};
+static const MapType<std::string, ROUGHSST_MODEL> RoughSST_Model_Map = {
+  MakePair("NONE", ROUGHSST_MODEL::NONE)
+  MakePair("WILCOX1998", ROUGHSST_MODEL::WILCOX1998)
+  MakePair("WILCOX2006", ROUGHSST_MODEL::WILCOX2006)
+  MakePair("LIMITER_KNOPP", ROUGHSST_MODEL::LIMITER_KNOPP)
+  MakePair("LIMITER_AUPOIX", ROUGHSST_MODEL::LIMITER_AUPOIX)
+};
+
+/*!
+ * \brief Structure containing parsed SST rough-wall boundary conditions options.
+ */
+struct ROUGH_SST_ParsedOptions {
+  ROUGHSST_MODEL version = ROUGHSST_MODEL::NONE;         /*!< \brief KWBC base model. */
+  bool wilcox1998 = false;                           /*!< \brief Use Wilcox boundary conditions for rough walls (1998). */
+  bool wilcox2006 = false;                          /*!< \brief Use Wilcox boundary conditions for rough walls (2006). */
+  bool limiter_knopp = false;                        /*!< \brief Use Knopp eddy viscosity limiter. */
+  bool limiter_aupoix = false;                      /*!< \brief Use Aupoix eddy viscosity limiter. */
+};
+
+/*!
+ * \brief Function to parse SST rough-wall boundary conditions options.
+ * \param[in] ROUGHSST_Options - Selected SST rough-wall boundary conditions option from config.
+ * \param[in] nROUGHSST_Options - Number of options selected.
+ * \param[in] rank - MPI rank.
+ * \return Struct with SST options.
+ */
+inline ROUGH_SST_ParsedOptions ParseROUGHSSTOptions(ROUGHSST_MODEL sstbcs_option) {
+    ROUGH_SST_ParsedOptions opts;
+    opts.version = sstbcs_option;
+
+    switch(sstbcs_option) {
+        case ROUGHSST_MODEL::WILCOX1998: opts.wilcox1998 = true; break;
+        case ROUGHSST_MODEL::WILCOX2006: opts.wilcox2006 = true; break;
+        case ROUGHSST_MODEL::LIMITER_KNOPP: opts.limiter_knopp = true; break;
+        case ROUGHSST_MODEL::LIMITER_AUPOIX: opts.limiter_aupoix = true; break;
+        default: break;
+    }
+
+    return opts;
+}
+
+
 /*!
  * \brief SA Options
  */

Common/src/CConfig.cpp

@@ -1140,6 +1140,8 @@ void CConfig::SetConfig_Options() {
   /*!\brief SST_OPTIONS \n DESCRIPTION: Specify SA turbulence model options/corrections. \n Options: see \link SA_Options_Map \endlink \n DEFAULT: NONE \ingroup Config*/
   addEnumListOption("SA_OPTIONS", nSA_Options, SA_Options, SA_Options_Map);
 
+  /*!\brief ROUGHSST_OPTIONS \n DESCRIPTION: Specify type of boundary condition for rough walls for SST turbulence model. \n Options: see \link ROUGHSST_Options_Map \endlink \n DEFAULT: wilcox1998 \ingroup Config*/
+  addEnumOption("KIND_ROUGHSST_MODEL", Kind_RoughSST_Model, RoughSST_Model_Map, ROUGHSST_MODEL::NONE);
   /*!\brief KIND_TRANS_MODEL \n DESCRIPTION: Specify transition model OPTIONS: see \link Trans_Model_Map \endlink \n DEFAULT: NONE \ingroup Config*/
   addEnumOption("KIND_TRANS_MODEL", Kind_Trans_Model, Trans_Model_Map, TURB_TRANS_MODEL::NONE);
   /*!\brief SST_OPTIONS \n DESCRIPTION: Specify LM transition model options/correlations. \n Options: see \link LM_Options_Map \endlink \n DEFAULT: NONE \ingroup Config*/
@@ -3598,6 +3600,7 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
   /*--- Postprocess SST_OPTIONS into structure. ---*/
   if (Kind_Turb_Model == TURB_MODEL::SST) {
     sstParsedOptions = ParseSSTOptions(SST_Options, nSST_Options, rank);
+    roughsstParsedOptions = ParseROUGHSSTOptions(Kind_RoughSST_Model);
   } else if (Kind_Turb_Model == TURB_MODEL::SA) {
     saParsedOptions = ParseSAOptions(SA_Options, nSA_Options, rank);
   }

SU2_CFD/include/solvers/CTurbSSTSolver.hpp

@@ -39,6 +39,7 @@ class CTurbSSTSolver final : public CTurbSolver {
 private:
   su2double constants[11] = {0.0}; /*!< \brief Constants for the model. */
   SST_ParsedOptions sstParsedOptions;
+  ROUGH_SST_ParsedOptions roughsstParsedOptions;
 
   /*!
    * \brief Compute nu tilde from the wall functions.

SU2_CFD/src/solvers/CTurbSSTSolver.cpp

@@ -425,9 +425,12 @@ void CTurbSSTSolver::BC_HeatFlux_Wall(CGeometry *geometry, CSolver **solver_cont
 
     /*--- Check if the node belongs to the domain (i.e, not a halo node) ---*/
     if (geometry->nodes->GetDomain(iPoint)) {
+      const auto options = config->GetROUGHSSTParsedOptions();
 
-      if (rough_wall) {
+      /*--- distance to closest neighbor ---*/
+      su2double wall_dist = geometry->vertex[val_marker][iVertex]->GetNearestNeighborDistance();
 
+      if (rough_wall) {
         /*--- Set wall values ---*/
         su2double density = solver_container[FLOW_SOL]->GetNodes()->GetDensity(iPoint);
         su2double laminar_viscosity = solver_container[FLOW_SOL]->GetNodes()->GetLaminarViscosity(iPoint);
@@ -441,33 +444,54 @@ void CTurbSSTSolver::BC_HeatFlux_Wall(CGeometry *geometry, CSolver **solver_cont
         su2double kPlus = FrictionVel*Roughness_Height*density/laminar_viscosity;
 
         su2double S_R= 0.0;
+        su2double solution[2] = {};
+        /*--- Modify the omega and k to account for a rough wall. ---*/
+
         /*--- Reference 1 original Wilcox (1998) ---*/
-        /*if (kPlus <= 25)
+        if (options.wilcox1998) {
+          if (kPlus <= 25)
             S_R = (50/(kPlus+EPS))*(50/(kPlus+EPS));
           else
-            S_R = 100/(kPlus+EPS);*/
-
+          S_R = 100/(kPlus+EPS);
+          
+          solution[0] = 0.0;
+          solution[1] = FrictionVel*FrictionVel*S_R/(laminar_viscosity/density);          
+        } 
+        else if (options.wilcox2006) {
         /*--- Reference 2 from D.C. Wilcox Turbulence Modeling for CFD (2006) ---*/
-        if (kPlus <= 5)
-          S_R = (200/(kPlus+EPS))*(200/(kPlus+EPS));
-        else
-          S_R = 100/(kPlus+EPS) + ((200/(kPlus+EPS))*(200/(kPlus+EPS)) - 100/(kPlus+EPS))*exp(5-kPlus);
-
-        /*--- Modify the omega to account for a rough wall. ---*/
-        su2double solution[2];
-        solution[0] = 0.0;
-        solution[1] = FrictionVel*FrictionVel*S_R/(laminar_viscosity/density);
-
+          if (kPlus <= 5)
+            S_R = pow(200/(kPlus+EPS),2);
+          else
+            S_R = 100/(kPlus+EPS) + (pow(200/(kPlus+EPS),2) - 100/(kPlus+EPS))*exp(5-kPlus);
+          
+          solution[0] = 0.0;
+          solution[1] = FrictionVel*FrictionVel*S_R/(laminar_viscosity/density);
+        } 
+        /*--- Knopp eddy viscosity limiter ---*/
+        else if (options.limiter_knopp) {
+          su2double d0 = 0.03*Roughness_Height*min(1.0, pow((kPlus + EPS )/30.0, 2.0/3.0))*min(1.0, pow((kPlus + EPS)/45.0, 0.25))*min(1.0, pow((kPlus + EPS) /60, 0.25));
+          solution[0] = (FrictionVel*FrictionVel / sqrt(constants[6]))*min(1.0, kPlus / 90.0);
+
+          const su2double kappa = config->GetwallModel_Kappa();
+          su2double beta_1 = constants[4];
+          solution[1] = min( FrictionVel/(sqrt(constants[6])*d0*kappa), 60.0*laminar_viscosity/(density*beta_1*pow(wall_dist,2))); 
+        }
+        /*--- Aupoix eddy viscosity limiter ---*/ 
+        else if (options.limiter_aupoix) {
+          su2double k0Plus = ( 1.0 /sqrt( constants[6])) * tanh((log10((kPlus +EPS ) / 30.0) + 1.0 - 1.0*tanh( (kPlus + EPS) / 125.0))*tanh((kPlus + EPS) / 125.0));
+          su2double kwallPlus = max(0.0, k0Plus);
+          su2double kwall = kwallPlus*FrictionVel*FrictionVel;
+      
+          su2double omegawallPlus = (300.0 / pow(kPlus + EPS, 2.0)) * pow(tanh(15.0 / (4.0*kPlus)), -1.0) + (191.0 / (kPlus + EPS))*(1.0 - exp(-kPlus / 250.0));
+
+          solution[0] = kwall;
+          solution[1] = omegawallPlus*FrictionVel*FrictionVel*density/laminar_viscosity;
+        }
         /*--- Set the solution values and zero the residual ---*/
         nodes->SetSolution_Old(iPoint,solution);
         nodes->SetSolution(iPoint,solution);
         LinSysRes.SetBlock_Zero(iPoint);
-
       } else { // smooth wall
-
-        /*--- distance to closest neighbor ---*/
-        su2double wall_dist = geometry->vertex[val_marker][iVertex]->GetNearestNeighborDistance();
-
         /*--- Set wall values ---*/
         su2double density = solver_container[FLOW_SOL]->GetNodes()->GetDensity(iPoint);
         su2double laminar_viscosity = solver_container[FLOW_SOL]->GetNodes()->GetLaminarViscosity(iPoint);

TestCases/parallel_regression.py

@@ -408,6 +408,22 @@ def main():
     turb_flatplate_CC_Sarkar.test_vals = [-1.195053, 2.089306, 1.529063, 5.164703, -3.700917, 8.162921]
     test_list.append(turb_flatplate_CC_Sarkar)
 
+    # FLAT PLATE, ROUGHNESS BC KNOPP SST
+    turb_flatplate_sst_roughBCKnopp           = TestCase('turb_sst_flatplate_roughBCKnopp')
+    turb_flatplate_sst_roughBCKnopp.cfg_dir   = "rans/flatplate/roughness/bc_knopp"
+    turb_flatplate_sst_roughBCKnopp.cfg_file  = "turb_SST_flatplate_roughBCKnopp.cfg"
+    turb_flatplate_sst_roughBCKnopp.test_iter = 10
+    turb_flatplate_sst_roughBCKnopp.test_vals = [-5.058634, -2.460850, -2.847064, 0.447200, -2.595042, 1.497149, -0.188079, 0.004571]
+    test_list.append(turb_flatplate_sst_roughBCKnopp)
+
+    # FLAT PLATE, ROUGHNESS BC AUPOIX SST
+    turb_flatplate_sst_roughBCAupoix           = TestCase('turb_sst_flatplate_roughBCAupoix')
+    turb_flatplate_sst_roughBCAupoix.cfg_dir   = "rans/flatplate/roughness/bc_aupoix"
+    turb_flatplate_sst_roughBCAupoix.cfg_file  = "turb_SST_flatplate_roughBCAupoix.cfg"
+    turb_flatplate_sst_roughBCAupoix.test_iter = 10
+    turb_flatplate_sst_roughBCAupoix.test_vals = [-5.278097, -2.297701, -2.883899, 0.228298, -1.375945, 3.209449, -0.188736, 0.007201]
+    test_list.append(turb_flatplate_sst_roughBCAupoix)
+
     # ONERA M6 Wing
     turb_oneram6           = TestCase('turb_oneram6')
     turb_oneram6.cfg_dir   = "rans/oneram6"

TestCases/rans/flatplate/roughness/bc_aupoix/turb_SST_flatplate_roughBCAupoix.cfg

@@ -0,0 +1,109 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%                                                                                    %
+% SU2 configuration file                                                             %
+% Case description: Turbulent flow over rough flat plate with zero pressure gradient %
+% Author: Thomas D. Economon                                                         %
+% Institution: Stanford University                                                   %
+% Date: 2011.11.10                                                                   %
+% File Version 8.3.0 "Harrier"                                                       %
+%                                                                                    %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------%
+%
+SOLVER= RANS
+KIND_TURB_MODEL= SST
+WALL_ROUGHNESS = (wall, 400e-6)
+KIND_ROUGHSST_MODEL = LIMITER_AUPOIX
+MATH_PROBLEM= DIRECT
+RESTART_SOL= YES
+READ_BINARY_RESTART= YES
+
+% ----------- COMPRESSIBLE AND INCOMPRESSIBLE FREE-STREAM DEFINITION ----------%
+MACH_NUMBER= 0.2
+AOA= 0.0
+SIDESLIP_ANGLE= 0.0
+FREESTREAM_TEMPERATURE= 300.0
+REYNOLDS_NUMBER= 5000000.0
+REYNOLDS_LENGTH= 1.0
+
+% ---------------------- REFERENCE VALUE DEFINITION ---------------------------%
+%
+REF_ORIGIN_MOMENT_X = 0.25
+REF_ORIGIN_MOMENT_Y = 0.00
+REF_ORIGIN_MOMENT_Z = 0.00
+REF_LENGTH= 1.0
+REF_AREA= 2.0
+
+% -------------------- BOUNDARY CONDITION DEFINITION --------------------------%
+%
+MARKER_HEATFLUX= ( wall, 0.0 )
+MARKER_INLET= ( inlet, 302.4, 118309.784, 1.0, 0.0, 0.0 )
+MARKER_OUTLET= ( outlet, 115056.0, farfield, 115056.0 )
+MARKER_SYM= ( symmetry )
+MARKER_PLOTTING= ( wall )
+MARKER_MONITORING= ( wall )
+
+% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------%
+%
+NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES
+CFL_NUMBER= 100.0
+CFL_ADAPT= YES
+CFL_ADAPT_PARAM= ( 0.1, 2.0, 100.0, 1e5 )
+RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 )
+ITER= 99999
+
+% ----------------------- SLOPE LIMITER DEFINITION ----------------------------%
+%
+VENKAT_LIMITER_COEFF= 0.1
+ADJ_SHARP_LIMITER_COEFF= 3.0
+REF_SHARP_EDGES= 3.0
+SENS_REMOVE_SHARP= NO
+
+% -------------------------- MULTIGRID PARAMETERS -----------------------------%
+%
+MGLEVEL= 0
+% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------%
+%
+CONV_NUM_METHOD_FLOW= ROE
+%
+MUSCL_FLOW= YES
+SLOPE_LIMITER_FLOW= NONE
+JST_SENSOR_COEFF= ( 0.5, 0.02 )
+TIME_DISCRE_FLOW= EULER_IMPLICIT
+
+% -------------------- TURBULENT NUMERICAL METHOD DEFINITION ------------------%
+%
+
+CONV_NUM_METHOD_TURB= SCALAR_UPWIND
+MUSCL_TURB= NO
+SLOPE_LIMITER_TURB= VENKATAKRISHNAN
+TIME_DISCRE_TURB= EULER_IMPLICIT
+
+% --------------------------- CONVERGENCE PARAMETERS --------------------------%
+%
+CONV_FIELD= RMS_DENSITY
+CONV_RESIDUAL_MINVAL= -14
+CONV_STARTITER= 10
+CONV_CAUCHY_ELEMS= 100
+CONV_CAUCHY_EPS= 1E-6
+%
+
+% ------------------------- INPUT/OUTPUT INFORMATION --------------------------%
+%
+MESH_FILENAME= mesh_flatplate_turb_137x97
+MESH_FORMAT= SU2
+MESH_OUT_FILENAME= mesh_out
+SOLUTION_FILENAME= restart_flow
+SOLUTION_ADJ_FILENAME= restart_adj
+TABULAR_FORMAT= CSV
+CONV_FILENAME= history
+RESTART_FILENAME= restart_flow
+RESTART_ADJ_FILENAME= restart_adj
+VOLUME_FILENAME= flow
+VOLUME_ADJ_FILENAME= adjoint
+GRAD_OBJFUNC_FILENAME= of_grad.dat
+SURFACE_FILENAME= surface_flow
+SURFACE_ADJ_FILENAME= surface_adjoint
+OUTPUT_WRT_FREQ= 1000
+SCREEN_OUTPUT= (INNER_ITER, WALL_TIME, RMS_RES, LIFT, DRAG)