[WIP] Inconsistencies and improvements to SST model

Common/include/CConfig.hpp

@@ -865,7 +865,7 @@ class CConfig {
   nMu_Temperature_Ref,           /*!< \brief Number of species reference temperature for Sutherland model. */
   nMu_S,                         /*!< \brief Number of species reference S for Sutherland model. */
   nThermal_Conductivity_Constant,/*!< \brief Number of species constant thermal conductivity. */
-  nPrandtl_Lam,                  /*!< \brief Number of species laminar Prandtl number. */
+  nPrandtl_Lam,                  /*!< \brief Prandtl number. */
   nPrandtl_Turb,                 /*!< \brief Number of species turbulent Prandtl number. */
   nConstant_Lewis_Number;       /*!< \brief Number of species Lewis Number. */
   su2double Diffusivity_Constant;   /*!< \brief Constant mass diffusivity for scalar transport.  */
@@ -877,6 +877,7 @@ class CConfig {
   array<su2double, N_POLY_COEFFS> MuPolyCoefficientsND{{0.0}};  /*!< \brief Definition of the non-dimensional temperature polynomial coefficients for viscosity. */
   array<su2double, N_POLY_COEFFS> KtPolyCoefficientsND{{0.0}};  /*!< \brief Definition of the non-dimensional temperature polynomial coefficients for thermal conductivity. */
   su2double TurbIntensityAndViscRatioFreeStream[2]; /*!< \brief Freestream turbulent intensity and viscosity ratio for turbulence and transition models. */
+
   su2double Energy_FreeStream,     /*!< \brief Free-stream total energy of the fluid.  */
   ModVel_FreeStream,               /*!< \brief Magnitude of the free-stream velocity of the fluid.  */
   ModVel_FreeStreamND,             /*!< \brief Non-dimensional magnitude of the free-stream velocity of the fluid.  */
@@ -1176,6 +1177,8 @@ class CConfig {
   nHistoryOutput, nVolumeOutput;  /*!< \brief Number of variables printed to the history file. */
   bool Multizone_Residual;        /*!< \brief Determines if memory should be allocated for the multizone residual. */
   SST_ParsedOptions sstParsedOptions; /*!< \brief Additional parameters for the SST turbulence model. */
+  su2double prodLimConst;
+  su2double LDomain;
   SA_ParsedOptions saParsedOptions;   /*!< \brief Additional parameters for the SA turbulence model. */
   LM_ParsedOptions lmParsedOptions;   /*!< \brief Additional parameters for the LM transition model. */
   su2double uq_delta_b;         /*!< \brief Parameter used to perturb eigenvalues of Reynolds Stress Matrix */
@@ -9908,6 +9911,9 @@ class CConfig {
    */
   SST_ParsedOptions GetSSTParsedOptions() const { return sstParsedOptions; }
 
+  su2double GetProdLimConst() const { return prodLimConst; }
+  su2double GetLDomain() const { return LDomain; }
+
   /*!
    * \brief Get parsed SA option data structure.
    * \return SA option data structure.

Common/include/option_structure.hpp

@@ -994,21 +994,26 @@ enum class SST_OPTIONS {
   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. */
   DLL,         /*!< \brief Menter k-w SST model with dimensionless lower limit clipping of turbulence variables. */
+  FULLPROD,    /*!< \brief Menter k-w SST model with full production term. */
+  PRODLIM,     /*!< \brief Menter k-w SST model with user-defined production limiter constant. */
+  NEWBC,       /*!< \brief Menter k-w SST model with new boundary conditions. */
 };
 static const MapType<std::string, SST_OPTIONS> SST_Options_Map = {
   MakePair("NONE", SST_OPTIONS::NONE)
   MakePair("V1994m", SST_OPTIONS::V1994m)
   MakePair("V2003m", SST_OPTIONS::V2003m)
-  /// TODO: For now we do not support "unmodified" versions of SST.
-  //MakePair("V1994", SST_OPTIONS::V1994)
-  //MakePair("V2003", SST_OPTIONS::V2003)
+  MakePair("V1994", SST_OPTIONS::V1994)
+  MakePair("V2003", SST_OPTIONS::V2003)
   MakePair("SUSTAINING", SST_OPTIONS::SUST)
   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)
   MakePair("DIMENSIONLESS_LIMIT", SST_OPTIONS::DLL)
+  MakePair("FULLPROD", SST_OPTIONS::FULLPROD)
+  MakePair("PRODLIM", SST_OPTIONS::PRODLIM)
+  MakePair("NEWBC", SST_OPTIONS::NEWBC)
 };
 
 /*!
@@ -1023,6 +1028,9 @@ struct SST_ParsedOptions {
   bool compWilcox = false;                    /*!< \brief Bool for compressibility correction of Wilcox. */
   bool compSarkar = false;                    /*!< \brief Bool for compressibility correction of Sarkar. */
   bool dll = false;                           /*!< \brief Bool dimensionless lower limit. */
+  bool fullProd = false;                      /*!< \brief Bool for full production term. */
+  bool prodLim = false;                       /*!< \brief Bool for user-defined production limiter constant. */
+  bool newBC = false;                         /*!< \brief Bool for new boundary conditions. */
 };
 
 /*!
@@ -1040,6 +1048,10 @@ inline SST_ParsedOptions ParseSSTOptions(const SST_OPTIONS *SST_Options, unsigne
     return std::find(SST_Options, sst_options_end, option) != sst_options_end;
   };
 
+  const bool found_fullProd = IsPresent(SST_OPTIONS::FULLPROD);
+  const bool found_prodLim = IsPresent(SST_OPTIONS::PRODLIM);
+  const bool found_newBC = IsPresent(SST_OPTIONS::NEWBC);
+
   const bool found_1994 = IsPresent(SST_OPTIONS::V1994);
   const bool found_2003 = IsPresent(SST_OPTIONS::V2003);
   const bool found_1994m = IsPresent(SST_OPTIONS::V1994m);
@@ -1097,6 +1109,9 @@ inline SST_ParsedOptions ParseSSTOptions(const SST_OPTIONS *SST_Options, unsigne
   SSTParsedOptions.compSarkar = sst_compSarkar;
   SSTParsedOptions.dll = sst_dll;
 
+  SSTParsedOptions.fullProd = found_fullProd;
+  SSTParsedOptions.prodLim = found_prodLim;
+  SSTParsedOptions.newBC = found_newBC;
   return SSTParsedOptions;
 }
 

Common/src/CConfig.cpp

@@ -1119,6 +1119,9 @@ 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);
 
+  addDoubleOption("PROD_LIM_CONST", prodLimConst, 20.0);
+  addDoubleOption("L_DOMAIN", LDomain, 1.0);
+
   /*!\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*/
@@ -6275,6 +6278,9 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) {
             else cout << "\nusing default hard coded lower limit clipping";
 
             cout << "." << endl;
+
+            if (sstParsedOptions.prodLim) cout << "Changing the value of the TKE production limiter constant to " << prodLimConst << endl;
+
             break;
         }
         switch (Kind_Trans_Model) {

SU2_CFD/include/numerics/CNumerics.hpp

@@ -190,6 +190,8 @@ class CNumerics {
 
   bool bounded_scalar = false;    /*!< \brief Flag for bounded scalar problem */
 
+  su2double ProdDistr[5];
+
 public:
   /*!
    * \brief Return type used in some "ComputeResidual" overloads to give a
@@ -1605,6 +1607,8 @@ class CNumerics {
    * \return is_bounded_scalar : scalar solver uses bounded scalar convective transport
    */
   inline bool GetBoundedScalar() const { return bounded_scalar;}
+
+  inline su2double GetProdDest(int index) const {return ProdDistr[index];}
 };
 
 /*!

SU2_CFD/include/numerics/turbulent/turb_sources.hpp

@@ -656,6 +656,8 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
   /*--- Ambient values for SST-SUST. ---*/
   const su2double kAmb, omegaAmb;
 
+  su2double Pk, Dk, Pw, Dw;
+
   su2double F1_i, F2_i, CDkw_i;
   su2double Residual[2];
   su2double* Jacobian_i[2];
@@ -869,15 +871,20 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
       const su2double prod_limit = prod_lim_const * beta_star * Density_i * ScalarVar_i[1] * ScalarVar_i[0];
 
       su2double P = Eddy_Viscosity_i * pow(P_Base, 2);
+      if (sstParsedOptions.fullProd) P -= Eddy_Viscosity_i * diverg*diverg * 2.0/3.0;
+      if (!sstParsedOptions.modified) P -= Density_i * ScalarVar_i[0] * diverg * 2.0/3.0;
+
+      su2double PLim = 0.0;
+      if ( P > prod_limit ) PLim = 1.0;
       su2double pk = max(0.0, min(P, prod_limit));
 
-      const auto& eddy_visc_var = sstParsedOptions.version == SST_OPTIONS::V1994 ? VorticityMag : StrainMag_i;
+      const su2double eddy_visc_var = sstParsedOptions.version == SST_OPTIONS::V1994 ? VorticityMag : StrainMag_i;
       const su2double zeta = max(ScalarVar_i[1], eddy_visc_var * F2_i / a1);
 
       /*--- Production limiter only for V2003, recompute for V1994. ---*/
       su2double pw;
       if (sstParsedOptions.version == SST_OPTIONS::V1994) {
-        pw = alfa_blended * Density_i * pow(P_Base, 2);
+        pw = alfa_blended * Density_i * P / Eddy_Viscosity_i;
       } else {
         pw = (alfa_blended * Density_i / Eddy_Viscosity_i) * pk;
       }
@@ -912,17 +919,26 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
       }
 
       /*--- Add the production terms to the residuals. ---*/
