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 bcd2e9c06328..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); @@ -1200,13 +1208,10 @@ inline SST_ParsedOptions ParseSSTOptions(const SST_OPTIONS *SST_Options, unsigne SSTParsedOptions.production = SST_OPTIONS::UQ; } - // Parse compressibility options + // Parse compressibility options, stored in compWilcox and compSarkar so that they can be combined with a + // production modifier if (sst_compWilcox && sst_compSarkar) { SU2_MPI::Error("Please select only one compressibility correction (COMPRESSIBILITY-WILCOX or COMPRESSIBILITY-SARKAR).", CURRENT_FUNCTION); - } else if (sst_compWilcox) { - SSTParsedOptions.production = SST_OPTIONS::COMP_Wilcox; - } else if (sst_compSarkar) { - SSTParsedOptions.production = SST_OPTIONS::COMP_Sarkar; } SSTParsedOptions.sust = sst_sust; @@ -1216,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 fd4207b0c5fe..35af64b17b8a 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*/ @@ -4974,6 +4981,14 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i Kind_Solver == MAIN_SOLVER::FEM_EULER) Kind_Turb_Model = TURB_MODEL::NONE; + /*--- SST has no engine or actuator-disk boundary conditions: the faces would get no turbulence flux at all. + Checked after the turbulence model of Euler zones is cleared (multizone). ---*/ + if (Kind_Turb_Model == TURB_MODEL::SST && + (nMarker_EngineInflow + nMarker_EngineExhaust + nMarker_ActDiskInlet + nMarker_ActDiskOutlet) > 0) { + SU2_MPI::Error("MARKER_ENGINE_INFLOW, MARKER_ENGINE_EXHAUST and MARKER_ACTDISK are not supported with the SST model.", + CURRENT_FUNCTION); + } + Kappa_2nd_Flow = jst_coeff[0]; Kappa_4th_Flow = jst_coeff[1]; Kappa_2nd_AdjFlow = jst_adj_coeff[0]; @@ -6676,16 +6691,12 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { cout << "\nperturbing the Reynold's Stress Matrix towards " << eig_val_comp << " component turbulence"; if (uq_permute) cout << " (permuting eigenvectors)"; break; - case SST_OPTIONS::COMP_Wilcox: - cout << " with compressibility correction of Wilcox"; - break; - case SST_OPTIONS::COMP_Sarkar: - cout << " with compressibility correction of Sarkar"; - break; default: cout << " with no production modification"; break; } + if (sstParsedOptions.compWilcox) cout << ", with compressibility correction of Wilcox"; + if (sstParsedOptions.compSarkar) cout << ", with compressibility correction of Sarkar"; if (sstParsedOptions.dll){ cout << "\nusing non dimensional lower limits relative to infinity values clipping by Coefficients:" ; diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index ec511a3fa1db..b26e911c3542 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -1146,10 +1146,11 @@ class CNumerics { * \param[in] val_normal - Normal vector, the norm of the vector is the area of the face. * \param[in] val_scale - Scale of the projection. * \param[out] val_Proj_Jac_tensor - Pointer to the projected inviscid Jacobian. + * \param[in] val_tke - Turbulent kinetic energy contained in the total energy (SST), held fixed. */ void GetInviscidProjJac(const su2double *val_velocity, const su2double *val_energy, const su2double *val_normal, su2double val_scale, - su2double **val_Proj_Jac_tensor) const; + su2double **val_Proj_Jac_tensor, su2double val_tke = 0.0) const; /*! * \brief Compute the projection of the inviscid Jacobian matrices (incompressible). @@ -1222,11 +1223,12 @@ class CNumerics { * \param[in] val_normal - Normal vector, the norm of the vector is the area of the face. * \param[in] val_scale - Scale of the projection. * \param[out] val_Proj_Jac_tensor - Pointer to the projected inviscid Jacobian. + * \param[in] val_tke - Turbulent kinetic energy contained in the total energy (SST), held fixed. */ void GetInviscidProjJac(const su2double *val_velocity, const su2double *val_enthalphy, const su2double *val_chi, const su2double *val_kappa, const su2double *val_normal, su2double val_scale, - su2double **val_Proj_Jac_tensor) const; + su2double **val_Proj_Jac_tensor, su2double val_tke = 0.0) const; /*! * \brief Mapping between primitives variables P and conservatives variables C. @@ -1253,7 +1255,7 @@ class CNumerics { void GetPMatrix(const su2double *val_density, const su2double *val_velocity, const su2double *val_soundspeed, const su2double *val_enthalpy, const su2double *val_chi, const su2double *val_kappa, - const su2double *val_normal, su2double **val_p_tensor) const; + const su2double *val_normal, su2double **val_p_tensor, su2double val_tke = 0.0) const; /*! * \brief Computation of the matrix P, this matrix diagonalize the conservative Jacobians in @@ -1267,12 +1269,13 @@ class CNumerics { template void GetPMatrix(const su2double& density, const su2double* velocity, const su2double& soundspeed, const su2double* normal, - Matrix& p_tensor) const { + Matrix& p_tensor, su2double tke = 0.0) const { + /*--- With SST the total energy contains k (tke): it is added to the kinetic energy of the energy row. ---*/ const su2double rhooc = density / soundspeed; const su2double rhoxc = density * soundspeed; if (nDim == 2) { - const su2double ke = 0.5 * GeometryToolbox::SquaredNorm(2, velocity); + const su2double ke = 0.5 * GeometryToolbox::SquaredNorm(2, velocity) + tke; const su2double projvel = GeometryToolbox::DotProduct(2, velocity, normal); p_tensor[0][0] = 1.0; @@ -1295,7 +1298,7 @@ class CNumerics { p_tensor[3][2] = 0.5 * (ke * rhooc + density * projvel + rhoxc / Gamma_Minus_One); p_tensor[3][3] = 0.5 * (ke * rhooc - density * projvel + rhoxc / Gamma_Minus_One); } else { - const su2double ke = 0.5 * GeometryToolbox::SquaredNorm(3, velocity); + const su2double ke = 0.5 * GeometryToolbox::SquaredNorm(3, velocity) + tke; const su2double projvel = GeometryToolbox::DotProduct(3, velocity, normal); p_tensor[0][0] = normal[0]; @@ -1430,7 +1433,7 @@ class CNumerics { void GetPMatrix_inv(su2double **val_invp_tensor, const su2double *val_density, const su2double *val_velocity, const su2double *val_soundspeed, const su2double *val_chi, const su2double *val_kappa, - const su2double *val_normal) const; + const su2double *val_normal, su2double val_tke = 0.0) const; /*! * \brief Computation of the matrix P^{-1}, this matrix diagonalize the conservative Jacobians @@ -1444,7 +1447,8 @@ class CNumerics { template void GetPMatrix_inv(const su2double& density, const su2double* velocity, const su2double& soundspeed, const su2double* normal, - Matrix& inv_p_tensor) const { + Matrix& inv_p_tensor, su2double tke = 0.0) const { + /*--- With SST the total energy contains k (tke), dp/drho = (gamma-1) (|u|^2/2 - k). ---*/ const su2double rhoxc = density * soundspeed; const su2double c2 = pow(soundspeed, 2); const su2double gm1 = Gamma_Minus_One; @@ -1455,7 +1459,7 @@ class CNumerics { if (nDim == 3) { const su2double k2orho = normal[2] / density; - const su2double ke = 0.5 * GeometryToolbox::SquaredNorm(3, velocity); + const su2double ke = 0.5 * GeometryToolbox::SquaredNorm(3, velocity) - tke; const su2double projvel_o_rho = GeometryToolbox::DotProduct(3, velocity, normal) / density; inv_p_tensor[0][0] = normal[0] + k1orho * velocity[2] - k2orho * velocity[1] - normal[0] * gm1_o_c2 * ke; @@ -1488,7 +1492,7 @@ class CNumerics { inv_p_tensor[4][3] = -k2orho - gm1_o_rhoxc * velocity[2]; inv_p_tensor[4][4] = gm1_o_rhoxc; } else { - const su2double ke = 0.5 * GeometryToolbox::SquaredNorm(2, velocity); + const su2double ke = 0.5 * GeometryToolbox::SquaredNorm(2, velocity) - tke; const su2double projvel_o_rho = GeometryToolbox::DotProduct(2, velocity, normal) / density; inv_p_tensor[0][0] = 1 - gm1_o_c2 * ke; diff --git a/SU2_CFD/include/numerics/flow/convection/roe.hpp b/SU2_CFD/include/numerics/flow/convection/roe.hpp index 01c0874f8e98..8664c60e49e1 100644 --- a/SU2_CFD/include/numerics/flow/convection/roe.hpp +++ b/SU2_CFD/include/numerics/flow/convection/roe.hpp @@ -45,6 +45,7 @@ class CUpwRoeBase_Flow : public CNumerics { su2double *ProjFlux_j = nullptr, *Conservatives_j = nullptr; su2double **P_Tensor = nullptr, **invP_Tensor = nullptr; su2double RoeDensity, RoeEnthalpy, RoeSoundSpeed, ProjVelocity, RoeSoundSpeed2, kappa; + su2double RoeTke = 0.0; /*!< \brief Roe-averaged k (SST, contained in the total energy). */ su2double* Flux = nullptr; /*!< \brief The flux accross the face. */ su2double** Jacobian_i = nullptr; /*!< \brief The Jacobian w.r.t. point i after computation. */ diff --git a/SU2_CFD/include/numerics/flow/flow_sources.hpp b/SU2_CFD/include/numerics/flow/flow_sources.hpp index 83c3811326cd..cf06eb753d8a 100644 --- a/SU2_CFD/include/numerics/flow/flow_sources.hpp +++ b/SU2_CFD/include/numerics/flow/flow_sources.hpp @@ -72,7 +72,9 @@ class CSourceBase_Flow : public CNumerics { */ class CSourceAxisymmetric_Flow : public CSourceBase_Flow { protected: - bool implicit, viscous, rans; + 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/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp index 9f72b6e37cc4..405be31e7f86 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp @@ -938,39 +938,50 @@ class CSourcePieceWise_TurbSST final : public CNumerics { P_Base = sqrt(StrainMag_i*VorticityMag); break; - case SST_OPTIONS::COMP_Wilcox: - P_Base = StrainMag_i; - if (Mt >= 0.25) { - zetaFMt = 2.0 * (Mt * Mt - 0.25 * 0.25); - } - break; - - case SST_OPTIONS::COMP_Sarkar: - P_Base = StrainMag_i; - if (Mt >= 0.25) { - zetaFMt = 0.5 * (Mt * Mt); - } - break; - default: /*--- Base production term for SST-1994 and SST-2003 ---*/ P_Base = StrainMag_i; break; } + /*--- Compressibility corrections, independent of the production modifier. ---*/ + if (sstParsedOptions.compWilcox && Mt >= 0.25) { + zetaFMt = 2.0 * (Mt * Mt - 0.25 * 0.25); + } + if (sstParsedOptions.compSarkar && Mt >= 0.25) { + zetaFMt = 0.5 * (Mt * Mt); + } + /*--- 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; } @@ -988,7 +999,7 @@ class CSourcePieceWise_TurbSST final : public CNumerics { pw = max(pw, sust_w); } - if (sstParsedOptions.production == SST_OPTIONS::COMP_Sarkar) { + if (sstParsedOptions.compSarkar) { const su2double Dilatation_Sarkar = -0.15 * pk * Mt + 0.2 * beta_star * (1.0 +zetaFMt) * Density_i * ScalarVar_i[1] * ScalarVar_i[0] * Mt * Mt; pk += Dilatation_Sarkar; } @@ -996,7 +1007,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 +1034,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); diff --git a/SU2_CFD/include/numerics_simd/flow/convection/centered.hpp b/SU2_CFD/include/numerics_simd/flow/convection/centered.hpp index faaf56912351..193b5ec7624f 100644 --- a/SU2_CFD/include/numerics_simd/flow/convection/centered.hpp +++ b/SU2_CFD/include/numerics_simd/flow/convection/centered.hpp @@ -52,6 +52,7 @@ class CCenteredBase : public Base { const su2double fixFactor; const bool dynamicGrid; const su2double stretchParam = 0.3; + const CVariable* const tkeVars; /*!< \brief SST turbulence variables (k is part of the total energy), else nullptr. */ /*! * \brief Constructor, store some constants and forward args to base. @@ -60,7 +61,8 @@ class CCenteredBase : public Base { CCenteredBase(const CConfig& config, Ts&... args) : Base(config, args...), gamma(config.GetGamma()), fixFactor(config.GetCent_Jac_Fix_Factor()), - dynamicGrid(config.GetDynamic_Grid()) { + dynamicGrid(config.GetDynamic_Grid()), + tkeVars(config.GetKind_Turb_Model() == TURB_MODEL::SST ? findTurbVars(args...) : nullptr) { } /*! @@ -140,8 +142,14 @@ class CCenteredBase : public Base { MatrixDbl jac_i, jac_j; if (implicit) { - jac_i = inviscidProjJac(gamma, V.i.velocity(), U.i.energy(), normal, 0.5); - jac_j = inviscidProjJac(gamma, V.j.velocity(), U.j.energy(), normal, 0.5); + /*--- With SST the total energy contains k, held fixed in the Jacobian. ---*/ + Double tke_i = 0.0, tke_j = 0.0; + if (tkeVars) { + tke_i = gatherVariables(iPoint, tkeVars->GetSolution()); + tke_j = gatherVariables(jPoint, tkeVars->GetSolution()); + } + jac_i = inviscidProjJac(gamma, V.i.velocity(), U.i.energy(), normal, 0.5, tke_i); + jac_j = inviscidProjJac(gamma, V.j.velocity(), U.j.energy(), normal, 0.5, tke_j); } /*--- Grid motion. ---*/ @@ -283,6 +291,7 @@ class CJSTmatScheme : public CCenteredBase,Decorator> { using Base::nVar; using Base::gamma; using Base::fixFactor; + using Base::tkeVars; const su2double kappa2; const su2double kappa4; const su2double entropyFix; @@ -350,11 +359,24 @@ class CJSTmatScheme : public CCenteredBase,Decorator> { const auto unitProjVel = dot(avgV.velocity(), unitNormal); + /*--- With SST the total energy contains k, averaged between the cells. ---*/ + Double tke_i = 0.0, tke_j = 0.0; + if (tkeVars) { + tke_i = gatherVariables(iPoint, tkeVars->GetSolution()); + tke_j = gatherVariables(jPoint, tkeVars->GetSolution()); + } + const Double avgTke = 0.5 * (tke_i + tke_j); + + /*--- Delta(rho k) = avg(rho) Delta k + avg(k) Delta rho. The part with Delta k is not a pressure jump: remove it + from the 2nd order dissipation before the projection and advect it with the contact wave. The 4th order term + keeps it (the Laplacian of k is not available). ---*/ + const Double rhoDeltaTke = eps2 * avgV.density() * (tke_i - tke_j); + scalarDissip(nVar-1) -= rhoDeltaTke; auto pMat = pMatrix(gamma, avgV.density(), avgV.velocity(), - unitProjVel, avgV.speedSound(), unitNormal); + unitProjVel, avgV.speedSound(), unitNormal, avgTke); auto pMatInv = pMatrixInv(gamma, avgV.density(), avgV.velocity(), - unitProjVel, avgV.speedSound(), unitNormal); + unitProjVel, avgV.speedSound(), unitNormal, avgTke); /*--- Compute limited absolute eigenvalues (times area). ---*/ @@ -394,6 +416,7 @@ class CJSTmatScheme : public CCenteredBase,Decorator> { } } } + flux(nVar-1) += lambda(0) * rhoDeltaTke; } }; diff --git a/SU2_CFD/include/numerics_simd/flow/convection/common.hpp b/SU2_CFD/include/numerics_simd/flow/convection/common.hpp index 800167d4b0f5..2d09e62d6b3b 100644 --- a/SU2_CFD/include/numerics_simd/flow/convection/common.hpp +++ b/SU2_CFD/include/numerics_simd/flow/convection/common.hpp @@ -134,9 +134,11 @@ FORCEINLINE CPair reconstructPrimitives(const Int& iEdge, */ template FORCEINLINE MatrixDbl pMatrix(const Double& gamma, const Double& density, const RandomAccessIterator& velocity, - const Double& projVel, const Double& speedSound, const VectorDbl& normal) { + const Double& projVel, const Double& speedSound, const VectorDbl& normal, + const Double& tke = Double(0.0)) { MatrixDbl pMat; - const Double vel2 = 0.5*squaredNorm(velocity); + /*--- With SST the total energy contains k (tke). ---*/ + const Double vel2 = 0.5*squaredNorm(velocity) + tke; if (nDim == 2) { pMat(0,0) = 1.0; @@ -197,11 +199,13 @@ FORCEINLINE MatrixDbl pMatrix(const Double& gamma, const Double& density template FORCEINLINE MatrixDbl pMatrixInv(const Double& gamma, const Double& density, const RandomAccessIterator& velocity, const Double& projVel, - const Double& speedSound, const VectorDbl& normal) { + const Double& speedSound, const VectorDbl& normal, + const Double& tke = Double(0.0)) { MatrixDbl pMatInv; const Double c2 = pow(speedSound,2); - const Double vel2 = 0.5*squaredNorm(velocity); + /*--- With SST the total energy contains k (tke), dp/drho = (gamma-1) (|u|^2/2 - k). ---*/ + const Double vel2 = 0.5*squaredNorm(velocity) - tke; const Double oneOnRho = 1 / density; if (nDim == 2) { @@ -273,19 +277,30 @@ FORCEINLINE VectorDbl inviscidProjFlux(const PrimVarType& V, return flux; } +/*! + * \brief Turbulence variables among the constructor arguments of a scheme (nullptr if there are none). + */ +inline const CVariable* findTurbVars() { return nullptr; } +template +const CVariable* findTurbVars(T& first, Ts&... rest) { + if constexpr (std::is_convertible::value) return first; + else return findTurbVars(rest...); +} + /*! * \brief Jacobian of the convective flux (compressible flow, ideal gas). + * \note With SST the total energy contains k (tke), which is held fixed. */ template FORCEINLINE MatrixDbl inviscidProjJac(const Double& gamma, RandomAccessIterator velocity, const Double& energy, const VectorDbl& normal, - const Double& scale) { + const Double& scale, const Double& tke = Double(0.0)) { MatrixDbl jac; Double projVel = dot(velocity, normal); Double gamma_m_1 = gamma-1; - Double phi = 0.5*gamma_m_1*squaredNorm(velocity); - Double a1 = gamma*energy - phi; + Double phi = gamma_m_1*(0.5*squaredNorm(velocity) - tke); // dp/drho, k held fixed + Double a1 = gamma*energy - gamma_m_1*(0.5*squaredNorm(velocity) + tke); // total enthalpy jac(0,0) = 0.0; for (size_t iDim = 0; iDim < nDim; ++iDim) { diff --git a/SU2_CFD/include/numerics_simd/flow/convection/upwind.hpp b/SU2_CFD/include/numerics_simd/flow/convection/upwind.hpp index 38ec02c902b5..f8f0e43e983a 100644 --- a/SU2_CFD/include/numerics_simd/flow/convection/upwind.hpp +++ b/SU2_CFD/include/numerics_simd/flow/convection/upwind.hpp @@ -61,6 +61,7 @@ class CUpwindBase : public Base { const bool muscl; const su2double umusclKappa; const LIMITER typeLimiter; + const CVariable* const tkeVars; /*!< \brief SST turbulence variables (k is part of the total energy), else nullptr. */ /*! * \brief Constructor, store some constants and forward args to base. @@ -73,7 +74,8 @@ class CUpwindBase : public Base { dynamicGrid(config.GetDynamic_Grid()), muscl(finestGrid && config.GetMUSCL_Flow()), umusclKappa(config.GetMUSCL_Kappa_Flow()), - typeLimiter(config.GetKind_SlopeLimit_Flow()) { + typeLimiter(config.GetKind_SlopeLimit_Flow()), + tkeVars(config.GetKind_Turb_Model() == TURB_MODEL::SST ? findTurbVars(args...) : nullptr) { } public: @@ -178,6 +180,7 @@ class CRoeScheme : public CUpwindBase, Decorator> { using Base::gamma; using Base::gasConst; using Base::dynamicGrid; + using Base::tkeVars; public: /*! @@ -190,6 +193,7 @@ class CRoeScheme : public CUpwindBase, Decorator> { typeDissip(static_cast(config.GetKind_RoeLowDiss())) { } + /*! * \brief Updates flux and Jacobians with standard Roe dissipation. * \note "Ts" is here just in case other schemes in the family need extra args. @@ -210,9 +214,14 @@ class CRoeScheme : public CUpwindBase, Decorator> { const CEulerVariable& solution, const CGeometry& geometry, Ts&...) const { - /*--- Roe averaged variables. ---*/ + /*--- Roe averaged variables. With SST the total enthalpy of the cells contains k. ---*/ - auto roeAvg = roeAveragedVariables(gamma, V, unitNormal); + Double tke_i = 0.0, tke_j = 0.0; + if (tkeVars) { + tke_i = gatherVariables(iPoint, tkeVars->GetSolution()); + tke_j = gatherVariables(jPoint, tkeVars->GetSolution()); + } + auto roeAvg = roeAveragedVariables(gamma, V, unitNormal, tke_i, tke_j); /*--- Grid motion. ---*/ @@ -249,8 +258,8 @@ class CRoeScheme : public CUpwindBase, Decorator> { flux(iVar) = 0.5 * (flux_i(iVar) + flux_j(iVar)); } if (implicit) { - jac_i = inviscidProjJac(gamma, V.i.velocity(), U.i.energy(), normal, kappa); - jac_j = inviscidProjJac(gamma, V.j.velocity(), U.j.energy(), normal, kappa); + jac_i = inviscidProjJac(gamma, V.i.velocity(), U.i.energy(), normal, kappa, tke_i); + jac_j = inviscidProjJac(gamma, V.j.velocity(), U.j.energy(), normal, kappa, tke_j); } /*--- Correct for grid motion. ---*/ @@ -270,12 +279,12 @@ class CRoeScheme : public CUpwindBase, Decorator> { /*--- P tensor. ---*/ auto pMat = pMatrix(gamma, roeAvg.density, roeAvg.velocity, - roeAvg.projVel, roeAvg.speedSound, unitNormal); + roeAvg.projVel, roeAvg.speedSound, unitNormal, roeAvg.tke); /*--- Inverse P tensor. ---*/ auto pMatInv = pMatrixInv(gamma, roeAvg.density, roeAvg.velocity, - roeAvg.projVel, roeAvg.speedSound, unitNormal); + roeAvg.projVel, roeAvg.speedSound, unitNormal, roeAvg.tke); /*--- Diference between conservative variables at jPoint and iPoint. ---*/ @@ -284,6 +293,11 @@ class CRoeScheme : public CUpwindBase, Decorator> { deltaU(iVar) = U.j.all(iVar) - U.i.all(iVar); } + /*--- With SST rho*E contains rho*k, and Delta(rho k) = rho_roe Delta k + k_roe Delta rho. The part with Delta k + is not a pressure jump: remove it before the projection and advect it with the contact wave. ---*/ + const Double rhoDeltaTke = roeAvg.density * (tke_j - tke_i); + deltaU(nVar-1) -= rhoDeltaTke; + /*--- Dissipation terms. ---*/ Double dissipation = roeDissipation(iPoint, jPoint, typeDissip, solution); @@ -309,6 +323,7 @@ class CRoeScheme : public CUpwindBase, Decorator> { } } } + flux(nVar-1) -= (1-kappa) * area * dissipation * lambda(0) * rhoDeltaTke; } }; @@ -330,6 +345,7 @@ class CMSWScheme : public CUpwindBase, Decorator> { using Base::gamma; using Base::gasConst; using Base::dynamicGrid; + using Base::tkeVars; public: /*! @@ -397,10 +413,16 @@ class CMSWScheme : public CUpwindBase, Decorator> { lambda(nDim) = fmax(projVel_i + soundSpeed_i, 0); lambda(nDim+1) = fmax(projVel_i - soundSpeed_i, 0); + /*--- With SST the total energy contains k, taken from the cells. ---*/ + Double tke_i = 0.0, tke_j = 0.0; + if (tkeVars) { + tke_i = gatherVariables(iPoint, tkeVars->GetSolution()); + tke_j = gatherVariables(jPoint, tkeVars->GetSolution()); + } auto pMat = pMatrix(gamma, Vweighted.i.density(), Vweighted.i.velocity(), - projVel_i, soundSpeed_i, unitNormal); + projVel_i, soundSpeed_i, unitNormal, tke_i); auto pMatInv = pMatrixInv(gamma, Vweighted.i.density(), Vweighted.i.velocity(), - projVel_i, soundSpeed_i, unitNormal); + projVel_i, soundSpeed_i, unitNormal, tke_i); auto updateFlux = [&](const auto& u, auto& jac) { for (size_t iVar = 0; iVar < nVar; ++iVar) { @@ -429,9 +451,9 @@ class CMSWScheme : public CUpwindBase, Decorator> { lambda(nDim+1) = fmin(projVel_j - soundSpeed_j, 0); pMat = pMatrix(gamma, Vweighted.j.density(), Vweighted.j.velocity(), - projVel_j, soundSpeed_j, unitNormal); + projVel_j, soundSpeed_j, unitNormal, tke_j); pMatInv = pMatrixInv(gamma, Vweighted.j.density(), Vweighted.j.velocity(), - projVel_j, soundSpeed_j, unitNormal); + projVel_j, soundSpeed_j, unitNormal, tke_j); updateFlux(U.j.all, jac_j); } }; 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/include/numerics_simd/flow/variables.hpp b/SU2_CFD/include/numerics_simd/flow/variables.hpp index 4406a06abf7c..abeabf049ed1 100644 --- a/SU2_CFD/include/numerics_simd/flow/variables.hpp +++ b/SU2_CFD/include/numerics_simd/flow/variables.hpp @@ -107,15 +107,19 @@ struct CRoeVariables { Double enthalpy; Double speedSound; Double projVel; + Double tke; /*!< \brief Roe-averaged k (SST, part of the total enthalpy). */ }; /*! * \brief Compute Roe-averaged variables from pair of primitive variables. + * \note With SST the total enthalpy contains k (tke_i, tke_j), which is not part of the speed of sound. */ template FORCEINLINE CRoeVariables roeAveragedVariables(const Double& gamma, const CPair& V, - const VectorDbl& normal) { + const VectorDbl& normal, + const Double& tke_i = Double(0.0), + const Double& tke_j = Double(0.0)) { CRoeVariables roeAvg; Double R = sqrt(V.j.density() / V.i.density()); Double D = 1 / (R+1); @@ -124,7 +128,8 @@ FORCEINLINE CRoeVariables roeAveragedVariables(const Double& gamma, roeAvg.velocity(iDim) = (R*V.j.velocity(iDim) + V.i.velocity(iDim)) * D; } roeAvg.enthalpy = (R*V.j.enthalpy() + V.i.enthalpy()) * D; - roeAvg.speedSound = sqrt((gamma-1) * (roeAvg.enthalpy - 0.5*squaredNorm(roeAvg.velocity))); + roeAvg.tke = (R*tke_j + tke_i) * D; + roeAvg.speedSound = sqrt((gamma-1) * (roeAvg.enthalpy - 0.5*squaredNorm(roeAvg.velocity) - roeAvg.tke)); roeAvg.projVel = dot(roeAvg.velocity, normal); return roeAvg; } diff --git a/SU2_CFD/include/solvers/CEulerSolver.hpp b/SU2_CFD/include/solvers/CEulerSolver.hpp index 345be6d3d3a4..ed4b5db05748 100644 --- a/SU2_CFD/include/solvers/CEulerSolver.hpp +++ b/SU2_CFD/include/solvers/CEulerSolver.hpp @@ -419,11 +419,13 @@ class CEulerSolver : public CFVMFlowSolverBase CUpwAUSMPLUS_SLAU_Base_Flow::ComputeResidual(const CCo AD::SetPreaccIn(Normal, nDim); AD::SetPreaccIn(V_i, nDim+4); AD::SetPreaccIn(V_j, nDim+4); + AD::SetPreaccIn(turb_ke_i); AD::SetPreaccIn(turb_ke_j); // SST: k in the total energy /*--- Variables for the general form and primitives for mass flux and pressure calculation. ---*/ /*--- F_{1/2} = ||A|| ( 0.5 * mdot * (psi_i+psi_j) - 0.5 * |mdot| * (psi_i-psi_j) + N * pf ) ---*/ @@ -413,8 +416,9 @@ void CUpwAUSMPLUSUP_Flow::ComputeMassAndPressureFluxes(const CConfig* config, su /*--- Compute interface speed of sound (aF) ---*/ - su2double astarL = sqrt(2.0*(Gamma-1.0)/(Gamma+1.0)*Enthalpy_i); - su2double astarR = sqrt(2.0*(Gamma-1.0)/(Gamma+1.0)*Enthalpy_j); + /*--- With SST the total enthalpy contains k, which is not part of the critical speed of sound. ---*/ + su2double astarL = sqrt(2.0*(Gamma-1.0)/(Gamma+1.0)*(Enthalpy_i-turb_ke_i)); + su2double astarR = sqrt(2.0*(Gamma-1.0)/(Gamma+1.0)*(Enthalpy_j-turb_ke_j)); su2double ahatL = astarL*astarL/max(astarL, ProjVelocity_i); su2double ahatR = astarR*astarR/max(astarR,-ProjVelocity_j); @@ -584,7 +588,7 @@ void CUpwAUSMPLUSUP_Flow::ComputeMassAndPressureFluxes(const CConfig* config, su astar_b *= 2.0*tmp; Vn_i_b -= tmp*tmp * aF_b; } - H_i_b = sqrt(0.5*(Gamma-1.0)/((Gamma+1.0)*Enthalpy_i)) * astar_b; + H_i_b = sqrt(0.5*(Gamma-1.0)/((Gamma+1.0)*(Enthalpy_i-turb_ke_i))) * astar_b; H_j_b = 0.0; } else { @@ -593,7 +597,7 @@ void CUpwAUSMPLUSUP_Flow::ComputeMassAndPressureFluxes(const CConfig* config, su astar_b *= 2.0*tmp; Vn_j_b += tmp*tmp * aF_b; } - H_j_b = sqrt(0.5*(Gamma-1.0)/((Gamma+1.0)*Enthalpy_j)) * astar_b; + H_j_b = sqrt(0.5*(Gamma-1.0)/((Gamma+1.0)*(Enthalpy_j-turb_ke_j))) * astar_b; H_i_b = 0.0; } @@ -639,8 +643,9 @@ void CUpwAUSMPLUSUP2_Flow::ComputeMassAndPressureFluxes(const CConfig* config, s /*--- Compute interface speed of sound (aF) ---*/ - su2double astarL = sqrt(2.0*(Gamma-1.0)/(Gamma+1.0)*Enthalpy_i); - su2double astarR = sqrt(2.0*(Gamma-1.0)/(Gamma+1.0)*Enthalpy_j); + /*--- With SST the total enthalpy contains k, which is not part of the critical speed of sound. ---*/ + su2double astarL = sqrt(2.0*(Gamma-1.0)/(Gamma+1.0)*(Enthalpy_i-turb_ke_i)); + su2double astarR = sqrt(2.0*(Gamma-1.0)/(Gamma+1.0)*(Enthalpy_j-turb_ke_j)); su2double ahatL = astarL*astarL/max(astarL, ProjVelocity_i); su2double ahatR = astarR*astarR/max(astarR,-ProjVelocity_j); @@ -731,11 +736,12 @@ void CUpwSLAU_Flow::ComputeMassAndPressureFluxes(const CConfig* config, su2doubl sq_velj += Velocity_j[iDim]*Velocity_j[iDim]; } + /*--- With SST the total energy contains k, which is not part of the speed of sound. ---*/ su2double Energy_i = Enthalpy_i - Pressure_i/Density_i; - SoundSpeed_i = sqrt(fabs(Gamma*Gamma_Minus_One*(Energy_i-0.5*sq_veli))); + SoundSpeed_i = sqrt(fabs(Gamma*Gamma_Minus_One*(Energy_i-0.5*sq_veli-turb_ke_i))); su2double Energy_j = Enthalpy_j - Pressure_j/Density_j; - SoundSpeed_j = sqrt(fabs(Gamma*Gamma_Minus_One*(Energy_j-0.5*sq_velj))); + SoundSpeed_j = sqrt(fabs(Gamma*Gamma_Minus_One*(Energy_j-0.5*sq_velj-turb_ke_j))); /*--- Compute interface speed of sound (aF), and left/right Mach number ---*/ @@ -842,6 +848,7 @@ CNumerics::ResidualType<> CUpwAUSM_Flow::ComputeResidual(const CConfig* config) AD::SetPreaccIn(Normal, nDim); AD::SetPreaccIn(V_i, nDim+4); AD::SetPreaccIn(V_j, nDim+4); + AD::SetPreaccIn(turb_ke_i); AD::SetPreaccIn(turb_ke_j); // SST: k in the total energy /*--- Face area (norm or the normal vector) ---*/ Area = GeometryToolbox::Norm(nDim, Normal); @@ -860,7 +867,7 @@ CNumerics::ResidualType<> CUpwAUSM_Flow::ComputeResidual(const CConfig* config) Density_i = V_i[nDim+2]; Enthalpy_i = V_i[nDim+3]; Energy_i = Enthalpy_i - Pressure_i/Density_i; - SoundSpeed_i = sqrt(fabs(Gamma*Gamma_Minus_One*(Energy_i-0.5*sq_vel))); + SoundSpeed_i = sqrt(fabs(Gamma*Gamma_Minus_One*(Energy_i-0.5*sq_vel-turb_ke_i))); // SST: E contains k /*--- Primitive variables at point j ---*/ sq_vel = 0.0; @@ -872,7 +879,7 @@ CNumerics::ResidualType<> CUpwAUSM_Flow::ComputeResidual(const CConfig* config) Density_j = V_j[nDim+2]; Enthalpy_j = V_j[nDim+3]; Energy_j = Enthalpy_j - Pressure_j/Density_j; - SoundSpeed_j = sqrt(fabs(Gamma*Gamma_Minus_One*(Energy_j-0.5*sq_vel))); + SoundSpeed_j = sqrt(fabs(Gamma*Gamma_Minus_One*(Energy_j-0.5*sq_vel-turb_ke_j))); // SST: E contains k /*--- Projected velocities ---*/ ProjVelocity_i = 0.0; ProjVelocity_j = 0.0; @@ -925,10 +932,11 @@ CNumerics::ResidualType<> CUpwAUSM_Flow::ComputeResidual(const CConfig* config) sq_vel += RoeVelocity[iDim]*RoeVelocity[iDim]; } RoeEnthalpy = (R*Enthalpy_j+Enthalpy_i)/(R+1); - RoeSoundSpeed = sqrt(fabs((Gamma-1)*(RoeEnthalpy-0.5*sq_vel))); + /*--- With SST the total enthalpy contains k, which is not part of the speed of sound. ---*/ + RoeSoundSpeed = sqrt(fabs((Gamma-1)*(RoeEnthalpy-0.5*sq_vel-(R*turb_ke_j+turb_ke_i)/(R+1)))); /*--- Compute P and Lambda (do it with the Normal) ---*/ - GetPMatrix(RoeDensity, RoeVelocity, RoeSoundSpeed, UnitNormal, P_Tensor); + GetPMatrix(RoeDensity, RoeVelocity, RoeSoundSpeed, UnitNormal, P_Tensor, (R*turb_ke_j+turb_ke_i)/(R+1)); ProjVelocity = 0.0; ProjVelocity_i = 0.0; ProjVelocity_j = 0.0; for (iDim = 0; iDim < nDim; iDim++) { @@ -944,11 +952,11 @@ CNumerics::ResidualType<> CUpwAUSM_Flow::ComputeResidual(const CConfig* config) Lambda[nVar-1] = ProjVelocity - RoeSoundSpeed; /*--- Compute inverse P ---*/ - GetPMatrix_inv(RoeDensity, RoeVelocity, RoeSoundSpeed, UnitNormal, invP_Tensor); + GetPMatrix_inv(RoeDensity, RoeVelocity, RoeSoundSpeed, UnitNormal, invP_Tensor, (R*turb_ke_j+turb_ke_i)/(R+1)); /*--- Jacobias of the inviscid flux, scale = 0.5 because val_residual ~ 0.5*(fc_i+fc_j)*Normal ---*/ - GetInviscidProjJac(Velocity_i, &Energy_i, Normal, 0.5, Jacobian_i); - GetInviscidProjJac(Velocity_j, &Energy_j, Normal, 0.5, Jacobian_j); + GetInviscidProjJac(Velocity_i, &Energy_i, Normal, 0.5, Jacobian_i, turb_ke_i); + GetInviscidProjJac(Velocity_j, &Energy_j, Normal, 0.5, Jacobian_j, turb_ke_j); /*--- Roe's Flux approximation ---*/ for (iVar = 0; iVar < nVar; iVar++) { diff --git a/SU2_CFD/src/numerics/flow/convection/fvs.cpp b/SU2_CFD/src/numerics/flow/convection/fvs.cpp index 7655b793c851..83684cb34b7a 100644 --- a/SU2_CFD/src/numerics/flow/convection/fvs.cpp +++ b/SU2_CFD/src/numerics/flow/convection/fvs.cpp @@ -45,6 +45,8 @@ CNumerics::ResidualType<> CUpwMSW_Flow::ComputeResidual(const CConfig* config) { AD::StartPreacc(); AD::SetPreaccIn(V_i, nDim + 4); AD::SetPreaccIn(V_j, nDim + 4); + AD::SetPreaccIn(turb_ke_i); // SST: k in the total energy + AD::SetPreaccIn(turb_ke_j); AD::SetPreaccIn(Sensor_i, Sensor_j); AD::SetPreaccIn(Normal, nDim); if (dynamic_grid) { @@ -71,8 +73,9 @@ CNumerics::ResidualType<> CUpwMSW_Flow::ComputeResidual(const CConfig* config) { /*--- Recompute the speed of sound because it is not MUSCL-reconstructed. ---*/ const su2double sqvel_i = GeometryToolbox::SquaredNorm(nDim, V_i + 1); const su2double sqvel_j = GeometryToolbox::SquaredNorm(nDim, V_j + 1); - const su2double c_i = sqrt(fmax((Gamma - 1) * (H_i - 0.5 * sqvel_i), EPS)); - const su2double c_j = sqrt(fmax((Gamma - 1) * (H_j - 0.5 * sqvel_j), EPS)); + /*--- With SST the total enthalpy contains k, which is not part of the speed of sound. ---*/ + const su2double c_i = sqrt(fmax((Gamma - 1) * (H_i - 0.5 * sqvel_i - turb_ke_i), EPS)); + const su2double c_j = sqrt(fmax((Gamma - 1) * (H_j - 0.5 * sqvel_j - turb_ke_j), EPS)); /*--- Recompute conservatives ---*/ @@ -109,6 +112,8 @@ CNumerics::ResidualType<> CUpwMSW_Flow::ComputeResidual(const CConfig* config) { } Vst_i[nDim + 4] = onemw * c_i + w * c_j; Vst_j[nDim + 4] = onemw * c_j + w * c_i; + const su2double tkest_i = onemw * turb_ke_i + w * turb_ke_j; + const su2double tkest_j = onemw * turb_ke_j + w * turb_ke_i; su2double Velst_i[MAXNDIM] = {}, Velst_j[MAXNDIM] = {}; su2double ProjVelst_i{}, ProjVelst_j{}; @@ -132,8 +137,8 @@ CNumerics::ResidualType<> CUpwMSW_Flow::ComputeResidual(const CConfig* config) { /*--- Compute projected P, invP, and Lambda ---*/ su2double P_Tensor[MAXNVAR][MAXNVAR], invP_Tensor[MAXNVAR][MAXNVAR]; - GetPMatrix(Vst_i[nDim + 2], Velst_i, Vst_i[nDim + 4], UnitNormal, P_Tensor); - GetPMatrix_inv(Vst_i[nDim + 2], Velst_i, Vst_i[nDim + 4], UnitNormal, invP_Tensor); + GetPMatrix(Vst_i[nDim + 2], Velst_i, Vst_i[nDim + 4], UnitNormal, P_Tensor, tkest_i); + GetPMatrix_inv(Vst_i[nDim + 2], Velst_i, Vst_i[nDim + 4], UnitNormal, invP_Tensor, tkest_i); /*--- Projected flux (f+) at i ---*/ @@ -167,8 +172,8 @@ CNumerics::ResidualType<> CUpwMSW_Flow::ComputeResidual(const CConfig* config) { /*--- Compute projected P, invP, and Lambda ---*/ - GetPMatrix(Vst_j[nDim + 2], Velst_j, Vst_j[nDim + 4], UnitNormal, P_Tensor); - GetPMatrix_inv(Vst_j[nDim + 2], Velst_j, Vst_j[nDim + 4], UnitNormal, invP_Tensor); + GetPMatrix(Vst_j[nDim + 2], Velst_j, Vst_j[nDim + 4], UnitNormal, P_Tensor, tkest_j); + GetPMatrix_inv(Vst_j[nDim + 2], Velst_j, Vst_j[nDim + 4], UnitNormal, invP_Tensor, tkest_j); /*--- Projected flux (f-) ---*/ diff --git a/SU2_CFD/src/numerics/flow/convection/hllc.cpp b/SU2_CFD/src/numerics/flow/convection/hllc.cpp index ef6a2a5a41e9..4bcacbd92d9a 100644 --- a/SU2_CFD/src/numerics/flow/convection/hllc.cpp +++ b/SU2_CFD/src/numerics/flow/convection/hllc.cpp @@ -120,8 +120,9 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config) Energy_i = Enthalpy_i - Pressure_i / Density_i; Energy_j = Enthalpy_j - Pressure_j / Density_j; - SoundSpeed_i = sqrt((Enthalpy_i - 0.5 * sq_vel_i) * Gamma_Minus_One); - SoundSpeed_j = sqrt((Enthalpy_j - 0.5 * sq_vel_j) * Gamma_Minus_One); + /*--- With SST the total enthalpy contains k, which is not part of the speed of sound. ---*/ + SoundSpeed_i = sqrt((Enthalpy_i - 0.5 * sq_vel_i - turb_ke_i) * Gamma_Minus_One); + SoundSpeed_j = sqrt((Enthalpy_j - 0.5 * sq_vel_j - turb_ke_j) * Gamma_Minus_One); /*--- Projected velocities ---*/ @@ -161,7 +162,8 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config) /*--- Roe-averaged speed of sound ---*/ - RoeSoundSpeed = sqrt( Gamma_Minus_One * ( RoeEnthalpy - 0.5 * sq_velRoe ) ) - ProjInterfaceVel; + const su2double RoeTke = ( sqrt(Density_j) * turb_ke_j + sqrt(Density_i) * turb_ke_i ) / Rrho; + RoeSoundSpeed = sqrt( Gamma_Minus_One * ( RoeEnthalpy - 0.5 * sq_velRoe - RoeTke ) ) - ProjInterfaceVel; /*--- Speed of sound at L and R ---*/ @@ -254,7 +256,7 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config) for (jVar = 0; jVar < nVar; jVar++) Jacobian_j[iVar][jVar] = 0; - GetInviscidProjJac(Velocity_i, &Energy_i, UnitNormal, 1.0, Jacobian_i); + GetInviscidProjJac(Velocity_i, &Energy_i, UnitNormal, 1.0, Jacobian_i, turb_ke_i); } else { @@ -270,7 +272,7 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config) /*--- Computing pressure derivatives d/dU_L (PI) ---*/ - dPI_dU[0] = 0.5 * Gamma_Minus_One * sq_vel_i; + dPI_dU[0] = Gamma_Minus_One * (0.5 * sq_vel_i - turb_ke_i); // k (SST) held fixed for (iDim = 0; iDim < nDim; iDim++) dPI_dU[iDim+1] = - Gamma_Minus_One * Velocity_i[iDim]; dPI_dU[nVar-1] = Gamma_Minus_One; @@ -389,7 +391,7 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config) for (jVar = 0; jVar < nVar; jVar++) Jacobian_i[iVar][jVar] = 0; - GetInviscidProjJac(Velocity_j, &Energy_j, UnitNormal, 1.0, Jacobian_j); + GetInviscidProjJac(Velocity_j, &Energy_j, UnitNormal, 1.0, Jacobian_j, turb_ke_j); } else { @@ -448,7 +450,7 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config) /*--- Computing pressure derivatives d/dU_R (PI) ---*/ - dPI_dU[0] = 0.5 * Gamma_Minus_One * sq_vel_j; + dPI_dU[0] = Gamma_Minus_One * (0.5 * sq_vel_j - turb_ke_j); // k (SST) held fixed for (iDim = 0; iDim < nDim; iDim++) dPI_dU[iDim+1] = - Gamma_Minus_One * Velocity_j[iDim]; dPI_dU[nVar-1] = Gamma_Minus_One; @@ -623,8 +625,9 @@ CNumerics::ResidualType<> CUpwGeneralHLLC_Flow::ComputeResidual(const CConfig* c } Energy_i = Enthalpy_i - Pressure_i / Density_i; - StaticEnthalpy_i = Enthalpy_i - 0.5 * sq_vel_i; - StaticEnergy_i = Energy_i - 0.5 * sq_vel_i; + /*--- With SST the total energy and enthalpy contain k. ---*/ + StaticEnthalpy_i = Enthalpy_i - 0.5 * sq_vel_i - turb_ke_i; + StaticEnergy_i = Energy_i - 0.5 * sq_vel_i - turb_ke_i; Kappa_i = S_i[1] / Density_i; Chi_i = S_i[0] - Kappa_i * StaticEnergy_i; @@ -632,8 +635,8 @@ CNumerics::ResidualType<> CUpwGeneralHLLC_Flow::ComputeResidual(const CConfig* c Energy_j = Enthalpy_j - Pressure_j / Density_j; - StaticEnthalpy_j = Enthalpy_j - 0.5 * sq_vel_j; - StaticEnergy_j = Energy_j - 0.5 * sq_vel_j; + StaticEnthalpy_j = Enthalpy_j - 0.5 * sq_vel_j - turb_ke_j; + StaticEnergy_j = Energy_j - 0.5 * sq_vel_j - turb_ke_j; Kappa_j = S_j[1] / Density_j; Chi_j = S_j[0] - Kappa_j * StaticEnergy_j; @@ -691,7 +694,8 @@ CNumerics::ResidualType<> CUpwGeneralHLLC_Flow::ComputeResidual(const CConfig* c /*--- Roe-averaged speed of sound ---*/ //RoeSoundSpeed2 = RoeChi + RoeKappa * ( RoeEnthalpy - 0.5 * sq_velRoe ); - RoeSoundSpeed = sqrt( RoeChi + RoeKappa * ( RoeEnthalpy - 0.5 * sq_velRoe ) ) - ProjInterfaceVel; + RoeSoundSpeed = sqrt( RoeChi + RoeKappa * ( RoeEnthalpy - 0.5 * sq_velRoe + - ( sqrt(Density_j) * turb_ke_j + sqrt(Density_i) * turb_ke_i ) / Rrho ) ) - ProjInterfaceVel; /*--- Speed of sound at L and R ---*/ @@ -784,7 +788,7 @@ CNumerics::ResidualType<> CUpwGeneralHLLC_Flow::ComputeResidual(const CConfig* c Jacobian_j[iVar][jVar] = 0; - GetInviscidProjJac(Velocity_i, &Enthalpy_i, &Chi_i, &Kappa_i, UnitNormal, 1.0, Jacobian_i); + GetInviscidProjJac(Velocity_i, &Enthalpy_i, &Chi_i, &Kappa_i, UnitNormal, 1.0, Jacobian_i, turb_ke_i); } else { @@ -800,7 +804,7 @@ CNumerics::ResidualType<> CUpwGeneralHLLC_Flow::ComputeResidual(const CConfig* c /*--- Computing pressure derivatives d/dU_L (PI) ---*/ - dPI_dU[0] = Chi_i - 0.5 * Kappa_i * sq_vel_i; + dPI_dU[0] = Chi_i - 0.5 * Kappa_i * sq_vel_i - Kappa_i * turb_ke_i; // k (SST) held fixed for (iDim = 0; iDim < nDim; iDim++) dPI_dU[iDim+1] = - Kappa_i * Velocity_i[iDim]; dPI_dU[nVar-1] = Kappa_i; @@ -874,7 +878,7 @@ CNumerics::ResidualType<> CUpwGeneralHLLC_Flow::ComputeResidual(const CConfig* c /*--- Computing pressure derivatives d/dU_R (PI) ---*/ - dPI_dU[0] = Chi_j - 0.5 * Kappa_j * sq_vel_j; + dPI_dU[0] = Chi_j - 0.5 * Kappa_j * sq_vel_j - Kappa_j * turb_ke_j; // k (SST) held fixed for (iDim = 0; iDim < nDim; iDim++) dPI_dU[iDim+1] = - Kappa_j * Velocity_j[iDim]; dPI_dU[nVar-1] = Kappa_j; @@ -928,7 +932,7 @@ CNumerics::ResidualType<> CUpwGeneralHLLC_Flow::ComputeResidual(const CConfig* c for (jVar = 0; jVar < nVar; jVar++) Jacobian_i[iVar][jVar] = 0; - GetInviscidProjJac(Velocity_j, &Enthalpy_j, &Chi_j, &Kappa_j, UnitNormal, 1.0, Jacobian_j); + GetInviscidProjJac(Velocity_j, &Enthalpy_j, &Chi_j, &Kappa_j, UnitNormal, 1.0, Jacobian_j, turb_ke_j); } else { @@ -944,7 +948,7 @@ CNumerics::ResidualType<> CUpwGeneralHLLC_Flow::ComputeResidual(const CConfig* c /*--- Computing pressure derivatives d/dU_L (PI) ---*/ - dPI_dU[0] = Chi_i - 0.5 * Kappa_i * sq_vel_i; + dPI_dU[0] = Chi_i - 0.5 * Kappa_i * sq_vel_i - Kappa_i * turb_ke_i; // k (SST) held fixed for (iDim = 0; iDim < nDim; iDim++) dPI_dU[iDim+1] = - Kappa_i * Velocity_i[iDim]; dPI_dU[nVar-1] = Kappa_i; @@ -995,7 +999,7 @@ CNumerics::ResidualType<> CUpwGeneralHLLC_Flow::ComputeResidual(const CConfig* c /*--- Computing pressure derivatives d/dU_R (PI) ---*/ - dPI_dU[0] = Chi_j - 0.5 * Kappa_j * sq_vel_j; + dPI_dU[0] = Chi_j - 0.5 * Kappa_j * sq_vel_j - Kappa_j * turb_ke_j; // k (SST) held fixed for (iDim = 0; iDim < nDim; iDim++) dPI_dU[iDim+1] = - Kappa_j * Velocity_j[iDim]; dPI_dU[nVar-1] = Kappa_j; diff --git a/SU2_CFD/src/numerics/flow/convection/roe.cpp b/SU2_CFD/src/numerics/flow/convection/roe.cpp index 73114576d409..9fffc3435651 100644 --- a/SU2_CFD/src/numerics/flow/convection/roe.cpp +++ b/SU2_CFD/src/numerics/flow/convection/roe.cpp @@ -101,6 +101,7 @@ CNumerics::ResidualType<> CUpwRoeBase_Flow::ComputeResidual(const CConfig* confi AD::StartPreacc(); AD::SetPreaccIn(V_i, nDim+4); AD::SetPreaccIn(V_j, nDim+4); AD::SetPreaccIn(Normal, nDim); + AD::SetPreaccIn(turb_ke_i); AD::SetPreaccIn(turb_ke_j); // SST: k in the total energy if (dynamic_grid) { AD::SetPreaccIn(GridVel_i, nDim); AD::SetPreaccIn(GridVel_j, nDim); } @@ -146,7 +147,9 @@ CNumerics::ResidualType<> CUpwRoeBase_Flow::ComputeResidual(const CConfig* confi sq_vel += RoeVelocity[iDim]*RoeVelocity[iDim]; } RoeEnthalpy = (R*Enthalpy_j+Enthalpy_i)/(R+1); - RoeSoundSpeed2 = (Gamma-1)*(RoeEnthalpy-0.5*sq_vel); + /*--- With SST the total enthalpy contains k, which is not part of the speed of sound. ---*/ + RoeTke = (R*turb_ke_j+turb_ke_i)/(R+1); + RoeSoundSpeed2 = (Gamma-1)*(RoeEnthalpy-0.5*sq_vel-RoeTke); /*--- Negative RoeSoundSpeed^2, the jump variables is too large, clear fluxes and exit. ---*/ @@ -170,7 +173,7 @@ CNumerics::ResidualType<> CUpwRoeBase_Flow::ComputeResidual(const CConfig* confi /*--- P tensor ---*/ - GetPMatrix(RoeDensity, RoeVelocity, RoeSoundSpeed, UnitNormal, P_Tensor); + GetPMatrix(RoeDensity, RoeVelocity, RoeSoundSpeed, UnitNormal, P_Tensor, RoeTke); /*--- Projected velocity adjusted for mesh motion ---*/ @@ -222,8 +225,8 @@ CNumerics::ResidualType<> CUpwRoeBase_Flow::ComputeResidual(const CConfig* confi Flux[iVar] = 0.5*(ProjFlux_i[iVar]+ProjFlux_j[iVar]); if (implicit) { - GetInviscidProjJac(Velocity_i, &Energy_i, Normal, 0.5, Jacobian_i); - GetInviscidProjJac(Velocity_j, &Energy_j, Normal, 0.5, Jacobian_j); + GetInviscidProjJac(Velocity_i, &Energy_i, Normal, 0.5, Jacobian_i, turb_ke_i); + GetInviscidProjJac(Velocity_j, &Energy_j, Normal, 0.5, Jacobian_j, turb_ke_j); } /*--- Finalize in children class ---*/ @@ -261,12 +264,17 @@ void CUpwRoe_Flow::FinalizeResidual(su2double *val_residual, su2double **val_Jac unsigned short iVar, jVar, kVar; /*--- Compute inverse P tensor ---*/ - GetPMatrix_inv(RoeDensity, RoeVelocity, RoeSoundSpeed, UnitNormal, invP_Tensor); + GetPMatrix_inv(RoeDensity, RoeVelocity, RoeSoundSpeed, UnitNormal, invP_Tensor, RoeTke); /*--- Diference between conservative variables at jPoint and iPoint ---*/ for (iVar = 0; iVar < nVar; iVar++) Diff_U[iVar] = Conservatives_j[iVar]-Conservatives_i[iVar]; + /*--- With SST rho*E contains rho*k, and Delta(rho k) = RoeDensity Delta k + RoeTke Delta rho. The part with + Delta k is not a pressure jump: remove it before the projection and advect it with the contact wave. ---*/ + const su2double rhoDeltaTke = RoeDensity*(turb_ke_j-turb_ke_i); + Diff_U[nVar-1] -= rhoDeltaTke; + /*--- Low dissipation formulation ---*/ if (roe_low_dissipation) Dissipation_ij = GetRoe_Dissipation(Dissipation_i, Dissipation_j, Sensor_i, Sensor_j, config); @@ -291,6 +299,7 @@ void CUpwRoe_Flow::FinalizeResidual(su2double *val_residual, su2double **val_Jac } } } + val_residual[nVar-1] -= (1.0-kappa)*Lambda[0]*rhoDeltaTke*Area*Dissipation_ij; } @@ -347,11 +356,14 @@ void CUpwL2Roe_Flow::FinalizeResidual(su2double *val_residual, su2double **val_J for (kVar = 0; kVar < nVar; kVar++) val_residual[iVar] -= (1.0-kappa)*Lambda[kVar]*delta_wave[kVar]*P_Tensor[iVar][kVar]*Area; + /*--- With SST, the jump of k in rho*E is advected with the contact wave. ---*/ + val_residual[nVar-1] -= (1.0-kappa)*Lambda[0]*RoeDensity*(turb_ke_j-turb_ke_i)*Area; + if (!implicit) return; /*--- If implicit use the Jacobians of the standard Roe scheme as an approximation ---*/ - GetPMatrix_inv(RoeDensity, RoeVelocity, RoeSoundSpeed, UnitNormal, invP_Tensor); + GetPMatrix_inv(RoeDensity, RoeVelocity, RoeSoundSpeed, UnitNormal, invP_Tensor, RoeTke); for (iVar = 0; iVar < nVar; iVar++) { for (jVar = 0; jVar < nVar; jVar++) { @@ -420,11 +432,14 @@ void CUpwLMRoe_Flow::FinalizeResidual(su2double *val_residual, su2double **val_J for (kVar = 0; kVar < nVar; kVar++) val_residual[iVar] -= (1.0-kappa)*Lambda[kVar]*delta_wave[kVar]*P_Tensor[iVar][kVar]*Area; + /*--- With SST, the jump of k in rho*E is advected with the contact wave. ---*/ + val_residual[nVar-1] -= (1.0-kappa)*Lambda[0]*RoeDensity*(turb_ke_j-turb_ke_i)*Area; + if (!implicit) return; /*--- If implicit use the Jacobians of the standard Roe scheme as an approximation ---*/ - GetPMatrix_inv(RoeDensity, RoeVelocity, RoeSoundSpeed, UnitNormal, invP_Tensor); + GetPMatrix_inv(RoeDensity, RoeVelocity, RoeSoundSpeed, UnitNormal, invP_Tensor, RoeTke); for (iVar = 0; iVar < nVar; iVar++) { for (jVar = 0; jVar < nVar; jVar++) { @@ -562,7 +577,8 @@ CNumerics::ResidualType<> CUpwTurkel_Flow::ComputeResidual(const CConfig* config sq_vel += RoeVelocity[iDim]*RoeVelocity[iDim]; } RoeEnthalpy = (R*Enthalpy_j+Enthalpy_i)/(R+1); - RoeSoundSpeed = sqrt(fabs((Gamma-1)*(RoeEnthalpy-0.5*sq_vel))); + /*--- With SST the total enthalpy contains k, which is not part of the speed of sound. ---*/ + RoeSoundSpeed = sqrt(fabs((Gamma-1)*(RoeEnthalpy-0.5*sq_vel-(R*turb_ke_j+turb_ke_i)/(R+1)))); RoePressure = RoeDensity/Gamma*RoeSoundSpeed*RoeSoundSpeed; /*--- Compute ProjFlux_i ---*/ @@ -629,8 +645,8 @@ CNumerics::ResidualType<> CUpwTurkel_Flow::ComputeResidual(const CConfig* config if (implicit) { /*--- Jacobians of the inviscid flux, scaled by 0.5 because Flux ~ 0.5*(fc_i+fc_j)*Normal ---*/ - GetInviscidProjJac(Velocity_i, &Energy_i, Normal, 0.5, Jacobian_i); - GetInviscidProjJac(Velocity_j, &Energy_j, Normal, 0.5, Jacobian_j); + GetInviscidProjJac(Velocity_i, &Energy_i, Normal, 0.5, Jacobian_i, turb_ke_i); + GetInviscidProjJac(Velocity_j, &Energy_j, Normal, 0.5, Jacobian_j, turb_ke_j); } for (iVar = 0; iVar < nVar; iVar ++) { @@ -741,6 +757,7 @@ CNumerics::ResidualType<> CUpwGeneralRoe_Flow::ComputeResidual(const CConfig* co AD::StartPreacc(); AD::SetPreaccIn(V_i, nDim+4); AD::SetPreaccIn(V_j, nDim+4); AD::SetPreaccIn(Normal, nDim); + AD::SetPreaccIn(turb_ke_i); AD::SetPreaccIn(turb_ke_j); // SST: k in the total energy AD::SetPreaccIn(S_i, 2); AD::SetPreaccIn(S_j, 2); if (dynamic_grid) { AD::SetPreaccIn(GridVel_i, nDim); AD::SetPreaccIn(GridVel_j, nDim); @@ -768,7 +785,7 @@ CNumerics::ResidualType<> CUpwGeneralRoe_Flow::ComputeResidual(const CConfig* co Density_i = V_i[nDim+2]; Enthalpy_i = V_i[nDim+3]; Energy_i = Enthalpy_i - Pressure_i/Density_i; - StaticEnthalpy_i = Enthalpy_i - 0.5*Velocity2_i; + StaticEnthalpy_i = Enthalpy_i - 0.5*Velocity2_i - turb_ke_i; // SST: H contains k StaticEnergy_i = StaticEnthalpy_i - Pressure_i/Density_i; Kappa_i = S_i[1]/Density_i; @@ -788,7 +805,7 @@ CNumerics::ResidualType<> CUpwGeneralRoe_Flow::ComputeResidual(const CConfig* co Enthalpy_j = V_j[nDim+3]; Energy_j = Enthalpy_j - Pressure_j/Density_j; - StaticEnthalpy_j = Enthalpy_j - 0.5*Velocity2_j; + StaticEnthalpy_j = Enthalpy_j - 0.5*Velocity2_j - turb_ke_j; // SST: H contains k StaticEnergy_j = StaticEnthalpy_j - Pressure_j/Density_j; Kappa_j = S_j[1]/Density_j; @@ -831,7 +848,8 @@ CNumerics::ResidualType<> CUpwGeneralRoe_Flow::ComputeResidual(const CConfig* co /*--- Compute P and Lambda (do it with the Normal) ---*/ - GetPMatrix(&RoeDensity, RoeVelocity, &RoeSoundSpeed, &RoeEnthalpy, &RoeChi, &RoeKappa, UnitNormal, P_Tensor); + GetPMatrix(&RoeDensity, RoeVelocity, &RoeSoundSpeed, &RoeEnthalpy, &RoeChi, &RoeKappa, UnitNormal, P_Tensor, + (R*turb_ke_j+turb_ke_i)/(R+1)); ProjVelocity = 0.0; ProjVelocity_i = 0.0; ProjVelocity_j = 0.0; for (iDim = 0; iDim < nDim; iDim++) { @@ -914,6 +932,8 @@ CNumerics::ResidualType<> CUpwGeneralRoe_Flow::ComputeResidual(const CConfig* co for (jVar = 0; jVar < nVar; jVar++) Flux[iVar] -= 0.5*Lambda[jVar]*delta_wave[jVar]*P_Tensor[iVar][jVar]*Area; } + /*--- With SST, the jump of k in rho*E is advected with the contact wave. ---*/ + Flux[nVar-1] -= 0.5*Lambda[0]*RoeDensity*(turb_ke_j-turb_ke_i)*Area; /*--- Flux contribution due to grid motion ---*/ if (dynamic_grid) { @@ -929,20 +949,26 @@ CNumerics::ResidualType<> CUpwGeneralRoe_Flow::ComputeResidual(const CConfig* co /*--- Compute inverse P ---*/ - GetPMatrix_inv(invP_Tensor, &RoeDensity, RoeVelocity, &RoeSoundSpeed, &RoeChi , &RoeKappa, UnitNormal); + GetPMatrix_inv(invP_Tensor, &RoeDensity, RoeVelocity, &RoeSoundSpeed, &RoeChi , &RoeKappa, UnitNormal, + (R*turb_ke_j+turb_ke_i)/(R+1)); /*--- Jacobians of the inviscid flux, scaled by 0.5 because val_resconv ~ 0.5*(fc_i+fc_j)*Normal ---*/ - GetInviscidProjJac(Velocity_i, &Enthalpy_i, &Chi_i, &Kappa_i, Normal, 0.5, Jacobian_i); + GetInviscidProjJac(Velocity_i, &Enthalpy_i, &Chi_i, &Kappa_i, Normal, 0.5, Jacobian_i, turb_ke_i); - GetInviscidProjJac(Velocity_j, &Enthalpy_j, &Chi_j, &Kappa_j, Normal, 0.5, Jacobian_j); + GetInviscidProjJac(Velocity_j, &Enthalpy_j, &Chi_j, &Kappa_j, Normal, 0.5, Jacobian_j, turb_ke_j); /*--- Diference variables iPoint and jPoint ---*/ for (iVar = 0; iVar < nVar; iVar++) Diff_U[iVar] = U_j[iVar]-U_i[iVar]; + /*--- With SST, the part of the jump of rho*E due to the jump of k is not a pressure jump: remove it before the + projection and advect it with the contact wave. ---*/ + const su2double rhoDeltaTke = RoeDensity*(turb_ke_j-turb_ke_i); + Diff_U[nVar-1] -= rhoDeltaTke; + /*--- Roe's Flux approximation ---*/ for (iVar = 0; iVar < nVar; iVar++) { Flux[iVar] = 0.5*(ProjFlux_i[iVar]+ProjFlux_j[iVar]); @@ -959,6 +985,7 @@ CNumerics::ResidualType<> CUpwGeneralRoe_Flow::ComputeResidual(const CConfig* co Jacobian_j[iVar][jVar] -= (1.0-kappa)*Proj_ModJac_Tensor_ij*Area; } } + Flux[nVar-1] -= (1.0-kappa)*Lambda[0]*rhoDeltaTke*Area; /*--- Jacobian contributions due to grid motion ---*/ if (dynamic_grid) { @@ -1017,6 +1044,6 @@ void CUpwGeneralRoe_Flow::ComputeRoeAverage() { // // } - RoeSoundSpeed2 = RoeChi + RoeKappa*(RoeEnthalpy-0.5*sq_vel); + RoeSoundSpeed2 = RoeChi + RoeKappa*(RoeEnthalpy-0.5*sq_vel-(R*turb_ke_j+turb_ke_i)/(R+1)); // SST: H contains k } 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 8cf9642b4251..129d34f1462d 100644 --- a/SU2_CFD/src/numerics/flow/flow_sources.cpp +++ b/SU2_CFD/src/numerics/flow/flow_sources.cpp @@ -54,7 +54,8 @@ CSourceAxisymmetric_Flow::CSourceAxisymmetric_Flow(unsigned short val_nDim, unsi implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); viscous = config->GetViscous(); - rans = (config->GetKind_Turb_Model() != TURB_MODEL::NONE); + tkeInEnergy = (config->GetKind_Turb_Model() == TURB_MODEL::SST); + tkeInStress = tkeInEnergy && !config->GetSSTParsedOptions().modified; } @@ -73,7 +74,9 @@ CNumerics::ResidualType<> CSourceAxisymmetric_Flow::ComputeResidual(const CConfi sq_vel += Velocity_i *Velocity_i; } - Pressure_i = Gamma_Minus_One*U_i[0]*(U_i[nDim+1]/U_i[0]-0.5*sq_vel); + /*--- With SST the total energy contains k. ---*/ + const su2double tke = tkeInEnergy ? turb_ke_i : 0.0; + Pressure_i = Gamma_Minus_One*U_i[0]*(U_i[nDim+1]/U_i[0]-0.5*sq_vel-tke); Enthalpy_i = (U_i[nDim+1] + Pressure_i) / U_i[0]; residual[0] = yinv*Volume*U_i[2]; @@ -103,6 +106,7 @@ CNumerics::ResidualType<> CSourceAxisymmetric_Flow::ComputeResidual(const CConfi jacobian[3][1] = -(Gamma-1)*U_i[2]*U_i[1]/(U_i[0]*U_i[0]); jacobian[3][2] = Gamma*U_i[3]/U_i[0] - 1/2*(Gamma-1)*( (U_i[1]*U_i[1]+U_i[2]*U_i[2])/(U_i[0]*U_i[0]) + 2*U_i[2]*U_i[2]/(U_i[0]*U_i[0]) ); jacobian[3][3] = Gamma*U_i[2]/U_i[0]; + jacobian[3][2] -= Gamma_Minus_One*tke; // k is a variable of the turbulence solver for (iVar=0; iVar < nVar; iVar++) for (jVar=0; jVar < nVar; jVar++) @@ -135,7 +139,8 @@ CNumerics::ResidualType<> CSourceAxisymmetric_Flow::ComputeResidual(const CConfi void CSourceAxisymmetric_Flow::ResidualDiffusion(){ - if (!rans){ turb_ke_i = 0.0; } + /*--- -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]; @@ -155,7 +160,8 @@ void CSourceAxisymmetric_Flow::ResidualDiffusion(){ -TWO3*AuxVar_Grad_i[0][1]); 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+U_i[0]*turb_ke_i)) + -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/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 32c265e55656..de07a569f44a 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -65,6 +65,7 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config, (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND); const bool time_stepping = (config->GetTime_Marching() == TIME_MARCHING::TIME_STEPPING); const bool adjoint = config->GetContinuous_Adjoint() || config->GetDiscrete_Adjoint(); + const bool tkeNeeded = (rans && config->GetKind_Turb_Model() == TURB_MODEL::SST); int Unst_RestartIter = 0; unsigned long iPoint, iMarker, counter_local = 0, counter_global = 0; @@ -246,6 +247,7 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config, Density_Inf = config->GetDensity_FreeStreamND(); Energy_Inf = config->GetEnergy_FreeStreamND(); Mach_Inf = config->GetMach(); + const su2double TKE_Inf = config->GetTke_FreeStreamND(); /*--- Initialize the secondary values for direct derivative approximations ---*/ @@ -315,6 +317,7 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config, Velocity2 += pow(nodes->GetSolution(iPoint,iDim+1)/Density,2); StaticEnergy= nodes->GetEnergy(iPoint) - 0.5*Velocity2; + if (tkeNeeded) StaticEnergy -= TKE_Inf; GetFluidModel()->SetTDState_rhoe(Density, StaticEnergy); Pressure= GetFluidModel()->GetPressure(); @@ -1009,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 { @@ -1019,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. ---*/ @@ -1086,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) { @@ -1843,6 +1870,9 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain const bool muscl = (config->GetMUSCL_Flow() && (iMesh == MESH_0)); const bool limiter = (config->GetKind_SlopeLimit_Flow() != LIMITER::NONE); + /*--- SST: the total enthalpy of the cells contains k. ---*/ + const bool tkeInEnergy = (config->GetKind_Turb_Model() == TURB_MODEL::SST); + const CVariable* turbNodes = tkeInEnergy ? solver_container[TURB_SOL]->GetNodes() : nullptr; const bool van_albada = (config->GetKind_SlopeLimit_Flow() == LIMITER::VAN_ALBADA_EDGE); const su2double kappa = config->GetMUSCL_Kappa_Flow(); @@ -1883,6 +1913,13 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain numerics->SetNormal(geometry->edges->GetNormal(iEdge)); + su2double tke_i = 0.0, tke_j = 0.0; + if (tkeInEnergy) { + tke_i = turbNodes->GetSolution(iPoint, 0); + tke_j = turbNodes->GetSolution(jPoint, 0); + numerics->SetTurbKineticEnergy(tke_i, tke_j); + } + auto Coord_i = geometry->nodes->GetCoord(iPoint); auto Coord_j = geometry->nodes->GetCoord(jPoint); @@ -1949,15 +1986,17 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain /*--- Recompute the reconstructed quantities in a thermodynamically consistent way. ---*/ + /*--- With SST the total enthalpy contains k, taken from the cells (as the vectorized reconstruction does). ---*/ + if (!ideal_gas || low_mach_corr) { - ComputeConsistentExtrapolation(GetFluidModel(), nDim, Primitive_i, Secondary_i); - ComputeConsistentExtrapolation(GetFluidModel(), nDim, Primitive_j, Secondary_j); + ComputeConsistentExtrapolation(GetFluidModel(), nDim, Primitive_i, Secondary_i, tke_i); + ComputeConsistentExtrapolation(GetFluidModel(), nDim, Primitive_j, Secondary_j, tke_j); } /*--- Low-Mach number correction. ---*/ if (low_mach_corr) { - LowMachPrimitiveCorrection(GetFluidModel(), nDim, Primitive_i, Primitive_j); + LowMachPrimitiveCorrection(GetFluidModel(), nDim, Primitive_i, Primitive_j, tke_i, tke_j); } /*--- Check for non-physical solutions after reconstruction. If found, use the @@ -1976,7 +2015,14 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain } su2double RoeEnthalpy = (R * Primitive_j[prim_idx.Enthalpy()] + Primitive_i[prim_idx.Enthalpy()]) / (R+1); - const bool neg_sound_speed = ((Gamma-1)*(RoeEnthalpy-0.5*sq_vel) < 0.0); + const su2double RoeTke = (R * tke_j + tke_i) / (R+1); // SST: the total enthalpy contains k + /*--- Speed of sound of the Roe average and of each side (schemes such as HLLC use both). ---*/ + auto negSoundSpeed2 = [&](const su2double* prim, su2double tke) { + const su2double vel2 = GeometryToolbox::SquaredNorm(nDim, &prim[prim_idx.Velocity()]); + return (Gamma-1)*(prim[prim_idx.Enthalpy()] - 0.5*vel2 - tke) < 0.0; + }; + const bool neg_sound_speed = ((Gamma-1)*(RoeEnthalpy-0.5*sq_vel-RoeTke) < 0.0) || + negSoundSpeed2(Primitive_i, tke_i) || negSoundSpeed2(Primitive_j, tke_j); bool bad_recon = neg_sound_speed || neg_pres_or_rho_i || neg_pres_or_rho_j; bad_recon = nodes->UpdateNonPhysicalEdgeCounter(iEdge, bad_recon); counter_local += bad_recon; @@ -2036,7 +2082,7 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain } void CEulerSolver::ComputeConsistentExtrapolation(CFluidModel *fluidModel, unsigned short nDim, - su2double *primitive, su2double *secondary) { + su2double *primitive, su2double *secondary, su2double tke) { SU2_ZONE_SCOPED const CEulerVariable::CIndices prim_idx(nDim, 0); const su2double density = primitive[prim_idx.Density()]; @@ -2046,7 +2092,7 @@ void CEulerSolver::ComputeConsistentExtrapolation(CFluidModel *fluidModel, unsig fluidModel->SetTDState_Prho(pressure, density); primitive[prim_idx.Temperature()] = fluidModel->GetTemperature(); - primitive[prim_idx.Enthalpy()] = fluidModel->GetStaticEnergy() + pressure / density + 0.5*velocity2; + primitive[prim_idx.Enthalpy()] = fluidModel->GetStaticEnergy() + pressure / density + 0.5*velocity2 + tke; primitive[prim_idx.SoundSpeed()] = fluidModel->GetSoundSpeed(); secondary[0] = fluidModel->GetdPdrho_e(); secondary[1] = fluidModel->GetdPde_rho(); @@ -2054,7 +2100,8 @@ void CEulerSolver::ComputeConsistentExtrapolation(CFluidModel *fluidModel, unsig } void CEulerSolver::LowMachPrimitiveCorrection(CFluidModel *fluidModel, unsigned short nDim, - su2double *primitive_i, su2double *primitive_j) { + su2double *primitive_i, su2double *primitive_j, + su2double tke_i, su2double tke_j) { SU2_ZONE_SCOPED unsigned short iDim; @@ -2085,10 +2132,10 @@ void CEulerSolver::LowMachPrimitiveCorrection(CFluidModel *fluidModel, unsigned } fluidModel->SetEnergy_Prho(primitive_i[nDim+1], primitive_i[nDim+2]); - primitive_i[nDim+3]= fluidModel->GetStaticEnergy() + primitive_i[nDim+1]/primitive_i[nDim+2] + 0.5*velocity2_i; + primitive_i[nDim+3]= fluidModel->GetStaticEnergy() + primitive_i[nDim+1]/primitive_i[nDim+2] + 0.5*velocity2_i + tke_i; fluidModel->SetEnergy_Prho(primitive_j[nDim+1], primitive_j[nDim+2]); - primitive_j[nDim+3]= fluidModel->GetStaticEnergy() + primitive_j[nDim+1]/primitive_j[nDim+2] + 0.5*velocity2_j; + primitive_j[nDim+3]= fluidModel->GetStaticEnergy() + primitive_j[nDim+1]/primitive_j[nDim+2] + 0.5*velocity2_j + tke_j; } @@ -2205,7 +2252,7 @@ void CEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_contain numerics->SetAuxVarGrad(nodes->GetAuxVarGradient(iPoint), nullptr); /*--- Set turbulence kinetic energy ---*/ - if (rans){ + if (config->GetKind_Turb_Model() == TURB_MODEL::SST){ CVariable* turbNodes = solver_container[TURB_SOL]->GetNodes(); numerics->SetTurbKineticEnergy(turbNodes->GetSolution(iPoint,0), turbNodes->GetSolution(iPoint,0)); } @@ -5016,7 +5063,9 @@ void CEulerSolver::BC_Far_Field(CGeometry *geometry, CSolver **solver_container, bool implicit = config->GetKind_TimeIntScheme() == EULER_IMPLICIT; bool viscous = config->GetViscous(); - bool tkeNeeded = config->GetKind_Turb_Model() == TURB_MODEL::SST; + bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST); + CVariable* turbNodes = nullptr; + if (tkeNeeded) turbNodes = solver_container[TURB_SOL]->GetNodes(); auto *Normal = new su2double[nDim]; @@ -5166,7 +5215,8 @@ void CEulerSolver::BC_Far_Field(CGeometry *geometry, CSolver **solver_container, } Pressure = Density*SoundSpeed*SoundSpeed/Gamma; Energy = Pressure/(Gamma_Minus_One*Density) + 0.5*Velocity2; - if (tkeNeeded) Energy += GetTke_Inf(); + /*--- k of the boundary node, not the free-stream value, which decays right away in the domain (#1851). ---*/ + if (tkeNeeded) Energy += turbNodes->GetSolution(iPoint,0); /*--- Store new primitive state for computing the flux. ---*/ @@ -5183,6 +5233,16 @@ void CEulerSolver::BC_Far_Field(CGeometry *geometry, CSolver **solver_container, conv_numerics->SetPrimitive(V_domain, V_infty); + /*--- SST: both states contain the k of the node in their total enthalpy. ---*/ + + if (config->GetKind_Turb_Model() == TURB_MODEL::SST) { + + const su2double tke = solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint, 0); + + conv_numerics->SetTurbKineticEnergy(tke, tke); + + } + if (dynamic_grid) { conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), geometry->nodes->GetGridVel(iPoint)); @@ -5263,6 +5323,8 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, tkeNeeded = config->GetKind_Turb_Model() == TURB_MODEL::SST, ideal_gas = config->GetKind_FluidModel() == STANDARD_AIR || config->GetKind_FluidModel() == IDEAL_GAS; + CVariable* turbNodes = nullptr; + if (tkeNeeded) turbNodes = solver_container[TURB_SOL]->GetNodes(); su2double **P_Tensor = new su2double*[nVar], **invP_Tensor = new su2double*[nVar]; @@ -5308,7 +5370,8 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, const auto Density_i = nodes->GetDensity(iPoint); const auto Energy_i = nodes->GetEnergy(iPoint); - const su2double StaticEnergy_i = Energy_i - 0.5*Velocity2_i; + su2double StaticEnergy_i = Energy_i - 0.5*Velocity2_i; + if (tkeNeeded) StaticEnergy_i -= turbNodes->GetSolution(iPoint, 0); GetFluidModel()->SetTDState_rhoe(Density_i, StaticEnergy_i); @@ -5362,7 +5425,6 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, Density_e = GetFluidModel()->GetDensity(); StaticEnergy_e = GetFluidModel()->GetStaticEnergy(); Energy_e = StaticEnergy_e + 0.5 * Velocity2_e; - if (tkeNeeded) Energy_e += GetTke_Inf(); break; case STATIC_SUPERSONIC_INFLOW_PT: @@ -5388,7 +5450,6 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, Density_e = GetFluidModel()->GetDensity(); StaticEnergy_e = GetFluidModel()->GetStaticEnergy(); Energy_e = StaticEnergy_e + 0.5 * Velocity2_e; - if (tkeNeeded) Energy_e += GetTke_Inf(); break; case STATIC_SUPERSONIC_INFLOW_PD: @@ -5415,7 +5476,6 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, Density_e = GetFluidModel()->GetDensity(); StaticEnergy_e = GetFluidModel()->GetStaticEnergy(); Energy_e = StaticEnergy_e + 0.5 * Velocity2_e; - if (tkeNeeded) Energy_e += GetTke_Inf(); break; case DENSITY_VELOCITY: @@ -5431,7 +5491,9 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, for (auto iDim = 0u; iDim < nDim; iDim++) Velocity_e[iDim] = VelMag_e*Flow_Dir[iDim]; + /*--- The total energy of the node already contains k, which is added below to every case. ---*/ Energy_e = Energy_i; + if (tkeNeeded) Energy_e -= turbNodes->GetSolution(iPoint, 0); break; case STATIC_PRESSURE: @@ -5455,11 +5517,18 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, break; } + /*--- k of the boundary node, not the free-stream value, which decays right away in the domain (#1851). + It is also removed from the static energy of the boundary state. ---*/ + const su2double Tke_e = tkeNeeded ? turbNodes->GetSolution(iPoint, 0) : 0.0; + Energy_e += Tke_e; + /*--- Compute P (matrix of right eigenvectors) ---*/ - conv_numerics->GetPMatrix(&Density_i, Velocity_i, &SoundSpeed_i, &Enthalpy_i, &Chi_i, &Kappa_i, UnitNormal, P_Tensor); + conv_numerics->GetPMatrix(&Density_i, Velocity_i, &SoundSpeed_i, &Enthalpy_i, &Chi_i, &Kappa_i, UnitNormal, P_Tensor, + Tke_e); /*--- Compute inverse P (matrix of left eigenvectors)---*/ - conv_numerics->GetPMatrix_inv(invP_Tensor, &Density_i, Velocity_i, &SoundSpeed_i, &Chi_i, &Kappa_i, UnitNormal); + conv_numerics->GetPMatrix_inv(invP_Tensor, &Density_i, Velocity_i, &SoundSpeed_i, &Chi_i, &Kappa_i, UnitNormal, + Tke_e); /*--- eigenvalues contribution due to grid motion ---*/ su2double ProjVelocity_i = GeometryToolbox::DotProduct(nDim, Velocity_i, UnitNormal); @@ -5508,7 +5577,7 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, const auto Velocity2_b = GeometryToolbox::SquaredNorm(nDim, Velocity_b); const su2double Energy_b = u_b[nVar-1]/Density_b; - const su2double StaticEnergy_b = Energy_b - 0.5*Velocity2_b; + const su2double StaticEnergy_b = Energy_b - 0.5*Velocity2_b - Tke_e; GetFluidModel()->SetTDState_rhoe(Density_b, StaticEnergy_b); /*--- Store number of Newton iterations at BC ---*/ @@ -5562,7 +5631,7 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, } /*--- Compute flux Jacobian in state b ---*/ - conv_numerics->GetInviscidProjJac(Velocity_b, &Enthalpy_b, &Chi_b, &Kappa_b, Normal, 1.0, Jacobian_b); + conv_numerics->GetInviscidProjJac(Velocity_b, &Enthalpy_b, &Chi_b, &Kappa_b, Normal, 1.0, Jacobian_b, Tke_e); /*--- Jacobian contribution due to grid motion ---*/ if (dynamic_grid){ @@ -5652,7 +5721,7 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, /*--- Turbulent kinetic energy ---*/ - if (config->GetKind_Turb_Model() == TURB_MODEL::SST) + if (tkeNeeded) visc_numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0), solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0)); @@ -5708,6 +5777,8 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST); const bool ideal_gas = config->GetKind_FluidModel() == STANDARD_AIR || config->GetKind_FluidModel() == IDEAL_GAS; + CVariable* turbNodes = nullptr; + if (tkeNeeded) turbNodes = solver_container[TURB_SOL]->GetNodes(); su2double *Normal, *turboNormal, *UnitNormal, *FlowDirMix, FlowDirMixMag, *turboVelocity; Normal = new su2double[nDim]; @@ -5781,6 +5852,7 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain Energy_i = nodes->GetEnergy(iPoint); StaticEnergy_i = Energy_i - 0.5*Velocity2_i; + if (tkeNeeded) StaticEnergy_i -= turbNodes->GetSolution(iPoint, 0); GetFluidModel()->SetTDState_rhoe(Density_i, StaticEnergy_i); @@ -5844,7 +5916,7 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain Density_e = GetFluidModel()->GetDensity(); StaticEnergy_e = GetFluidModel()->GetStaticEnergy(); Energy_e = StaticEnergy_e + 0.5 * Velocity2_e; - if (tkeNeeded) Energy_e += GetTke_Inf(); + if (tkeNeeded) Energy_e += turbNodes->GetSolution(iPoint, 0); break; case MIXING_IN: @@ -5879,7 +5951,7 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain Density_e = GetFluidModel()->GetDensity(); StaticEnergy_e = GetFluidModel()->GetStaticEnergy(); Energy_e = StaticEnergy_e + 0.5 * Velocity2_e; - // if (tkeNeeded) Energy_e += GetTke_Inf(); + if (tkeNeeded) Energy_e += turbNodes->GetSolution(iPoint, 0); break; @@ -5897,6 +5969,7 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain Velocity2_e += Velocity_e[iDim]*Velocity_e[iDim]; } Energy_e = GetFluidModel()->GetStaticEnergy() + 0.5*Velocity2_e; + if (tkeNeeded) Energy_e += turbNodes->GetSolution(iPoint, 0); break; case STATIC_PRESSURE: @@ -5914,6 +5987,7 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain Velocity2_e += Velocity_e[iDim]*Velocity_e[iDim]; } Energy_e = GetFluidModel()->GetStaticEnergy() + 0.5*Velocity2_e; + if (tkeNeeded) Energy_e += turbNodes->GetSolution(iPoint, 0); break; @@ -5931,6 +6005,7 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain Velocity2_e += Velocity_e[iDim]*Velocity_e[iDim]; } Energy_e = GetFluidModel()->GetStaticEnergy() + 0.5*Velocity2_e; + if (tkeNeeded) Energy_e += turbNodes->GetSolution(iPoint, 0); break; default: @@ -5939,10 +6014,13 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain } /*--- Compute P (matrix of right eigenvectors) ---*/ - conv_numerics->GetPMatrix(&Density_i, Velocity_i, &SoundSpeed_i, &Enthalpy_i, &Chi_i, &Kappa_i, UnitNormal, P_Tensor); + const su2double Tke_i = tkeNeeded ? turbNodes->GetSolution(iPoint, 0) : 0.0; // SST: k in the total energy + conv_numerics->GetPMatrix(&Density_i, Velocity_i, &SoundSpeed_i, &Enthalpy_i, &Chi_i, &Kappa_i, UnitNormal, P_Tensor, + Tke_i); /*--- Compute inverse P (matrix of left eigenvectors)---*/ - conv_numerics->GetPMatrix_inv(invP_Tensor, &Density_i, Velocity_i, &SoundSpeed_i, &Chi_i, &Kappa_i, UnitNormal); + conv_numerics->GetPMatrix_inv(invP_Tensor, &Density_i, Velocity_i, &SoundSpeed_i, &Chi_i, &Kappa_i, UnitNormal, + Tke_i); /*--- eigenvalues contribution due to grid motion ---*/ if (dynamic_grid){ @@ -6006,7 +6084,7 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain Velocity2_b += Velocity_b[iDim]*Velocity_b[iDim]; } Energy_b = u_b[nVar-1]/Density_b; - StaticEnergy_b = Energy_b - 0.5*Velocity2_b; + StaticEnergy_b = Energy_b - 0.5*Velocity2_b - Tke_i; // SST: k in the total energy GetFluidModel()->SetTDState_rhoe(Density_b, StaticEnergy_b); Pressure_b = GetFluidModel()->GetPressure(); Temperature_b = GetFluidModel()->GetTemperature(); @@ -6062,7 +6140,7 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain } /*--- Compute flux Jacobian in state b ---*/ - conv_numerics->GetInviscidProjJac(Velocity_b, &Enthalpy_b, &Chi_b, &Kappa_b, Normal, 1.0, Jacobian_b); + conv_numerics->GetInviscidProjJac(Velocity_b, &Enthalpy_b, &Chi_b, &Kappa_b, Normal, 1.0, Jacobian_b, Tke_i); /*--- Jacobian contribution due to grid motion ---*/ if (dynamic_grid) @@ -6764,6 +6842,9 @@ void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container, CNu Energy_i = nodes->GetEnergy(iPoint); StaticEnergy_i = Energy_i - 0.5*Velocity2_i; + /*--- SST: the total energy contains k. ---*/ + if (config->GetKind_Turb_Model() == TURB_MODEL::SST) + StaticEnergy_i -= solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint, 0); GetFluidModel()->SetTDState_rhoe(Density_i, StaticEnergy_i); @@ -7020,6 +7101,9 @@ void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container, CNu Energy_b = GetFluidModel()->GetStaticEnergy() + 0.5*Velocity2_b; Temperature_b= GetFluidModel()->GetTemperature(); Enthalpy_b = Energy_b + Pressure_b/Density_b; + /*--- SST: the boundary state carries the k of the node in its total enthalpy, as the interior one. ---*/ + if (config->GetKind_Turb_Model() == TURB_MODEL::SST) + Enthalpy_b += solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint, 0); /*--- Primitive variables, using the derived quantities ---*/ V_boundary[0] = Temperature_b; @@ -7037,6 +7121,16 @@ void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container, CNu /*--- Set various quantities in the solver class ---*/ conv_numerics->SetPrimitive(V_domain, V_boundary); + + /*--- SST: both states contain the k of the node in their total enthalpy. ---*/ + + if (config->GetKind_Turb_Model() == TURB_MODEL::SST) { + + const su2double tke = solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint, 0); + + conv_numerics->SetTurbKineticEnergy(tke, tke); + + } conv_numerics->SetSecondary(S_domain, S_boundary); @@ -7260,6 +7354,8 @@ void CEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, Velocity2 += Velocity[iDim]*Velocity[iDim]; } Energy = V_domain[nDim+3] - V_domain[nDim+1]/V_domain[nDim+2]; + /*--- SST: the total enthalpy of the domain contains k, the imposed total enthalpy does not. ---*/ + if (tkeNeeded) Energy -= solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0); Pressure = V_domain[nDim+1]; H_Total = (Gamma*Gas_Constant/Gamma_Minus_One)*T_Total; SoundSpeed2 = Gamma*Pressure/Density; @@ -7330,7 +7426,6 @@ void CEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, /*--- Using pressure, density, & velocity, compute the energy ---*/ Energy = Pressure/(Density*Gamma_Minus_One) + 0.5*Velocity2; - if (tkeNeeded) Energy += GetTke_Inf(); /*--- Primitive variables, using the derived quantities ---*/ @@ -7341,6 +7436,9 @@ void CEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, V_inlet[nDim+2] = Density; V_inlet[nDim+3] = Energy + Pressure/Density; + /*--- k of the boundary node, not the free-stream value, which decays right away in the domain (#1851). ---*/ + if (tkeNeeded) V_inlet[nDim+3] += solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0); + break; } /*--- Mass flow has been specified at the inlet. ---*/ @@ -7402,7 +7500,6 @@ void CEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, /*--- Energy for the fictitious inlet state ---*/ Energy = Pressure/(Density*Gamma_Minus_One) + 0.5*Vel_Mag*Vel_Mag; - if (tkeNeeded) Energy += GetTke_Inf(); /*--- Primitive variables, using the derived quantities ---*/ Temperature = Pressure / ( Gas_Constant * Density); @@ -7413,6 +7510,9 @@ void CEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, V_inlet[nDim+2] = Density; V_inlet[nDim+3] = Energy + Pressure/Density; + /*--- k of the boundary node, not the free-stream value, which decays right away in the domain (#1851). ---*/ + if (tkeNeeded) V_inlet[nDim+3] += solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0); + break; } default: @@ -7424,6 +7524,16 @@ void CEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, conv_numerics->SetPrimitive(V_domain, V_inlet); + /*--- SST: both states contain the k of the node in their total enthalpy. ---*/ + + if (config->GetKind_Turb_Model() == TURB_MODEL::SST) { + + const su2double tke = solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint, 0); + + conv_numerics->SetTurbKineticEnergy(tke, tke); + + } + if (dynamic_grid) conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), geometry->nodes->GetGridVel(iPoint)); @@ -7461,6 +7571,9 @@ void CEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, bool gravity = (config->GetGravityForce()); bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST); + CVariable* turbNodes = nullptr; + if (tkeNeeded) turbNodes = solver_container[TURB_SOL]->GetNodes(); + auto *Normal = new su2double[nDim]; /*--- Loop over all the vertices on this boundary marker ---*/ @@ -7542,7 +7655,7 @@ void CEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, Velocity2 += Velocity[iDim]*Velocity[iDim]; } Energy = P_Exit/(Density*Gamma_Minus_One) + 0.5*Velocity2; - if (tkeNeeded) Energy += GetTke_Inf(); + if (tkeNeeded) Energy += turbNodes->GetSolution(iPoint,0); /*--- Conservative variables, using the derived quantities ---*/ V_outlet[0] = Pressure / ( Gas_Constant * Density); @@ -7556,6 +7669,11 @@ void CEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, /*--- Set various quantities in the solver class ---*/ conv_numerics->SetPrimitive(V_domain, V_outlet); + /*--- SST: both states contain the k of the node in their total enthalpy. ---*/ + if (config->GetKind_Turb_Model() == TURB_MODEL::SST) { + const su2double tke = solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint, 0); + conv_numerics->SetTurbKineticEnergy(tke, tke); + } if (dynamic_grid) conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), geometry->nodes->GetGridVel(iPoint)); @@ -7653,13 +7771,7 @@ void CEulerSolver::BC_Supersonic_Inlet(CGeometry *geometry, CSolver **solver_con /*--- Compute the energy from the specified state. ---*/ const su2double Velocity2 = GeometryToolbox::SquaredNorm(int(MAXNDIM), Velocity); - su2double Energy = Pressure / (Density * Gamma_Minus_One) + 0.5 * Velocity2; - if (tkeNeeded) { - const su2double* Turb_Properties = config->GetInlet_TurbVal(Marker_Tag); - const su2double Intensity = Turb_Properties[0]; - const su2double Tke = 3.0 / 2.0 * (Velocity2 * pow(Intensity, 2)); - Energy += Tke; - } + const su2double Energy = Pressure / (Density * Gamma_Minus_One) + 0.5 * Velocity2; /*--- Primitive variables, using the derived quantities. ---*/ @@ -7671,6 +7783,9 @@ void CEulerSolver::BC_Supersonic_Inlet(CGeometry *geometry, CSolver **solver_con for (unsigned short iDim = 0; iDim < nDim; iDim++) V_inlet[iDim+prim_idx.Velocity()] = Velocity[iDim]; + /*--- k of the boundary node, not the free-stream value, which decays right away in the domain (#1851). ---*/ + if (tkeNeeded) V_inlet[prim_idx.Enthalpy()] += solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0); + /*--- Current solution at this boundary node. ---*/ const auto* V_domain = nodes->GetPrimitive(iPoint); @@ -7685,6 +7800,11 @@ void CEulerSolver::BC_Supersonic_Inlet(CGeometry *geometry, CSolver **solver_con conv_numerics->SetNormal(Normal); conv_numerics->SetPrimitive(V_domain, V_inlet); + /*--- SST: both states contain the k of the node in their total enthalpy. ---*/ + if (config->GetKind_Turb_Model() == TURB_MODEL::SST) { + const su2double tke = solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint, 0); + conv_numerics->SetTurbKineticEnergy(tke, tke); + } if (dynamic_grid) conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), @@ -7766,6 +7886,11 @@ void CEulerSolver::BC_Supersonic_Outlet(CGeometry *geometry, CSolver **solver_co conv_numerics->SetNormal(Normal); conv_numerics->SetPrimitive(V_domain, V_outlet); + /*--- SST: both states contain the k of the node in their total enthalpy. ---*/ + if (config->GetKind_Turb_Model() == TURB_MODEL::SST) { + const su2double tke = solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint, 0); + conv_numerics->SetTurbKineticEnergy(tke, tke); + } if (dynamic_grid) conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), @@ -7982,7 +8107,7 @@ void CEulerSolver::BC_Engine_Inflow(CGeometry *geometry, CSolver **solver_contai } Energy = Inflow_Pressure/(Density*Gamma_Minus_One) + 0.5*Velocity2; - if (tkeNeeded) Energy += GetTke_Inf(); + if (tkeNeeded) Energy += solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0); // node k, see #1851 /*--- Conservative variables, using the derived quantities ---*/ @@ -7998,6 +8123,11 @@ void CEulerSolver::BC_Engine_Inflow(CGeometry *geometry, CSolver **solver_contai conv_numerics->SetNormal(Normal); conv_numerics->SetPrimitive(V_domain, V_inflow); + /*--- SST: both states contain the k of the node in their total enthalpy. ---*/ + if (config->GetKind_Turb_Model() == TURB_MODEL::SST) { + const su2double tke = solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint, 0); + conv_numerics->SetTurbKineticEnergy(tke, tke); + } /*--- Set grid movement ---*/ @@ -8147,6 +8277,8 @@ void CEulerSolver::BC_Engine_Exhaust(CGeometry *geometry, CSolver **solver_conta Velocity2 += Velocity[iDim]*Velocity[iDim]; } Energy = V_domain[nDim+3] - V_domain[nDim+1]/V_domain[nDim+2]; + /*--- SST: the total enthalpy of the domain contains k, the imposed total enthalpy does not. ---*/ + if (tkeNeeded) Energy -= solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0); Pressure = V_domain[nDim+1]; H_Exhaust = (Gamma*Gas_Constant/Gamma_Minus_One)*Exhaust_Temperature; SoundSpeed2 = Gamma*Pressure/Density; @@ -8219,7 +8351,7 @@ void CEulerSolver::BC_Engine_Exhaust(CGeometry *geometry, CSolver **solver_conta /*--- Using pressure, density, & velocity, compute the energy ---*/ Energy = Pressure/(Density*Gamma_Minus_One) + 0.5*Velocity2; - if (tkeNeeded) Energy += GetTke_Inf(); + if (tkeNeeded) Energy += solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0); // node k, see #1851 /*--- Primitive variables, using the derived quantities ---*/ @@ -8250,6 +8382,11 @@ void CEulerSolver::BC_Engine_Exhaust(CGeometry *geometry, CSolver **solver_conta conv_numerics->SetNormal(Normal); conv_numerics->SetPrimitive(V_domain, V_exhaust); + /*--- SST: both states contain the k of the node in their total enthalpy. ---*/ + if (config->GetKind_Turb_Model() == TURB_MODEL::SST) { + const su2double tke = solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint, 0); + conv_numerics->SetTurbKineticEnergy(tke, tke); + } /*--- Set grid movement ---*/ @@ -8511,7 +8648,7 @@ void CEulerSolver::BC_ActDisk(CGeometry *geometry, CSolver **solver_container, C Velocity2 += Velocity[iDim]*Velocity[iDim]; } Energy = Pressure/(Density*Gamma_Minus_One) + 0.5*Velocity2; - if (tkeNeeded) Energy += GetTke_Inf(); + if (tkeNeeded) Energy += solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0); // node k, see #1851 /*--- Conservative variables, using the derived quantities ---*/ @@ -8523,6 +8660,11 @@ void CEulerSolver::BC_ActDisk(CGeometry *geometry, CSolver **solver_container, C V_inlet[nDim+3] = Energy + Pressure/Density; V_inlet[nDim+4] = SoundSpeed; conv_numerics->SetPrimitive(V_domain, V_inlet); + /*--- SST: both states contain the k of the node in their total enthalpy. ---*/ + if (config->GetKind_Turb_Model() == TURB_MODEL::SST) { + const su2double tke = solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint, 0); + conv_numerics->SetTurbKineticEnergy(tke, tke); + } } @@ -8598,6 +8740,8 @@ void CEulerSolver::BC_ActDisk(CGeometry *geometry, CSolver **solver_container, C Velocity2 += Velocity[iDim]*Velocity[iDim]; } Energy = V_domain[nDim+3] - V_domain[nDim+1]/V_domain[nDim+2]; + /*--- SST: the total enthalpy of the domain contains k, the imposed total enthalpy does not. ---*/ + if (tkeNeeded) Energy -= solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0); Pressure = V_domain[nDim+1]; H_Total = (Gamma*Gas_Constant/Gamma_Minus_One)*T_Total; SoundSpeed2 = Gamma*Pressure/Density; @@ -8666,7 +8810,7 @@ void CEulerSolver::BC_ActDisk(CGeometry *geometry, CSolver **solver_container, C /*--- Using pressure, density, & velocity, compute the energy ---*/ Energy = Pressure/(Density*Gamma_Minus_One) + 0.5*Velocity2; - if (tkeNeeded) Energy += GetTke_Inf(); + if (tkeNeeded) Energy += solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0); // node k, see #1851 /*--- Primitive variables, using the derived quantities ---*/ @@ -8678,6 +8822,11 @@ void CEulerSolver::BC_ActDisk(CGeometry *geometry, CSolver **solver_container, C V_outlet[nDim+3] = Energy + Pressure/Density; V_outlet[nDim+4] = sqrt(SoundSpeed2); conv_numerics->SetPrimitive(V_domain, V_outlet); + /*--- SST: both states contain the k of the node in their total enthalpy. ---*/ + if (config->GetKind_Turb_Model() == TURB_MODEL::SST) { + const su2double tke = solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint, 0); + conv_numerics->SetTurbKineticEnergy(tke, tke); + } } @@ -8914,7 +9063,7 @@ void CEulerSolver::BC_ActDisk_VariableLoad(CGeometry *geometry, CSolver **solver Velocity2 += Velocity[iDim]*Velocity[iDim]; } Energy = Pressure/(Density*Gamma_Minus_One) + 0.5*Velocity2; - if (tkeNeeded) Energy += GetTke_Inf(); + if (tkeNeeded) Energy += solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0); // node k, see #1851 /*--- Conservative variables, using the derived quantities ---*/ @@ -8925,6 +9074,11 @@ void CEulerSolver::BC_ActDisk_VariableLoad(CGeometry *geometry, CSolver **solver V_inlet[nDim+3] = Energy + Pressure/Density; V_inlet[nDim+4] = SoundSpeed; conv_numerics->SetPrimitive(V_domain, V_inlet); + /*--- SST: both states contain the k of the node in their total enthalpy. ---*/ + if (config->GetKind_Turb_Model() == TURB_MODEL::SST) { + const su2double tke = solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint, 0); + conv_numerics->SetTurbKineticEnergy(tke, tke); + } } else { /*--- Acoustic Riemann invariant extrapolation form the interior domain. ---*/ @@ -8962,7 +9116,7 @@ void CEulerSolver::BC_ActDisk_VariableLoad(CGeometry *geometry, CSolver **solver /*--- Computation of the enthalpy, total energy, temperature and speed of sound. ---*/ H_out = H_in/Density_in + Fa/Density_out; Energy = H_out - Pressure_out/Density_out; - if (tkeNeeded) Energy += GetTke_Inf(); + if (tkeNeeded) Energy += solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0); // node k, see #1851 Temperature_out = (Energy-0.5*Velocity2/(pow(Density_out,2)))*(Gamma_Minus_One/Gas_Constant); SoS_out = sqrt(Gamma*Gas_Constant*Temperature_out); @@ -8976,6 +9130,11 @@ void CEulerSolver::BC_ActDisk_VariableLoad(CGeometry *geometry, CSolver **solver V_outlet[nDim+3] = H_out; V_outlet[nDim+4] = SoS_out; conv_numerics->SetPrimitive(V_domain, V_outlet); + /*--- SST: both states contain the k of the node in their total enthalpy. ---*/ + if (config->GetKind_Turb_Model() == TURB_MODEL::SST) { + const su2double tke = solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint, 0); + conv_numerics->SetTurbKineticEnergy(tke, tke); + } } /*--- Grid Movement (NOT TESTED!)---*/ diff --git a/SU2_CFD/src/solvers/CFEM_DG_EulerSolver.cpp b/SU2_CFD/src/solvers/CFEM_DG_EulerSolver.cpp index 0f5db8c0f8a7..96955c3b6555 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(); // Should it be non-dimensionalized for the Reynolds length? + + 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 512a99a97dc6..8feab8bc55be 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..41252b6d8c13 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(); // Should it be non-dimensionalized for the Reynolds length? + + 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/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); diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp index d6f98ecdee6b..9b52d9edc2f1 100644 --- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp @@ -124,6 +124,25 @@ 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, or (e.g. for MACH_NUMBER= 0 with grid motion) SST_OPTIONS= SUST with\n" + "positive SST_SUST_TKE_AMB and SST_SUST_OMEGA_AMB.", CURRENT_FUNCTION); + } + Solution_Inf[0] = kine_Inf; Solution_Inf[1] = omega_Inf; @@ -251,7 +270,7 @@ void CTurbSSTSolver::Postprocessing(CGeometry *geometry, CSolver **solver_contai const su2double kine = nodes->GetSolution(iPoint,0); const su2double omega = nodes->GetSolution(iPoint,1); - const auto& eddy_visc_var = sstParsedOptions.version == SST_OPTIONS::V1994 ? VorticityMag : StrainMag; + const su2double eddy_visc_var = sstParsedOptions.version == SST_OPTIONS::V1994 ? VorticityMag : StrainMag; const su2double muT = max(0.0, rho * a1 * kine / max(a1 * omega, eddy_visc_var * F2)); nodes->SetmuT(iPoint, muT); @@ -513,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 ---*/ @@ -580,7 +605,6 @@ void CTurbSSTSolver::SetTurbVars_WF(CGeometry *geometry, CSolver **solver_contai su2double solution[MAXNVAR] = {k, omega}; nodes->SetSolution_Old(iPoint_Neighbor,solution); - nodes->SetSolution(iPoint,solution); LinSysRes.SetBlock_Zero(iPoint_Neighbor); @@ -641,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]); @@ -777,15 +806,15 @@ void CTurbSSTSolver::BC_Inlet_Turbo(CGeometry *geometry, CSolver **solver_contai su2double rho = flowSolver->GetAverageDensity(val_marker, iSpan); su2double pressure = flowSolver->GetAveragePressure(val_marker, iSpan); - su2double kine = flowSolver->GetAverageKine(val_marker, iSpan); FluidModel->SetTDState_Prho(pressure, rho); su2double muLam = FluidModel->GetLaminarViscosity(); su2double VelMag2 = GeometryToolbox::SquaredNorm(nDim, flowSolver->GetAverageTurboVelocity(val_marker, iSpan)); + /*--- Imposed k from the intensity, omega from the imposed k and the viscosity ratio. ---*/ su2double kine_b = 3.0/2.0*(VelMag2*Intensity*Intensity); - su2double omega_b = rho*kine/(muLam*viscRatio); + su2double omega_b = rho*kine_b/(muLam*viscRatio); const su2double solution_j[] = {kine_b, omega_b}; 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 ---*/ diff --git a/TestCases/axisymmetric_rans/air_nozzle/air_nozzle_density_velocity.cfg b/TestCases/axisymmetric_rans/air_nozzle/air_nozzle_density_velocity.cfg new file mode 100644 index 000000000000..208ad6b63159 --- /dev/null +++ b/TestCases/axisymmetric_rans/air_nozzle/air_nozzle_density_velocity.cfg @@ -0,0 +1,100 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Axisymmetric supersonic converging-diverging air nozzle, % +% Riemann DENSITY_VELOCITY inlet (k in the boundary energy) % +% Author: Florian Dittmann % +% Date: 2021.12.02 % +% File Version 8.5.0 "Harrier" % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +% +SOLVER= RANS +KIND_TURB_MODEL= SST +RESTART_SOL= YES +AXISYMMETRIC= YES + +% -------------------- COMPRESSIBLE FREE-STREAM DEFINITION --------------------% +% +MACH_NUMBER= 1E-9 +INIT_OPTION= TD_CONDITIONS +FREESTREAM_OPTION= TEMPERATURE_FS +FREESTREAM_PRESSURE= 1400000 +FREESTREAM_TEMPERATURE= 373.15 +REF_DIMENSIONALIZATION= DIMENSIONAL + +% ---- IDEAL GAS, POLYTROPIC, VAN DER WAALS AND PENG ROBINSON CONSTANTS -------% +% +FLUID_MODEL= STANDARD_AIR + +% --------------------------- VISCOSITY MODEL ---------------------------------% +% +VISCOSITY_MODEL= CONSTANT_VISCOSITY +MU_CONSTANT= 1.716E-5 + +% --------------------------- THERMAL CONDUCTIVITY MODEL ----------------------% +% +CONDUCTIVITY_MODEL= CONSTANT_PRANDTL +PRANDTL_LAM= 0.72 +PRANDTL_TURB= 0.90 + +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +MARKER_HEATFLUX= ( WALL, 0.0 ) +MARKER_SYM= ( SYMMETRY ) +MARKER_RIEMANN= ( INFLOW, DENSITY_VELOCITY, 13.07, 100.0, 1.0, 0.0, 0.0, \ + OUTFLOW, STATIC_PRESSURE, 100000.0, 0.0, 0.0, 0.0, 0.0 ) +MARKER_MONITORING = (WALL) +MARKER_ANALYZE= (INFLOW) +% +% High free-stream turbulence, so that k is not negligible in the energy of the inlet state +FREESTREAM_TURBULENCEINTENSITY= 0.2 + +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +NUM_METHOD_GRAD= GREEN_GAUSS +CFL_NUMBER= 1000.0 +CFL_ADAPT= NO +MAX_DELTA_TIME= 1E6 +OBJECTIVE_FUNCTION= DRAG + +% ----------- SLOPE LIMITER AND DISSIPATION SENSOR DEFINITION -----------------% +% +MUSCL_FLOW= YES +SLOPE_LIMITER_FLOW= NONE + +% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% +% +LINEAR_SOLVER= FGMRES +LINEAR_SOLVER_PREC= ILU +LINEAR_SOLVER_ILU_FILL_IN= 0 +LINEAR_SOLVER_ERROR= 0.01 +LINEAR_SOLVER_ITER= 10 + +% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% +% +CONV_NUM_METHOD_FLOW= ROE +ENTROPY_FIX_COEFF= 0.1 +TIME_DISCRE_FLOW= EULER_IMPLICIT + +% -------------------- TURBULENT NUMERICAL METHOD DEFINITION ------------------% +% +CONV_NUM_METHOD_TURB= SCALAR_UPWIND +TIME_DISCRE_TURB= EULER_IMPLICIT +CFL_REDUCTION_TURB= 1.0 + +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +ITER= 3000 +CONV_RESIDUAL_MINVAL= -12 +CONV_STARTITER= 10 + +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +MESH_FILENAME= nozzle.su2 +SOLUTION_FILENAME= solution_flow +RESTART_FILENAME= restart_flow +OUTPUT_WRT_FREQ= 1000 +SCREEN_OUTPUT= (INNER_ITER, RMS_DENSITY, RMS_ENERGY, RMS_TKE, RMS_DISSIPATION, SURFACE_STATIC_TEMPERATURE) diff --git a/TestCases/hybrid_regression.py b/TestCases/hybrid_regression.py index 149f4a921607..b88f0988cfe9 100644 --- a/TestCases/hybrid_regression.py +++ b/TestCases/hybrid_regression.py @@ -149,7 +149,7 @@ def main(): poiseuille_profile.cfg_dir = "navierstokes/poiseuille" poiseuille_profile.cfg_file = "profile_poiseuille.cfg" poiseuille_profile.test_iter = 10 - poiseuille_profile.test_vals = [-12.004080, -7.638116, -0.000000, 2.089953] + poiseuille_profile.test_vals = [-12.004082, -7.637491, -0.000000, 2.089953] poiseuille_profile.test_vals_aarch64 = [-12.004276, -7.636719, -0.000000, 2.089953] test_list.append(poiseuille_profile) @@ -178,7 +178,7 @@ def main(): rae2822_sst.cfg_dir = "rans/rae2822" rae2822_sst.cfg_file = "turb_SST_RAE2822.cfg" rae2822_sst.test_iter = 20 - rae2822_sst.test_vals = [-1.501958, 5.889330, 0.635453, 0.021771, 100.000000] + rae2822_sst.test_vals = [-1.479160, 5.898617, 0.678431, 0.024642, 100.000000] test_list.append(rae2822_sst) # RAE2822 SST_SUST @@ -186,7 +186,7 @@ def main(): rae2822_sst_sust.cfg_dir = "rans/rae2822" rae2822_sst_sust.cfg_file = "turb_SST_SUST_RAE2822.cfg" rae2822_sst_sust.test_iter = 20 - rae2822_sst_sust.test_vals = [-2.465534, 5.844290, 0.496763, 0.041047] + rae2822_sst_sust.test_vals = [-2.942880, 5.054180, 0.509065, 0.043399] test_list.append(rae2822_sst_sust) # Flat plate @@ -194,7 +194,7 @@ def main(): turb_flatplate.cfg_dir = "rans/flatplate" turb_flatplate.cfg_file = "turb_SA_flatplate.cfg" turb_flatplate.test_iter = 20 - turb_flatplate.test_vals = [-0.187373, 0.003723, 10.000000, -1.443902] + turb_flatplate.test_vals = [-0.187373, 0.003723, 10.000000, -1.444696] test_list.append(turb_flatplate) # ONERA M6 Wing @@ -210,7 +210,7 @@ def main(): turb_naca0012_sa.cfg_dir = "rans/naca0012" turb_naca0012_sa.cfg_file = "turb_NACA0012_sa.cfg" turb_naca0012_sa.test_iter = 5 - turb_naca0012_sa.test_vals = [-12.038075, -16.332088, 1.080346, 0.018385, 20.000000, -2.873507, 0.000000, -14.250270, 0.000000] + turb_naca0012_sa.test_vals = [-12.038059, -16.332088, 1.080346, 0.018385, 20.000000, -2.873343, 0.000000, -14.250270, 0.000000] turb_naca0012_sa.test_vals_aarch64 = [-12.038091, -16.332090, 1.080346, 0.018385, 20.000000, -2.873236, 0.000000, -14.250271, 0.000000] test_list.append(turb_naca0012_sa) @@ -219,7 +219,7 @@ def main(): turb_naca0012_sst.cfg_dir = "rans/naca0012" turb_naca0012_sst.cfg_file = "turb_NACA0012_sst.cfg" turb_naca0012_sst.test_iter = 10 - turb_naca0012_sst.test_vals = [-12.094026, -15.250728, -5.906323, 1.070413, 0.015775, -2.854334, 0.000000] + turb_naca0012_sst.test_vals = [-5.948944, -10.295370, -3.739562, 1.069316, 0.015855, -3.048239, 0.000000] turb_naca0012_sst.test_vals_aarch64 = [-12.075928, -15.246732, -5.861249, 1.070036, 0.015841, -2.835263, 0] test_list.append(turb_naca0012_sst) @@ -228,7 +228,7 @@ def main(): turb_naca0012_sst_sust.cfg_dir = "rans/naca0012" turb_naca0012_sst_sust.cfg_file = "turb_NACA0012_sst_sust.cfg" turb_naca0012_sst_sust.test_iter = 10 - turb_naca0012_sst_sust.test_vals = [-12.080919, -14.837175, -5.732906, 1.000893, 0.019109, -2.119961] + turb_naca0012_sst_sust.test_vals = [-7.626467, -9.698545, -2.145340, 1.006188, 0.019328, -1.353713] turb_naca0012_sst_sust.test_vals_aarch64 = [-12.073210, -14.836724, -5.732627, 1.000050, 0.019144, -2.629689] test_list.append(turb_naca0012_sst_sust) @@ -237,7 +237,7 @@ def main(): turb_naca0012_sst_fixedvalues.cfg_dir = "rans/naca0012" turb_naca0012_sst_fixedvalues.cfg_file = "turb_NACA0012_sst_fixedvalues.cfg" turb_naca0012_sst_fixedvalues.test_iter = 10 - turb_naca0012_sst_fixedvalues.test_vals = [-5.192391, -10.448219, 0.773965, 1.022535, 0.040529, -2.383381] + turb_naca0012_sst_fixedvalues.test_vals = [-5.192392, -10.448218, 0.773965, 1.022535, 0.040529, -2.383336] test_list.append(turb_naca0012_sst_fixedvalues) # NACA0012 (SST, explicit Euler for flow and turbulence equations) @@ -245,7 +245,7 @@ def main(): turb_naca0012_sst_expliciteuler.cfg_dir = "rans/naca0012" turb_naca0012_sst_expliciteuler.cfg_file = "turb_NACA0012_sst_expliciteuler.cfg" turb_naca0012_sst_expliciteuler.test_iter = 10 - turb_naca0012_sst_expliciteuler.test_vals = [-3.532365, -3.157224, 3.743381, 1.124798, 0.501715, -16.000000] + turb_naca0012_sst_expliciteuler.test_vals = [-3.408564, -3.157224, 3.743383, 1.124798, 0.501715, -16.000000] test_list.append(turb_naca0012_sst_expliciteuler) # PROPELLER @@ -265,7 +265,7 @@ def main(): axi_rans_air_nozzle_restart.cfg_dir = "axisymmetric_rans/air_nozzle" axi_rans_air_nozzle_restart.cfg_file = "air_nozzle_restart.cfg" axi_rans_air_nozzle_restart.test_iter = 10 - axi_rans_air_nozzle_restart.test_vals = [-11.083068, -5.374686, -8.880093, -4.073548, 0.000000] + axi_rans_air_nozzle_restart.test_vals = [-4.875434, 0.608496, -4.040017, 0.185569, 0.000000] axi_rans_air_nozzle_restart.test_vals_aarch64 = [-14.140441, -9.154674, -10.886121, -5.806594, 0.000000] test_list.append(axi_rans_air_nozzle_restart) @@ -279,7 +279,7 @@ def main(): turb_naca0012_sst_restart_mg.cfg_file = "turb_NACA0012_sst_multigrid_restart.cfg" turb_naca0012_sst_restart_mg.test_iter = 20 turb_naca0012_sst_restart_mg.ntest_vals = 5 - turb_naca0012_sst_restart_mg.test_vals = [-6.551665, -5.057150, 0.830240, -0.008811, 0.078187] + turb_naca0012_sst_restart_mg.test_vals = [-3.689602, -5.056964, 0.830559, -0.008732, 0.078126] test_list.append(turb_naca0012_sst_restart_mg) ############################# @@ -291,7 +291,7 @@ def main(): turb_naca0012_1c.cfg_dir = "rans_uq/naca0012" turb_naca0012_1c.cfg_file = "turb_NACA0012_uq_1c.cfg" turb_naca0012_1c.test_iter = 10 - turb_naca0012_1c.test_vals = [-4.978722, 1.343530, 0.443909, -0.029125] + turb_naca0012_1c.test_vals = [-4.964342, 1.345656, 0.446527, -0.028331] turb_naca0012_1c.test_vals_aarch64 = [-4.976620, 1.345983, 0.433171, -0.033685] test_list.append(turb_naca0012_1c) @@ -300,7 +300,7 @@ def main(): turb_naca0012_2c.cfg_dir = "rans_uq/naca0012" turb_naca0012_2c.cfg_file = "turb_NACA0012_uq_2c.cfg" turb_naca0012_2c.test_iter = 10 - turb_naca0012_2c.test_vals = [-5.482844, 1.260870, 0.404376, -0.040361] + turb_naca0012_2c.test_vals = [-5.482844, 1.260888, 0.404947, -0.040154] turb_naca0012_2c.test_vals_aarch64 = [-5.485484, 1.263406, 0.411442, -0.040859] test_list.append(turb_naca0012_2c) @@ -309,7 +309,7 @@ def main(): turb_naca0012_3c.cfg_dir = "rans_uq/naca0012" turb_naca0012_3c.cfg_file = "turb_NACA0012_uq_3c.cfg" turb_naca0012_3c.test_iter = 10 - turb_naca0012_3c.test_vals = [-5.583730, 1.228732, 0.381967, -0.046233] + turb_naca0012_3c.test_vals = [-5.583731, 1.228717, 0.379824, -0.046992] turb_naca0012_3c.test_vals_aarch64 = [-5.583737, 1.232005, 0.390258, -0.046305] test_list.append(turb_naca0012_3c) @@ -318,7 +318,7 @@ def main(): turb_naca0012_p1c1.cfg_dir = "rans_uq/naca0012" turb_naca0012_p1c1.cfg_file = "turb_NACA0012_uq_p1c1.cfg" turb_naca0012_p1c1.test_iter = 10 - turb_naca0012_p1c1.test_vals = [-5.134181, 1.283464, 0.553971, 0.011916] + turb_naca0012_p1c1.test_vals = [-5.167299, 1.279420, 0.549536, 0.011043] turb_naca0012_p1c1.test_vals_aarch64 = [-5.114189, 1.285037, 0.406851, -0.043003] test_list.append(turb_naca0012_p1c1) @@ -327,7 +327,7 @@ def main(): turb_naca0012_p1c2.cfg_dir = "rans_uq/naca0012" turb_naca0012_p1c2.cfg_file = "turb_NACA0012_uq_p1c2.cfg" turb_naca0012_p1c2.test_iter = 10 - turb_naca0012_p1c2.test_vals = [-5.553988, 1.234031, 0.424168, -0.033501] + turb_naca0012_p1c2.test_vals = [-5.553265, 1.234113, 0.424941, -0.033285] turb_naca0012_p1c2.test_vals_aarch64 = [-5.548245, 1.236384, 0.381821, -0.050337] test_list.append(turb_naca0012_p1c2) @@ -460,7 +460,7 @@ def main(): inc_turb_naca0012_sst_sust.cfg_dir = "incomp_rans/naca0012" inc_turb_naca0012_sst_sust.cfg_file = "naca0012_SST_SUST.cfg" inc_turb_naca0012_sst_sust.test_iter = 20 - inc_turb_naca0012_sst_sust.test_vals = [-7.170018, 0.332638, 0.000002, 0.312117] + inc_turb_naca0012_sst_sust.test_vals = [-9.784873, 3.354565, 0.000008, 0.309470] test_list.append(inc_turb_naca0012_sst_sust) # Weakly coupled heat equation diff --git a/TestCases/hybrid_regression_AD.py b/TestCases/hybrid_regression_AD.py index b2086915f9df..f5371bbe5200 100644 --- a/TestCases/hybrid_regression_AD.py +++ b/TestCases/hybrid_regression_AD.py @@ -87,7 +87,7 @@ def main(): discadj_rans_naca0012_sst.cfg_dir = "disc_adj_rans/naca0012" discadj_rans_naca0012_sst.cfg_file = "turb_NACA0012_sst.cfg" discadj_rans_naca0012_sst.test_iter = 10 - discadj_rans_naca0012_sst.test_vals = [-2.201555, -0.175213, 3.045200, -0.041846] + discadj_rans_naca0012_sst.test_vals = [-2.201612, -0.175090, 3.046200, -0.041861] discadj_rans_naca0012_sst.test_vals_aarch64 = [-2.201855, -0.172443, 3.043400, -0.041820] test_list.append(discadj_rans_naca0012_sst) diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 36766da22528..0548c75e0134 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -438,7 +438,7 @@ def main(): rae2822_sst.cfg_dir = "rans/rae2822" rae2822_sst.cfg_file = "turb_SST_RAE2822.cfg" rae2822_sst.test_iter = 20 - rae2822_sst.test_vals = [-1.495621, 5.875629, 0.638353, 0.020410, 100.000000] + rae2822_sst.test_vals = [-1.505224, 5.876281, 0.639383, 0.020780, 100.000000] test_list.append(rae2822_sst) # RAE2822 SST_SUST @@ -446,7 +446,7 @@ def main(): rae2822_sst_sust.cfg_dir = "rans/rae2822" rae2822_sst_sust.cfg_file = "turb_SST_SUST_RAE2822.cfg" rae2822_sst_sust.test_iter = 20 - rae2822_sst_sust.test_vals = [-2.465087, 5.851605, 0.496009, 0.041100] + rae2822_sst_sust.test_vals = [-2.939568, 5.056410, 0.509111, 0.042674] test_list.append(rae2822_sst_sust) # Flat plate @@ -486,7 +486,7 @@ def main(): turb_flatplate_sst_roughBCKnopp.cfg_dir = "rans/flatplate/roughness/bc_knopp" turb_flatplate_sst_roughBCKnopp.cfg_file = "turb_SST_flatplate_roughBCKnopp.cfg" turb_flatplate_sst_roughBCKnopp.test_iter = 10 - turb_flatplate_sst_roughBCKnopp.test_vals = [-5.058634, -2.460850, -2.847064, 0.447200, -2.595042, 1.497149, -0.188079, 0.004571] + turb_flatplate_sst_roughBCKnopp.test_vals = [-4.843998, -2.247185, -2.637470, 0.662708, -2.565774, 1.507195, -0.188018, 0.004571] test_list.append(turb_flatplate_sst_roughBCKnopp) # FLAT PLATE, ROUGHNESS BC AUPOIX SST @@ -494,7 +494,7 @@ def main(): turb_flatplate_sst_roughBCAupoix.cfg_dir = "rans/flatplate/roughness/bc_aupoix" turb_flatplate_sst_roughBCAupoix.cfg_file = "turb_SST_flatplate_roughBCAupoix.cfg" turb_flatplate_sst_roughBCAupoix.test_iter = 10 - turb_flatplate_sst_roughBCAupoix.test_vals = [-5.278728, -2.302982, -2.890419, 0.227585, -1.393786, 3.192601, -0.188729, 0.0071821] + turb_flatplate_sst_roughBCAupoix.test_vals = [-5.069751, -2.238097, -2.717921, 0.438657, -1.393740, 3.192608, -0.188736, 0.007182] test_list.append(turb_flatplate_sst_roughBCAupoix) # ONERA M6 Wing @@ -540,7 +540,7 @@ def main(): turb_naca0012_sst.cfg_dir = "rans/naca0012" turb_naca0012_sst.cfg_file = "turb_NACA0012_sst.cfg" turb_naca0012_sst.test_iter = 10 - turb_naca0012_sst.test_vals = [-12.094692, -15.251093, -5.906365, 1.070413, 0.015775, -2.376043, 0.000000] + turb_naca0012_sst.test_vals = [-5.949126, -10.294721, -3.763091, 1.069311, 0.015852, -3.462062, 0.000000] turb_naca0012_sst.test_vals_aarch64 = [-12.075620, -15.246688, -5.861276, 1.070036, 0.015841, -1.991001, 0.000000] turb_naca0012_sst.timeout = 3200 test_list.append(turb_naca0012_sst) @@ -550,7 +550,7 @@ def main(): turb_naca0012_sst_sust.cfg_dir = "rans/naca0012" turb_naca0012_sst_sust.cfg_file = "turb_NACA0012_sst_sust.cfg" turb_naca0012_sst_sust.test_iter = 10 - turb_naca0012_sst_sust.test_vals = [-12.081966, -14.837177, -5.733436, 1.000893, 0.019109, -2.241123] + turb_naca0012_sst_sust.test_vals = [-7.620992, -9.772595, -2.145273, 1.006174, 0.019314, -1.562485] turb_naca0012_sst_sust.test_vals_aarch64 = [-12.073964, -14.836726, -5.732390, 1.000050, 0.019144, -2.229074] turb_naca0012_sst_sust.timeout = 3200 test_list.append(turb_naca0012_sst_sust) @@ -560,7 +560,7 @@ def main(): turb_naca0012_sst_2003_Vm.cfg_dir = "rans/naca0012" turb_naca0012_sst_2003_Vm.cfg_file = "turb_NACA0012_sst_2003-Vm.cfg" turb_naca0012_sst_2003_Vm.test_iter = 10 - turb_naca0012_sst_2003_Vm.test_vals = [-10.168739, -3.668709, 1.060563, 0.019147, -2.282800] + turb_naca0012_sst_2003_Vm.test_vals = [-10.082677, -3.584114, 1.059491, 0.019212, -3.091163] turb_naca0012_sst_2003_Vm.timeout = 3200 test_list.append(turb_naca0012_sst_2003_Vm) @@ -569,7 +569,7 @@ def main(): turb_naca0012_sst_1994_KLm.cfg_dir = "rans/naca0012" turb_naca0012_sst_1994_KLm.cfg_file = "turb_NACA0012_sst_1994-KLm.cfg" turb_naca0012_sst_1994_KLm.test_iter = 10 - turb_naca0012_sst_1994_KLm.test_vals = [-10.390871, -3.701584, 1.062164, 0.019076, -2.426123] + turb_naca0012_sst_1994_KLm.test_vals = [-10.202368, -3.596058, 1.061073, 0.019145, -3.189038] turb_naca0012_sst_1994_KLm.timeout = 3200 test_list.append(turb_naca0012_sst_1994_KLm) @@ -587,7 +587,7 @@ def main(): turb_naca0012_sst_expliciteuler.cfg_dir = "rans/naca0012" turb_naca0012_sst_expliciteuler.cfg_file = "turb_NACA0012_sst_expliciteuler.cfg" turb_naca0012_sst_expliciteuler.test_iter = 10 - turb_naca0012_sst_expliciteuler.test_vals = [-3.532365, -3.157224, 3.743381, 1.124798, 0.501715, -16.000000] + turb_naca0012_sst_expliciteuler.test_vals = [-3.408564, -3.157224, 3.743383, 1.124798, 0.501715, -16.000000] turb_naca0012_sst_expliciteuler.timeout = 3200 test_list.append(turb_naca0012_sst_expliciteuler) @@ -619,11 +619,19 @@ def main(): axi_rans_air_nozzle_restart.cfg_dir = "axisymmetric_rans/air_nozzle" axi_rans_air_nozzle_restart.cfg_file = "air_nozzle_restart.cfg" axi_rans_air_nozzle_restart.test_iter = 10 - axi_rans_air_nozzle_restart.test_vals = [-11.056236, -5.334116, -8.842310, -4.067917, 0.000000] + axi_rans_air_nozzle_restart.test_vals = [-4.875217, 0.608350, -4.006405, 0.187713, 0.000000] axi_rans_air_nozzle_restart.test_vals_aarch64 = [-14.143310, -9.163287, -10.858232, -5.787715, 0.000000] axi_rans_air_nozzle_restart.tol = 0.0001 test_list.append(axi_rans_air_nozzle_restart) + # Riemann DENSITY_VELOCITY inlet with SST: k of the node in the energy of the boundary state + axi_rans_air_nozzle_density_velocity = TestCase('axi_rans_air_nozzle_density_velocity') + axi_rans_air_nozzle_density_velocity.cfg_dir = "axisymmetric_rans/air_nozzle" + axi_rans_air_nozzle_density_velocity.cfg_file = "air_nozzle_density_velocity.cfg" + axi_rans_air_nozzle_density_velocity.test_iter = 10 + axi_rans_air_nozzle_density_velocity.test_vals = [-1.377739, 4.590218, 0.124013, 4.863032, 698.28, 698.28] + test_list.append(axi_rans_air_nozzle_density_velocity) + ################################# ## Compressible RANS Restart ### ################################# @@ -634,7 +642,7 @@ def main(): turb_naca0012_sst_restart_mg.cfg_file = "turb_NACA0012_sst_multigrid_restart.cfg" turb_naca0012_sst_restart_mg.test_iter = 20 turb_naca0012_sst_restart_mg.ntest_vals = 5 - turb_naca0012_sst_restart_mg.test_vals = [-6.560190, -5.057149, 0.830238, -0.008812, 0.078167] + turb_naca0012_sst_restart_mg.test_vals = [-3.692233, -5.056998, 0.830505, -0.008795, 0.078083] turb_naca0012_sst_restart_mg.timeout = 3200 turb_naca0012_sst_restart_mg.tol = 0.000001 test_list.append(turb_naca0012_sst_restart_mg) @@ -797,7 +805,7 @@ def main(): inc_turb_naca0012_sst_sust.cfg_dir = "incomp_rans/naca0012" inc_turb_naca0012_sst_sust.cfg_file = "naca0012_SST_SUST.cfg" inc_turb_naca0012_sst_sust.test_iter = 20 - inc_turb_naca0012_sst_sust.test_vals = [-7.169837, 0.332730, -0.000001, 0.312131] + inc_turb_naca0012_sst_sust.test_vals = [-9.784857, 3.354566, -0.000002, 0.309484] test_list.append(inc_turb_naca0012_sst_sust) # Flat plate, pressure-based @@ -1041,7 +1049,7 @@ def main(): turb_naca0012_1c.cfg_dir = "rans_uq/naca0012" turb_naca0012_1c.cfg_file = "turb_NACA0012_uq_1c.cfg" turb_naca0012_1c.test_iter = 10 - turb_naca0012_1c.test_vals = [-4.979554, 1.344291, 0.666805, 0.010247] + turb_naca0012_1c.test_vals = [-4.965784, 1.346107, 0.672718, 0.011889] turb_naca0012_1c.test_vals_aarch64 = [-4.981036, 1.345868, 0.673232, 0.010091] test_list.append(turb_naca0012_1c) @@ -1050,7 +1058,7 @@ def main(): turb_naca0012_2c.cfg_dir = "rans_uq/naca0012" turb_naca0012_2c.cfg_file = "turb_NACA0012_uq_2c.cfg" turb_naca0012_2c.test_iter = 10 - turb_naca0012_2c.test_vals = [-5.482691, 1.262918, 0.529450, -0.023537] + turb_naca0012_2c.test_vals = [-5.482675, 1.262934, 0.535624, -0.021772] turb_naca0012_2c.test_vals_aarch64 = [-5.484365, 1.264701, 0.501741, -0.033109] test_list.append(turb_naca0012_2c) @@ -1059,7 +1067,7 @@ def main(): turb_naca0012_3c.cfg_dir = "rans_uq/naca0012" turb_naca0012_3c.cfg_file = "turb_NACA0012_uq_3c.cfg" turb_naca0012_3c.test_iter = 10 - turb_naca0012_3c.test_vals = [-5.583617, 1.229594, 0.463384, -0.035444] + turb_naca0012_3c.test_vals = [-5.583618, 1.229581, 0.465402, -0.035057] test_list.append(turb_naca0012_3c) # NACA0012 p1c1 @@ -1067,7 +1075,7 @@ def main(): turb_naca0012_p1c1.cfg_dir = "rans_uq/naca0012" turb_naca0012_p1c1.cfg_file = "turb_NACA0012_uq_p1c1.cfg" turb_naca0012_p1c1.test_iter = 10 - turb_naca0012_p1c1.test_vals = [-5.129555, 1.283892, 0.808074, 0.047784] + turb_naca0012_p1c1.test_vals = [-5.164490, 1.279381, 0.810659, 0.048539] turb_naca0012_p1c1.test_vals_aarch64 = [-5.122100, 1.284478, 0.608744, -0.008593] test_list.append(turb_naca0012_p1c1) @@ -1076,7 +1084,7 @@ def main(): turb_naca0012_p1c2.cfg_dir = "rans_uq/naca0012" turb_naca0012_p1c2.cfg_file = "turb_NACA0012_uq_p1c2.cfg" turb_naca0012_p1c2.test_iter = 10 - turb_naca0012_p1c2.test_vals = [-5.554045, 1.234446, 0.603731, -0.008717] + turb_naca0012_p1c2.test_vals = [-5.553369, 1.234366, 0.608821, -0.007429] test_list.append(turb_naca0012_p1c2) ###################################### @@ -1137,7 +1145,7 @@ def main(): square_cylinder.cfg_dir = "unsteady/square_cylinder" square_cylinder.cfg_file = "turb_square.cfg" square_cylinder.test_iter = 3 - square_cylinder.test_vals = [-1.175979, 0.062203, 1.399352, 2.219226, 1.399298, 2.217472, 0.000000] + square_cylinder.test_vals = [-1.175978, 0.062220, 1.399352, 2.219226, 1.399298, 2.217472, 0.000000] square_cylinder.unsteady = True test_list.append(square_cylinder) @@ -1231,7 +1239,7 @@ def main(): datadriven_fluidModel.cfg_dir = "nicf/datadriven" datadriven_fluidModel.cfg_file = "datadriven_nozzle.cfg" datadriven_fluidModel.test_iter = 50 - datadriven_fluidModel.test_vals = [-4.879421, -2.353319, -2.344004, 0.390613, -1.559800, 1.322761] + datadriven_fluidModel.test_vals = [-4.788556, -1.721315, -2.318534, 0.759228, -1.558870, 1.321951] test_list.append(datadriven_fluidModel) ###################################### @@ -1360,7 +1368,7 @@ def main(): bars_SST_2D.cfg_dir = "sliding_interface/bars_SST_2D" bars_SST_2D.cfg_file = "bars.cfg" bars_SST_2D.test_iter = 13 - bars_SST_2D.test_vals = [13.000000, -0.457912, -1.541047] + bars_SST_2D.test_vals = [13.000000, -0.457937, -1.540952] bars_SST_2D.multizone = True test_list.append(bars_SST_2D) @@ -1560,7 +1568,7 @@ def main(): pywrapper_turb_naca0012_sst.cfg_dir = "rans/naca0012" pywrapper_turb_naca0012_sst.cfg_file = "turb_NACA0012_sst.cfg" pywrapper_turb_naca0012_sst.test_iter = 10 - pywrapper_turb_naca0012_sst.test_vals = [-12.094692, -15.251093, -5.906365, 1.070413, 0.015775, -2.376043, 0.000000] + pywrapper_turb_naca0012_sst.test_vals = [-5.949126, -10.294721, -3.763091, 1.069311, 0.015852, -3.462062, 0.000000] pywrapper_turb_naca0012_sst.test_vals_aarch64 = [-12.075620, -15.246688, -5.861276, 1.070036, 0.015841, -1.991001, 0.000000] pywrapper_turb_naca0012_sst.command = TestCase.Command("mpirun -np 2", "SU2_CFD.py", "--parallel -f") pywrapper_turb_naca0012_sst.timeout = 3200 @@ -1571,7 +1579,7 @@ def main(): pywrapper_square_cylinder.cfg_dir = "unsteady/square_cylinder" pywrapper_square_cylinder.cfg_file = "turb_square.cfg" pywrapper_square_cylinder.test_iter = 10 - pywrapper_square_cylinder.test_vals = [-1.178521, -0.349772, 1.401402, 2.359105, 1.401729, 2.300974, 0.000000] + pywrapper_square_cylinder.test_vals = [-1.178520, -0.349721, 1.401402, 2.359105, 1.401729, 2.300974, 0.000000] pywrapper_square_cylinder.command = TestCase.Command("mpirun -np 2", "SU2_CFD.py", "--parallel -f") pywrapper_square_cylinder.unsteady = True test_list.append(pywrapper_square_cylinder) @@ -1611,7 +1619,7 @@ def main(): pywrapper_unsteadyFSI.cfg_dir = "py_wrapper/dyn_fsi" pywrapper_unsteadyFSI.cfg_file = "config.cfg" pywrapper_unsteadyFSI.test_iter = 4 - pywrapper_unsteadyFSI.test_vals = [0, 31, 5, 47, -1.756677, -2.828286, -7.638545, -6.863930, 0.000156] + pywrapper_unsteadyFSI.test_vals = [0.000000, 31.000000, 5.000000, 48.000000, -1.756677, -2.828286, -7.638545, -6.863930, 0.000156] pywrapper_unsteadyFSI.command = TestCase.Command("mpirun -np 2", "python", "run.py") pywrapper_unsteadyFSI.unsteady = True pywrapper_unsteadyFSI.multizone = True @@ -1622,7 +1630,7 @@ def main(): pywrapper_unsteadyCHT.cfg_dir = "py_wrapper/flatPlate_unsteady_CHT" pywrapper_unsteadyCHT.cfg_file = "unsteady_CHT_FlatPlate_Conf.cfg" pywrapper_unsteadyCHT.test_iter = 5 - pywrapper_unsteadyCHT.test_vals = [-1.614168, 2.260203, -0.008065, 0.201068] + pywrapper_unsteadyCHT.test_vals = [-1.614169, 2.260208, -0.008127, 0.201931] pywrapper_unsteadyCHT.command = TestCase.Command("mpirun -np 2", "python", "launch_unsteady_CHT_FlatPlate.py --parallel -f") pywrapper_unsteadyCHT.unsteady = True test_list.append(pywrapper_unsteadyCHT) diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py index 0b5d8724f62b..24cd6087573a 100644 --- a/TestCases/parallel_regression_AD.py +++ b/TestCases/parallel_regression_AD.py @@ -91,7 +91,7 @@ def main(): discadj_rans_naca0012_sst.cfg_dir = "disc_adj_rans/naca0012" discadj_rans_naca0012_sst.cfg_file = "turb_NACA0012_sst.cfg" discadj_rans_naca0012_sst.test_iter = 10 - discadj_rans_naca0012_sst.test_vals = [-2.226820, -0.256327, -2.299900, -0.003171] + discadj_rans_naca0012_sst.test_vals = [-2.226859, -0.256181, -2.298600, -0.003183] discadj_rans_naca0012_sst.test_vals_aarch64 = [-2.226852, -0.253100, -2.299300, -0.003164] test_list.append(discadj_rans_naca0012_sst) @@ -150,7 +150,7 @@ def main(): discadj_axisymmetric_rans_nozzle.cfg_dir = "axisymmetric_rans/air_nozzle" discadj_axisymmetric_rans_nozzle.cfg_file = "air_nozzle_restart.cfg" discadj_axisymmetric_rans_nozzle.test_iter = 10 - discadj_axisymmetric_rans_nozzle.test_vals = [9.909657, 5.078045, 7.129068, 2.490955] + discadj_axisymmetric_rans_nozzle.test_vals = [9.909640, 5.078072, 7.132311, 2.490930] discadj_axisymmetric_rans_nozzle.no_restart = True test_list.append(discadj_axisymmetric_rans_nozzle) diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index ab1de7fe1431..7c54aec582c5 100755 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -208,7 +208,7 @@ def main(): poiseuille_profile.cfg_dir = "navierstokes/poiseuille" poiseuille_profile.cfg_file = "profile_poiseuille.cfg" poiseuille_profile.test_iter = 10 - poiseuille_profile.test_vals = [-12.003115, -7.626023, -0.000000, 2.089953] + poiseuille_profile.test_vals = [-12.003113, -7.626005, -0.000000, 2.089953] poiseuille_profile.test_vals_aarch64 = [-12.009012, -7.262299, -0.000000, 2.089953] #last 4 columns test_list.append(poiseuille_profile) @@ -236,7 +236,7 @@ def main(): rae2822_sst.cfg_dir = "rans/rae2822" rae2822_sst.cfg_file = "turb_SST_RAE2822.cfg" rae2822_sst.test_iter = 20 - rae2822_sst.test_vals = [-1.745363, -1.484660, 5.883796, 0.579910, 0.017299, 100.000000] + rae2822_sst.test_vals = [-1.682306, -1.473692, 5.871878, 0.593725, 0.020638, 100.000000] test_list.append(rae2822_sst) # RAE2822 SST_SUST @@ -244,7 +244,7 @@ def main(): rae2822_sst_sust.cfg_dir = "rans/rae2822" rae2822_sst_sust.cfg_file = "turb_SST_SUST_RAE2822.cfg" rae2822_sst_sust.test_iter = 20 - rae2822_sst_sust.test_vals = [-2.464319, 5.852571, 0.490309, 0.042521] + rae2822_sst_sust.test_vals = [-2.941125, 5.058633, 0.507642, 0.044280] test_list.append(rae2822_sst_sust) # Flat plate @@ -260,7 +260,7 @@ def main(): turb_wallfunction_flatplate_sst.cfg_dir = "wallfunctions/flatplate/compressible_SST" turb_wallfunction_flatplate_sst.cfg_file = "turb_SST_flatplate.cfg" turb_wallfunction_flatplate_sst.test_iter = 10 - turb_wallfunction_flatplate_sst.test_vals = [-4.036300, -1.916264, -1.821790, 1.443878, -1.579554, 1.521535, 10.000000, -2.351460, 0.030011, 0.002409] + turb_wallfunction_flatplate_sst.test_vals = [-4.037523, -1.925464, -1.823711, 1.442524, -1.579599, 1.519008, 10.000000, -2.403318, 0.030004, 0.002409] test_list.append(turb_wallfunction_flatplate_sst) # FLAT PLATE, ROUGHNESS BC WILCOX2006 SST @@ -268,7 +268,7 @@ def main(): turb_flatplate_sst_roughBCWilcox2006.cfg_dir = "rans/flatplate/roughness/bc_wilcox2006" turb_flatplate_sst_roughBCWilcox2006.cfg_file = "turb_SST_flatplate_roughBCWilcox2006.cfg" turb_flatplate_sst_roughBCWilcox2006.test_iter = 10 - turb_flatplate_sst_roughBCWilcox2006.test_vals = [-5.117900, -2.534224, -2.904279, 0.381710, -3.100346, 1.180161, -0.188797, 0.004029] + turb_flatplate_sst_roughBCWilcox2006.test_vals = [-5.235113, -2.562052, -2.931476, 0.265993, -3.126054, 1.259616, -0.188832, 0.004029] test_list.append(turb_flatplate_sst_roughBCWilcox2006) # FLAT PLATE, WALL FUNCTIONS, COMPRESSIBLE SA @@ -303,7 +303,7 @@ def main(): turb_naca0012_sst.cfg_dir = "rans/naca0012" turb_naca0012_sst.cfg_file = "turb_NACA0012_sst.cfg" turb_naca0012_sst.test_iter = 10 - turb_naca0012_sst.test_vals = [-12.094431, -15.251082, -5.906366, 1.070413, 0.015775, -3.178469, 0.000000] + turb_naca0012_sst.test_vals = [-5.948762, -10.293853, -3.762064, 1.069291, 0.015853, -3.348638, 0.000000] turb_naca0012_sst.test_vals_aarch64 = [-12.076068, -15.246740, -5.861280, 1.070036, 0.015841, -3.297854, 0.000000] turb_naca0012_sst.timeout = 3200 test_list.append(turb_naca0012_sst) @@ -326,7 +326,7 @@ def main(): turb_naca0012_sst_2003m.cfg_dir = "rans/naca0012" turb_naca0012_sst_2003m.cfg_file = "turb_NACA0012_sst_2003m.cfg" turb_naca0012_sst_2003m.test_iter = 10 - turb_naca0012_sst_2003m.test_vals = [-7.129622, -10.410728, -3.699440, 1.061671, 0.019047, -2.435494, 0.000000] + turb_naca0012_sst_2003m.test_vals = [-5.953752, -10.211152, -3.594485, 1.060530, 0.019121, -3.188773, 0.000000] turb_naca0012_sst_2003m.timeout = 3200 test_list.append(turb_naca0012_sst_2003m) @@ -335,7 +335,7 @@ def main(): turb_naca0012_sst_sust_restart.cfg_dir = "rans/naca0012" turb_naca0012_sst_sust_restart.cfg_file = "turb_NACA0012_sst_sust.cfg" turb_naca0012_sst_sust_restart.test_iter = 10 - turb_naca0012_sst_sust_restart.test_vals = [-12.080496, -14.837169, -5.733461, 1.000893, 0.019109, -2.634008] + turb_naca0012_sst_sust_restart.test_vals = [-7.626856, -9.773425, -2.145353, 1.006177, 0.019310, -1.684024] turb_naca0012_sst_sust_restart.test_vals_aarch64 = [-12.074189, -14.836725, -5.732398, 1.000050, 0.019144, -3.315560] turb_naca0012_sst_sust_restart.timeout = 3200 test_list.append(turb_naca0012_sst_sust_restart) @@ -367,17 +367,25 @@ def main(): axi_rans_air_nozzle_restart.cfg_dir = "axisymmetric_rans/air_nozzle" axi_rans_air_nozzle_restart.cfg_file = "air_nozzle_restart.cfg" axi_rans_air_nozzle_restart.test_iter = 10 - axi_rans_air_nozzle_restart.test_vals = [-11.054281, -5.328905, -8.835591, -4.056830, 0.000000] + axi_rans_air_nozzle_restart.test_vals = [-4.875503, 0.608561, -4.031599, 0.232810, 0.000000] axi_rans_air_nozzle_restart.test_vals_aarch64 = [-14.143715, -9.170705, -10.848554, -5.776746, 0.000000] axi_rans_air_nozzle_restart.tol = 0.0001 test_list.append(axi_rans_air_nozzle_restart) + # Riemann DENSITY_VELOCITY inlet with SST: k of the node in the energy of the boundary state + axi_rans_air_nozzle_density_velocity = TestCase('axi_rans_air_nozzle_density_velocity') + axi_rans_air_nozzle_density_velocity.cfg_dir = "axisymmetric_rans/air_nozzle" + axi_rans_air_nozzle_density_velocity.cfg_file = "air_nozzle_density_velocity.cfg" + axi_rans_air_nozzle_density_velocity.test_iter = 10 + axi_rans_air_nozzle_density_velocity.test_vals = [-1.372885, 4.597063, 0.099806, 4.871497, 697.47, 697.47] + test_list.append(axi_rans_air_nozzle_density_velocity) + # Axisymmetric air nozzle species axi_rans_air_nozzle_species = TestCase('axi_rans_air_nozzle_species') axi_rans_air_nozzle_species.cfg_dir = "axisymmetric_rans/air_nozzle" axi_rans_air_nozzle_species.cfg_file = "air_nozzle_species.cfg" axi_rans_air_nozzle_species.test_iter = 10 - axi_rans_air_nozzle_species.test_vals = [-1.690665, 3.882506, -2.928702, 5.760933, -3.144560, 0.000000] + axi_rans_air_nozzle_species.test_vals = [-1.690622, 3.882548, -2.928703, 5.760932, -3.146519, 0.000000] axi_rans_air_nozzle_species.tol = 0.0001 test_list.append(axi_rans_air_nozzle_species) @@ -391,7 +399,7 @@ def main(): turb_naca0012_sst_restart_mg.cfg_file = "turb_NACA0012_sst_multigrid_restart.cfg" turb_naca0012_sst_restart_mg.test_iter = 50 turb_naca0012_sst_restart_mg.ntest_vals = 5 - turb_naca0012_sst_restart_mg.test_vals = [-6.570974, -5.081421, 0.810883, -0.008830, 0.077967] + turb_naca0012_sst_restart_mg.test_vals = [-3.649967, -5.081630, 0.810585, -0.009208, 0.078506] turb_naca0012_sst_restart_mg.timeout = 3200 turb_naca0012_sst_restart_mg.tol = 0.000001 test_list.append(turb_naca0012_sst_restart_mg) @@ -508,7 +516,7 @@ def main(): inc_turb_naca0012_sst_sust.cfg_dir = "incomp_rans/naca0012" inc_turb_naca0012_sst_sust.cfg_file = "naca0012_SST_SUST.cfg" inc_turb_naca0012_sst_sust.test_iter = 20 - inc_turb_naca0012_sst_sust.test_vals = [-7.169704, 0.332779, 0.000021, 0.312114] + inc_turb_naca0012_sst_sust.test_vals = [-9.784804, 3.354566, 0.000023, 0.309468] test_list.append(inc_turb_naca0012_sst_sust) # Flat plate, pressure-based @@ -744,7 +752,7 @@ def main(): turb_naca0012_1c.cfg_dir = "rans_uq/naca0012" turb_naca0012_1c.cfg_file = "turb_NACA0012_uq_1c.cfg" turb_naca0012_1c.test_iter = 10 - turb_naca0012_1c.test_vals = [-4.981181, 1.343583, 0.597121, 0.017223] + turb_naca0012_1c.test_vals = [-4.963256, 1.346299, 0.604691, 0.019761] turb_naca0012_1c.test_vals_aarch64 = [-4.992791, 1.342873, 0.557941, 0.003269] test_list.append(turb_naca0012_1c) @@ -753,7 +761,7 @@ def main(): turb_naca0012_2c.cfg_dir = "rans_uq/naca0012" turb_naca0012_2c.cfg_file = "turb_NACA0012_uq_2c.cfg" turb_naca0012_2c.test_iter = 10 - turb_naca0012_2c.test_vals = [-5.482907, 1.262063, 0.459671, -0.026410] + turb_naca0012_2c.test_vals = [-5.482902, 1.262062, 0.464303, -0.024986] test_list.append(turb_naca0012_2c) # NACA0012 3c @@ -761,7 +769,7 @@ def main(): turb_naca0012_3c.cfg_dir = "rans_uq/naca0012" turb_naca0012_3c.cfg_file = "turb_NACA0012_uq_3c.cfg" turb_naca0012_3c.test_iter = 10 - turb_naca0012_3c.test_vals = [-5.583768, 1.229824, 0.426258, -0.033150] + turb_naca0012_3c.test_vals = [-5.583769, 1.229839, 0.426520, -0.033213] test_list.append(turb_naca0012_3c) # NACA0012 p1c1 @@ -769,7 +777,7 @@ def main(): turb_naca0012_p1c1.cfg_dir = "rans_uq/naca0012" turb_naca0012_p1c1.cfg_file = "turb_NACA0012_uq_p1c1.cfg" turb_naca0012_p1c1.test_iter = 10 - turb_naca0012_p1c1.test_vals = [-5.127701, 1.284875, 0.779301, 0.086111] + turb_naca0012_p1c1.test_vals = [-5.163879, 1.280584, 0.794994, 0.091481] turb_naca0012_p1c1.test_vals_aarch64 = [-5.119942, 1.283920, 0.486264, -0.021518] test_list.append(turb_naca0012_p1c1) @@ -778,7 +786,7 @@ def main(): turb_naca0012_p1c2.cfg_dir = "rans_uq/naca0012" turb_naca0012_p1c2.cfg_file = "turb_NACA0012_uq_p1c2.cfg" turb_naca0012_p1c2.test_iter = 10 - turb_naca0012_p1c2.test_vals = [-5.553987, 1.235071, 0.521197, -0.004487] + turb_naca0012_p1c2.test_vals = [-5.553241, 1.235093, 0.577763, 0.018264] test_list.append(turb_naca0012_p1c2) ###################################### @@ -838,7 +846,7 @@ def main(): square_cylinder.cfg_dir = "unsteady/square_cylinder" square_cylinder.cfg_file = "turb_square.cfg" square_cylinder.test_iter = 3 - square_cylinder.test_vals = [-2.560678, -1.175979, 0.062202, 1.399351, 2.219234, 1.399297, 2.217479, 0.000000] + square_cylinder.test_vals = [-2.560698, -1.175978, 0.062220, 1.399351, 2.219234, 1.399297, 2.217479, 0.000000] square_cylinder.unsteady = True test_list.append(square_cylinder) @@ -1078,7 +1086,7 @@ def main(): bars_SST_2D.cfg_dir = "sliding_interface/bars_SST_2D" bars_SST_2D.cfg_file = "bars.cfg" bars_SST_2D.test_iter = 13 - bars_SST_2D.test_vals = [13.000000, -0.456143, -1.541051] + bars_SST_2D.test_vals = [13.000000, -0.456168, -1.540956] bars_SST_2D.multizone = True test_list.append(bars_SST_2D) @@ -1156,7 +1164,7 @@ def main(): fsi_cht.cfg_dir = "fea_fsi/stat_fsi" fsi_cht.cfg_file = "config.cfg" fsi_cht.test_iter = 20 - fsi_cht.test_vals = [5.000000, -5.077002, -5.379450, -9.247804, -9.320014, -9.185034, 608.350000, -0.012973, 0.000000, 30.000000] + fsi_cht.test_vals = [5.000000, -5.077006, -5.379442, -9.247807, -9.319231, -9.184766, 608.350000, -0.012973, 0.000000, 30.000000] fsi_cht.multizone = True test_list.append(fsi_cht) @@ -1645,7 +1653,7 @@ def main(): pywrapper_turb_naca0012_sst.cfg_dir = "rans/naca0012" pywrapper_turb_naca0012_sst.cfg_file = "turb_NACA0012_sst.cfg" pywrapper_turb_naca0012_sst.test_iter = 10 - pywrapper_turb_naca0012_sst.test_vals = [-12.094431, -15.251082, -5.906366, 1.070413, 0.015775, -3.178469, 0.000000] + pywrapper_turb_naca0012_sst.test_vals = [-5.948762, -10.293853, -3.762064, 1.069291, 0.015853, -3.348638, 0.000000] pywrapper_turb_naca0012_sst.test_vals_aarch64 = [-12.076068, -15.246740, -5.861280, 1.070036, 0.015841, -3.297854, 0.000000] pywrapper_turb_naca0012_sst.command = TestCase.Command(exec = "SU2_CFD.py", param = "-f") pywrapper_turb_naca0012_sst.timeout = 3200 @@ -1659,7 +1667,7 @@ def main(): pywrapper_square_cylinder.cfg_dir = "unsteady/square_cylinder" pywrapper_square_cylinder.cfg_file = "turb_square.cfg" pywrapper_square_cylinder.test_iter = 3 - pywrapper_square_cylinder.test_vals = [-2.560678, -1.175979, 0.062202, 1.399351, 2.219234, 1.399297, 2.217479, 0.000000] + pywrapper_square_cylinder.test_vals = [-2.560698, -1.175978, 0.062220, 1.399351, 2.219234, 1.399297, 2.217479, 0.000000] pywrapper_square_cylinder.command = TestCase.Command(exec = "SU2_CFD.py", param = "-f") pywrapper_square_cylinder.timeout = 1600 pywrapper_square_cylinder.tol = 0.00001 diff --git a/TestCases/serial_regression_AD.py b/TestCases/serial_regression_AD.py index 8a0277294406..ef2b21b9915d 100644 --- a/TestCases/serial_regression_AD.py +++ b/TestCases/serial_regression_AD.py @@ -96,7 +96,7 @@ def main(): discadj_rans_naca0012_sst.cfg_dir = "disc_adj_rans/naca0012" discadj_rans_naca0012_sst.cfg_file = "turb_NACA0012_sst.cfg" discadj_rans_naca0012_sst.test_iter = 10 - discadj_rans_naca0012_sst.test_vals = [-2.096075, -0.181237, 0.350380, -0.023435] + discadj_rans_naca0012_sst.test_vals = [-2.096080, -0.181131, 0.352000, -0.023451] test_list.append(discadj_rans_naca0012_sst) ####################################### diff --git a/TestCases/tutorials.py b/TestCases/tutorials.py index d59aad17bf96..cab0d9c23f48 100644 --- a/TestCases/tutorials.py +++ b/TestCases/tutorials.py @@ -249,7 +249,7 @@ def main(): tutorial_trans_flatplate_T3A.cfg_dir = "../Tutorials/compressible_flow/Transitional_Flat_Plate/Langtry_and_Menter/T3A" tutorial_trans_flatplate_T3A.cfg_file = "transitional_LM_model_ConfigFile.cfg" tutorial_trans_flatplate_T3A.test_iter = 20 - tutorial_trans_flatplate_T3A.test_vals = [-5.790137, -2.054834, -3.894659, -0.255074, -1.747220, 5.119341, -3.493237, 0.393262] + tutorial_trans_flatplate_T3A.test_vals = [-5.790251, -2.054841, -3.895394, -0.255140, -1.750418, 5.119342, -3.493255, 0.392596] tutorial_trans_flatplate_T3A.test_vals_aarch64 = [-5.808996, -2.070606, -3.969765, -0.277943, -1.953289, 1.708472, -3.514943, 0.357411] tutorial_trans_flatplate_T3A.no_restart = True test_list.append(tutorial_trans_flatplate_T3A) @@ -259,7 +259,7 @@ def main(): tutorial_trans_flatplate_T3Am.cfg_dir = "../Tutorials/compressible_flow/Transitional_Flat_Plate/Langtry_and_Menter/T3A-" tutorial_trans_flatplate_T3Am.cfg_file = "transitional_LM_model_ConfigFile.cfg" tutorial_trans_flatplate_T3Am.test_iter = 20 - tutorial_trans_flatplate_T3Am.test_vals = [-5.587389, -1.700868, -3.093935, -0.102834, -3.750523, 3.287643, -2.394575, 1.119623] + tutorial_trans_flatplate_T3Am.test_vals = [-5.585365, -1.700868, -3.089065, -0.101095, -3.750523, 3.287642, -2.394575, 1.119623] tutorial_trans_flatplate_T3Am.test_vals_aarch64 = [-5.540938, -1.681627, -2.878831, -0.058224, -3.695533, 3.413628, -2.385345, 1.103633] tutorial_trans_flatplate_T3Am.no_restart = True test_list.append(tutorial_trans_flatplate_T3Am) diff --git a/TestCases/vandv.py b/TestCases/vandv.py index 516f51cc5ebd..0c9c9723d3b3 100644 --- a/TestCases/vandv.py +++ b/TestCases/vandv.py @@ -63,7 +63,7 @@ def main(): flatplate_sst1994m.cfg_dir = "vandv/rans/flatplate" flatplate_sst1994m.cfg_file = "turb_flatplate_sst.cfg" flatplate_sst1994m.test_iter = 5 - flatplate_sst1994m.test_vals = [-13.040943, -10.136971, -10.942764, -7.985181, -10.323857, -4.732487, 0.002801] + flatplate_sst1994m.test_vals = [-6.758709, -3.881664, -4.245031, -1.294380, -4.474846, -0.677798, 0.002801] flatplate_sst1994m.test_vals_aarch64 = [-13.021715, -9.534786, -10.401912, -7.501836, -9.750800, -4.850665, 0.002807] test_list.append(flatplate_sst1994m) @@ -72,7 +72,7 @@ def main(): bump_sst1994m.cfg_dir = "vandv/rans/bump_in_channel" bump_sst1994m.cfg_file = "turb_bump_sst.cfg" bump_sst1994m.test_iter = 5 - bump_sst1994m.test_vals = [-11.928292, -10.095796, -9.512953, -6.445652, -11.774088, -6.988752, 0.004931] + bump_sst1994m.test_vals = [-8.014868, -5.138354, -6.180378, -2.558011, -6.570741, -2.603570, 0.004931] bump_sst1994m.test_vals_aarch64 = [-13.042689, -10.812982, -10.604523, -7.655547, -10.816257, -5.308083, 0.004911] test_list.append(bump_sst1994m) @@ -91,7 +91,7 @@ def main(): swbli_sst.cfg_dir = "vandv/rans/swbli" swbli_sst.cfg_file = "config_sst.cfg" swbli_sst.test_iter = 5 - swbli_sst.test_vals = [-11.319578, -10.641523, -11.224600, -10.150215, -11.407538, -2.637660, 0.001816, -1.839484, -3.514593, 11.136000] + swbli_sst.test_vals = [-10.620096, -9.467013, -10.489091, -9.162510, -10.530142, -2.612181, 0.001816, -1.568660, -3.540055, 8.615300] test_list.append(swbli_sst) # DSMA661 - SA @@ -108,7 +108,7 @@ def main(): dsma661_sst.cfg_dir = "vandv/rans/dsma661" dsma661_sst.cfg_file = "dsma661_sst_config.cfg" dsma661_sst.test_iter = 5 - dsma661_sst.test_vals = [-11.025153, -8.156995, -9.057021, -5.947228, -10.650874, -7.884423, 0.155882, 0.023344] + dsma661_sst.test_vals = [-2.616553, 0.009594, -0.244123, 2.840416, -2.180298, 1.413440, 0.156003, 0.023336] dsma661_sst.test_vals_aarch64 = [-10.977195, -8.403731, -8.747068, -5.808899, -10.522786, -7.369851, 0.155875, 0.023353] test_list.append(dsma661_sst) diff --git a/UnitTests/SU2_CFD/numerics/ausm_tke_tests.cpp b/UnitTests/SU2_CFD/numerics/ausm_tke_tests.cpp new file mode 100644 index 000000000000..1a8decc8a0a0 --- /dev/null +++ b/UnitTests/SU2_CFD/numerics/ausm_tke_tests.cpp @@ -0,0 +1,173 @@ +/*! + * \file ausm_tke_tests.cpp + * \brief SLAU, SLAU2 and AUSM with the turbulent kinetic energy of SST in the total energy. + * \version 8.5.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#include "catch.hpp" +#include +#include +#include +#include +#include +#include "../../../SU2_CFD/include/numerics/flow/convection/ausm_slau.hpp" +#include "../../../SU2_CFD/include/numerics/flow/convection/fvs.hpp" + +namespace { + +constexpr su2double gamma = 1.4; + +std::unique_ptr MakeConfig(bool implicit = false) { + std::stringstream options; + options << "SOLVER= EULER\nCONV_NUM_METHOD_FLOW= SLAU\nMACH_NUMBER= 0.5\n"; + options << "TIME_DISCRE_FLOW= " + << (implicit ? "EULER_IMPLICIT\nUSE_ACCURATE_FLUX_JACOBIANS= YES\n" : "EULER_EXPLICIT\n"); + return std::make_unique(options, SU2_COMPONENT::SU2_CFD, false); +} + +/*--- Primitive variables (T, u, v, [w], p, rho, h, c) of an ideal gas whose total energy contains k. ---*/ +void Primitives(unsigned short nDim, su2double rho, const su2double* vel, su2double p, su2double k, su2double* V) { + su2double q2 = 0.0; + for (unsigned short iDim = 0; iDim < nDim; iDim++) { + V[iDim + 1] = vel[iDim]; + q2 += vel[iDim] * vel[iDim]; + } + V[0] = p / (rho * 287.058); + V[nDim + 1] = p; + V[nDim + 2] = rho; + V[nDim + 3] = gamma / (gamma - 1) * p / rho + 0.5 * q2 + k; + V[nDim + 4] = sqrt(gamma * p / rho); +} + +/*--- k changes neither the pressure nor the speed of sound: with the same density, velocity and pressure, the mass + and momentum fluxes must not depend on k, and the energy flux changes only by the advected k. ---*/ +void CheckTkeIndependence(CNumerics& numerics, const CConfig* config, unsigned short nDim, bool uniformTke = false) { + const su2double vel_i[3] = {0.8, 0.1, -0.05}, vel_j[3] = {0.6, -0.1, 0.05}, n[3] = {2.0 / 7, 3.0 / 7, 6.0 / 7}; + const su2double n2[2] = {0.6, 0.8}; + su2double Vi[10] = {0.0}, Vj[10] = {0.0}; + numerics.SetNormal(nDim == 2 ? n2 : n); + + su2double flux[2][5] = {{0.0}}; + /*--- MSW mixes the two states, so k must be uniform for the fluxes to be independent of it. ---*/ + const su2double k_i[2] = {0.0, 0.3}, k_j[2] = {0.0, uniformTke ? 0.3 : 0.2}; + for (int c = 0; c < 2; c++) { + Primitives(nDim, 1.2, vel_i, 1.1, k_i[c], Vi); + Primitives(nDim, 0.9, vel_j, 0.8, k_j[c], Vj); + numerics.SetPrimitive(Vi, Vj); + numerics.SetTurbKineticEnergy(k_i[c], k_j[c]); + const auto res = numerics.ComputeResidual(config); + for (unsigned short iVar = 0; iVar < nDim + 2; iVar++) flux[c][iVar] = res.residual[iVar]; + } + CAPTURE(nDim); + for (unsigned short iVar = 0; iVar < nDim + 1; iVar++) CHECK(flux[1][iVar] == Approx(flux[0][iVar]).epsilon(1e-12)); + /*--- Energy: the extra flux is the mass flux times the upwind k (the flow is from i to j here). ---*/ + CHECK(flux[1][nDim + 1] - flux[0][nDim + 1] == Approx(flux[0][0] * k_i[1]).epsilon(1e-10)); +} + +/*--- Largest difference between the (accurate) Jacobians and central differences of the flux with respect to the + conservative variables, k held fixed, relative to the largest entry. ---*/ +su2double JacobianError(CNumerics& numerics, const CConfig* config, su2double k) { + constexpr unsigned short nDim = 2, nVar = 4; + using State = std::array; + const su2double normal[2] = {0.6, 0.8}; + numerics.SetNormal(normal); + auto evaluate = [&](const State& Ui, const State& Uj, su2double* flux, su2double(*jac_i)[nVar], + su2double(*jac_j)[nVar]) { + su2double Vi[10] = {0.0}, Vj[10] = {0.0}; + for (int side = 0; side < 2; side++) { + const State& U = side ? Uj : Ui; + const su2double vel[2] = {U[1] / U[0], U[2] / U[0]}; + const su2double p = (gamma - 1) * (U[3] - 0.5 * (U[1] * U[1] + U[2] * U[2]) / U[0] - U[0] * k); + Primitives(nDim, U[0], vel, p, k, side ? Vj : Vi); + } + numerics.SetPrimitive(Vi, Vj); + numerics.SetTurbKineticEnergy(k, k); + const auto res = numerics.ComputeResidual(config); + for (unsigned short iVar = 0; iVar < nVar; iVar++) { + flux[iVar] = res.residual[iVar]; + if (jac_i) + for (unsigned short jVar = 0; jVar < nVar; jVar++) { + jac_i[iVar][jVar] = res.jacobian_i[iVar][jVar]; + jac_j[iVar][jVar] = res.jacobian_j[iVar][jVar]; + } + } + }; + const State Ui = {1.2, 0.72, 0.24, 1.1 / (gamma - 1) + 0.5 * (0.72 * 0.72 + 0.24 * 0.24) / 1.2 + 1.2 * k}; + const State Uj = {0.9, 0.36, -0.18, 0.8 / (gamma - 1) + 0.5 * (0.36 * 0.36 + 0.18 * 0.18) / 0.9 + 0.9 * k}; + su2double flux[nVar], jac_i[nVar][nVar], jac_j[nVar][nVar]; + evaluate(Ui, Uj, flux, jac_i, jac_j); + su2double maxErr = 0.0, maxJac = 0.0; + for (int side = 0; side < 2; side++) { + for (unsigned short jVar = 0; jVar < nVar; jVar++) { + State Up = side ? Uj : Ui, Um = Up; + const su2double h = 1e-6; + Up[jVar] += h; + Um[jVar] -= h; + su2double fp[nVar], fm[nVar]; + if (side) { + evaluate(Ui, Up, fp, nullptr, nullptr); + evaluate(Ui, Um, fm, nullptr, nullptr); + } else { + evaluate(Up, Uj, fp, nullptr, nullptr); + evaluate(Um, Uj, fm, nullptr, nullptr); + } + for (unsigned short iVar = 0; iVar < nVar; iVar++) { + const su2double jac = side ? jac_j[iVar][jVar] : jac_i[iVar][jVar]; + maxErr = std::max(maxErr, std::abs(jac - (fp[iVar] - fm[iVar]) / (2 * h))); + maxJac = std::max(maxJac, std::abs(jac)); + } + } + } + return maxErr / maxJac; +} + +} // namespace + +TEST_CASE("Accurate SLAU and AUSM+up Jacobians hold k fixed", "[AUSM][SST]") { + auto config = MakeConfig(true); + for (const su2double k : {0.0, 0.2}) { + CAPTURE(k); + CUpwSLAU_Flow slau(2, 4, config.get(), false); + CHECK(JacobianError(slau, config.get(), k) < 1e-4); + CUpwAUSMPLUSUP_Flow ausmup(2, 4, config.get()); + CHECK(JacobianError(ausmup, config.get(), k) < 1e-4); + } +} + +TEST_CASE("SLAU, SLAU2, AUSM, AUSM+up, AUSM+up2 and MSW do not use k in the speed of sound", "[AUSM][SST]") { + auto config = MakeConfig(); + for (const unsigned short nDim : {2, 3}) { + CUpwSLAU_Flow slau(nDim, nDim + 2, config.get(), false); + CheckTkeIndependence(slau, config.get(), nDim); + CUpwSLAU2_Flow slau2(nDim, nDim + 2, config.get(), false); + CheckTkeIndependence(slau2, config.get(), nDim); + CUpwAUSM_Flow ausm(nDim, nDim + 2, config.get()); + CheckTkeIndependence(ausm, config.get(), nDim); + CUpwAUSMPLUSUP_Flow ausmup(nDim, nDim + 2, config.get()); + CheckTkeIndependence(ausmup, config.get(), nDim); + CUpwAUSMPLUSUP2_Flow ausmup2(nDim, nDim + 2, config.get()); + CheckTkeIndependence(ausmup2, config.get(), nDim); + CUpwMSW_Flow msw(nDim, nDim + 2, config.get()); + CheckTkeIndependence(msw, config.get(), nDim, true); + } +} diff --git a/UnitTests/SU2_CFD/numerics/roe_tke_tests.cpp b/UnitTests/SU2_CFD/numerics/roe_tke_tests.cpp new file mode 100644 index 000000000000..eb2a97061d24 --- /dev/null +++ b/UnitTests/SU2_CFD/numerics/roe_tke_tests.cpp @@ -0,0 +1,169 @@ +/*! + * \file roe_tke_tests.cpp + * \brief Roe eigenvector matrices with the turbulent kinetic energy of SST in the total energy. + * \version 8.5.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#include "catch.hpp" +#include +#include +#include +#include "../../../SU2_CFD/include/numerics/flow/convection/roe.hpp" + +namespace { + +constexpr su2double gamma = 1.4; + +std::unique_ptr MakeConfig(bool implicit = true) { + std::stringstream options; + options << "SOLVER= EULER\nCONV_NUM_METHOD_FLOW= ROE\nENTROPY_FIX_COEFF= 0.0\nROE_KAPPA= 0.5\n"; + options << "TIME_DISCRE_FLOW= " << (implicit ? "EULER_IMPLICIT" : "EULER_EXPLICIT") << "\n"; + return std::make_unique(options, SU2_COMPONENT::SU2_CFD, false); +} + +/*--- Primitive variables (T, u, v, [w], p, rho, h, c) of an ideal gas whose total energy contains k. ---*/ +void Primitives(unsigned short nDim, su2double rho, const su2double* vel, su2double p, su2double k, su2double* V) { + su2double q2 = 0.0; + for (unsigned short iDim = 0; iDim < nDim; iDim++) { + V[iDim + 1] = vel[iDim]; + q2 += vel[iDim] * vel[iDim]; + } + V[0] = p / (rho * 287.058); + V[nDim + 1] = p; + V[nDim + 2] = rho; + V[nDim + 3] = gamma / (gamma - 1) * p / rho + 0.5 * q2 + k; + V[nDim + 4] = sqrt(gamma * p / rho); +} + +/*--- Secondary variables (dp/drho_e, dp/de_rho) of the ideal gas, for the general-gas schemes. ---*/ +void Secondary(su2double rho, su2double p, su2double* S) { + S[0] = p / rho; + S[1] = (gamma - 1) * rho; +} + +/*--- Common constructor for the schemes of the family. ---*/ +template +std::unique_ptr MakeScheme(unsigned short nDim, unsigned short nVar, const CConfig* config) { + return std::make_unique(nDim, nVar, config); +} +template <> +std::unique_ptr MakeScheme(unsigned short nDim, unsigned short nVar, const CConfig* config) { + return std::make_unique(nDim, nVar, config, false); +} + +/*--- Roe flux between two states with the same pressure and velocity and different density and k. With the + entropy fix off this is a contact discontinuity, whose exact (upwind) flux is the flux of the upstream state. ---*/ +template +void CheckContactWithTkeJump(const char* name, unsigned short nDim, const su2double* vel, bool implicit) { + auto config = MakeConfig(implicit); + const unsigned short nVar = nDim + 2; + auto numerics = MakeScheme(nDim, nVar, config.get()); + + const su2double p = 1.0, rho_i = 1.4, rho_j = 0.7, k_i = 0.1, k_j = 0.3; + const su2double n2[2] = {0.6, 0.8}, n3[3] = {2.0 / 7, 3.0 / 7, 6.0 / 7}; + const su2double* n = (nDim == 2) ? n2 : n3; + su2double Vi[10] = {0.0}, Vj[10] = {0.0}, Si[2], Sj[2]; + Primitives(nDim, rho_i, vel, p, k_i, Vi); + Primitives(nDim, rho_j, vel, p, k_j, Vj); + Secondary(rho_i, p, Si); + Secondary(rho_j, p, Sj); + numerics->SetPrimitive(Vi, Vj); + numerics->SetSecondary(Si, Sj); + numerics->SetNormal(n); + numerics->SetTurbKineticEnergy(k_i, k_j); + const auto res = numerics->ComputeResidual(config.get()); + + su2double un = 0.0; + for (unsigned short iDim = 0; iDim < nDim; iDim++) un += vel[iDim] * n[iDim]; + const su2double rho_up = (un >= 0.0) ? rho_i : rho_j; + const su2double* V_up = (un >= 0.0) ? Vi : Vj; + + CAPTURE(name, nDim, un, implicit); + CHECK(res.residual[0] == Approx(rho_up * un).margin(1e-12)); + for (unsigned short iDim = 0; iDim < nDim; iDim++) + CHECK(res.residual[iDim + 1] == Approx(rho_up * vel[iDim] * un + p * n[iDim]).margin(1e-12)); + CHECK(res.residual[nVar - 1] == Approx(rho_up * V_up[nDim + 3] * un).margin(1e-12)); +} + +template +void CheckContactsWithTkeJump(const char* name, bool implicit = true, unsigned short maxDim = 3) { + const su2double zero[3] = {0.0, 0.0, 0.0}, moving[3] = {0.3, -0.1, 0.2}, reverse[3] = {-0.3, 0.1, -0.2}; + for (unsigned short nDim = 2; nDim <= maxDim; nDim++) { + for (const auto* vel : {zero, moving, reverse}) CheckContactWithTkeJump(name, nDim, vel, implicit); + } +} + +} // namespace + +TEST_CASE("Roe schemes resolve a contact with a jump of k exactly", "[Roe][SST]") { + CheckContactsWithTkeJump("Roe"); + CheckContactsWithTkeJump("general Roe", true); + /*--- 2D only: in 3D, the wave strengths of these schemes do not match the columns of the P matrix (not related + to k, fixed separately). ---*/ + CheckContactsWithTkeJump("L2Roe", true, 2); + CheckContactsWithTkeJump("LMRoe", true, 2); + CheckContactsWithTkeJump("general Roe", false, 2); +} + +TEST_CASE("Roe eigenvector matrices are inverse with k in the total energy", "[Roe][SST]") { + auto config = MakeConfig(); + for (const unsigned short nDim : {2, 3}) { + const unsigned short nVar = nDim + 2; + CUpwRoe_Flow numerics(nDim, nVar, config.get(), false); + const su2double velocity[3] = {0.4, -0.3, 0.2}, normal[3] = {2.0 / 7, 3.0 / 7, 6.0 / 7}; + const su2double n2[3] = {0.6, 0.8, 0.0}; + su2double P[5][5], Pinv[5][5]; + numerics.GetPMatrix(1.1, velocity, 1.2, nDim == 3 ? normal : n2, P, 0.3); + numerics.GetPMatrix_inv(1.1, velocity, 1.2, nDim == 3 ? normal : n2, Pinv, 0.3); + for (unsigned short i = 0; i < nVar; i++) { + for (unsigned short j = 0; j < nVar; j++) { + su2double prod = 0.0; + for (unsigned short l = 0; l < nVar; l++) prod += P[i][l] * Pinv[l][j]; + CAPTURE(nDim, i, j); + CHECK(prod == Approx(i == j ? 1.0 : 0.0).margin(1e-12)); + } + } + } +} + +TEST_CASE("Roe does not dissipate a stationary contact with k in the total energy", "[Roe][SST]") { + /*--- Zero velocity, same pressure and k on both sides, different density: a contact discontinuity at rest. + The Roe flux must be the pressure only (no mass or energy flux) when the entropy fix is off. ---*/ + auto config = MakeConfig(); + constexpr unsigned short nDim = 2, nVar = nDim + 2; + CUpwRoe_Flow numerics(nDim, nVar, config.get(), false); + + const su2double zero[2] = {0.0, 0.0}, p = 1.0, k = 0.2; + su2double Vi[10] = {0.0}, Vj[10] = {0.0}, normal[2] = {0.6, 0.8}; + Primitives(nDim, 1.4, zero, p, k, Vi); + Primitives(nDim, 0.7, zero, p, k, Vj); + numerics.SetPrimitive(Vi, Vj); + numerics.SetNormal(normal); + numerics.SetTurbKineticEnergy(k, k); + const auto res = numerics.ComputeResidual(config.get()); + + CHECK(res.residual[0] == Approx(0.0).margin(1e-12)); + CHECK(res.residual[1] == Approx(p * normal[0]).epsilon(1e-12)); + CHECK(res.residual[2] == Approx(p * normal[1]).epsilon(1e-12)); + CHECK(res.residual[3] == Approx(0.0).margin(1e-12)); +} diff --git a/UnitTests/meson.build b/UnitTests/meson.build index fc99bc101b3a..08e556861451 100644 --- a/UnitTests/meson.build +++ b/UnitTests/meson.build @@ -15,6 +15,8 @@ su2_cfd_tests = files(['Common/CConfig_tests.cpp', 'Common/containers/CLookupTable_tests.cpp', 'Common/toolboxes/multilayer_perceptron/CLookUp_ANN_tests.cpp', 'SU2_CFD/numerics/CNumerics_tests.cpp', + 'SU2_CFD/numerics/roe_tke_tests.cpp', + 'SU2_CFD/numerics/ausm_tke_tests.cpp', 'SU2_CFD/fluid/CFluidModel_tests.cpp', 'SU2_CFD/gradients.cpp', 'SU2_CFD/solvers/CNEMOEulerSolver_tests.cpp', diff --git a/config_template.cfg b/config_template.cfg index 88533d33c689..938e9817af7e 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,