Boundary conditions for SST for rough walls

Author: marcosch27

Common/include/CConfig.hpp

@@ -749,9 +749,11 @@ 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_OPTIONS *ROUGHSST_Options;         /*!< \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. */
+  unsigned short nROUGHSST_Options;       /*!< \brief Number of rough-wall boundary conditions for SST turbulence model options specified. */ 
   WALL_FUNCTIONS  *Kind_WallFunctions;        /*!< \brief The kind of wall function to use for the corresponding markers. */
   unsigned short  **IntInfo_WallFunctions;    /*!< \brief Additional integer information for the wall function markers. */
   su2double       **DoubleInfo_WallFunctions; /*!< \brief Additional double information for the wall function markers. */
@@ -1200,6 +1202,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 */
@@ -9938,6 +9941,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,59 @@ inline SST_ParsedOptions ParseSSTOptions(const SST_OPTIONS *SST_Options, unsigne
   return SSTParsedOptions;
 }
 
+/*!
+ * \brief SST rough-wall boundary conditions Options
+ */
+enum class ROUGHSST_OPTIONS {
+  NONE,                 /*!< \brief No option / default. */
+  WILCOX1998,           /*!< \brief Wilcox 1998 boundary conditions for rough walls  / default. */
+  WILCOX2006,           /*!< \brief Wilcox 2006 boundary conditions for rough walls. */
+  LIMITER_KNOPP,        /*!< \brief Knopp eddy viscosity limiter. */
+  LIMITER_AUPOIX,       /*!< \brief Aupoix eddy viscosity limiter. */
+};
+static const MapType<std::string, ROUGHSST_OPTIONS> ROUGHSST_Options_Map = {
+  MakePair("NONE", ROUGHSST_OPTIONS::NONE)
+  MakePair("WILCOX1998", ROUGHSST_OPTIONS::WILCOX1998)
+  MakePair("WILCOX2006", ROUGHSST_OPTIONS::WILCOX2006)
+  MakePair("LIMITER_KNOPP", ROUGHSST_OPTIONS::LIMITER_KNOPP)
+  MakePair("LIMITER_AUPOIX", ROUGHSST_OPTIONS::LIMITER_AUPOIX)
+};
+
+/*!
+ * \brief Structure containing parsed SST rough-wall boundary conditions options.
+ */
+struct ROUGH_SST_ParsedOptions {
+  ROUGHSST_OPTIONS version = ROUGHSST_OPTIONS::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 SA options.
+ */
+inline ROUGH_SST_ParsedOptions ParseROUGHSSTOptions(const ROUGHSST_OPTIONS *ROUGHSST_Options, unsigned short nROUGHSST_Options, int rank) {
+  ROUGH_SST_ParsedOptions ROUGHSSTParsedOptions;
+
+  auto IsPresent = [&](ROUGHSST_OPTIONS option) {
+    const auto roughsst_options_end = ROUGHSST_Options + nROUGHSST_Options;
+    return std::find(ROUGHSST_Options, roughsst_options_end, option) != roughsst_options_end;
+  };
+
+  ROUGHSSTParsedOptions.wilcox2006 = IsPresent(ROUGHSST_OPTIONS::WILCOX2006);
+  ROUGHSSTParsedOptions.limiter_knopp = IsPresent(ROUGHSST_OPTIONS::LIMITER_KNOPP);
+  ROUGHSSTParsedOptions.limiter_aupoix = IsPresent(ROUGHSST_OPTIONS::LIMITER_AUPOIX);
+  
+
+  return ROUGHSSTParsedOptions;
+}
+
+
 /*!
  * \brief SA Options
  */

Common/src/CConfig.cpp

@@ -1129,7 +1129,9 @@ void CConfig::SetConfig_Options() {
   addEnumListOption("SST_OPTIONS", nSST_Options, SST_Options, SST_Options_Map);
   /*!\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*/
+  addEnumListOption("ROUGHSST_OPTIONS", nROUGHSST_Options, ROUGHSST_Options, ROUGHSST_Options_Map);
   /*!\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*/
@@ -3555,6 +3557,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(ROUGHSST_Options, nROUGHSST_Options, rank);  
   } 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

@@ -426,6 +426,9 @@ 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)) {
 
+      /*--- distance to closest neighbor ---*/
+      su2double wall_dist = geometry->vertex[val_marker][iVertex]->GetNearestNeighborDistance();
+
       if (rough_wall) {
 
         /*--- Set wall values ---*/
@@ -441,33 +444,61 @@ 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 (roughsstParsedOptions.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 (roughsstParsedOptions.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 = (200/(kPlus+EPS))*(200/(kPlus+EPS));
+          else
+            S_R = 100/(kPlus+EPS) + ((200/(kPlus+EPS))*(200/(kPlus+EPS)) - 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 (roughsstParsedOptions.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 (roughsstParsedOptions.limiter_aupoix) {
+          
+          su2double k0Plus = ( 1.0 /sqrt( constants[6])) * tanh(((log((kPlus +EPS ) / 30.0) / log(10.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);