+      Pk = pk;
+      Pw = pw;
 
       Residual[0] += pk * Volume;
       Residual[1] += pw * Volume;
 
       /*--- Add the dissipation  terms to the residuals.---*/
+      Dk = dk;
+      Dw = dw;
 
       Residual[0] -= dk * Volume;
       Residual[1] -= dw * Volume;
 
-      /*--- Cross diffusion ---*/
+      ProdDistr[0] = Pk;
+      ProdDistr[1] = Dk;
+      ProdDistr[2] = Pw;
+      ProdDistr[3] = Dw;
+      ProdDistr[4] = PLim;
 
+      /*--- Cross diffusion ---*/
       Residual[1] += (1.0 - F1_i) * CDkw_i * Volume;
 
       /*--- Contribution due to 2D axisymmetric formulation ---*/
@@ -932,6 +948,7 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
       /*--- Implicit part ---*/
 
       Jacobian_i[0][0] = -beta_star * ScalarVar_i[1] * Volume * (1.0 + zetaFMt);
+      if (!sstParsedOptions.modified) Jacobian_i[0][0] -= diverg * Volume*2.0/3.0;
       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 * (1.0 - 0.09/beta_blended * zetaFMt);

SU2_CFD/include/numerics_simd/flow/convection/centered.hpp

