Merge pull request #2670 from su2code/pedro/smooth_the_roughness

smooth the SST roughness implementation

Author: pcarruscag

Common/include/CConfig.hpp

@@ -1227,7 +1227,6 @@ 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 */
@@ -10176,10 +10175,9 @@ 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.
+   * \brief Get rough-wall boundary conditions for SST.
    */
-  ROUGH_SST_ParsedOptions GetROUGHSSTParsedOptions() const { return roughsstParsedOptions; }
+  ROUGHSST_MODEL GetKindRoughSSTModel() const { return Kind_RoughSST_Model; }
 
   /*!
    * \brief Get parsed option data structure for data-driven fluid model.

Common/include/option_structure.hpp

@@ -1105,54 +1105,18 @@ inline SST_ParsedOptions ParseSSTOptions(const SST_OPTIONS *SST_Options, unsigne
  * \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

@@ -1141,7 +1141,7 @@ void CConfig::SetConfig_Options() {
   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);
+  addEnumOption("KIND_ROUGHSST_MODEL", Kind_RoughSST_Model, RoughSST_Model_Map, ROUGHSST_MODEL::WILCOX1998);
   /*!\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*/
@@ -3600,7 +3600,6 @@ 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,7 +39,6 @@ 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.
@@ -54,13 +53,13 @@ class CTurbSSTSolver final : public CTurbSolver {
                      CSolver **solver_container,
                      const CConfig *config,
                      unsigned short val_marker);
-                     
+
   /*!
    * \brief Compute a suitable under-relaxation parameter to limit the change in the solution variables over
    * a nonlinear iteration for stability.
    * \param[in] config - Definition of the particular problem.
    */
-  void ComputeUnderRelaxationFactor(const CConfig *config);
+  void ComputeUnderRelaxationFactor(const CConfig *config) override;
 
 public:
   /*!

SU2_CFD/src/solvers/CTurbSSTSolver.cpp

@@ -405,11 +405,10 @@ void CTurbSSTSolver::BC_HeatFlux_Wall(CGeometry *geometry, CSolver **solver_cont
 
   const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT);
 
-  bool rough_wall = false;
   string Marker_Tag = config->GetMarker_All_TagBound(val_marker);
   WALL_TYPE WallType; su2double Roughness_Height;
   tie(WallType, Roughness_Height) = config->GetWallRoughnessProperties(Marker_Tag);
-  if (WallType == WALL_TYPE::ROUGH) rough_wall = true;
+  const bool rough_wall = WallType == WALL_TYPE::ROUGH && Roughness_Height > 0;
 
   /*--- Evaluate nu tilde at the closest point to the surface using the wall functions. ---*/
 
@@ -424,94 +423,90 @@ void CTurbSSTSolver::BC_HeatFlux_Wall(CGeometry *geometry, CSolver **solver_cont
     const auto iPoint = geometry->vertex[val_marker][iVertex]->GetNode();
 
     /*--- 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 (!geometry->nodes->GetDomain(iPoint)) continue;
+
+    /*--- distance to closest neighbor ---*/
+    su2double wall_dist = geometry->vertex[val_marker][iVertex]->GetNearestNeighborDistance();
 
-      /*--- distance to closest neighbor ---*/
-      su2double wall_dist = geometry->vertex[val_marker][iVertex]->GetNearestNeighborDistance();
+    su2double solution[MAXNVAR];
 
-      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);
-        su2double WallShearStress = solver_container[FLOW_SOL]->GetWallShearStress(val_marker, iVertex);
+    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);
+      su2double WallShearStress = solver_container[FLOW_SOL]->GetWallShearStress(val_marker, iVertex);
 
-        /*--- Compute non-dimensional velocity ---*/
-        su2double FrictionVel = sqrt(fabs(WallShearStress)/density);
+      /*--- Compute non-dimensional velocity ---*/
+      su2double FrictionVel = sqrt(fabs(WallShearStress)/density);
 
-        /*--- Compute roughness in wall units. ---*/
-        //su2double Roughness_Height = config->GetWall_RoughnessHeight(Marker_Tag);
-        su2double kPlus = FrictionVel*Roughness_Height*density/laminar_viscosity;
+      /*--- Compute roughness in wall units. ---*/
+      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. ---*/
+      /*--- Modify the omega and k to account for a rough wall. ---*/
 
