Align the SST model with the NASA TMR definitions and fix the boundary energy (#1851)

Common/include/CConfig.hpp

@@ -1279,6 +1279,9 @@ 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 LDomain;           /*!< \brief Approximate length of the domain, for the far-field omega of TMRBC. */
+  su2double sstSustTkeAmb;     /*!< \brief Ambient k of the SST sustaining terms (dimensional), <= 0 for the default. */
+  su2double sstSustOmegaAmb;   /*!< \brief Ambient omega of the SST sustaining terms (dimensional), <= 0 for the default. */
   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 */
@@ -10360,6 +10363,27 @@ class CConfig {
    */
   SST_ParsedOptions GetSSTParsedOptions() const { return sstParsedOptions; }
 
+  su2double GetLDomain() const { return LDomain; }
+
+  /*!
+   * \brief Ambient (free-stream) k of the SST sustaining terms, dimensional.
+   * \note Default of Spalart and Rumsey (AIAA J 45(10), 2007), as in the NASA TMR SST-sust: 1e-6 U^2.
+   * \param[in] velMag - Free-stream velocity magnitude (dimensional).
+   */
+  su2double GetSSTSust_TkeAmb(su2double velMag) const {
+    return sstSustTkeAmb > 0.0 ? sstSustTkeAmb : 1e-6 * velMag * velMag;
+  }
+
+  /*!
+   * \brief Ambient (free-stream) omega of the SST sustaining terms, dimensional.
+   * \note Default of Spalart and Rumsey (AIAA J 45(10), 2007), as in the NASA TMR SST-sust: 5 U / L, with L the
+   *       defining length of the problem, taken as REYNOLDS_LENGTH.
+   * \param[in] velMag - Free-stream velocity magnitude (dimensional).
+   */
+  su2double GetSSTSust_OmegaAmb(su2double velMag) const {
+    return sstSustOmegaAmb > 0.0 ? sstSustOmegaAmb : 5.0 * velMag / Length_Reynolds;
+  }
+
   /*!
    * \brief Get parsed SA option data structure.
    * \return SA option data structure.

Common/include/option_structure.hpp

@@ -1101,33 +1101,36 @@ inline TURB_FAMILY TurbModelFamily(TURB_MODEL model) {
  * \brief SST Options
  */
 enum class SST_OPTIONS {
-  NONE,        /*!< \brief No SST Turb model. */
-  V1994,       /*!< \brief 1994 Menter k-w SST model. */
-  V2003,       /*!< \brief 2003 Menter k-w SST model. */
-  V1994m,      /*!< \brief 1994m Menter k-w SST model. */
-  V2003m,      /*!< \brief 2003m Menter k-w SST model. */
-  SUST,        /*!< \brief Menter k-w SST model with sustaining terms. */
-  V,           /*!< \brief Menter k-w SST model with vorticity production terms. */
-  KL,          /*!< \brief Menter k-w SST model with Kato-Launder production terms. */
-  UQ,          /*!< \brief Menter k-w SST model with uncertainty quantification modifications. */
-  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. */
+  NONE,          /*!< \brief No SST Turb model. */
+  V1994,         /*!< \brief 1994 Menter k-w SST model. */
+  V2003,         /*!< \brief 2003 Menter k-w SST model. */
+  V1994m,        /*!< \brief 1994m Menter k-w SST model. */
+  V2003m,        /*!< \brief 2003m Menter k-w SST model. */
+  SUST,          /*!< \brief Menter k-w SST model with sustaining terms. */
+  V,             /*!< \brief Menter k-w SST model with vorticity production terms. */
+  KL,            /*!< \brief Menter k-w SST model with Kato-Launder production terms. */
+  UQ,            /*!< \brief Menter k-w SST model with uncertainty quantification modifications. */
+  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. */
+  WALL_OMEGA_LIMIT, /*!< \brief Clip the omega wall value to the upper limit of omega. */
+  TMRBC,         /*!< \brief Far-field omega = 10 U / L_DOMAIN, original reference of the NASA TMR SST page. */
 };
 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("TMRBC", SST_OPTIONS::TMRBC)
+  MakePair("WALL_OMEGA_LIMIT", SST_OPTIONS::WALL_OMEGA_LIMIT)
 };
 
 /*!
@@ -1142,6 +1145,8 @@ 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 wallOmegaLimit = false;                /*!< \brief Bool for clipping the omega wall value (WALL_OMEGA_LIMIT). */
+  bool tmrBC = false;                         /*!< \brief Bool for the far-field values of the NASA TMR (TMRBC). */
 };
 
 /*!
@@ -1159,6 +1164,9 @@ 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_tmrBC = IsPresent(SST_OPTIONS::TMRBC);
+  const bool found_wallOmegaLimit = IsPresent(SST_OPTIONS::WALL_OMEGA_LIMIT);
+
   const bool found_1994 = IsPresent(SST_OPTIONS::V1994);
   const bool found_2003 = IsPresent(SST_OPTIONS::V2003);
   const bool found_1994m = IsPresent(SST_OPTIONS::V1994m);
@@ -1216,6 +1224,8 @@ inline SST_ParsedOptions ParseSSTOptions(const SST_OPTIONS *SST_Options, unsigne
   SSTParsedOptions.compSarkar = sst_compSarkar;
   SSTParsedOptions.dll = sst_dll;
 
+  SSTParsedOptions.tmrBC = found_tmrBC;
+  SSTParsedOptions.wallOmegaLimit = found_wallOmegaLimit;
   return SSTParsedOptions;
 }
 

Common/src/CConfig.cpp

@@ -1202,6 +1202,13 @@ 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 L_DOMAIN \n DESCRIPTION: Approximate length of the computational domain, for the far-field omega of SST_OPTIONS= TMRBC (NASA TMR). \ingroup Config*/
+  addDoubleOption("L_DOMAIN", LDomain, 1.0);
+  /*!\brief SST_SUST_TKE_AMB \n DESCRIPTION: Ambient k of the SST sustaining terms (m^2/s^2), <= 0 for 1e-6 U^2. \ingroup Config*/
+  addDoubleOption("SST_SUST_TKE_AMB", sstSustTkeAmb, 0.0);
+  /*!\brief SST_SUST_OMEGA_AMB \n DESCRIPTION: Ambient omega of the SST sustaining terms (1/s), <= 0 for 5 U / REYNOLDS_LENGTH. \ingroup Config*/
+  addDoubleOption("SST_SUST_OMEGA_AMB", sstSustOmegaAmb, 0.0);
+
   /*!\brief KIND_INCOMP_SYSTEM \n DESCRIPTION: Incomp type \n OPTIONS: see \link Incomp_Map \endlink DEFAULT: NONE \ingroup Config*/
   addEnumOption("KIND_INCOMP_SYSTEM", Kind_Incomp_System, Incomp_Map, INCOMP_SYSTEM::DENSITY_BASED);
   /*!\brief KIND_PB_ITER \n  DESCRIPTION: Kind_PBIter \n OPTIONS: see \link PBIter_Map \endlink \ingroup Config*/

