Add SST Limiter Output

Author: marcosch27

SU2_CFD/include/variables/CTurbSSTVariable.hpp

@@ -41,9 +41,11 @@ class CTurbSSTVariable final : public CTurbVariable {
   su2double sigma_om2;
   su2double beta_star;
   su2double prod_lim_const;
+  su2double a1;
   VectorType F1;
   VectorType F2;    /*!< \brief Menter blending function for blending of k-w and k-eps. */
   VectorType CDkw;  /*!< \brief Cross-diffusion. */
+  VectorType SST_Limiter; 
   SST_ParsedOptions sstParsedOptions;
 public:
   /*!
@@ -71,7 +73,7 @@ class CTurbSSTVariable final : public CTurbVariable {
    * \param[in] val_dist - Value of the distance to the wall.
    * \param[in] val_density - Value of the density.
    */
-  void SetBlendingFunc(unsigned long iPoint, su2double val_viscosity, su2double val_dist, su2double val_density, TURB_TRANS_MODEL trans_model) override;
+  void SetBlendingFunc(unsigned long iPoint, su2double val_viscosity, su2double val_dist, su2double val_density, su2double vorticity_mag, TURB_TRANS_MODEL trans_model) override;
 
   /*!
    * \brief Get the first blending function.
@@ -87,4 +89,9 @@ class CTurbSSTVariable final : public CTurbVariable {
    * \brief Get the value of the cross diffusion of tke and omega.
    */
   inline su2double GetCrossDiff(unsigned long iPoint) const override { return CDkw(iPoint); }
+  /*!
+   * \brief Get the value of the SST Limiter.
+   */
+  inline su2double GetSSTLimiter(unsigned long iPoint) const { return SST_Limiter[iPoint]; }
+
 };

SU2_CFD/include/variables/CVariable.hpp

@@ -1656,7 +1656,7 @@ class CVariable {
    * \param[in] val_density - Value of the density.
    * \param[in] val_dist - Value of the distance to the wall.
    */
-  inline virtual void SetBlendingFunc(unsigned long iPoint, su2double val_viscosity, su2double val_dist, su2double val_density, TURB_TRANS_MODEL trans_model) {}
+  inline virtual void SetBlendingFunc(unsigned long iPoint, su2double val_viscosity, su2double val_dist, su2double val_density, su2double vorticity_mag, TURB_TRANS_MODEL trans_model) {}
 
   /*!
    * \brief Get the first blending function of the SST model.
@@ -1673,6 +1673,17 @@ class CVariable {
    */
   inline virtual su2double GetCrossDiff(unsigned long iPoint) const { return 0.0; }
 
+  /*!
+   * \brief Get the value of the SST Limiter.
+   */
+  inline virtual su2double GetSSTLimiter(unsigned long iPoint) const { return 0.0; }
+
+  
+  /*!
+   * \brief Set the value of the SST Limiter.
+   */
+  inline virtual su2double SetSSTLimiter(unsigned long iPoint) const { return 0.0; }
+
   /*!
    * \brief Get the value of the eddy viscosity.
    * \return the value of the eddy viscosity.

SU2_CFD/src/output/CFlowOutput.cpp

@@ -1422,6 +1422,10 @@ void CFlowOutput::SetVolumeOutputFieldsScalarPrimitive(const CConfig* config) {
 
   if (config->GetKind_Turb_Model() != TURB_MODEL::NONE) {
     AddVolumeOutput("EDDY_VISCOSITY", "Eddy_Viscosity", "PRIMITIVE", "Turbulent eddy viscosity");
+    /*--- SST Limiter ---*/
+    if (config->GetKind_Turb_Model() == TURB_MODEL::SST) {
+      AddVolumeOutput("SST_LIMITER", "SST_Limiter", "PRIMITIVE", "SST viscosity limiter");
+    }
   }
 
 }
@@ -1559,6 +1563,11 @@ void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* con
     SetVolumeOutputValue("EDDY_VISCOSITY", iPoint, Node_Flow->GetEddyViscosity(iPoint));
     SetVolumeOutputValue("TURB_DELTA_TIME", iPoint, Node_Turb->GetDelta_Time(iPoint));
     SetVolumeOutputValue("TURB_CFL", iPoint, Node_Turb->GetLocalCFL(iPoint));