@@ -53,6 +53,8 @@ class CCenteredBase : public Base {
   const bool dynamicGrid;
   const su2double stretchParam = 0.3;
 
+  using Base::turbVars;
+
   /*!
    * \brief Constructor, store some constants and forward args to base.
    */
@@ -96,6 +98,8 @@ class CCenteredBase : public Base {
     const bool implicit = (config.GetKind_TimeIntScheme() == EULER_IMPLICIT);
     const auto& solution = static_cast<const CEulerVariable&>(solution_);
 
+    const bool tkeNeeded = config.GetKind_Turb_Model() == TURB_MODEL::SST;
+
     const auto iPoint = geometry.edges->GetNode(iEdge,0);
     const auto jPoint = geometry.edges->GetNode(iEdge,1);
 
@@ -119,6 +123,13 @@ class CCenteredBase : public Base {
       avgV.all(iVar) = 0.5 * (V.i.all(iVar) + V.j.all(iVar));
     }
 
+    if (tkeNeeded) {
+      V.i.allTurb = gatherVariables<1>(iPoint, turbVars->GetSolution());
+      V.j.allTurb = gatherVariables<1>(jPoint, turbVars->GetSolution());
+
+      avgV.allTurb(0) = 0.5*(V.i.allTurb(0)+V.j.allTurb(0));
+    }
+
     /*--- Compute conservative variables. ---*/
 
     CPair<CCompressibleConservatives<nDim> > U;

SU2_CFD/include/numerics_simd/flow/convection/common.hpp

@@ -113,7 +113,8 @@ FORCEINLINE CPair<ReconVarType> reconstructPrimitives(Int iEdge, Int iPoint, Int
                                                       bool muscl, LIMITER limiterType,
                                                       const CPair<PrimVarType>& V1st,
                                                       const VectorDbl<nDim>& vector_ij,
-                                                      const VariableType& solution) {
+                                                      const VariableType& solution,
+                                                      const bool tkeNeeded) {
   static_assert(ReconVarType::nVar <= PrimVarType::nVar,"");
 
   const auto& gradients = solution.GetGradient_Reconstruction();
@@ -125,6 +126,11 @@ FORCEINLINE CPair<ReconVarType> reconstructPrimitives(Int iEdge, Int iPoint, Int
     V.i.all(iVar) = V1st.i.all(iVar);
     V.j.all(iVar) = V1st.j.all(iVar);
   }
+  // Only first order for turbulence
+  if (tkeNeeded) {
+    V.i.allTurb(0) = V1st.i.allTurb(0);
+    V.j.allTurb(0) = V1st.j.allTurb(0);
+  }
 
   if (muscl) {
     /*--- Recompute density and enthalpy instead of reconstructing. ---*/

SU2_CFD/include/numerics_simd/flow/convection/roe.hpp

@@ -63,6 +63,8 @@ class CRoeBase : public Base {
   const bool muscl;
   const LIMITER typeLimiter;
 
+  using Base::turbVars;
+
   /*!
    * \brief Constructor, store some constants and forward args to base.
    */
@@ -98,6 +100,8 @@ class CRoeBase : public Base {
     const bool implicit = (config.GetKind_TimeIntScheme() == EULER_IMPLICIT);
     const auto& solution = static_cast<const CEulerVariable&>(solution_);
 
+    const bool tkeNeeded = config.GetKind_Turb_Model() == TURB_MODEL::SST;
+
     const auto iPoint = geometry.edges->GetNode(iEdge,0);
     const auto jPoint = geometry.edges->GetNode(iEdge,1);
 
@@ -118,8 +122,13 @@ class CRoeBase : public Base {
     V1st.i.all = gatherVariables<nPrimVar>(iPoint, solution.GetPrimitive());
     V1st.j.all = gatherVariables<nPrimVar>(jPoint, solution.GetPrimitive());
 
+    if (tkeNeeded) {
+      V1st.i.allTurb = gatherVariables(iPoint, turbVars->GetSolution());
+      V1st.j.allTurb = gatherVariables(jPoint, turbVars->GetSolution());
+    }
+
     auto V = reconstructPrimitives<CCompressiblePrimitives<nDim,nPrimVarGrad> >(
-        iEdge, iPoint, jPoint, gamma, gasConst, muscl, typeLimiter, V1st, vector_ij, solution);
+        iEdge, iPoint, jPoint, gamma, gasConst, muscl, typeLimiter, V1st, vector_ij, solution, tkeNeeded);
 
     /*--- Compute conservative variables. ---*/
 
@@ -231,6 +240,8 @@ class CRoeScheme : public CRoeBase<CRoeScheme<Decorator>,Decorator> {
   using Base::kappa;
   const ENUM_ROELOWDISS typeDissip;
 
+  using Base::turbVars;
+
 public:
   /*!
    * \brief Constructor, store some constants and forward to base.

SU2_CFD/include/numerics_simd/flow/diffusion/common.hpp

@@ -71,13 +71,13 @@ FORCEINLINE void correctGradient(const PrimitiveType& V,
  */
 template<size_t nVar, size_t nDim>
 FORCEINLINE MatrixDbl<nDim> stressTensor(Double viscosity,
-                                         const MatrixDbl<nVar,nDim>& grad) {
+                                         const MatrixDbl<nVar,nDim>& grad, Double density=0.0, Double tke=0.0) {
   /*--- Hydrostatic term. ---*/
   Double velDiv = 0.0;
   for (size_t iDim = 0; iDim < nDim; ++iDim) {
     velDiv += grad(iDim+1,iDim);
   }
-  Double pTerm = 2.0/3.0 * viscosity * velDiv;
+  Double pTerm = 2.0/3.0 * (viscosity * velDiv + density * tke);
 
   MatrixDbl<nDim> tau;
   for (size_t iDim = 0; iDim < nDim; ++iDim) {

SU2_CFD/include/numerics_simd/flow/diffusion/viscous_fluxes.hpp

@@ -50,6 +50,7 @@ class CNoViscousFlux : public CNumericsSIMD {
 protected:
   static constexpr size_t nDim = NDIM;
   static constexpr size_t nPrimVar = 0;
+  const CVariable* turbVars = nullptr;
 
   template<class... Ts>
   CNoViscousFlux(Ts&...) {}
@@ -140,6 +141,8 @@ class CCompressibleViscousFluxBase : public CNumericsSIMD {
     const auto& solution = static_cast<const CNSVariable&>(solution_);
     const auto& gradient = solution.GetGradient_Primitive();
 
+    const bool tkeNeeded = (config.GetKind_Turb_Model() == TURB_MODEL::SST) && !(config.GetSSTParsedOptions().modified);
+
     /*--- Compute distance and handle zero without "ifs" by making it large. ---*/
 
     auto dist2_ij = squaredNorm(vector_ij);
@@ -152,12 +155,15 @@ class CCompressibleViscousFluxBase : public CNumericsSIMD {
     if(correct) correctGradient(V, vector_ij, dist2_ij, avgGrad);
 
     /*--- Stress and heat flux tensors. ---*/
-
-    auto tau = stressTensor(avgV.laminarVisc() + (uq? Double(0.0) : avgV.eddyVisc()), avgGrad);
+    Double turb_ke = 0.0;
+    Double density = avgV.density();
+    if (tkeNeeded)  turb_ke = 0.5*(gatherVariables(iPoint, turbVars->GetSolution()) +
+                                   gatherVariables(jPoint, turbVars->GetSolution()));
+    auto tau = stressTensor(avgV.laminarVisc() + (uq? Double(0.0) : avgV.eddyVisc()), avgGrad, density, turb_ke);
     if(useSA_QCR) addQCR(avgGrad, tau);
     if(uq) {
-      Double turb_ke = 0.5*(gatherVariables(iPoint, turbVars->GetSolution()) +
-                            gatherVariables(jPoint, turbVars->GetSolution()));
+      turb_ke = 0.5*(gatherVariables(iPoint, turbVars->GetSolution()) +
+                     gatherVariables(jPoint, turbVars->GetSolution()));
       addPerturbedRSM(avgV, avgGrad, turb_ke, tau,
                       uq_eigval_comp, uq_permute, uq_delta_b, uq_urlx);
     }
@@ -263,6 +269,7 @@ class CCompressibleViscousFlux : public CCompressibleViscousFluxBase<NDIM, CComp
   using Base::prandtlLam;
   using Base::prandtlTurb;
   using Base::cp;
+  using Base::turbVars;
 
   /*!
    * \brief Constructor, initialize constants and booleans.
@@ -319,6 +326,8 @@ class CGeneralCompressibleViscousFlux : public CCompressibleViscousFluxBase<NDIM
   static constexpr size_t nSecVar = 4;
   using Base = CCompressibleViscousFluxBase<NDIM, CGeneralCompressibleViscousFlux<NDIM> >;
   using Base::prandtlTurb;
+  using Base::turbVars;
+;
 
   /*!
    * \brief Constructor, initialize constants and booleans.