SU2_CFD/include/numerics/turbulent/turb_sources.hpp

@@ -961,16 +961,33 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
       /*--- Production limiter. ---*/
       const su2double prod_limit = prod_lim_const * beta_star * Density_i * ScalarVar_i[1] * ScalarVar_i[0];
 
+      /*--- The modified (m) versions use P = mu_t S^2 (NASA TMR). The standard versions use the exact production,
+       P = tau_ij du_i/dx_j = mu_t (S^2 - 2/3 div(u)^2) - 2/3 rho k div(u); with the vorticity (V) and Kato-Launder (KL)
+       forms of mu_t S^2 only the -2/3 rho k div(u) term is added. ---*/
+      const bool strainProduction = sstParsedOptions.production != SST_OPTIONS::V &&
+                                    sstParsedOptions.production != SST_OPTIONS::KL;
       su2double P = Eddy_Viscosity_i * pow(P_Base, 2);
+      if (!sstParsedOptions.modified) {
+        if (strainProduction) P -= Eddy_Viscosity_i * diverg * diverg * 2.0/3.0;
+        P -= Density_i * ScalarVar_i[0] * diverg * 2.0/3.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);
+        /*--- gamma/nu_t * P, with P/mu_t expanded so that it is defined where mu_t = 0. The last term,
+         * rho k / mu_t, is bounded since mu_t is proportional to k (it vanishes with k). ---*/
+        su2double P_over_muT = pow(P_Base, 2);
+        if (!sstParsedOptions.modified) {
+          if (strainProduction) P_over_muT -= diverg * diverg * 2.0/3.0;
+          P_over_muT -= Density_i * ScalarVar_i[0] * diverg * 2.0/3.0 / max(Eddy_Viscosity_i, EPS);
+        }
+        pw = alfa_blended * Density_i * P_over_muT;
       } else {
         pw = (alfa_blended * Density_i / Eddy_Viscosity_i) * pk;
       }
