From 482d18a93090457f60e03c0eec1af49adc6b3e5c Mon Sep 17 00:00:00 2001 From: rois1995 Date: Wed, 30 Sep 2026 14:48:19 +0200 Subject: [PATCH 1/5] Enable the standard SST versions, add the TMRBC, WALL_OMEGA_LIMIT and ambient value options - Accept V1994 and V2003 in SST_OPTIONS (the standard versions of the NASA TMR, with the exact production and 2/3 rho k in the stress tensor), and add TMRBC (far-field omega of the original reference of the NASA TMR) and WALL_OMEGA_LIMIT (clip the omega wall value to the upper limit of omega). - New options L_DOMAIN (approximate length of the domain, for TMRBC) and SST_SUST_TKE_AMB / SST_SUST_OMEGA_AMB (ambient values of the sustaining terms, <= 0 for the defaults of Spalart and Rumsey). - Document the versions, the new options and the range of mu_t/mu of the NASA TMR in config_template.cfg. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6 --- Common/include/CConfig.hpp | 24 +++++++++++++++++ Common/include/option_structure.hpp | 40 ++++++++++++++++++----------- Common/src/CConfig.cpp | 7 +++++ config_template.cfg | 20 +++++++++++++-- 4 files changed, 74 insertions(+), 17 deletions(-) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index bf5bbd0684cb..4a2128362e5b 100644 --- a/Common/include/CConfig.hpp +++ b/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. diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index b8222da53784..73df22810ef7 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -1101,26 +1101,27 @@ 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 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) @@ -1128,6 +1129,8 @@ static const MapType SST_Options_Map = { 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); @@ -1213,6 +1221,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; } diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 852909cff588..d32b309f20d6 100644 --- a/Common/src/CConfig.cpp +++ b/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*/ diff --git a/config_template.cfg b/config_template.cfg index 5c7132467cc4..ca035b77ea35 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -21,9 +21,24 @@ SOLVER= EULER % Specify turbulence model (NONE, SA, SST) KIND_TURB_MODEL= NONE % -% Specify versions/corrections of the SST model (V2003m, V1994m, VORTICITY, KATO_LAUNDER, UQ, SUSTAINING, COMPRESSIBILITY-WILCOX, COMPRESSIBILITY-SARKAR, DIMENSIONLESS_LIMIT) +% Specify versions/corrections of the SST model (V2003m, V1994m, V2003, V1994, VORTICITY, KATO_LAUNDER, UQ, SUSTAINING, +% COMPRESSIBILITY-WILCOX, COMPRESSIBILITY-SARKAR, DIMENSIONLESS_LIMIT, TMRBC, WALL_OMEGA_LIMIT). The "m" versions ignore 2/3 rho k in the +% stress tensor and use P = mu_t S^2 (NASA TMR); V2003 and V1994 are the standard versions with the exact production. +% TMRBC sets the far-field (and inlet) omega as in the original reference of the NASA TMR, omega = 10 U / L_DOMAIN, +% with k from FREESTREAM_TURB2LAMVISCRATIO (or the inlet viscosity ratio). +% WALL_OMEGA_LIMIT clips the omega wall value 60 nu / (beta_1 d^2) to the upper limit of omega (off by default). SST_OPTIONS= NONE % +% Approximate length of the computational domain for TMRBC (m) +L_DOMAIN= 1.0 +% +% Ambient values of the SST sustaining terms (SST_OPTIONS= SUSTAINING), also used as free-stream +% values of k (m^2/s^2) and omega (1/s). Values <= 0 select the defaults of Spalart and Rumsey +% (AIAA J 45(10), 2007) used by the NASA TMR SST-sust: k = 1e-6 U^2 (Tu = 0.08165%) and +% omega = 5 U / L, with U the free-stream velocity and L = REYNOLDS_LENGTH. +SST_SUST_TKE_AMB= 0.0 +SST_SUST_OMEGA_AMB= 0.0 +% % Specify versions/corrections of the SA model (NEGATIVE, EDWARDS, WITHFT2, QCR2000, COMPRESSIBILITY, ROTATION, BCM, EXPERIMENTAL) SA_OPTIONS= NONE % @@ -277,7 +292,8 @@ TURB_FIXED_VALUES= NO % normalized far-field velocity is less than this parameter. TURB_FIXED_VALUES_DOMAIN= -1.0 % -% Free-stream ratio between turbulent and laminar viscosity +% Free-stream ratio between turbulent and laminar viscosity. For the SST model the NASA TMR (original reference +% of Menter) recommends far-field values giving mu_t/mu between 1e-5 and 1e-2. FREESTREAM_TURB2LAMVISCRATIO= 10.0 % % Compressible flow non-dimensionalization (DIMENSIONAL, FREESTREAM_PRESS_EQ_ONE, From fa241d72c7c7427fd8afed7fb9ce53ff19e41783 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Wed, 30 Sep 2026 14:48:54 +0200 Subject: [PATCH 2/5] Use the exact production for the standard SST versions, fix the CDkw clip per version - V1994 and V2003 use the exact production of the NASA TMR, P = mu_t (S^2 - 2/3 div(u)^2) - 2/3 rho k div(u) (with VORTICITY and KATO-LAUNDER only the k term is added); the m versions keep P = mu_t S^2. - V1994 evaluates gamma/nu_t P with P/mu_t expanded, guarded by max(mu_t, EPS) where mu_t = 0. - Add the derivative of -2/3 rho k div(u) to the Jacobian diagonal only where it increases it (div(u) > 0). - Set the lower bound of CD_kw from the version (1e-20 for 1994, 1e-10 for 2003, as in the NASA TMR) instead of from the production limiter constant, which happens to have the same exponent. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6 --- .../numerics/turbulent/turb_sources.hpp | 22 ++++++++++++++++++- SU2_CFD/src/variables/CTurbSSTVariable.cpp | 4 +++- 2 files changed, 24 insertions(+), 2 deletions(-) diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp index d1cfd2f09054..405be31e7f86 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp @@ -955,7 +955,17 @@ 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 su2double eddy_visc_var = sstParsedOptions.version == SST_OPTIONS::V1994 ? VorticityMag : StrainMag_i; @@ -964,7 +974,14 @@ class CSourcePieceWise_TurbSST final : public CNumerics { /*--- 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; } @@ -1017,6 +1034,9 @@ 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 - beta_star/beta_blended * zetaFMt); diff --git a/SU2_CFD/src/variables/CTurbSSTVariable.cpp b/SU2_CFD/src/variables/CTurbSSTVariable.cpp index a0a08669c766..6fc2d48ce147 100644 --- a/SU2_CFD/src/variables/CTurbSSTVariable.cpp +++ b/SU2_CFD/src/variables/CTurbSSTVariable.cpp @@ -69,7 +69,9 @@ void CTurbSSTVariable::SetBlendingFunc(unsigned long iPoint, su2double val_visco for (unsigned long iDim = 0; iDim < nDim; iDim++) CDkw(iPoint) += Gradient(iPoint,0,iDim)*Gradient(iPoint,1,iDim); CDkw(iPoint) *= 2.0*val_density*sigma_om2/Solution(iPoint,1); - CDkw(iPoint) = max(CDkw(iPoint), pow(10.0, -prod_lim_const)); + su2double exponent = 20.0; + if (sstParsedOptions.version == SST_OPTIONS::V2003) exponent = 10.0; + CDkw(iPoint) = max(CDkw(iPoint), pow(10.0, -exponent)); /*--- F1 ---*/ From 274cecbec47706c02098883b99dde3d57cc40a9a Mon Sep 17 00:00:00 2001 From: rois1995 Date: Wed, 30 Sep 2026 14:49:49 +0200 Subject: [PATCH 3/5] Include 2/3 rho k in the stress tensor for the standard SST versions The standard (non-m) versions of the NASA TMR keep -2/3 rho k in the Boussinesq stress tensor, which the m versions (and #797) ignore. Add it, only for V1994 and V2003, in: - the viscous fluxes (CAvgGrad_Base::SetStressTensor and the vectorized viscous fluxes; with UQ it is already part of the perturbed Reynolds stress), - the viscous part of the axisymmetric source (radial normal stress in the energy equation), - the wall shear stress of the wall functions (CNSSolver::SetTau_Wall_WF), - the viscous sensitivities of the continuous adjoint (CAdjNSSolver::Viscous_Sensitivity). The other models keep calling the stress tensor without k, so that their results do not change. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6 --- SU2_CFD/include/numerics/flow/flow_sources.hpp | 1 + .../numerics_simd/flow/diffusion/viscous_fluxes.hpp | 11 +++++++++++ SU2_CFD/src/numerics/flow/flow_diffusion.cpp | 11 +++++++++-- SU2_CFD/src/numerics/flow/flow_sources.cpp | 7 ++++++- SU2_CFD/src/solvers/CAdjNSSolver.cpp | 5 +++++ SU2_CFD/src/solvers/CNSSolver.cpp | 10 +++++++++- 6 files changed, 41 insertions(+), 4 deletions(-) diff --git a/SU2_CFD/include/numerics/flow/flow_sources.hpp b/SU2_CFD/include/numerics/flow/flow_sources.hpp index ca0613a0ebfc..cf06eb753d8a 100644 --- a/SU2_CFD/include/numerics/flow/flow_sources.hpp +++ b/SU2_CFD/include/numerics/flow/flow_sources.hpp @@ -74,6 +74,7 @@ class CSourceAxisymmetric_Flow : public CSourceBase_Flow { protected: bool implicit, viscous; bool tkeInEnergy; /*!< \brief SST: k is part of the total energy (turb_ke_i is set by the solver). */ + bool tkeInStress; /*!< \brief Standard (non-m) SST versions: -2/3 rho k is part of the stress tensor. */ su2double yinv{0.0}; /*! diff --git a/SU2_CFD/include/numerics_simd/flow/diffusion/viscous_fluxes.hpp b/SU2_CFD/include/numerics_simd/flow/diffusion/viscous_fluxes.hpp index 1a8745b10bc2..57466b20068b 100644 --- a/SU2_CFD/include/numerics_simd/flow/diffusion/viscous_fluxes.hpp +++ b/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())); diff --git a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp index 66945b2f700a..6ac65167c99b 100644 --- a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp +++ b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp @@ -131,6 +131,7 @@ 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++) @@ -138,8 +139,14 @@ void CAvgGrad_Base::SetStressTensor(const su2double *val_primvar, 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 ---*/ diff --git a/SU2_CFD/src/numerics/flow/flow_sources.cpp b/SU2_CFD/src/numerics/flow/flow_sources.cpp index 44e5cf579dcf..8c2b004de543 100644 --- a/SU2_CFD/src/numerics/flow/flow_sources.cpp +++ b/SU2_CFD/src/numerics/flow/flow_sources.cpp @@ -55,6 +55,7 @@ CSourceAxisymmetric_Flow::CSourceAxisymmetric_Flow(unsigned short val_nDim, unsi implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); viscous = config->GetViscous(); tkeInEnergy = (config->GetKind_Turb_Model() == TURB_MODEL::SST); + tkeInStress = tkeInEnergy && !config->GetSSTParsedOptions().modified; } @@ -164,7 +165,10 @@ CNumerics::ResidualType<> CSourceAxisymmetric_Flow::ComputeResidual(const CConfi void CSourceAxisymmetric_Flow::ResidualDiffusion(){ - su2double laminar_viscosity_i = V_i[nDim+5]; + /*--- -2/3 rho k of the radial normal stress (only where k is part of the stress tensor). ---*/ + const su2double tke = tkeInStress ? turb_ke_i : 0.0; + + su2double laminar_viscosity_i = V_i[nDim+5]; su2double eddy_viscosity_i = V_i[nDim+6]; su2double thermal_conductivity_i = V_i[nDim+7]; su2double heat_capacity_cp_i = V_i[nDim+8]; @@ -183,6 +187,7 @@ void CSourceAxisymmetric_Flow::ResidualDiffusion(){ residual[3] -= Volume*(yinv*(total_viscosity_i*(u*(PrimVar_Grad_i[2][0]+PrimVar_Grad_i[1][1]) +v*TWO3*(2*PrimVar_Grad_i[2][1]-PrimVar_Grad_i[1][0] -v*yinv)) + -v*TWO3*U_i[0]*tke +total_conductivity_i*PrimVar_Grad_i[0][1]) -TWO3*(AuxVar_Grad_i[1][1]+AuxVar_Grad_i[2][0])); } diff --git a/SU2_CFD/src/solvers/CAdjNSSolver.cpp b/SU2_CFD/src/solvers/CAdjNSSolver.cpp index 3403522f109d..4c3c628b65b9 100644 --- a/SU2_CFD/src/solvers/CAdjNSSolver.cpp +++ b/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) ---*/ diff --git a/SU2_CFD/src/solvers/CNSSolver.cpp b/SU2_CFD/src/solvers/CNSSolver.cpp index 643cb9a4f6c9..860a6153d58d 100644 --- a/SU2_CFD/src/solvers/CNSSolver.cpp +++ b/SU2_CFD/src/solvers/CNSSolver.cpp @@ -810,6 +810,9 @@ void CNSSolver::SetTau_Wall_WF(CGeometry *geometry, CSolver **solver_container, const unsigned short max_iter = config->GetwallModel_MaxIter(); const su2double relax = config->GetwallModel_RelFac(); + const bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST && !(config->GetSSTParsedOptions().modified)); + const auto* Node_Turb = (tkeNeeded) ? solver_container[TURB_SOL]->GetNodes() : nullptr; + /*--- Compute the recovery factor * use Molecular (Laminar) Prandtl number (see Nichols & Nelson, nomenclature ) ---*/ @@ -909,7 +912,12 @@ void CNSSolver::SetTau_Wall_WF(CGeometry *geometry, CSolver **solver_container, const su2double Lam_Visc_Wall = nodes->GetLaminarViscosity(iPoint); su2double Eddy_Visc_Wall = nodes->GetEddyViscosity(iPoint); - CNumerics::ComputeStressTensor(nDim, tau, nodes->GetVelocityGradient(iPoint), Lam_Visc_Wall); + if (tkeNeeded) { + CNumerics::ComputeStressTensor(nDim, tau, nodes->GetVelocityGradient(iPoint), Lam_Visc_Wall, + nodes->GetDensity(iPoint), Node_Turb->GetSolution(iPoint, 0)); + } else { + CNumerics::ComputeStressTensor(nDim, tau, nodes->GetVelocityGradient(iPoint), Lam_Visc_Wall); + } su2double TauTangent[MAXNDIM] = {0.0}; GeometryToolbox::TangentProjection(nDim, tau, UnitNormal, TauTangent); From dc104a2a429d4eaff0778b949def0cb4ce358d03 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Wed, 30 Sep 2026 14:50:33 +0200 Subject: [PATCH 4/5] Add the far-field values of the NASA TMR, the SST-sust ambient values and WALL_OMEGA_LIMIT - SST_OPTIONS= TMRBC: far-field (and inlet) omega = 10 U / L_DOMAIN, the original reference of the NASA TMR SST page, with k from FREESTREAM_TURB2LAMVISCRATIO (or the inlet viscosity ratio), in the non-dimensionalization of all flow solvers (also in the k of the free-stream energy), in the far-field state of CTurbSSTSolver and in BC_Inlet. - SST_OPTIONS= SUSTAINING: the ambient values of the sustaining terms, SST_SUST_TKE_AMB and SST_SUST_OMEGA_AMB (defaults 1e-6 U^2 and 5 U / REYNOLDS_LENGTH of Spalart and Rumsey, as the NASA TMR SST-sust), are also the free-stream values; mention them in the error for non-positive free-stream values. - SST_OPTIONS= WALL_OMEGA_LIMIT: clip the omega wall value 60 nu / (beta_1 d^2), which grows as 1/d^2 on fine wall grids, to the upper limit of omega used in the rest of the domain. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6 --- SU2_CFD/src/solvers/CEulerSolver.cpp | 40 ++++++++++++++++----- SU2_CFD/src/solvers/CFEM_DG_EulerSolver.cpp | 24 ++++++++++--- SU2_CFD/src/solvers/CIncEulerSolver.cpp | 24 ++++++++++--- SU2_CFD/src/solvers/CNEMOEulerSolver.cpp | 24 ++++++++++--- SU2_CFD/src/solvers/CTurbSSTSolver.cpp | 27 ++++++++++++-- 5 files changed, 113 insertions(+), 26 deletions(-) diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 8d490cb2c9d1..d584f80fc9f3 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -1012,6 +1012,17 @@ void CEulerSolver::SetNondimensionalization(CConfig *config, unsigned short iMes Tke_FreeStream = 3.0/2.0*(ModVel_FreeStream*ModVel_FreeStream*config->GetTurbulenceIntensity_FreeStream()*config->GetTurbulenceIntensity_FreeStream()); + if (config->GetSSTParsedOptions().tmrBC) { + su2double Omega_Freestream = 10 * ModVel_FreeStream / config->GetLDomain(); + Tke_FreeStream = Omega_Freestream*(Viscosity_FreeStream*config->GetTurb2LamViscRatio_FreeStream())/Density_FreeStream; + } else if (config->GetSSTParsedOptions().sust) { + Tke_FreeStream = config->GetSSTSust_TkeAmb(ModVel_FreeStream); + } + + /*-- Compute the freestream energy. ---*/ + + if (tkeNeeded) Energy_FreeStream += Tke_FreeStream; + } else { @@ -1022,9 +1033,8 @@ void CEulerSolver::SetNondimensionalization(CConfig *config, unsigned short iMes } - /*-- Compute the freestream energy. ---*/ - - if (tkeNeeded) { Energy_FreeStream += Tke_FreeStream; }; config->SetEnergy_FreeStream(Energy_FreeStream); + /*-- Set the freestream energy. ---*/ + config->SetEnergy_FreeStream(Energy_FreeStream); /*--- Compute non dimensional quantities. By definition, Lref is one because we have converted the grid to meters. ---*/ @@ -1089,15 +1099,29 @@ void CEulerSolver::SetNondimensionalization(CConfig *config, unsigned short iMes config->SetSpecificHeatCp_FreeStreamND(SpecificHeat_Cp_FreeStreamND); Tke_FreeStream = 3.0/2.0*(ModVel_FreeStream*ModVel_FreeStream*config->GetTurbulenceIntensity_FreeStream()*config->GetTurbulenceIntensity_FreeStream()); - config->SetTke_FreeStream(Tke_FreeStream); - Tke_FreeStreamND = 3.0/2.0*(ModVel_FreeStreamND*ModVel_FreeStreamND*config->GetTurbulenceIntensity_FreeStream()*config->GetTurbulenceIntensity_FreeStream()); - config->SetTke_FreeStreamND(Tke_FreeStreamND); Omega_FreeStream = Density_FreeStream*Tke_FreeStream/max(Viscosity_FreeStream*config->GetTurb2LamViscRatio_FreeStream(), EPS); - config->SetOmega_FreeStream(Omega_FreeStream); - Omega_FreeStreamND = Density_FreeStreamND*Tke_FreeStreamND/max(Viscosity_FreeStreamND*config->GetTurb2LamViscRatio_FreeStream(), EPS); + + if (config->GetSSTParsedOptions().tmrBC) { + Omega_FreeStream = 10 * ModVel_FreeStream / config->GetLDomain(); + Omega_FreeStreamND = 10 * ModVel_FreeStreamND / config->GetLDomain(); // the reference length is 1 + + Tke_FreeStream = Omega_FreeStream*(Viscosity_FreeStream*config->GetTurb2LamViscRatio_FreeStream())/Density_FreeStream; + Tke_FreeStreamND = Omega_FreeStreamND*(Viscosity_FreeStreamND*config->GetTurb2LamViscRatio_FreeStream())/Density_FreeStreamND; + } else if (config->GetSSTParsedOptions().sust) { + /*--- Ambient values of the sustaining terms, also used as free-stream values. ---*/ + Omega_FreeStream = config->GetSSTSust_OmegaAmb(ModVel_FreeStream); + Omega_FreeStreamND = Omega_FreeStream / Omega_Ref; + Tke_FreeStream = config->GetSSTSust_TkeAmb(ModVel_FreeStream); + Tke_FreeStreamND = Tke_FreeStream / pow(Velocity_Ref, 2); + } + + config->SetTke_FreeStream(Tke_FreeStream); + config->SetTke_FreeStreamND(Tke_FreeStreamND); + + config->SetOmega_FreeStream(Omega_FreeStream); config->SetOmega_FreeStreamND(Omega_FreeStreamND); if (config->GetTurbulenceIntensity_FreeStream() *100 <= 1.3) { diff --git a/SU2_CFD/src/solvers/CFEM_DG_EulerSolver.cpp b/SU2_CFD/src/solvers/CFEM_DG_EulerSolver.cpp index 0f5db8c0f8a7..a32eba8e1e5e 100644 --- a/SU2_CFD/src/solvers/CFEM_DG_EulerSolver.cpp +++ b/SU2_CFD/src/solvers/CFEM_DG_EulerSolver.cpp @@ -1076,15 +1076,29 @@ void CFEM_DG_EulerSolver::SetNondimensionalization(CConfig *config, config->SetSpecificHeatCp_FreeStreamND(SpecificHeat_Cp_FreeStreamND); Tke_FreeStream = 3.0/2.0*(ModVel_FreeStream*ModVel_FreeStream*config->GetTurbulenceIntensity_FreeStream()*config->GetTurbulenceIntensity_FreeStream()); - config->SetTke_FreeStream(Tke_FreeStream); - Tke_FreeStreamND = 3.0/2.0*(ModVel_FreeStreamND*ModVel_FreeStreamND*config->GetTurbulenceIntensity_FreeStream()*config->GetTurbulenceIntensity_FreeStream()); - config->SetTke_FreeStreamND(Tke_FreeStreamND); Omega_FreeStream = Density_FreeStream*Tke_FreeStream/max(Viscosity_FreeStream*config->GetTurb2LamViscRatio_FreeStream(), EPS); - config->SetOmega_FreeStream(Omega_FreeStream); - Omega_FreeStreamND = Density_FreeStreamND*Tke_FreeStreamND/max(Viscosity_FreeStreamND*config->GetTurb2LamViscRatio_FreeStream(), EPS); + + if (config->GetSSTParsedOptions().tmrBC) { + Omega_FreeStream = 10 * ModVel_FreeStream / config->GetLDomain(); + Omega_FreeStreamND = 10 * ModVel_FreeStreamND / config->GetLDomain(); // the reference length is 1 + + Tke_FreeStream = Omega_FreeStream*(Viscosity_FreeStream*config->GetTurb2LamViscRatio_FreeStream())/Density_FreeStream; + Tke_FreeStreamND = Omega_FreeStreamND*(Viscosity_FreeStreamND*config->GetTurb2LamViscRatio_FreeStream())/Density_FreeStreamND; + } else if (config->GetSSTParsedOptions().sust) { + /*--- Ambient values of the sustaining terms, also used as free-stream values. ---*/ + Omega_FreeStream = config->GetSSTSust_OmegaAmb(ModVel_FreeStream); + Omega_FreeStreamND = Omega_FreeStream / Omega_Ref; + Tke_FreeStream = config->GetSSTSust_TkeAmb(ModVel_FreeStream); + Tke_FreeStreamND = Tke_FreeStream / pow(Velocity_Ref, 2); + } + + config->SetTke_FreeStream(Tke_FreeStream); + config->SetTke_FreeStreamND(Tke_FreeStreamND); + + config->SetOmega_FreeStream(Omega_FreeStream); config->SetOmega_FreeStreamND(Omega_FreeStreamND); /*--- Initialize the dimensionless Fluid Model that will be used to solve the dimensionless problem ---*/ diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index 8c6ee1719a01..ab968eb125af 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -494,15 +494,29 @@ void CIncEulerSolver::SetNondimensionalization(CConfig *config, unsigned short i config->SetSpecificHeatCp_FreeStreamND(SpecificHeat_Cp_FreeStreamND); Tke_FreeStream = 3.0/2.0*(ModVel_FreeStream*ModVel_FreeStream*config->GetTurbulenceIntensity_FreeStream()*config->GetTurbulenceIntensity_FreeStream()); - config->SetTke_FreeStream(Tke_FreeStream); - Tke_FreeStreamND = 3.0/2.0*(ModVel_FreeStreamND*ModVel_FreeStreamND*config->GetTurbulenceIntensity_FreeStream()*config->GetTurbulenceIntensity_FreeStream()); - config->SetTke_FreeStreamND(Tke_FreeStreamND); Omega_FreeStream = Density_FreeStream*Tke_FreeStream/max(Viscosity_FreeStream*config->GetTurb2LamViscRatio_FreeStream(), EPS); - config->SetOmega_FreeStream(Omega_FreeStream); - Omega_FreeStreamND = Density_FreeStreamND*Tke_FreeStreamND/max(Viscosity_FreeStreamND*config->GetTurb2LamViscRatio_FreeStream(), EPS); + + if (config->GetSSTParsedOptions().tmrBC) { + Omega_FreeStream = 10 * ModVel_FreeStream / config->GetLDomain(); + Omega_FreeStreamND = 10 * ModVel_FreeStreamND / config->GetLDomain(); // the reference length is 1 + + Tke_FreeStream = Omega_FreeStream*(Viscosity_FreeStream*config->GetTurb2LamViscRatio_FreeStream())/Density_FreeStream; + Tke_FreeStreamND = Omega_FreeStreamND*(Viscosity_FreeStreamND*config->GetTurb2LamViscRatio_FreeStream())/Density_FreeStreamND; + } else if (config->GetSSTParsedOptions().sust) { + /*--- Ambient values of the sustaining terms, also used as free-stream values. ---*/ + Omega_FreeStream = config->GetSSTSust_OmegaAmb(ModVel_FreeStream); + Omega_FreeStreamND = Omega_FreeStream / Omega_Ref; + Tke_FreeStream = config->GetSSTSust_TkeAmb(ModVel_FreeStream); + Tke_FreeStreamND = Tke_FreeStream / pow(Velocity_Ref, 2); + } + + config->SetTke_FreeStream(Tke_FreeStream); + config->SetTke_FreeStreamND(Tke_FreeStreamND); + + config->SetOmega_FreeStream(Omega_FreeStream); config->SetOmega_FreeStreamND(Omega_FreeStreamND); const su2double MassDiffusivityND = config->GetDiffusivity_Constant() / (Velocity_Ref * Length_Ref); diff --git a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp index 2a09a7a1f42e..7a0895731c08 100644 --- a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp @@ -1192,15 +1192,29 @@ void CNEMOEulerSolver::SetNondimensionalization(CConfig *config, unsigned short Viscosity_FreeStreamND = Viscosity_FreeStream / Viscosity_Ref; config->SetViscosity_FreeStreamND(Viscosity_FreeStreamND); Tke_FreeStream = 3.0/2.0*(ModVel_FreeStream*ModVel_FreeStream*config->GetTurbulenceIntensity_FreeStream()*config->GetTurbulenceIntensity_FreeStream()); - config->SetTke_FreeStream(Tke_FreeStream); - Tke_FreeStreamND = 3.0/2.0*(ModVel_FreeStreamND*ModVel_FreeStreamND*config->GetTurbulenceIntensity_FreeStream()*config->GetTurbulenceIntensity_FreeStream()); - config->SetTke_FreeStreamND(Tke_FreeStreamND); Omega_FreeStream = Density_FreeStream*Tke_FreeStream/max(Viscosity_FreeStream*config->GetTurb2LamViscRatio_FreeStream(), EPS); - config->SetOmega_FreeStream(Omega_FreeStream); - Omega_FreeStreamND = Density_FreeStreamND*Tke_FreeStreamND/max(Viscosity_FreeStreamND*config->GetTurb2LamViscRatio_FreeStream(), EPS); + + if (config->GetSSTParsedOptions().tmrBC) { + Omega_FreeStream = 10 * ModVel_FreeStream / config->GetLDomain(); + Omega_FreeStreamND = 10 * ModVel_FreeStreamND / config->GetLDomain(); // the reference length is 1 + + Tke_FreeStream = Omega_FreeStream*(Viscosity_FreeStream*config->GetTurb2LamViscRatio_FreeStream())/Density_FreeStream; + Tke_FreeStreamND = Omega_FreeStreamND*(Viscosity_FreeStreamND*config->GetTurb2LamViscRatio_FreeStream())/Density_FreeStreamND; + } else if (config->GetSSTParsedOptions().sust) { + /*--- Ambient values of the sustaining terms, also used as free-stream values. ---*/ + Omega_FreeStream = config->GetSSTSust_OmegaAmb(ModVel_FreeStream); + Omega_FreeStreamND = Omega_FreeStream / Omega_Ref; + Tke_FreeStream = config->GetSSTSust_TkeAmb(ModVel_FreeStream); + Tke_FreeStreamND = Tke_FreeStream / pow(Velocity_Ref, 2); + } + + config->SetTke_FreeStream(Tke_FreeStream); + config->SetTke_FreeStreamND(Tke_FreeStreamND); + + config->SetOmega_FreeStream(Omega_FreeStream); config->SetOmega_FreeStreamND(Omega_FreeStreamND); /*--- Initialize the dimensionless Fluid Model that will be used to solve the dimensionless problem ---*/ diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp index 1377defbaefc..a7ba084788bd 100644 --- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp @@ -124,13 +124,23 @@ CTurbSSTSolver::CTurbSSTSolver(CGeometry *geometry, CConfig *config, const CSolv su2double kine_Inf = 3.0/2.0*(VelMag2*Intensity*Intensity); su2double omega_Inf = rhoInf*kine_Inf/(muLamInf*viscRatio); + if (sstParsedOptions.tmrBC) { + omega_Inf = 10 * sqrt(VelMag2) / config->GetLDomain(); + kine_Inf = omega_Inf*(muLamInf*viscRatio)/rhoInf; + } else if (sstParsedOptions.sust) { + /*--- Ambient values of the sustaining terms, set by the flow solver (SST_SUST_TKE_AMB, SST_SUST_OMEGA_AMB). ---*/ + kine_Inf = config->GetTke_FreeStreamND(); + omega_Inf = config->GetOmega_FreeStreamND(); + } + /*--- A zero turbulence intensity, viscosity ratio or free-stream velocity gives k = 0 or omega = 0 (or infinite), and then mu_t = 0/0 in the whole field. ---*/ if (!(kine_Inf > 0.0 && omega_Inf > 0.0 && std::isfinite(SU2_TYPE::GetValue(kine_Inf)) && std::isfinite(SU2_TYPE::GetValue(omega_Inf)))) { SU2_MPI::Error("SST: the free-stream k and omega must be positive and finite.\n" "Set FREESTREAM_TURBULENCEINTENSITY > 0 and FREESTREAM_TURB2LAMVISCRATIO > 0 with a nonzero\n" - "free-stream velocity.", CURRENT_FUNCTION); + "free-stream velocity, or (e.g. for MACH_NUMBER= 0 with grid motion) SST_OPTIONS= SUSTAINING with\n" + "positive SST_SUST_TKE_AMB and SST_SUST_OMEGA_AMB.", CURRENT_FUNCTION); } Solution_Inf[0] = kine_Inf; @@ -522,6 +532,12 @@ void CTurbSSTSolver::BC_HeatFlux_Wall(CGeometry *geometry, CSolver **solver_cont su2double beta_1 = constants[4]; solution[0] = 0.0; solution[1] = 60.0*laminar_viscosity/(density*beta_1*pow(wall_dist,2)); + + /*--- Menter's wall value, omega_w = 10 * 6 nu / (beta_1 d^2) (AIAA J 32(8), 1994), grows as 1/d^2 with the + distance d of the first point off the wall, so on very fine wall grids it can exceed the upper limit used + to clip omega in the rest of the domain (upperlimit[1]). SST_OPTIONS= WALL_OMEGA_LIMIT clips it to the + same limit, so that the wall and the interior values stay consistent. Off by default. ---*/ + if (sstParsedOptions.wallOmegaLimit) solution[1] = min(solution[1], upperlimit[1]); } /*--- Set the solution values and zero the residual ---*/ @@ -649,8 +665,13 @@ void CTurbSSTSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, C const su2double viscRatio = Turb_Properties[1]; const su2double VelMag2 = GeometryToolbox::SquaredNorm(nDim, Velocity_Inlet); - Inlet_Vars[0] = 3.0 / 2.0 * (VelMag2 * pow(Intensity, 2)); - Inlet_Vars[1] = Density_Inlet * Inlet_Vars[0] / (Laminar_Viscosity_Inlet * viscRatio); + if (sstParsedOptions.tmrBC) { + Inlet_Vars[1] = 10 * sqrt(VelMag2) / config->GetLDomain(); + Inlet_Vars[0] = Inlet_Vars[1]*(Laminar_Viscosity_Inlet*viscRatio)/Density_Inlet; + } else { + Inlet_Vars[0] = 3.0 / 2.0 * (VelMag2 * pow(Intensity, 2)); + Inlet_Vars[1] = Density_Inlet * Inlet_Vars[0] / (Laminar_Viscosity_Inlet * viscRatio); + } } for (auto iVar = 0u; iVar < nVar; iVar++) ghostNodes->SetSolution(iVertex, iVar, Inlet_Vars[iVar]); From d8079252803090619b5b7bbb630c2ab5448ae91b Mon Sep 17 00:00:00 2001 From: rois1995 Date: Wed, 30 Sep 2026 15:39:14 +0200 Subject: [PATCH 5/5] Regression test of the standard SST-2003 with the NASA TMR far-field values Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6 --- TestCases/parallel_regression.py | 8 ++ .../rans/flatplate/turb_SST_flatplate_tmr.cfg | 114 ++++++++++++++++++ 2 files changed, 122 insertions(+) create mode 100644 TestCases/rans/flatplate/turb_SST_flatplate_tmr.cfg diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 8fff340267e9..e8fc08ed599d 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -457,6 +457,14 @@ def main(): turb_flatplate.test_vals = [-0.187388, 0.003749, 10.000000, -1.406182] test_list.append(turb_flatplate) + # Flat plate, standard SST-2003 with the far-field values of the NASA TMR + turb_flatplate_sst_tmr = TestCase('turb_flatplate_sst_tmr') + turb_flatplate_sst_tmr.cfg_dir = "rans/flatplate" + turb_flatplate_sst_tmr.cfg_file = "turb_SST_flatplate_tmr.cfg" + turb_flatplate_sst_tmr.test_iter = 19 + turb_flatplate_sst_tmr.test_vals = [-5.225275, 0.259404, -5.804219, 7.702339, -0.187496, 0.003596] + test_list.append(turb_flatplate_sst_tmr) + # Flat plate (compressible) with species inlet turb_flatplate_species = TestCase('turb_flatplate_species') turb_flatplate_species.cfg_dir = "rans/flatplate" diff --git a/TestCases/rans/flatplate/turb_SST_flatplate_tmr.cfg b/TestCases/rans/flatplate/turb_SST_flatplate_tmr.cfg new file mode 100644 index 000000000000..ae71283c5e19 --- /dev/null +++ b/TestCases/rans/flatplate/turb_SST_flatplate_tmr.cfg @@ -0,0 +1,114 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Turbulent flow over flat plate with zero pressure gradient % +% Author: Thomas D. Economon % +% Institution: Stanford University % +% Date: 2011.11.10 % +% File Version 8.5.0 "Harrier" % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +% +SOLVER= RANS +KIND_TURB_MODEL= SST +% Standard SST-2003 with the far-field values of the NASA TMR (omega = 10 U / L_DOMAIN) +SST_OPTIONS= V2003, TMRBC +L_DOMAIN= 2.33333 +MATH_PROBLEM= DIRECT +RESTART_SOL= NO + +% ----------- COMPRESSIBLE AND INCOMPRESSIBLE FREE-STREAM DEFINITION ----------% +% +MACH_NUMBER= 0.2 +AOA= 0.0 +SIDESLIP_ANGLE= 0.0 +FREESTREAM_TEMPERATURE= 300.0 +REYNOLDS_NUMBER= 5000000.0 +REYNOLDS_LENGTH= 1.0 +FREESTREAM_TURB2LAMVISCRATIO= 0.009 + +% ---------------------- REFERENCE VALUE DEFINITION ---------------------------% +% +REF_ORIGIN_MOMENT_X = 0.25 +REF_ORIGIN_MOMENT_Y = 0.00 +REF_ORIGIN_MOMENT_Z = 0.00 +REF_LENGTH= 1.0 +REF_AREA= 2.0 + +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +MARKER_HEATFLUX= ( wall, 0.0 ) +MARKER_INLET= ( inlet, 302.4, 118309.784, 1.0, 0.0, 0.0 ) +MARKER_OUTLET= ( outlet, 115056.0, farfield, 115056.0 ) +MARKER_SYM= ( symmetry ) +MARKER_PLOTTING= ( wall ) +MARKER_MONITORING= ( wall ) + +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES +CFL_NUMBER= 1000.0 +CFL_ADAPT= NO +CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.0, 100.0 ) +RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 ) +ITER= 20 + +% ----------------------- SLOPE LIMITER DEFINITION ----------------------------% +% +VENKAT_LIMITER_COEFF= 0.1 +ADJ_SHARP_LIMITER_COEFF= 3.0 +REF_SHARP_EDGES= 3.0 +SENS_REMOVE_SHARP= NO + +% -------------------------- MULTIGRID PARAMETERS -----------------------------% +% +MGLEVEL= 3 +MGCYCLE= V_CYCLE +MG_PRE_SMOOTH= ( 2, 2, 2, 2 ) +MG_POST_SMOOTH= ( 2, 2, 2, 2) +MG_CORRECTION_SMOOTH= ( 1, 1, 1, 1 ) +MG_DAMP_RESTRICTION= 0.5 +MG_DAMP_PROLONGATION= 0.5 +MG_IMPLICIT_LINES= YES + +% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% +% +CONV_NUM_METHOD_FLOW= ROE +MUSCL_FLOW= YES +SLOPE_LIMITER_FLOW= NONE +JST_SENSOR_COEFF= ( 0.5, 0.02 ) +TIME_DISCRE_FLOW= EULER_IMPLICIT + +% -------------------- TURBULENT NUMERICAL METHOD DEFINITION ------------------% +% +CONV_NUM_METHOD_TURB= SCALAR_UPWIND +MUSCL_TURB= NO +SLOPE_LIMITER_TURB= VENKATAKRISHNAN +TIME_DISCRE_TURB= EULER_IMPLICIT + +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +CONV_RESIDUAL_MINVAL= -12 +CONV_STARTITER= 10 +CONV_CAUCHY_ELEMS= 100 +CONV_CAUCHY_EPS= 1E-6 + +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +MESH_FILENAME= mesh_flatplate_turb_137x97.su2 +MESH_FORMAT= SU2 +MESH_OUT_FILENAME= mesh_out +SOLUTION_FILENAME= solution_flow +SOLUTION_ADJ_FILENAME= solution_adj +TABULAR_FORMAT= CSV +CONV_FILENAME= history_tmr +SCREEN_OUTPUT= (INNER_ITER, RMS_DENSITY, RMS_ENERGY, RMS_TKE, RMS_DISSIPATION, LIFT, DRAG) +RESTART_FILENAME= restart_flow +RESTART_ADJ_FILENAME= restart_adj +VOLUME_FILENAME= flow +VOLUME_ADJ_FILENAME= adjoint +GRAD_OBJFUNC_FILENAME= of_grad +SURFACE_FILENAME= surface_flow +SURFACE_ADJ_FILENAME= surface_adjoint +OUTPUT_WRT_FREQ= 1000