+    /*--- SST Limiter ---*/
+    if (config->GetKind_Turb_Model() == TURB_MODEL::SST) {
+        SetVolumeOutputValue("SST_LIMITER", iPoint, Node_Turb->GetSSTLimiter(iPoint));
+        su2double sst_lim = Node_Turb->GetSSTLimiter(iPoint);
+  }
   }
 
   if (config->GetSAParsedOptions().bc) {

SU2_CFD/src/solvers/CTurbSSTSolver.cpp

@@ -233,7 +233,7 @@ void CTurbSSTSolver::Postprocessing(CGeometry *geometry, CSolver **solver_contai
 
     const su2double VorticityMag = max(GeometryToolbox::Norm(3, flowNodes->GetVorticity(iPoint)), 1e-12);
     const su2double StrainMag = max(flowNodes->GetStrainMag(iPoint), 1e-12);
-    nodes->SetBlendingFunc(iPoint, mu, dist, rho, config->GetKind_Trans_Model());
+    nodes->SetBlendingFunc(iPoint, mu, dist, rho, VorticityMag, config->GetKind_Trans_Model());
 
     const su2double F2 = nodes->GetF2blending(iPoint);
 
@@ -425,12 +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();
 
       /*--- 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);
@@ -445,11 +445,10 @@ void CTurbSSTSolver::BC_HeatFlux_Wall(CGeometry *geometry, CSolver **solver_cont
 
         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 (roughsstParsedOptions.wilcox1998) {
+        if (options.wilcox1998) {
           if (kPlus <= 25)
             S_R = (50/(kPlus+EPS))*(50/(kPlus+EPS));
           else
@@ -458,7 +457,7 @@ void CTurbSSTSolver::BC_HeatFlux_Wall(CGeometry *geometry, CSolver **solver_cont
           solution[0] = 0.0;
           solution[1] = FrictionVel*FrictionVel*S_R/(laminar_viscosity/density);          
         } 
-        else if (roughsstParsedOptions.wilcox2006) {
+        else if (options.wilcox2006) {
         /*--- Reference 2 from D.C. Wilcox Turbulence Modeling for CFD (2006) ---*/
           if (kPlus <= 5)
             S_R = pow(200/(kPlus+EPS),2);
@@ -469,8 +468,7 @@ void CTurbSSTSolver::BC_HeatFlux_Wall(CGeometry *geometry, CSolver **solver_cont
           solution[1] = FrictionVel*FrictionVel*S_R/(laminar_viscosity/density);
         } 
         /*--- Knopp eddy viscosity limiter ---*/
-        else if (roughsstParsedOptions.limiter_knopp) {
-            
+        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);
 
@@ -479,8 +477,7 @@ void CTurbSSTSolver::BC_HeatFlux_Wall(CGeometry *geometry, CSolver **solver_cont
           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) {
-          
+        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;
@@ -489,16 +486,13 @@ void CTurbSSTSolver::BC_HeatFlux_Wall(CGeometry *geometry, CSolver **solver_cont
 
           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);

SU2_CFD/src/variables/CTurbSSTVariable.cpp

@@ -45,16 +45,18 @@ CTurbSSTVariable::CTurbSSTVariable(su2double kine, su2double omega, su2double mu
   sigma_om2 = constants[3];
   beta_star = constants[6];
   prod_lim_const = constants[10];
+  a1 = constants[7];
 
   F1.resize(nPoint) = su2double(1.0);
   F2.resize(nPoint) = su2double(0.0);
   CDkw.resize(nPoint) = su2double(0.0);
+  SST_Limiter.resize(nPoint, 1.0);
 
   muT.resize(nPoint) = mut;
 }
 
 void CTurbSSTVariable::SetBlendingFunc(unsigned long iPoint, su2double val_viscosity,
-                                       su2double val_dist, su2double val_density, TURB_TRANS_MODEL trans_model) {
+                                       su2double val_dist, su2double val_density, su2double vorticity_mag, TURB_TRANS_MODEL trans_model) {
   su2double arg2, arg2A, arg2B, arg1;
 
   AD::StartPreacc();
@@ -91,7 +93,13 @@ void CTurbSSTVariable::SetBlendingFunc(unsigned long iPoint, su2double val_visco
     F1(iPoint) = max(F1(iPoint), F3);
   }
 
-  AD::SetPreaccOut(F1(iPoint)); AD::SetPreaccOut(F2(iPoint)); AD::SetPreaccOut(CDkw(iPoint));
+  /*--- SST Limiter ---*/ 
+  
+  su2double limiter = a1*Solution(iPoint,1) / (F2(iPoint)*vorticity_mag + EPS);
+  limiter = min(1.0, limiter);
+  SST_Limiter(iPoint) = limiter;
+
+  AD::SetPreaccOut(F1(iPoint)); AD::SetPreaccOut(F2(iPoint)); AD::SetPreaccOut(CDkw(iPoint)); AD::SetPreaccOut(SST_Limiter(iPoint));
   AD::EndPreacc();
 
 }