@@ -996,7 +1013,7 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
       /*--- Dissipation ---*/
 
       su2double dk = beta_star * Density_i * ScalarVar_i[1] * ScalarVar_i[0] * (1.0 + zetaFMt);
-      su2double dw = beta_blended * Density_i * ScalarVar_i[1] * ScalarVar_i[1] * (1.0 - 0.09/beta_blended * zetaFMt);
+      su2double dw = beta_blended * Density_i * ScalarVar_i[1] * ScalarVar_i[1] * (1.0 - beta_star/beta_blended * zetaFMt);
 
       /*--- LM model coupling with production and dissipation term for k transport equation---*/
       if (config->GetKind_Trans_Model() == TURB_TRANS_MODEL::LM) {
@@ -1023,9 +1040,12 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
       /*--- Implicit part ---*/
 
       Jacobian_i[0][0] = -beta_star * ScalarVar_i[1] * Volume * (1.0 + zetaFMt);
+      /*--- Derivative of -2/3 rho k div(u), only where it adds to the diagonal (compression is left out
+       * of the Jacobian but kept in the residual, so that the linear system does not lose diagonal dominance). ---*/
+      if (!sstParsedOptions.modified) Jacobian_i[0][0] -= max(diverg, 0.0) * 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);
+      Jacobian_i[1][1] = -2.0 * beta_blended * ScalarVar_i[1] * Volume * (1.0 - beta_star/beta_blended * zetaFMt);
     }
 
     AD::SetPreaccOut(Residual, nVar);

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

