diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index 0c32f044288..ee47d4cacaf 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -877,6 +877,7 @@ class CConfig { array MuPolyCoefficientsND{{0.0}}; /*!< \brief Definition of the non-dimensional temperature polynomial coefficients for viscosity. */ array KtPolyCoefficientsND{{0.0}}; /*!< \brief Definition of the non-dimensional temperature polynomial coefficients for thermal conductivity. */ su2double TurbIntensityAndViscRatioFreeStream[2]; /*!< \brief Freestream turbulent intensity and viscosity ratio for turbulence and transition models. */ + su2double Energy_FreeStream, /*!< \brief Free-stream total energy of the fluid. */ ModVel_FreeStream, /*!< \brief Magnitude of the free-stream velocity of the fluid. */ ModVel_FreeStreamND, /*!< \brief Non-dimensional magnitude of the free-stream velocity of the fluid. */ @@ -1176,6 +1177,8 @@ class CConfig { nHistoryOutput, nVolumeOutput; /*!< \brief Number of variables printed to the history file. */ bool Multizone_Residual; /*!< \brief Determines if memory should be allocated for the multizone residual. */ SST_ParsedOptions sstParsedOptions; /*!< \brief Additional parameters for the SST turbulence model. */ + su2double prodLimConst; + su2double LDomain; SA_ParsedOptions saParsedOptions; /*!< \brief Additional parameters for the SA turbulence model. */ LM_ParsedOptions lmParsedOptions; /*!< \brief Additional parameters for the LM transition model. */ su2double uq_delta_b; /*!< \brief Parameter used to perturb eigenvalues of Reynolds Stress Matrix */ @@ -9908,6 +9911,9 @@ class CConfig { */ SST_ParsedOptions GetSSTParsedOptions() const { return sstParsedOptions; } + su2double GetProdLimConst() const { return prodLimConst; } + su2double GetLDomain() const { return LDomain; } + /*! * \brief Get parsed SA option data structure. * \return SA option data structure. diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index ddf570ab697..7e056ce703c 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -994,14 +994,16 @@ enum class SST_OPTIONS { COMP_Wilcox, /*!< \brief Menter k-w SST model with Compressibility correction of Wilcox. */ COMP_Sarkar, /*!< \brief Menter k-w SST model with Compressibility correction of Sarkar. */ DLL, /*!< \brief Menter k-w SST model with dimensionless lower limit clipping of turbulence variables. */ + FULLPROD, /*!< \brief Menter k-w SST model with full production term. */ + PRODLIM, /*!< \brief Menter k-w SST model with user-defined production limiter constant. */ + NEWBC, /*!< \brief Menter k-w SST model with new boundary conditions. */ }; static const MapType 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) @@ -1009,6 +1011,9 @@ 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("FULLPROD", SST_OPTIONS::FULLPROD) + MakePair("PRODLIM", SST_OPTIONS::PRODLIM) + MakePair("NEWBC", SST_OPTIONS::NEWBC) }; /*! @@ -1023,6 +1028,9 @@ struct SST_ParsedOptions { bool compWilcox = false; /*!< \brief Bool for compressibility correction of Wilcox. */ bool compSarkar = false; /*!< \brief Bool for compressibility correction of Sarkar. */ bool dll = false; /*!< \brief Bool dimensionless lower limit. */ + bool fullProd = false; /*!< \brief Bool for full production term. */ + bool prodLim = false; /*!< \brief Bool for user-defined production limiter constant. */ + bool newBC = false; /*!< \brief Bool for new boundary conditions. */ }; /*! @@ -1040,6 +1048,10 @@ inline SST_ParsedOptions ParseSSTOptions(const SST_OPTIONS *SST_Options, unsigne return std::find(SST_Options, sst_options_end, option) != sst_options_end; }; + const bool found_fullProd = IsPresent(SST_OPTIONS::FULLPROD); + const bool found_prodLim = IsPresent(SST_OPTIONS::PRODLIM); + const bool found_newBC = IsPresent(SST_OPTIONS::NEWBC); + const bool found_1994 = IsPresent(SST_OPTIONS::V1994); const bool found_2003 = IsPresent(SST_OPTIONS::V2003); const bool found_1994m = IsPresent(SST_OPTIONS::V1994m); @@ -1097,6 +1109,9 @@ inline SST_ParsedOptions ParseSSTOptions(const SST_OPTIONS *SST_Options, unsigne SSTParsedOptions.compSarkar = sst_compSarkar; SSTParsedOptions.dll = sst_dll; + SSTParsedOptions.fullProd = found_fullProd; + SSTParsedOptions.prodLim = found_prodLim; + SSTParsedOptions.newBC = found_newBC; return SSTParsedOptions; } diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 8a60579189c..ca37df523e9 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -1119,6 +1119,9 @@ void CConfig::SetConfig_Options() { /*!\brief SST_OPTIONS \n DESCRIPTION: Specify SA turbulence model options/corrections. \n Options: see \link SA_Options_Map \endlink \n DEFAULT: NONE \ingroup Config*/ addEnumListOption("SA_OPTIONS", nSA_Options, SA_Options, SA_Options_Map); + addDoubleOption("PROD_LIM_CONST", prodLimConst, 20.0); + addDoubleOption("L_DOMAIN", LDomain, 1.0); + /*!\brief KIND_TRANS_MODEL \n DESCRIPTION: Specify transition model OPTIONS: see \link Trans_Model_Map \endlink \n DEFAULT: NONE \ingroup Config*/ addEnumOption("KIND_TRANS_MODEL", Kind_Trans_Model, Trans_Model_Map, TURB_TRANS_MODEL::NONE); /*!\brief SST_OPTIONS \n DESCRIPTION: Specify LM transition model options/correlations. \n Options: see \link LM_Options_Map \endlink \n DEFAULT: NONE \ingroup Config*/ @@ -6275,6 +6278,9 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { else cout << "\nusing default hard coded lower limit clipping"; cout << "." << endl; + + if (sstParsedOptions.prodLim) cout << "Changing the value of the TKE production limiter constant to " << prodLimConst << endl; + break; } switch (Kind_Trans_Model) { diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index 535335947ca..bc456ea31fc 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -190,6 +190,8 @@ class CNumerics { bool bounded_scalar = false; /*!< \brief Flag for bounded scalar problem */ + su2double ProdDistr[5]; + public: /*! * \brief Return type used in some "ComputeResidual" overloads to give a @@ -1605,6 +1607,8 @@ class CNumerics { * \return is_bounded_scalar : scalar solver uses bounded scalar convective transport */ inline bool GetBoundedScalar() const { return bounded_scalar;} + + inline su2double GetProdDest(int index) const {return ProdDistr[index];} }; /*! diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp index c9464b1dfc3..42cabf926dd 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp @@ -656,6 +656,8 @@ class CSourcePieceWise_TurbSST final : public CNumerics { /*--- Ambient values for SST-SUST. ---*/ const su2double kAmb, omegaAmb; + su2double Pk, Dk, Pw, Dw; + su2double F1_i, F2_i, CDkw_i; su2double Residual[2]; su2double* Jacobian_i[2]; @@ -869,15 +871,20 @@ class CSourcePieceWise_TurbSST final : public CNumerics { const su2double prod_limit = prod_lim_const * beta_star * Density_i * ScalarVar_i[1] * ScalarVar_i[0]; su2double P = Eddy_Viscosity_i * pow(P_Base, 2); + if (sstParsedOptions.fullProd) P -= Eddy_Viscosity_i * diverg*diverg * 2.0/3.0; + if (!sstParsedOptions.modified) P -= Density_i * ScalarVar_i[0] * diverg * 2.0/3.0; + + su2double PLim = 0.0; + if ( P > prod_limit ) PLim = 1.0; su2double pk = max(0.0, min(P, prod_limit)); - const auto& eddy_visc_var = sstParsedOptions.version == SST_OPTIONS::V1994 ? VorticityMag : StrainMag_i; + const su2double eddy_visc_var = sstParsedOptions.version == SST_OPTIONS::V1994 ? VorticityMag : StrainMag_i; const su2double zeta = max(ScalarVar_i[1], eddy_visc_var * F2_i / a1); /*--- Production limiter only for V2003, recompute for V1994. ---*/ su2double pw; if (sstParsedOptions.version == SST_OPTIONS::V1994) { - pw = alfa_blended * Density_i * pow(P_Base, 2); + pw = alfa_blended * Density_i * P / Eddy_Viscosity_i; } else { pw = (alfa_blended * Density_i / Eddy_Viscosity_i) * pk; } @@ -912,17 +919,26 @@ class CSourcePieceWise_TurbSST final : public CNumerics { } /*--- Add the production terms to the residuals. ---*/ + Pk = pk; + Pw = pw; Residual[0] += pk * Volume; Residual[1] += pw * Volume; /*--- Add the dissipation terms to the residuals.---*/ + Dk = dk; + Dw = dw; Residual[0] -= dk * Volume; Residual[1] -= dw * Volume; - /*--- Cross diffusion ---*/ + ProdDistr[0] = Pk; + ProdDistr[1] = Dk; + ProdDistr[2] = Pw; + ProdDistr[3] = Dw; + ProdDistr[4] = PLim; + /*--- Cross diffusion ---*/ Residual[1] += (1.0 - F1_i) * CDkw_i * Volume; /*--- Contribution due to 2D axisymmetric formulation ---*/ @@ -932,6 +948,7 @@ class CSourcePieceWise_TurbSST final : public CNumerics { /*--- Implicit part ---*/ Jacobian_i[0][0] = -beta_star * ScalarVar_i[1] * Volume * (1.0 + zetaFMt); + if (!sstParsedOptions.modified) Jacobian_i[0][0] -= diverg * Volume*2.0/3.0; Jacobian_i[0][1] = -beta_star * ScalarVar_i[0] * Volume * (1.0 + zetaFMt); Jacobian_i[1][0] = 0.0; Jacobian_i[1][1] = -2.0 * beta_blended * ScalarVar_i[1] * Volume * (1.0 - 0.09/beta_blended * zetaFMt); diff --git a/SU2_CFD/include/numerics_simd/flow/convection/centered.hpp b/SU2_CFD/include/numerics_simd/flow/convection/centered.hpp index 4e73ab5b752..4cd37bfc3b3 100644 --- a/SU2_CFD/include/numerics_simd/flow/convection/centered.hpp +++ b/SU2_CFD/include/numerics_simd/flow/convection/centered.hpp @@ -53,6 +53,8 @@ class CCenteredBase : public Base { const bool dynamicGrid; const su2double stretchParam = 0.3; + using Base::turbVars; + /*! * \brief Constructor, store some constants and forward args to base. */ @@ -96,6 +98,8 @@ class CCenteredBase : public Base { const bool implicit = (config.GetKind_TimeIntScheme() == EULER_IMPLICIT); const auto& solution = static_cast(solution_); + const bool tkeNeeded = config.GetKind_Turb_Model() == TURB_MODEL::SST; + const auto iPoint = geometry.edges->GetNode(iEdge,0); const auto jPoint = geometry.edges->GetNode(iEdge,1); @@ -119,6 +123,13 @@ class CCenteredBase : public Base { avgV.all(iVar) = 0.5 * (V.i.all(iVar) + V.j.all(iVar)); } + if (tkeNeeded) { + V.i.allTurb = gatherVariables<1>(iPoint, turbVars->GetSolution()); + V.j.allTurb = gatherVariables<1>(jPoint, turbVars->GetSolution()); + + avgV.allTurb(0) = 0.5*(V.i.allTurb(0)+V.j.allTurb(0)); + } + /*--- Compute conservative variables. ---*/ CPair > U; diff --git a/SU2_CFD/include/numerics_simd/flow/convection/common.hpp b/SU2_CFD/include/numerics_simd/flow/convection/common.hpp index b5d3c39e3f7..da49bb36bbc 100644 --- a/SU2_CFD/include/numerics_simd/flow/convection/common.hpp +++ b/SU2_CFD/include/numerics_simd/flow/convection/common.hpp @@ -113,7 +113,8 @@ FORCEINLINE CPair reconstructPrimitives(Int iEdge, Int iPoint, Int bool muscl, LIMITER limiterType, const CPair& V1st, const VectorDbl& vector_ij, - const VariableType& solution) { + const VariableType& solution, + const bool tkeNeeded) { static_assert(ReconVarType::nVar <= PrimVarType::nVar,""); const auto& gradients = solution.GetGradient_Reconstruction(); @@ -125,6 +126,11 @@ FORCEINLINE CPair reconstructPrimitives(Int iEdge, Int iPoint, Int V.i.all(iVar) = V1st.i.all(iVar); V.j.all(iVar) = V1st.j.all(iVar); } + // Only first order for turbulence + if (tkeNeeded) { + V.i.allTurb(0) = V1st.i.allTurb(0); + V.j.allTurb(0) = V1st.j.allTurb(0); + } if (muscl) { /*--- Recompute density and enthalpy instead of reconstructing. ---*/ diff --git a/SU2_CFD/include/numerics_simd/flow/convection/roe.hpp b/SU2_CFD/include/numerics_simd/flow/convection/roe.hpp index 204738b1ef9..98be310304c 100644 --- a/SU2_CFD/include/numerics_simd/flow/convection/roe.hpp +++ b/SU2_CFD/include/numerics_simd/flow/convection/roe.hpp @@ -63,6 +63,8 @@ class CRoeBase : public Base { const bool muscl; const LIMITER typeLimiter; + using Base::turbVars; + /*! * \brief Constructor, store some constants and forward args to base. */ @@ -98,6 +100,8 @@ class CRoeBase : public Base { const bool implicit = (config.GetKind_TimeIntScheme() == EULER_IMPLICIT); const auto& solution = static_cast(solution_); + const bool tkeNeeded = config.GetKind_Turb_Model() == TURB_MODEL::SST; + const auto iPoint = geometry.edges->GetNode(iEdge,0); const auto jPoint = geometry.edges->GetNode(iEdge,1); @@ -118,8 +122,13 @@ class CRoeBase : public Base { V1st.i.all = gatherVariables(iPoint, solution.GetPrimitive()); V1st.j.all = gatherVariables(jPoint, solution.GetPrimitive()); + if (tkeNeeded) { + V1st.i.allTurb = gatherVariables(iPoint, turbVars->GetSolution()); + V1st.j.allTurb = gatherVariables(jPoint, turbVars->GetSolution()); + } + auto V = reconstructPrimitives >( - iEdge, iPoint, jPoint, gamma, gasConst, muscl, typeLimiter, V1st, vector_ij, solution); + iEdge, iPoint, jPoint, gamma, gasConst, muscl, typeLimiter, V1st, vector_ij, solution, tkeNeeded); /*--- Compute conservative variables. ---*/ @@ -231,6 +240,8 @@ class CRoeScheme : public CRoeBase,Decorator> { using Base::kappa; const ENUM_ROELOWDISS typeDissip; + using Base::turbVars; + public: /*! * \brief Constructor, store some constants and forward to base. diff --git a/SU2_CFD/include/numerics_simd/flow/diffusion/common.hpp b/SU2_CFD/include/numerics_simd/flow/diffusion/common.hpp index 279a623c0c7..555ad0c1390 100644 --- a/SU2_CFD/include/numerics_simd/flow/diffusion/common.hpp +++ b/SU2_CFD/include/numerics_simd/flow/diffusion/common.hpp @@ -71,13 +71,13 @@ FORCEINLINE void correctGradient(const PrimitiveType& V, */ template FORCEINLINE MatrixDbl stressTensor(Double viscosity, - const MatrixDbl& grad) { + const MatrixDbl& grad, Double density=0.0, Double tke=0.0) { /*--- Hydrostatic term. ---*/ Double velDiv = 0.0; for (size_t iDim = 0; iDim < nDim; ++iDim) { velDiv += grad(iDim+1,iDim); } - Double pTerm = 2.0/3.0 * viscosity * velDiv; + Double pTerm = 2.0/3.0 * (viscosity * velDiv + density * tke); MatrixDbl tau; for (size_t iDim = 0; iDim < nDim; ++iDim) { 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 73bfaac1537..90123fa4c62 100644 --- a/SU2_CFD/include/numerics_simd/flow/diffusion/viscous_fluxes.hpp +++ b/SU2_CFD/include/numerics_simd/flow/diffusion/viscous_fluxes.hpp @@ -50,6 +50,7 @@ class CNoViscousFlux : public CNumericsSIMD { protected: static constexpr size_t nDim = NDIM; static constexpr size_t nPrimVar = 0; + const CVariable* turbVars = nullptr; template CNoViscousFlux(Ts&...) {} @@ -140,6 +141,8 @@ class CCompressibleViscousFluxBase : public CNumericsSIMD { const auto& solution = static_cast(solution_); const auto& gradient = solution.GetGradient_Primitive(); + const bool tkeNeeded = (config.GetKind_Turb_Model() == TURB_MODEL::SST) && !(config.GetSSTParsedOptions().modified); + /*--- Compute distance and handle zero without "ifs" by making it large. ---*/ auto dist2_ij = squaredNorm(vector_ij); @@ -152,12 +155,15 @@ class CCompressibleViscousFluxBase : public CNumericsSIMD { if(correct) correctGradient(V, vector_ij, dist2_ij, avgGrad); /*--- Stress and heat flux tensors. ---*/ - - auto tau = stressTensor(avgV.laminarVisc() + (uq? Double(0.0) : avgV.eddyVisc()), avgGrad); + Double turb_ke = 0.0; + Double density = avgV.density(); + if (tkeNeeded) turb_ke = 0.5*(gatherVariables(iPoint, turbVars->GetSolution()) + + gatherVariables(jPoint, turbVars->GetSolution())); + auto tau = stressTensor(avgV.laminarVisc() + (uq? Double(0.0) : avgV.eddyVisc()), avgGrad, density, turb_ke); if(useSA_QCR) addQCR(avgGrad, tau); if(uq) { - Double turb_ke = 0.5*(gatherVariables(iPoint, turbVars->GetSolution()) + - gatherVariables(jPoint, turbVars->GetSolution())); + turb_ke = 0.5*(gatherVariables(iPoint, turbVars->GetSolution()) + + gatherVariables(jPoint, turbVars->GetSolution())); addPerturbedRSM(avgV, avgGrad, turb_ke, tau, uq_eigval_comp, uq_permute, uq_delta_b, uq_urlx); } @@ -263,6 +269,7 @@ class CCompressibleViscousFlux : public CCompressibleViscousFluxBase >; using Base::prandtlTurb; + using Base::turbVars; +; /*! * \brief Constructor, initialize constants and booleans. diff --git a/SU2_CFD/include/numerics_simd/flow/variables.hpp b/SU2_CFD/include/numerics_simd/flow/variables.hpp index 086adbf27dc..178edcd5017 100644 --- a/SU2_CFD/include/numerics_simd/flow/variables.hpp +++ b/SU2_CFD/include/numerics_simd/flow/variables.hpp @@ -50,12 +50,16 @@ struct CCompressiblePrimitives { FORCEINLINE const Double& velocity(size_t iDim) const { return all(iDim+1); } FORCEINLINE const Double* velocity() const { return &velocity(0); } + VectorDbl<1> allTurb; + FORCEINLINE const Double& tke() const { return allTurb(0); } + /*--- Un-reconstructed variables. ---*/ FORCEINLINE Double& speedSound() { return all(nDim+4); } FORCEINLINE Double& laminarVisc() { return all(nDim+5); } FORCEINLINE Double& eddyVisc() { return all(nDim+6); } FORCEINLINE Double& thermalCond() { return all(nDim+7); } FORCEINLINE Double& cp() { return all(nDim+8); } + FORCEINLINE Double& tke() { return all(nDim+9); } FORCEINLINE const Double& speedSound() const { return all(nDim+4); } FORCEINLINE const Double& laminarVisc() const { return all(nDim+5); } FORCEINLINE const Double& eddyVisc() const { return all(nDim+6); } @@ -70,17 +74,22 @@ template struct CCompressibleConservatives { static constexpr size_t nDim = nDim_; static constexpr size_t nVar = nDim+2; + static constexpr size_t nVarTurb = 1; VectorDbl all; + VectorDbl allTurb; FORCEINLINE Double& density() { return all(0); } FORCEINLINE Double& rhoEnergy() { return all(nDim+1); } FORCEINLINE Double& momentum(size_t iDim) { return all(iDim+1); } + FORCEINLINE Double& rhoTke() { return allTurb(0); } FORCEINLINE const Double& density() const { return all(0); } FORCEINLINE const Double& rhoEnergy() const { return all(nDim+1); } FORCEINLINE const Double& momentum(size_t iDim) const { return all(iDim+1); } + FORCEINLINE const Double& rhoTke() const { return allTurb(0); } FORCEINLINE Double energy() const { return rhoEnergy() / density(); } FORCEINLINE const Double* momentum() const { return &momentum(0); } + FORCEINLINE Double tke() const { return rhoTke() / density(); } }; /*! @@ -94,6 +103,7 @@ FORCEINLINE CCompressibleConservatives compressibleConservatives(const CCo U.momentum(iDim) = V.density() * V.velocity(iDim); } U.rhoEnergy() = V.density() * V.enthalpy() - V.pressure(); + U.rhoTke() = V.density() * V.tke(); return U; } @@ -105,6 +115,7 @@ struct CRoeVariables { Double density; VectorDbl velocity; Double enthalpy; + Double tke; Double speedSound; Double projVel; }; @@ -124,6 +135,7 @@ FORCEINLINE CRoeVariables roeAveragedVariables(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.tke = (R*V.j.tke() + V.i.tke()) * D; roeAvg.speedSound = sqrt((gamma-1) * (roeAvg.enthalpy - 0.5*squaredNorm(roeAvg.velocity))); roeAvg.projVel = dot(roeAvg.velocity, normal); return roeAvg; diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index 4797293529f..e7472025900 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -432,7 +432,7 @@ void CFVMFlowSolverBase::Viscous_Residual_impl(unsigned long iEdge, CGeome CNumerics *numerics, CConfig *config) { const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); - const bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST); + const bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST && !(config->GetSSTParsedOptions().modified)); CVariable* turbNodes = nullptr; if (tkeNeeded) turbNodes = solver_container[TURB_SOL]->GetNodes(); @@ -2369,6 +2369,8 @@ void CFVMFlowSolverBase::Friction_Forces(const CGeometry* geometr const su2double Prandtl_Lam = config->GetPrandtl_Lam(); const bool energy = config->GetEnergy_Equation(); const bool QCR = config->GetSAParsedOptions().qcr2000; + TURB_FAMILY TurbFamily = TurbModelFamily(config->GetKind_Turb_Model()); + const bool SSTm = config->GetSSTParsedOptions().modified; const bool axisymmetric = config->GetAxisymmetric(); const bool roughwall = (config->GetnRoughWall() > 0); const bool nemo = config->GetNEMOProblem(); @@ -2452,7 +2454,12 @@ void CFVMFlowSolverBase::Friction_Forces(const CGeometry* geometr } /*--- Evaluate Tau ---*/ - CNumerics::ComputeStressTensor(nDim, Tau, Grad_Vel, Viscosity); + su2double tke = 0.0; + // if (TurbFamily == TURB_FAMILY::KW) { + + // } + // Don't knowhow to retrieve the TKE yet + CNumerics::ComputeStressTensor(nDim, Tau, Grad_Vel, Viscosity, Density, tke); /*--- If necessary evaluate the QCR contribution to Tau ---*/ diff --git a/SU2_CFD/include/variables/CTurbVariable.hpp b/SU2_CFD/include/variables/CTurbVariable.hpp index 48f2087b2bf..85d2c0a446f 100644 --- a/SU2_CFD/include/variables/CTurbVariable.hpp +++ b/SU2_CFD/include/variables/CTurbVariable.hpp @@ -43,6 +43,11 @@ class CTurbVariable : public CScalarVariable { static constexpr size_t MAXNVAR = 2; VectorType turb_index; VectorType intermittency; /*!< \brief Value of the intermittency for the trans. model. */ + VectorType Pk; + VectorType Pw; + VectorType Dk; + VectorType Dw; + VectorType PkLim; /*! * \brief Constructor of the class. @@ -100,5 +105,38 @@ class CTurbVariable : public CScalarVariable { */ inline void SetIntermittency(unsigned long iPoint, su2double val_intermittency) final { intermittency(iPoint) = val_intermittency; } + /*! + * \brief Set the intermittency of the transition model. + * \param[in] iPoint - Point index. + * \param[in] val_intermittency - New value of the intermittency. + */ + inline void SetProdDestr(unsigned long iPoint, su2double* val_ProdDestr) final { + Pk(iPoint) = val_ProdDestr[0]; + Dk(iPoint) = val_ProdDestr[1]; + Pw(iPoint) = val_ProdDestr[2]; + Dw(iPoint) = val_ProdDestr[3]; + PkLim(iPoint) = val_ProdDestr[4]; + } + + /*! + * \brief Set the intermittency of the transition model. + * \param[in] iPoint - Point index. + * \param[in] val_intermittency - New value of the intermittency. + */ + inline su2double GetProdTKE(unsigned long iPoint) const final { + return Pk(iPoint); + } + inline su2double GetDestrTKE(unsigned long iPoint) const final { + return Dk(iPoint); + } + inline su2double GetProdW(unsigned long iPoint) const final { + return Pw(iPoint); + } + inline su2double GetDestrW(unsigned long iPoint) const final { + return Dw(iPoint); + } + inline su2double GetPkLim(unsigned long iPoint) const final { + return PkLim(iPoint); + } }; diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index bc855bc5dca..84c77e8d463 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -1588,6 +1588,12 @@ class CVariable { */ inline virtual su2double GetGradient_Primitive(unsigned long iPoint, unsigned long iVar, unsigned long iDim) const { return 0.0; } + /*! + * \brief A virtual member. + * \return Value of the primitive variables gradient. + */ + inline virtual CMatrixView GetGradient_Primitive(unsigned long iPoint, unsigned long iVar=0) { return nullptr; } + /*! * \brief Get the primitive variable gradients for all points. * \return Reference to primitive variable gradient. @@ -1609,12 +1615,6 @@ class CVariable { */ inline virtual su2double GetLimiter_Primitive(unsigned long iPoint, unsigned long iVar) const { return 0.0; } - /*! - * \brief A virtual member. - * \return Value of the primitive variables gradient. - */ - inline virtual CMatrixView GetGradient_Primitive(unsigned long iPoint, unsigned long iVar=0) { return nullptr; } - /*! * \brief A virtual member. * \return Value of the primitive variables gradient. @@ -2350,4 +2350,11 @@ class CVariable { inline virtual const su2double *GetScalarSources(unsigned long iPoint) const { return nullptr; } inline virtual const su2double *GetScalarLookups(unsigned long iPoint) const { return nullptr; } + + inline virtual void SetProdDestr(unsigned long iPoint, su2double* val_ProdDestr) { } + inline virtual su2double GetProdTKE(unsigned long iPoint) const { return 0.0; } + inline virtual su2double GetDestrTKE(unsigned long iPoint) const { return 0.0; } + inline virtual su2double GetProdW(unsigned long iPoint) const { return 0.0; } + inline virtual su2double GetDestrW(unsigned long iPoint) const { return 0.0; } + inline virtual su2double GetPkLim(unsigned long iPoint) const { return 0.0; } }; diff --git a/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp b/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp index cc35fdc9cbc..c66b29ca373 100644 --- a/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp +++ b/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp @@ -179,6 +179,8 @@ void CFlowTractionInterface::GetDonor_Variable(CSolver *flow_solution, CGeometry unsigned long Vertex_Flow, unsigned long Point_Struct) { const auto Point_Flow = flow_geometry->vertex[Marker_Flow][Vertex_Flow]->GetNode(); + TURB_FAMILY TurbFamily = TurbModelFamily(flow_config->GetKind_Turb_Model()); + const bool SSTm = flow_config->GetSSTParsedOptions().modified; // Get the normal at the vertex: this normal goes inside the fluid domain. const su2double* Normal_Flow = flow_geometry->vertex[Marker_Flow][Vertex_Flow]->GetNormal(); @@ -205,7 +207,12 @@ void CFlowTractionInterface::GetDonor_Variable(CSolver *flow_solution, CGeometry su2double Viscosity = flow_nodes->GetLaminarViscosity(Point_Flow); su2double tau[3][3]; - CNumerics::ComputeStressTensor(nVar, tau, flow_nodes->GetVelocityGradient(Point_Flow), Viscosity); + su2double TKE = 0.0, Density = 0.0; + if (TurbFamily == TURB_FAMILY::KW && !SSTm) { + TKE = 0.0; // Have to find out how to pass the tke value if needed + Density = 0.0; + } + CNumerics::ComputeStressTensor(nVar, tau, flow_nodes->GetVelocityGradient(Point_Flow), Viscosity, Density, TKE); for (auto iVar = 0u; iVar < nVar; iVar++) { for (auto jVar = 0u; jVar < nVar; jVar++) { // Viscous component in the tn vector --> Units of force (non-dimensional). diff --git a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp index a2cdb893533..d012092fbd5 100644 --- a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp +++ b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp @@ -132,7 +132,7 @@ void CAvgGrad_Base::SetStressTensor(const su2double *val_primvar, if (sstParsedOptions.uq) { // laminar part - ComputeStressTensor(nDim, tau, val_gradprimvar+1, val_laminar_viscosity); + ComputeStressTensor(nDim, tau, val_gradprimvar+1, val_laminar_viscosity, Density, val_turb_ke); // add turbulent part which was perturbed for (unsigned short iDim = 0 ; iDim < nDim; iDim++) for (unsigned short jDim = 0 ; jDim < nDim; jDim++) @@ -140,7 +140,7 @@ void CAvgGrad_Base::SetStressTensor(const su2double *val_primvar, } 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)); + ComputeStressTensor(nDim, tau, val_gradprimvar+1, total_viscosity, Density, val_turb_ke); } } @@ -370,6 +370,7 @@ CNumerics::ResidualType<> CAvgGrad_Flow::ComputeResidual(const CConfig* config) AD::SetPreaccIn(Normal, nDim); unsigned short iVar, jVar, iDim; + const bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST) && !(config->GetSSTParsedOptions().modified); /*--- Normalized normal vector ---*/ @@ -403,6 +404,7 @@ CNumerics::ResidualType<> CAvgGrad_Flow::ComputeResidual(const CConfig* config) Mean_Laminar_Viscosity = 0.5*(Laminar_Viscosity_i + Laminar_Viscosity_j); Mean_Eddy_Viscosity = 0.5*(Eddy_Viscosity_i + Eddy_Viscosity_j); Mean_turb_ke = 0.5*(turb_ke_i + turb_ke_j); + const su2double Mean_turb_ke_ForStressTensor = tkeNeeded ? Mean_turb_ke : 0.0; /*--- Mean gradient approximation ---*/ @@ -434,7 +436,7 @@ CNumerics::ResidualType<> CAvgGrad_Flow::ComputeResidual(const CConfig* config) /*--- Get projected flux tensor (viscous residual) ---*/ - SetStressTensor(Mean_PrimVar, Mean_GradPrimVar, Mean_turb_ke, + SetStressTensor(Mean_PrimVar, Mean_GradPrimVar, Mean_turb_ke_ForStressTensor, Mean_Laminar_Viscosity, Mean_Eddy_Viscosity); if (config->GetSAParsedOptions().qcr2000) AddQCR(nDim, &Mean_GradPrimVar[1], tau); if (Mean_TauWall > 0) AddTauWall(UnitNormal, Mean_TauWall); @@ -555,6 +557,7 @@ CNumerics::ResidualType<> CAvgGradInc_Flow::ComputeResidual(const CConfig* confi AD::SetPreaccIn(Normal, nDim); unsigned short iVar, jVar, iDim; + const bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST) && !(config->GetSSTParsedOptions().modified); /*--- Normalized normal vector ---*/ @@ -590,6 +593,7 @@ CNumerics::ResidualType<> CAvgGradInc_Flow::ComputeResidual(const CConfig* confi Mean_Eddy_Viscosity = 0.5*(Eddy_Viscosity_i + Eddy_Viscosity_j); Mean_turb_ke = 0.5*(turb_ke_i + turb_ke_j); Mean_Thermal_Conductivity = 0.5*(Thermal_Conductivity_i + Thermal_Conductivity_j); + const su2double Mean_turb_ke_ForStressTensor = tkeNeeded ? Mean_turb_ke : 0.0; /*--- Mean gradient approximation ---*/ @@ -618,7 +622,7 @@ CNumerics::ResidualType<> CAvgGradInc_Flow::ComputeResidual(const CConfig* confi } /*--- Get projected flux tensor (viscous residual) ---*/ - SetStressTensor(Mean_PrimVar, Mean_GradPrimVar, Mean_turb_ke, + SetStressTensor(Mean_PrimVar, Mean_GradPrimVar, Mean_turb_ke_ForStressTensor, Mean_Laminar_Viscosity, Mean_Eddy_Viscosity); if (config->GetSAParsedOptions().qcr2000) AddQCR(nDim, &Mean_GradPrimVar[1], tau); if (Mean_TauWall > 0) AddTauWall(UnitNormal, Mean_TauWall); @@ -873,6 +877,7 @@ CNumerics::ResidualType<> CGeneralAvgGrad_Flow::ComputeResidual(const CConfig* c AD::SetPreaccIn(Normal, nDim); unsigned short iVar, jVar, iDim; + const bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST) && !(config->GetSSTParsedOptions().modified); /*--- Normalized normal vector ---*/ @@ -918,6 +923,7 @@ CNumerics::ResidualType<> CGeneralAvgGrad_Flow::ComputeResidual(const CConfig* c Mean_turb_ke = 0.5*(turb_ke_i + turb_ke_j); Mean_Thermal_Conductivity = 0.5*(Thermal_Conductivity_i + Thermal_Conductivity_j); Mean_Cp = 0.5*(Cp_i + Cp_j); + const su2double Mean_turb_ke_ForStressTensor = tkeNeeded ? Mean_turb_ke : 0.0; /*--- Mean gradient approximation ---*/ @@ -949,7 +955,7 @@ CNumerics::ResidualType<> CGeneralAvgGrad_Flow::ComputeResidual(const CConfig* c /*--- Get projected flux tensor (viscous residual) ---*/ - SetStressTensor(Mean_PrimVar, Mean_GradPrimVar, Mean_turb_ke, + SetStressTensor(Mean_PrimVar, Mean_GradPrimVar, Mean_turb_ke_ForStressTensor, Mean_Laminar_Viscosity, Mean_Eddy_Viscosity); if (config->GetSAParsedOptions().qcr2000) AddQCR(nDim, &Mean_GradPrimVar[1], tau); if (Mean_TauWall > 0) AddTauWall(UnitNormal, Mean_TauWall); diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp index 3c651269b02..a4274758a62 100644 --- a/SU2_CFD/src/output/CFlowOutput.cpp +++ b/SU2_CFD/src/output/CFlowOutput.cpp @@ -1252,6 +1252,14 @@ void CFlowOutput::SetVolumeOutputFieldsScalarSolution(const CConfig* config){ case TURB_FAMILY::KW: AddVolumeOutput("TKE", "Turb_Kin_Energy", "SOLUTION", "Turbulent kinetic energy"); AddVolumeOutput("DISSIPATION", "Omega", "SOLUTION", "Rate of dissipation"); + AddVolumeOutput("PROD_TKE", "Prod_TKE", "SOLUTION", "Production of turbulent kinetic energy"); + AddVolumeOutput("DESTR_TKE", "Destr_TKE", "SOLUTION", "Destruction of turbulent kinetic energy"); + AddVolumeOutput("PROD_TKE_LIM", "Prod_TKE_Lim", "SOLUTION", "Check if production limiter has been used for TKE"); + AddVolumeOutput("PROD_W", "Prod_W", "SOLUTION", "Production of rate of dissipation"); + AddVolumeOutput("DESTR_W", "Destr_W", "SOLUTION", "Destruction of rate of dissipation"); + AddVolumeOutput("CDkw", "CDkw", "SOLUTION", "Cross-Diffusion term"); + AddVolumeOutput("F1", "F1", "SOLUTION", "F1 blending function"); + AddVolumeOutput("F2", "F2", "SOLUTION", "F2 blending function"); break; case TURB_FAMILY::NONE: @@ -1479,6 +1487,8 @@ void CFlowOutput::SetVolumeOutputFieldsScalarMisc(const CConfig* config) { AddVolumeOutput("VORTICITY", "Vorticity", "VORTEX_IDENTIFICATION", "Value of the vorticity"); } AddVolumeOutput("Q_CRITERION", "Q_Criterion", "VORTEX_IDENTIFICATION", "Value of the Q-Criterion"); + AddVolumeOutput("WALL_DISTANCE", "Wall_Distance", "MISC", "Wall distance value"); + AddVolumeOutput("STRAIN_MAG", "Strain_Magnitude", "MISC", "Strain magnitude value"); } // Timestep info @@ -1512,6 +1522,8 @@ void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* con SetVolumeOutputValue("VORTICITY", iPoint, Node_Flow->GetVorticity(iPoint)[2]); } SetVolumeOutputValue("Q_CRITERION", iPoint, GetQCriterion(Node_Flow->GetVelocityGradient(iPoint))); + SetVolumeOutputValue("WALL_DISTANCE", iPoint, Node_Geo->GetWall_Distance(iPoint)); + SetVolumeOutputValue("STRAIN_MAG", iPoint, Node_Turb->GetStrainMag(iPoint)); } const bool limiter = (config->GetKind_SlopeLimit_Turb() != LIMITER::NONE); @@ -1528,6 +1540,14 @@ void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* con case TURB_FAMILY::KW: SetVolumeOutputValue("TKE", iPoint, Node_Turb->GetSolution(iPoint, 0)); SetVolumeOutputValue("DISSIPATION", iPoint, Node_Turb->GetSolution(iPoint, 1)); + SetVolumeOutputValue("PROD_TKE", iPoint, Node_Turb->GetProdTKE(iPoint)); + SetVolumeOutputValue("DESTR_TKE", iPoint, Node_Turb->GetDestrTKE(iPoint)); + SetVolumeOutputValue("PROD_TKE_LIM", iPoint, Node_Turb->GetPkLim(iPoint)); + SetVolumeOutputValue("PROD_W", iPoint, Node_Turb->GetProdW(iPoint)); + SetVolumeOutputValue("DESTR_W", iPoint, Node_Turb->GetDestrW(iPoint)); + SetVolumeOutputValue("CDkw", iPoint, Node_Turb->GetCrossDiff(iPoint)); + SetVolumeOutputValue("F1", iPoint, Node_Turb->GetF1blending(iPoint)); + SetVolumeOutputValue("F2", iPoint, Node_Turb->GetF2blending(iPoint)); SetVolumeOutputValue("RES_TKE", iPoint, turb_solver->LinSysRes(iPoint, 0)); SetVolumeOutputValue("RES_DISSIPATION", iPoint, turb_solver->LinSysRes(iPoint, 1)); if (limiter) { diff --git a/SU2_CFD/src/solvers/CAdjNSSolver.cpp b/SU2_CFD/src/solvers/CAdjNSSolver.cpp index 54e6c717b27..e3904110ab4 100644 --- a/SU2_CFD/src/solvers/CAdjNSSolver.cpp +++ b/SU2_CFD/src/solvers/CAdjNSSolver.cpp @@ -592,6 +592,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]; @@ -762,6 +766,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 e89955b4f68..776a36ca273 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, const bool time_stepping = (config->GetTime_Marching() == TIME_MARCHING::TIME_STEPPING); const bool adjoint = config->GetContinuous_Adjoint() || config->GetDiscrete_Adjoint(); const bool centered = config->GetKind_ConvNumScheme_Flow() == SPACE_CENTERED; + bool tkeNeeded = (rans && config->GetKind_Turb_Model() == TURB_MODEL::SST && !(config->GetSSTParsedOptions().modified)); int Unst_RestartIter = 0; unsigned long iPoint, iMarker, counter_local = 0, counter_global = 0; @@ -247,6 +248,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 approxiations ---*/ @@ -316,6 +318,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(); @@ -328,7 +331,7 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config, Solution[0] = Density_Inf; for (iDim = 0; iDim < nDim; iDim++) Solution[iDim+1] = Velocity_Inf[iDim]*Density_Inf; - Solution[nDim+1] = Energy_Inf*Density_Inf; + Solution[nDim+1] = Energy_Inf*Density_Inf; nodes->SetSolution(iPoint,Solution); nodes->SetSolution_Old(iPoint,Solution); counter_local++; @@ -810,7 +813,8 @@ void CEulerSolver::SetNondimensionalization(CConfig *config, unsigned short iMes bool viscous = config->GetViscous(); bool gravity = config->GetGravityForce(); bool turbulent = (config->GetKind_Turb_Model() != TURB_MODEL::NONE); - bool tkeNeeded = (turbulent && config->GetKind_Turb_Model() == TURB_MODEL::SST); + bool tkeNeeded = (turbulent && config->GetKind_Turb_Model() == TURB_MODEL::SST && !(config->GetSSTParsedOptions().modified)); + // bool tkeNeeded = (turbulent && config->GetKind_Turb_Model() == TURB_MODEL::SST); bool free_stream_temp = (config->GetKind_FreeStreamOption() == FREESTREAM_OPTION::TEMPERATURE_FS); bool reynolds_init = (config->GetKind_InitOption() == REYNOLDS); bool aeroelastic = config->GetAeroelastic_Simulation(); @@ -975,6 +979,19 @@ 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().newBC) { + su2double Omega_Freestream = 10 * ModVel_FreeStream / config->GetLDomain(); + Tke_FreeStream = Omega_Freestream*(Viscosity_FreeStream*config->GetTurb2LamViscRatio_FreeStream())/Density_FreeStream; + } else if (config->GetSSTParsedOptions().sust) { + su2double Omega_Freestream = 5*ModVel_FreeStream/config->GetLength_Reynolds(); + /*--- not good ---*/ + Tke_FreeStream = pow(10.0, -6) * ModVel_FreeStream*ModVel_FreeStream; /*--- Equals a freestream turbulence intensity of 0.08165%. This should be the default one, not hard-coding the freestream TI. ---*/ + } + + /*-- Compute the freestream energy. ---*/ + + if (tkeNeeded) { Energy_FreeStream += Tke_FreeStream; }; + } else { @@ -985,9 +1002,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. ---*/ @@ -1048,15 +1064,30 @@ void CEulerSolver::SetNondimensionalization(CConfig *config, unsigned short iMes 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/(Viscosity_FreeStream*config->GetTurb2LamViscRatio_FreeStream()); - config->SetOmega_FreeStream(Omega_FreeStream); - Omega_FreeStreamND = Density_FreeStreamND*Tke_FreeStreamND/(Viscosity_FreeStreamND*config->GetTurb2LamViscRatio_FreeStream()); + + if (config->GetSSTParsedOptions().newBC) { + 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) { + Omega_FreeStream = 5*ModVel_FreeStream/config->GetLength_Reynolds(); + Omega_FreeStreamND = 5*ModVel_FreeStreamND / config->GetLength_Reynolds(); // Should it be non-dimensionalized by what? Dimensions are 1/s + + /*--- not good ---*/ + Tke_FreeStream = pow(10.0, -6) * ModVel_FreeStream*ModVel_FreeStream; /*--- Equals a freestream turbulence intensity of 0.08165%. This should be the default one, not hard-coding the freestream TI. ---*/ + Tke_FreeStreamND = pow(10.0, -6) * ModVel_FreeStreamND*ModVel_FreeStreamND; /*--- Equals a freestream turbulence intensity of 0.08165%. This should be the default one, not hard-coding the freestream TI. ---*/ + } + + 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) { @@ -4777,7 +4808,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 && !(config->GetSSTParsedOptions().modified)); + CVariable* turbNodes = nullptr; + if (tkeNeeded) turbNodes = solver_container[TURB_SOL]->GetNodes(); auto *Normal = new su2double[nDim]; @@ -4927,7 +4960,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(); + // if (tkeNeeded) Energy += GetTke_Inf(); + if (tkeNeeded) Energy += turbNodes->GetSolution(iPoint,0); /*--- Store new primitive state for computing the flux. ---*/ @@ -5020,7 +5054,9 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, const bool viscous = config->GetViscous(), implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT), gravity = config->GetGravityForce(), - tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST); + tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST && !(config->GetSSTParsedOptions().modified)); + CVariable* turbNodes = nullptr; + if (tkeNeeded) turbNodes = solver_container[TURB_SOL]->GetNodes(); su2double **P_Tensor = new su2double*[nVar], **invP_Tensor = new su2double*[nVar]; @@ -5066,7 +5102,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); @@ -5120,7 +5157,8 @@ 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(); + // if (tkeNeeded) Energy_e += GetTke_Inf(); + if (tkeNeeded) Energy_e += turbNodes->GetSolution(iPoint,0); break; case STATIC_SUPERSONIC_INFLOW_PT: @@ -5146,7 +5184,8 @@ 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(); + // if (tkeNeeded) Energy_e += GetTke_Inf(); + if (tkeNeeded) Energy_e += turbNodes->GetSolution(iPoint,0); break; case STATIC_SUPERSONIC_INFLOW_PD: @@ -5173,7 +5212,8 @@ 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(); + // if (tkeNeeded) Energy_e += GetTke_Inf(); + if (tkeNeeded) Energy_e += turbNodes->GetSolution(iPoint,0); break; case DENSITY_VELOCITY: @@ -5206,6 +5246,7 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, Velocity2_e = GeometryToolbox::SquaredNorm(nDim, Velocity_e); Energy_e = GetFluidModel()->GetStaticEnergy() + 0.5*Velocity2_e; + if (tkeNeeded) Energy_e += turbNodes->GetSolution(iPoint,0); break; default: @@ -5410,7 +5451,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)); @@ -5462,7 +5503,9 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain unsigned short nSpanWiseSections = geometry->GetnSpanWiseSections(config->GetMarker_All_TurbomachineryFlag(val_marker)); bool viscous = config->GetViscous(); bool gravity = (config->GetGravityForce()); - bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST); + bool tkeNeeded = ((config->GetKind_Turb_Model() == TURB_MODEL::SST) && !(config->GetSSTParsedOptions().modified)); + CVariable* turbNodes = nullptr; + if (tkeNeeded) turbNodes = solver_container[TURB_SOL]->GetNodes(); su2double *Normal, *turboNormal, *UnitNormal, *FlowDirMix, FlowDirMixMag, *turboVelocity; Normal = new su2double[nDim]; @@ -5536,6 +5579,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); @@ -5581,11 +5625,12 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain turboVelocity[iDim] = sqrt(Velocity2_e)*Flow_Dir[iDim]; ComputeBackVelocity(turboVelocity,turboNormal, Velocity_e, config->GetMarker_All_TurbomachineryFlag(val_marker),config->GetKind_TurboMachinery(iZone)); StaticEnthalpy_e = Enthalpy_e - 0.5 * Velocity2_e; + if(tkeNeeded) StaticEnthalpy_e -= turbNodes->GetSolution(iPoint, 0); GetFluidModel()->SetTDState_hs(StaticEnthalpy_e, Entropy_e); 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: @@ -5616,10 +5661,12 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain ComputeBackVelocity(turboVelocity,turboNormal, Velocity_e, config->GetMarker_All_TurbomachineryFlag(val_marker),config->GetKind_TurboMachinery(iZone)); StaticEnthalpy_e = Enthalpy_e - 0.5 * Velocity2_e; + if(tkeNeeded) StaticEnthalpy_e -= turbNodes->GetSolution(iPoint, 0); GetFluidModel()->SetTDState_hs(StaticEnthalpy_e, Entropy_e); Density_e = GetFluidModel()->GetDensity(); StaticEnergy_e = GetFluidModel()->GetStaticEnergy(); Energy_e = StaticEnergy_e + 0.5 * Velocity2_e; + if(tkeNeeded) Energy_e += turbNodes->GetSolution(iPoint, 0); // if (tkeNeeded) Energy_e += GetTke_Inf(); break; @@ -5638,6 +5685,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: @@ -5655,6 +5703,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; @@ -5672,6 +5721,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: @@ -6893,7 +6943,7 @@ void CEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, const su2double Gas_Constant = config->GetGas_ConstantND(); const auto Kind_Inlet = config->GetKind_Inlet(); const auto Marker_Tag = config->GetMarker_All_TagBound(val_marker); - const bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST); + const bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST && !(config->GetSSTParsedOptions().modified)); /*--- Loop over all the vertices on this boundary marker ---*/ @@ -7194,7 +7244,10 @@ void CEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, su2double Gas_Constant = config->GetGas_ConstantND(); string Marker_Tag = config->GetMarker_All_TagBound(val_marker); bool gravity = (config->GetGravityForce()); - bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST); + bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST && !(config->GetSSTParsedOptions().modified)); + + CVariable* turbNodes = nullptr; + if (tkeNeeded) turbNodes = solver_container[TURB_SOL]->GetNodes(); auto *Normal = new su2double[nDim]; @@ -7277,7 +7330,8 @@ 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 += GetTke_Inf(); + if (tkeNeeded) Energy += turbNodes->GetSolution(iPoint,0); /*--- Conservative variables, using the derived quantities ---*/ V_outlet[0] = Pressure / ( Gas_Constant * Density); @@ -7360,11 +7414,14 @@ void CEulerSolver::BC_Supersonic_Inlet(CGeometry *geometry, CSolver **solver_con const su2double Gas_Constant = config->GetGas_ConstantND(); const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); const auto Marker_Tag = config->GetMarker_All_TagBound(val_marker); - const bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST); + const bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST && !(config->GetSSTParsedOptions().modified)); + CVariable* turbNodes = nullptr; + if (tkeNeeded) turbNodes = solver_container[TURB_SOL]->GetNodes(); /*--- Supersonic inlet flow: there are no outgoing characteristics, so all flow variables can be imposed at the inlet. ---*/ + /*--- Loop over all the vertices on this boundary marker ---*/ SU2_OMP_FOR_DYN(OMP_MIN_SIZE) for (auto iVertex = 0ul; iVertex < geometry->nVertex[val_marker]; iVertex++) { const auto iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); @@ -7387,10 +7444,12 @@ 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) Energy += GetTke_Inf(); + const su2double Energy_woTKE = Pressure / (Density * Gamma_Minus_One) + 0.5 * Velocity2; + su2double Energy = Energy_woTKE; + + if (tkeNeeded) Energy = Energy_woTKE + turbNodes->GetSolution(iPoint,0); - /*--- Primitive variables, using the derived quantities. ---*/ + /*--- Allocate the value at the inlet ---*/ auto* V_inlet = GetCharacPrimVar(val_marker, iVertex); V_inlet[prim_idx.Temperature()] = Temperature; diff --git a/SU2_CFD/src/solvers/CFEM_DG_EulerSolver.cpp b/SU2_CFD/src/solvers/CFEM_DG_EulerSolver.cpp index 29c53116c84..15e4f9698cc 100644 --- a/SU2_CFD/src/solvers/CFEM_DG_EulerSolver.cpp +++ b/SU2_CFD/src/solvers/CFEM_DG_EulerSolver.cpp @@ -1055,15 +1055,30 @@ void CFEM_DG_EulerSolver::SetNondimensionalization(CConfig *config, 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()); + + Omega_FreeStream = Density_FreeStream*Tke_FreeStream/(Viscosity_FreeStream*config->GetTurb2LamViscRatio_FreeStream()); + Omega_FreeStreamND = Density_FreeStreamND*Tke_FreeStreamND/(Viscosity_FreeStreamND*config->GetTurb2LamViscRatio_FreeStream()); + + if (config->GetSSTParsedOptions().newBC) { + 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) { + Omega_FreeStream = 5*ModVel_FreeStream/config->GetLength_Reynolds(); + Omega_FreeStreamND = 5*ModVel_FreeStreamND / config->GetLength_Reynolds(); // Should it be non-dimensionalized by what? Dimensions are 1/s + + /*--- not good ---*/ + Tke_FreeStream = pow(10.0, -6) * ModVel_FreeStream*ModVel_FreeStream; /*--- Equals a freestream turbulence intensity of 0.08165%. This should be the default one, not hard-coding the freestream TI. ---*/ + Tke_FreeStreamND = pow(10.0, -6) * ModVel_FreeStreamND*ModVel_FreeStreamND; /*--- Equals a freestream turbulence intensity of 0.08165%. This should be the default one, not hard-coding the freestream TI. ---*/ + } + + config->SetTke_FreeStream(Tke_FreeStream); config->SetTke_FreeStreamND(Tke_FreeStreamND); - Omega_FreeStream = Density_FreeStream*Tke_FreeStream/max((Viscosity_FreeStream*config->GetTurb2LamViscRatio_FreeStream()), 1.e-25); config->SetOmega_FreeStream(Omega_FreeStream); - - Omega_FreeStreamND = Density_FreeStreamND*Tke_FreeStreamND/max((Viscosity_FreeStreamND*config->GetTurb2LamViscRatio_FreeStream()), 1.e-25); 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 6aaf09a5183..499c36cfd60 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -452,15 +452,30 @@ void CIncEulerSolver::SetNondimensionalization(CConfig *config, unsigned short i 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/(Viscosity_FreeStream*config->GetTurb2LamViscRatio_FreeStream()); - config->SetOmega_FreeStream(Omega_FreeStream); - Omega_FreeStreamND = Density_FreeStreamND*Tke_FreeStreamND/(Viscosity_FreeStreamND*config->GetTurb2LamViscRatio_FreeStream()); + + if (config->GetSSTParsedOptions().newBC) { + 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) { + Omega_FreeStream = 5*ModVel_FreeStream/config->GetLength_Reynolds(); + Omega_FreeStreamND = 5*ModVel_FreeStreamND / config->GetLength_Reynolds(); // Should it be non-dimensionalized by what? Dimensions are 1/s + + /*--- not good ---*/ + Tke_FreeStream = pow(10.0, -6) * ModVel_FreeStream*ModVel_FreeStream; /*--- Equals a freestream turbulence intensity of 0.08165%. This should be the default one, not hard-coding the freestream TI. ---*/ + Tke_FreeStreamND = pow(10.0, -6) * ModVel_FreeStreamND*ModVel_FreeStreamND; /*--- Equals a freestream turbulence intensity of 0.08165%. This should be the default one, not hard-coding the freestream TI. ---*/ + } + + 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/CIncNSSolver.cpp b/SU2_CFD/src/solvers/CIncNSSolver.cpp index 41cedcdb5f5..6e0f17ca0bf 100644 --- a/SU2_CFD/src/solvers/CIncNSSolver.cpp +++ b/SU2_CFD/src/solvers/CIncNSSolver.cpp @@ -287,7 +287,7 @@ unsigned long CIncNSSolver::SetPrimitive_Variables(CSolver **solver_container, c const TURB_MODEL turb_model = config->GetKind_Turb_Model(); const SPECIES_MODEL species_model = config->GetKind_Species_Model(); - bool tkeNeeded = (turb_model == TURB_MODEL::SST); + bool tkeNeeded = (turb_model == TURB_MODEL::SST && !(config->GetSSTParsedOptions().modified)); AD::StartNoSharedReading(); diff --git a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp index 1a41a04f66a..3d236080f47 100644 --- a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp @@ -1187,15 +1187,30 @@ 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/(Viscosity_FreeStream*config->GetTurb2LamViscRatio_FreeStream()); - config->SetOmega_FreeStream(Omega_FreeStream); - Omega_FreeStreamND = Density_FreeStreamND*Tke_FreeStreamND/(Viscosity_FreeStreamND*config->GetTurb2LamViscRatio_FreeStream()); + + if (config->GetSSTParsedOptions().newBC) { + 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) { + Omega_FreeStream = 5*ModVel_FreeStream/config->GetLength_Reynolds(); + Omega_FreeStreamND = 5*ModVel_FreeStreamND / config->GetLength_Reynolds(); // Should it be non-dimensionalized by what? Dimensions are 1/s + + /*--- not good ---*/ + Tke_FreeStream = pow(10.0, -6) * ModVel_FreeStream*ModVel_FreeStream; /*--- Equals a freestream turbulence intensity of 0.08165%. This should be the default one, not hard-coding the freestream TI. ---*/ + Tke_FreeStreamND = pow(10.0, -6) * ModVel_FreeStreamND*ModVel_FreeStreamND; /*--- Equals a freestream turbulence intensity of 0.08165%. This should be the default one, not hard-coding the freestream TI. ---*/ + } + + 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/CNEMONSSolver.cpp b/SU2_CFD/src/solvers/CNEMONSSolver.cpp index 01d63385739..227b2a5ae80 100644 --- a/SU2_CFD/src/solvers/CNEMONSSolver.cpp +++ b/SU2_CFD/src/solvers/CNEMONSSolver.cpp @@ -958,7 +958,7 @@ void CNEMONSSolver::BC_Smoluchowski_Maxwell(CGeometry *geometry, /*--- Compute wall shear stress (using the stress tensor) ---*/ su2double Tau[MAXNDIM][MAXNDIM] = {{0.0}}; - CNumerics::ComputeStressTensor(nDim, Tau, Grad_PrimVar+VEL_INDEX, Viscosity); + CNumerics::ComputeStressTensor(nDim, Tau, Grad_PrimVar+VEL_INDEX, Viscosity); // no need for TKE as it is on the wall su2double TauTangent[MAXNDIM] = {0.0}; GeometryToolbox::TangentProjection(nDim,Tau,UnitNormal,TauTangent); diff --git a/SU2_CFD/src/solvers/CNSSolver.cpp b/SU2_CFD/src/solvers/CNSSolver.cpp index 127ff3f5b08..36c900c9a43 100644 --- a/SU2_CFD/src/solvers/CNSSolver.cpp +++ b/SU2_CFD/src/solvers/CNSSolver.cpp @@ -133,7 +133,7 @@ unsigned long CNSSolver::SetPrimitive_Variables(CSolver **solver_container, cons unsigned long nonPhysicalPoints = 0; const TURB_MODEL turb_model = config->GetKind_Turb_Model(); - const bool tkeNeeded = (turb_model == TURB_MODEL::SST); + const bool tkeNeeded = (turb_model == TURB_MODEL::SST && !(config->GetSSTParsedOptions().modified)); AD::StartNoSharedReading(); @@ -799,6 +799,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 ) ---*/ @@ -904,7 +907,9 @@ 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); + su2double tke = 0.0; + if (tkeNeeded) tke = Node_Turb->GetSolution(iPoint, 0); + CNumerics::ComputeStressTensor(nDim, tau, nodes->GetVelocityGradient(iPoint), Lam_Visc_Wall, nodes->GetDensity(iPoint), tke); su2double TauTangent[MAXNDIM] = {0.0}; GeometryToolbox::TangentProjection(nDim, tau, UnitNormal, TauTangent); diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index e446c3f1133..6fff0f93408 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -4037,6 +4037,7 @@ void CSolver::ComputeVertexTractions(CGeometry *geometry, const CConfig *config) if (viscous_flow) { const su2double Viscosity = base_nodes->GetLaminarViscosity(iPoint); su2double Tau[3][3]; + // Here I do not care for modified version of SST since TKE at wall is 0 CNumerics::ComputeStressTensor(nDim, Tau, base_nodes->GetVelocityGradient(iPoint), Viscosity); for (unsigned short iDim = 0; iDim < nDim; iDim++) { auxForce[iDim] += GeometryToolbox::DotProduct(nDim, Tau[iDim], Normal); diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp index 02b91c8ab2c..def030f630c 100644 --- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp @@ -101,11 +101,13 @@ CTurbSSTSolver::CTurbSSTSolver(CGeometry *geometry, CConfig *config, unsigned sh constants[8] = constants[4]/constants[6] - constants[2]*0.41*0.41/sqrt(constants[6]); //alfa_1 constants[9] = constants[5]/constants[6] - constants[3]*0.41*0.41/sqrt(constants[6]); //alfa_2 constants[10] = 20.0; // production limiter constant + if (sstParsedOptions.prodLim) constants[10] = config->GetProdLimConst(); } else { /* SST-V2003 */ constants[8] = 5.0 / 9.0; //gamma_1 constants[9] = 0.44; //gamma_2 constants[10] = 10.0; // production limiter constant + if (sstParsedOptions.prodLim) constants[10] = config->GetProdLimConst(); } /*--- Far-field flow state quantities and initialization. ---*/ @@ -122,6 +124,15 @@ CTurbSSTSolver::CTurbSSTSolver(CGeometry *geometry, CConfig *config, unsigned sh su2double kine_Inf = 3.0/2.0*(VelMag2*Intensity*Intensity); su2double omega_Inf = rhoInf*kine_Inf/(muLamInf*viscRatio); + if (sstParsedOptions.newBC) { + omega_Inf = 10 * sqrt(VelMag2) / config->GetLDomain(); + kine_Inf = omega_Inf*(muLamInf*viscRatio)/rhoInf; + } else if (sstParsedOptions.sust) { + omega_Inf = 5*sqrt(VelMag2)/config->GetLength_Reynolds(); + /*--- not good ---*/ + kine_Inf = pow(10.0, -6) * VelMag2; /*--- Equals a freestream turbulence intensity of 0.08165%. This should be the default one, not hard-coding the freestream TI. ---*/ + } + Solution_Inf[0] = kine_Inf; Solution_Inf[1] = omega_Inf; @@ -240,7 +251,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); @@ -377,6 +388,14 @@ void CTurbSSTSolver::Source_Residual(CGeometry *geometry, CSolver **solver_conta nodes->SetIntermittency(iPoint, numerics->GetIntermittencyEff()); } + su2double ProdDestr[5]; + ProdDestr[0] = numerics->GetProdDest(0); + ProdDestr[1] = numerics->GetProdDest(1); + ProdDestr[2] = numerics->GetProdDest(2); + ProdDestr[3] = numerics->GetProdDest(3); + ProdDestr[4] = numerics->GetProdDest(4); + nodes->SetProdDestr(iPoint, ProdDestr); + /*--- Subtract residual and the Jacobian ---*/ LinSysRes.SubtractBlock(iPoint, residual); @@ -625,8 +644,14 @@ 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.newBC) { + 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); + } + } /*--- Set the turbulent variable states. Use free-stream SST diff --git a/SU2_CFD/src/variables/CTurbSSTVariable.cpp b/SU2_CFD/src/variables/CTurbSSTVariable.cpp index d1fd98c8a14..6da5333bdfa 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/SU2_CFD/src/variables/CTurbVariable.cpp b/SU2_CFD/src/variables/CTurbVariable.cpp index 46b686766a7..68083833de9 100644 --- a/SU2_CFD/src/variables/CTurbVariable.cpp +++ b/SU2_CFD/src/variables/CTurbVariable.cpp @@ -35,4 +35,12 @@ CTurbVariable::CTurbVariable(unsigned long npoint, unsigned long ndim, unsigned turb_index.resize(nPoint) = su2double(1.0); intermittency.resize(nPoint) = su2double(1.0); + if (TurbModelFamily(config->GetKind_Turb_Model()) == TURB_FAMILY::KW) { + Pk.resize(nPoint) = su2double(1.0); + Pw.resize(nPoint) = su2double(1.0); + Dk.resize(nPoint) = su2double(1.0); + Dw.resize(nPoint) = su2double(1.0); + PkLim.resize(nPoint) = su2double(1.0); + } + } diff --git a/externals/codi b/externals/codi index 762ba7698e3..c6b039e5c9e 160000 --- a/externals/codi +++ b/externals/codi @@ -1 +1 @@ -Subproject commit 762ba7698e3ceaa1b17aa299421ddd418f00b823 +Subproject commit c6b039e5c9edb7675f90ffc725f9dd8e66571264 diff --git a/externals/medi b/externals/medi index 7d550831e0e..ab3a7688f6d 160000 --- a/externals/medi +++ b/externals/medi @@ -1 +1 @@ -Subproject commit 7d550831e0e233a85b9d9af9c181d7ecb2929946 +Subproject commit ab3a7688f6d518f8d940eb61a341d89f51922ba4 diff --git a/externals/opdi b/externals/opdi index a6b9655c240..8c897988172 160000 --- a/externals/opdi +++ b/externals/opdi @@ -1 +1 @@ -Subproject commit a6b9655c240af2a35454a61727e5bbbbaa3a425f +Subproject commit 8c89798817253abb017d857a0ae7f0520187645c