-        /*--- Reference 1 original Wilcox (1998) ---*/
-        if (options.wilcox1998) {
+      switch (config->GetKindRoughSSTModel()) {
+        /*--- Reference 1 original Wilcox (1998). ---*/
+        case ROUGHSST_MODEL::WILCOX1998: {
+          su2double S_R = 0.0;
           if (kPlus <= 25)
-            S_R = (50/(kPlus+EPS))*(50/(kPlus+EPS));
+            S_R = pow(50/(kPlus+EPS), 2);
           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) {
+          solution[1] = FrictionVel*FrictionVel*S_R/(laminar_viscosity/density);
+        } break;
         /*--- Reference 2 from D.C. Wilcox Turbulence Modeling for CFD (2006) ---*/
+        case ROUGHSST_MODEL::WILCOX2006: {
+          su2double S_R = 0.0;
           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);
-        } 
+        } break;
         /*--- Knopp eddy viscosity limiter ---*/
-        else if (options.limiter_knopp) {
+        case ROUGHSST_MODEL::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) {
+          solution[1] = min( FrictionVel/(sqrt(constants[6])*d0*kappa), 60.0*laminar_viscosity/(density*beta_1*pow(wall_dist,2)));
+        } break;
+        /*--- Aupoix eddy viscosity limiter ---*/
+        case (ROUGHSST_MODEL::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
-        /*--- Set wall values ---*/
-        su2double density = solver_container[FLOW_SOL]->GetNodes()->GetDensity(iPoint);
-        su2double laminar_viscosity = solver_container[FLOW_SOL]->GetNodes()->GetLaminarViscosity(iPoint);
-
-        su2double beta_1 = constants[4];
-        su2double solution[MAXNVAR];
-        solution[0] = 0.0;
-        solution[1] = 60.0*laminar_viscosity/(density*beta_1*pow(wall_dist,2));
-
-        /*--- Set the solution values and zero the residual ---*/
-        nodes->SetSolution_Old(iPoint,solution);
-        nodes->SetSolution(iPoint,solution);
-        LinSysRes.SetBlock_Zero(iPoint);
+        } break;
       }
+    } else { // smooth wall
+      /*--- Set wall values ---*/
+      su2double density = solver_container[FLOW_SOL]->GetNodes()->GetDensity(iPoint);
+      su2double laminar_viscosity = solver_container[FLOW_SOL]->GetNodes()->GetLaminarViscosity(iPoint);
+
+      su2double beta_1 = constants[4];
+      solution[0] = 0.0;
+      solution[1] = 60.0*laminar_viscosity/(density*beta_1*pow(wall_dist,2));
+    }
 
-      if (implicit) {
-        /*--- Change rows of the Jacobian (includes 1 in the diagonal) ---*/
-        Jacobian.DeleteValsRowi(iPoint*nVar);
-        Jacobian.DeleteValsRowi(iPoint*nVar+1);
-      }
+    /*--- Set the solution values and zero the residual ---*/
+    nodes->SetSolution_Old(iPoint, solution);
+    nodes->SetSolution(iPoint, solution);
+    LinSysRes.SetBlock_Zero(iPoint);
+
+    if (implicit) {
+      /*--- Change rows of the Jacobian (includes 1 in the diagonal) ---*/
+      Jacobian.DeleteValsRowi(iPoint*nVar);
+      Jacobian.DeleteValsRowi(iPoint*nVar+1);
     }
   }
   END_SU2_OMP_FOR

config_template.cfg

@@ -1362,6 +1362,9 @@ MARKER_DISPLACEMENT= ( NONE )
 % This is a list of (string, double) each element corresponding to the MARKER defined in WALL_TYPE.
 WALL_ROUGHNESS = (wall1, ks1, wall2, ks2)
 %WALL_ROUGHNESS = (wall1, ks1, wall2, 0.0) %is also allowed
+%
+% Roughness model used for SST (WILCOX1998, WILCOX2006, LIMITER_KNOPP, LIMITER_AUPOIX).
+KIND_ROUGHSST_MODEL = WILCOX1998
 
 % ------------------------ WALL FUNCTION DEFINITION --------------------------%
 %