@@ -80,6 +80,7 @@ class CCompressibleViscousFluxBase : public CNumericsSIMD {
   const bool useSA_QCR;
   const bool wallFun;
   const bool uq;
+  const bool tkeInStress; /*!< \brief 2/3 rho k in the stress tensor (standard, non-m, SST versions). */
   const bool uq_permute;
   const size_t uq_eigval_comp;
   const su2double uq_delta_b;
@@ -100,6 +101,8 @@ class CCompressibleViscousFluxBase : public CNumericsSIMD {
     useSA_QCR(config.GetSAParsedOptions().qcr2000),
     wallFun(config.GetWall_Functions()),
     uq(config.GetSSTParsedOptions().uq),
+    tkeInStress(config.GetKind_Turb_Model() == TURB_MODEL::SST && !config.GetSSTParsedOptions().modified &&
+                !config.GetSSTParsedOptions().uq),
     uq_permute(config.GetUQ_Permute()),
     uq_eigval_comp(config.GetEig_Val_Comp()),
     uq_delta_b(config.GetUQ_Delta_B()),
@@ -152,6 +155,14 @@ class CCompressibleViscousFluxBase : public CNumericsSIMD {
     const Double eddyVisc = uq? Double(0.0) : avgV.eddyVisc();
     auto tau = stressTensor(avgV.laminarVisc() + eddyVisc, avgGrad);
     if(useSA_QCR) addQCR(avgGrad, tau, eddyVisc / (avgV.laminarVisc() + eddyVisc));
+    if(tkeInStress) {
+      /*--- 2/3 rho k term of the Boussinesq approximation, ignored by the modified (m) SST versions (with UQ it is
+       * part of the perturbed Reynolds stress). ---*/
+      const Double turb_ke = 0.5*(gatherVariables(iPoint, turbVars->GetSolution()) +
+                                  gatherVariables(jPoint, turbVars->GetSolution()));
+      const Double kTerm = 2.0/3.0 * avgV.density() * turb_ke;
+      for (size_t iDim = 0; iDim < nDim; ++iDim) tau(iDim,iDim) -= kTerm;
+    }
     if(uq) {
       Double turb_ke = 0.5*(gatherVariables(iPoint, turbVars->GetSolution()) +
                             gatherVariables(jPoint, turbVars->GetSolution()));

SU2_CFD/src/numerics/flow/flow_diffusion.cpp

@@ -131,15 +131,22 @@ void CAvgGrad_Base::SetStressTensor(const su2double *val_primvar,
 
   if (sstParsedOptions.uq) {
     // laminar part
+    // the 2/3 rho k term is already in the perturbed Reynolds stress
     ComputeStressTensor(nDim, tau, val_gradprimvar+1, val_laminar_viscosity);
     // add turbulent part which was perturbed
     for (unsigned short iDim = 0 ; iDim < nDim; iDim++)
       for (unsigned short jDim = 0 ; jDim < nDim; jDim++)
         tau[iDim][jDim] += (-Density) * MeanPerturbedRSM[iDim][jDim];
   } else {
     const su2double total_viscosity = val_laminar_viscosity + val_eddy_viscosity;
-    // turb_ke is not considered in the stress tensor, see #797
-    ComputeStressTensor(nDim, tau, val_gradprimvar+1, total_viscosity, Density, su2double(0.0));
+    /*--- 2/3 rho k only for the standard (non-m) SST versions, see #797. The other models call the function
+     * without it, so that their results do not change with how the compiler rounds a zero k term. ---*/
+    const bool tkeInStress = config->GetKind_Turb_Model() == TURB_MODEL::SST && !sstParsedOptions.modified;
+    if (tkeInStress) {
+      ComputeStressTensor(nDim, tau, val_gradprimvar+1, total_viscosity, Density, val_turb_ke);
+    } else {
+      ComputeStressTensor(nDim, tau, val_gradprimvar+1, total_viscosity, Density, su2double(0.0));
+    }
   }
 
   /* --- If the Stochastic Backscatter Model is active, add random contribution to stress tensor ---*/

SU2_CFD/src/solvers/CAdjNSSolver.cpp

@@ -598,6 +598,10 @@ void CAdjNSSolver::Viscous_Sensitivity(CGeometry *geometry, CSolver **solver_con
   dp_drv, dp_drw, dp_drE, dH_dr, dH_dru, dH_drv, dH_drw, dH_drE, H, D[3][3], Dd[3], Mach_Inf, eps, scale = 1.0,
   RefVel2, RefDensity, Mach2Vel, *Velocity_Inf, factor;
 
+  const bool SSTm = config->GetSSTParsedOptions().modified;
+  const bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST) && !SSTm;
+  const auto* Node_Turb = (tkeNeeded) ? solver_container[TURB_SOL]->GetNodes() : nullptr;
+
   auto *USens = new su2double[nVar];
   auto *UnitNormal = new su2double[nDim];
   auto *normal_grad_vel = new su2double[nDim];
@@ -766,6 +770,7 @@ void CAdjNSSolver::Viscous_Sensitivity(CGeometry *geometry, CSolver **solver_con
             /*--- Turbulent kinetic energy ---*/
             // turb_ke is not considered in the stress tensor, see #797
             val_turb_ke = 0.0;
+            if (tkeNeeded) val_turb_ke = Node_Turb->GetSolution(iPoint, 0);
             CNumerics::ComputeStressTensor(nDim, tau, PrimVar_Grad+1, Laminar_Viscosity, Density, val_turb_ke);
 
             /*--- Form normal_grad_gridvel = \partial_n (u_omega) ---*/