diff --git a/Common/include/geometry/dual_grid/CPoint.hpp b/Common/include/geometry/dual_grid/CPoint.hpp index c297b1b5df9..ae725fea258 100644 --- a/Common/include/geometry/dual_grid/CPoint.hpp +++ b/Common/include/geometry/dual_grid/CPoint.hpp @@ -108,10 +108,11 @@ class CPoint { su2vector ClosestWall_Elem; /*!< \brief Element index of closest wall element, for givenrank, zone and marker index. */ - su2activevector SharpEdge_Distance; /*!< \brief Distance to a sharp edge. */ - su2activevector Curvature; /*!< \brief Value of the surface curvature (SU2_GEO). */ - su2activevector MaxLength; /*!< \brief The maximum cell-center to cell-center length. */ - su2activevector RoughnessHeight; /*!< \brief Roughness of the nearest wall. */ + su2activevector SharpEdge_Distance; /*!< \brief Distance to a sharp edge. */ + su2activevector Curvature; /*!< \brief Value of the surface curvature (SU2_GEO). */ + su2activevector MaxLength; /*!< \brief The maximum cell-center to cell-center length. */ + su2activevector RoughnessHeight; /*!< \brief Roughness of the nearest wall. */ + su2activematrix Normals; /*!< \brief Normal of the nearest wall element. */ su2matrix AD_InputIndex; /*!< \brief Indices of Coord variables in the adjoint vector. */ su2matrix @@ -484,6 +485,27 @@ class CPoint { } inline void SetWall_Distance(unsigned long iPoint, su2double distance) { Wall_Distance(iPoint) = distance; } + /*! + * \brief Get the index of the closest wall element. + * \param[in] iPoint - Index of the point. + * \param[out] ClosestWall_Elem - ID of the closest element on a wall boundary. + */ + inline unsigned long GetClosestWall_Elem(unsigned long iPoint) {return ClosestWall_Elem(iPoint);} + + /*! + * \brief Get the marker of the closest wall marker. + * \param[in] iPoint - Index of the point. + * \param[out] ClosestWall_Marker - MarkerID of the closest wall boundary. + */ + inline unsigned long GetClosestWall_Marker(unsigned long iPoint) {return ClosestWall_Marker(iPoint);} + + /*! + * \brief Get the rank of the closest wall marker. + * \param[in] iPoint - Index of the point. + * \param[out] ClosestWall_Rank - RankID of the closest wall boundary. + */ + inline unsigned long GetClosestWall_Rank(unsigned long iPoint) {return ClosestWall_Rank(iPoint);} + /*! * \brief Get the value of the distance to the nearest wall. * \param[in] iPoint - Index of the point. @@ -506,6 +528,24 @@ class CPoint { */ inline su2double GetRoughnessHeight(unsigned long iPoint) const { return RoughnessHeight(iPoint); } + /*! + * \brief Set the value of the normal of the nearest wall element. + * \param[in] iPoint - Index of the point. + * \param[in] normal - Value of the normal. + */ + template + inline void SetNormal(unsigned long iPoint, Normals_type const&normal) { + for (unsigned long iDim = 0; iDim < nDim; iDim++) + Normals(iPoint,iDim) = normal[iDim]; + } + + /*! + * \brief Set the value of the normal of the nearest wall element. + * \param[in] iPoint - Index of the point. + * \return normal to the normal of the nearest wall element. + */ + inline su2double *GetNormal(unsigned long iPoint) { return Normals[iPoint]; } + /*! * \brief Set the value of the distance to a sharp edge. * \param[in] iPoint - Index of the point. @@ -906,4 +946,22 @@ class CPoint { } } } + + + /*! + * \brief Set wall normal according to stored closest wall information. + * \param[in] normals - Mapping [rank][zone][marker][element] -> normal + */ + template + void SetWallNormals(Normals_type const&normals){ + for (unsigned long iPoint=0; iPoint= 0) + SetNormal(iPoint, normals[rankID][zoneID][markerID][elementID]); + } + } + }; diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index d1f59abbd79..5f99b6c502b 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -1215,27 +1215,37 @@ static const MapType Trans_Model_Map = { * \brief LM Options */ enum class LM_OPTIONS { - NONE, /*!< \brief No option / default. */ - LM2015, /*!< \brief Cross-flow corrections. */ - MALAN, /*!< \brief Kind of transition correlation model (Malan). */ - SULUKSNA, /*!< \brief Kind of transition correlation model (Suluksna). */ - KRAUSE, /*!< \brief Kind of transition correlation model (Krause). */ - KRAUSE_HYPER, /*!< \brief Kind of transition correlation model (Krause hypersonic). */ - MEDIDA_BAEDER,/*!< \brief Kind of transition correlation model (Medida-Baeder). */ - MEDIDA, /*!< \brief Kind of transition correlation model (Medida). */ - MENTER_LANGTRY, /*!< \brief Kind of transition correlation model (Menter-Langtry). */ - DEFAULT /*!< \brief Kind of transition correlation model (Menter-Langtry if SST, MALAN if SA). */ + NONE, /*!< \brief No option / default. */ + CROSSFLOW, /*!< \brief Cross-flow corrections. */ + SLM, /*!< \brief Simplified version. */ + PRODLIM, /*!< \brief Add production term to Pk. */ + MALAN, /*!< \brief Kind of transition correlation model (Malan). */ + SULUKSNA, /*!< \brief Kind of transition correlation model (Suluksna). */ + KRAUSE, /*!< \brief Kind of transition correlation model (Krause). */ + KRAUSE_HYPER, /*!< \brief Kind of transition correlation model (Krause hypersonic). */ + MEDIDA_BAEDER, /*!< \brief Kind of transition correlation model (Medida-Baeder). */ + MEDIDA, /*!< \brief Kind of transition correlation model (Medida). */ + MENTER_LANGTRY, /*!< \brief Kind of transition correlation model (Menter-Langtry). */ + MENTER_SLM, /*!< \brief Kind of transition correlation model (Menter Simplified LM model). */ + CODER_SLM, /*!< \brief Kind of transition correlation model (Coder Simplified LM model). */ + MOD_EPPLER_SLM, /*!< \brief Kind of transition correlation model (Modified Eppler Simplified LM model). */ + DEFAULT /*!< \brief Kind of transition correlation model (Menter-Langtry if SST, MALAN if SA). */ }; static const MapType LM_Options_Map = { MakePair("NONE", LM_OPTIONS::NONE) - MakePair("LM2015", LM_OPTIONS::LM2015) + MakePair("CROSSFLOW", LM_OPTIONS::CROSSFLOW) + MakePair("SLM", LM_OPTIONS::SLM) + MakePair("PRODLIM", LM_OPTIONS::PRODLIM) MakePair("MALAN", LM_OPTIONS::MALAN) MakePair("SULUKSNA", LM_OPTIONS::SULUKSNA) MakePair("KRAUSE", LM_OPTIONS::KRAUSE) MakePair("KRAUSE_HYPER", LM_OPTIONS::KRAUSE_HYPER) MakePair("MEDIDA_BAEDER", LM_OPTIONS::MEDIDA_BAEDER) MakePair("MENTER_LANGTRY", LM_OPTIONS::MENTER_LANGTRY) + MakePair("MENTER_SLM", LM_OPTIONS::MENTER_SLM) + MakePair("CODER_SLM", LM_OPTIONS::CODER_SLM) + MakePair("MOD_EPPLER_SLM", LM_OPTIONS::MOD_EPPLER_SLM) MakePair("DEFAULT", LM_OPTIONS::DEFAULT) }; @@ -1253,13 +1263,26 @@ enum class TURB_TRANS_CORRELATION { DEFAULT /*!< \brief Kind of transition correlation model (Menter-Langtry if SST, MALAN if SA). */ }; +/*! + * \brief Types of transition correlations for Simplified LM model + */ +enum class TURB_TRANS_CORRELATION_SLM { + MENTER_SLM, /*!< \brief Kind of transition correlation model (Menter Simplified LM model). */ + CODER_SLM, /*!< \brief Kind of transition correlation model (Coder Simplified LM model). */ + MOD_EPPLER_SLM, /*!< \brief Kind of transition correlation model (Modified Eppler Simplified LM model). */ + DEFAULT /*!< \brief Kind of transition correlation model. */ +}; + /*! * \brief Structure containing parsed LM options. */ struct LM_ParsedOptions { LM_OPTIONS version = LM_OPTIONS::NONE; /*!< \brief LM base model. */ - bool LM2015 = false; /*!< \brief Use cross-flow corrections. */ + bool SLM = false; /*!< \brief Use simplified version. */ + bool ProdLim = false; /*!< \brief Add production term to Pk. */ + bool CrossFlow = false; /*!< \brief Use cross-flow corrections. */ TURB_TRANS_CORRELATION Correlation = TURB_TRANS_CORRELATION::DEFAULT; + TURB_TRANS_CORRELATION_SLM Correlation_SLM = TURB_TRANS_CORRELATION_SLM::DEFAULT; }; /*! @@ -1277,7 +1300,10 @@ inline LM_ParsedOptions ParseLMOptions(const LM_OPTIONS *LM_Options, unsigned sh return std::find(LM_Options, lm_options_end, option) != lm_options_end; }; - LMParsedOptions.LM2015 = IsPresent(LM_OPTIONS::LM2015); + LMParsedOptions.SLM = IsPresent(LM_OPTIONS::SLM); + LMParsedOptions.ProdLim = IsPresent(LM_OPTIONS::PRODLIM); + + LMParsedOptions.CrossFlow = IsPresent(LM_OPTIONS::CROSSFLOW); int NFoundCorrelations = 0; if (IsPresent(LM_OPTIONS::MALAN)) { @@ -1309,6 +1335,20 @@ inline LM_ParsedOptions ParseLMOptions(const LM_OPTIONS *LM_Options, unsigned sh NFoundCorrelations++; } + int NFoundCorrelations_SLM = 0; + if (IsPresent(LM_OPTIONS::MENTER_SLM)) { + LMParsedOptions.Correlation_SLM = TURB_TRANS_CORRELATION_SLM::MENTER_SLM; + NFoundCorrelations_SLM++; + } + if (IsPresent(LM_OPTIONS::CODER_SLM)) { + LMParsedOptions.Correlation_SLM = TURB_TRANS_CORRELATION_SLM::CODER_SLM; + NFoundCorrelations_SLM++; + } + if (IsPresent(LM_OPTIONS::MOD_EPPLER_SLM)) { + LMParsedOptions.Correlation_SLM = TURB_TRANS_CORRELATION_SLM::MOD_EPPLER_SLM; + NFoundCorrelations_SLM++; + } + if (NFoundCorrelations > 1) { SU2_MPI::Error("Two correlations selected for LM_OPTIONS. Please choose only one.", CURRENT_FUNCTION); } @@ -1321,6 +1361,10 @@ inline LM_ParsedOptions ParseLMOptions(const LM_OPTIONS *LM_Options, unsigned sh } } + if (LMParsedOptions.Correlation_SLM == TURB_TRANS_CORRELATION_SLM::DEFAULT){ + LMParsedOptions.Correlation_SLM = TURB_TRANS_CORRELATION_SLM::MENTER_SLM; + } + return LMParsedOptions; } diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 43f77fc9c52..f73ac5970fc 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -3495,9 +3495,9 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i if (Kind_Trans_Model == TURB_TRANS_MODEL::LM) { lmParsedOptions = ParseLMOptions(LM_Options, nLM_Options, rank, Kind_Turb_Model); - /*--- Check if problem is 2D and LM2015 has been selected ---*/ - if (lmParsedOptions.LM2015 && val_nDim == 2) { - SU2_MPI::Error("LM2015 is available only for 3D problems", CURRENT_FUNCTION); + /*--- Check if problem is 2D and CrossFlow has been selected ---*/ + if (lmParsedOptions.CrossFlow && val_nDim == 2) { + SU2_MPI::Error("Cross-flow corrections are available only for 3D problems", CURRENT_FUNCTION); } } @@ -6230,11 +6230,27 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { switch (Kind_Trans_Model) { case TURB_TRANS_MODEL::NONE: break; case TURB_TRANS_MODEL::LM: { - cout << "Transition model: Langtry and Menter's 4 equation model"; - if (lmParsedOptions.LM2015) { - cout << " w/ cross-flow corrections (2015)" << endl; + int NTurbEqs = 0; + switch (Kind_Turb_Model) { + case TURB_MODEL::SA: NTurbEqs = 1; break; + case TURB_MODEL::SST: NTurbEqs = 2; break; + case TURB_MODEL::NONE: SU2_MPI::Error("No turbulence model has been selected but LM transition model is active.", CURRENT_FUNCTION); break; + } + if (!lmParsedOptions.SLM) { + int NEquations = 2; + cout << "Transition model: Langtry and Menter's "<< NEquations+NTurbEqs <<" equation model"; } else { - cout << " (2009)" << endl; + int NEquations = 1; + cout << "Transition model: Simplified Langtry and Menter's "<< NEquations+NTurbEqs <<" equation model"; + } + if (lmParsedOptions.CrossFlow) { + cout << " w/ cross-flow corrections" << endl; + } else { + if (!lmParsedOptions.SLM) { + cout << " (2009)" << endl; + } else { + cout << " (2015)" << endl; + } } break; } @@ -6258,6 +6274,13 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { } break; } + cout << "Correlation Functions for Simplified LM model: "; + switch (lmParsedOptions.Correlation_SLM) { + case TURB_TRANS_CORRELATION_SLM::CODER_SLM: cout << "Coder et al. (2012)" << endl; break; + case TURB_TRANS_CORRELATION_SLM::MOD_EPPLER_SLM: cout << "Modified Eppler (from Coder et al. 2012)" << endl; break; + case TURB_TRANS_CORRELATION_SLM::MENTER_SLM: + case TURB_TRANS_CORRELATION_SLM::DEFAULT: cout << "Menter et al. (2015)" << endl; break; + } } cout << "Hybrid RANS/LES: "; switch (Kind_HybridRANSLES) { diff --git a/Common/src/geometry/CGeometry.cpp b/Common/src/geometry/CGeometry.cpp index b1c77372141..4c3a149389f 100644 --- a/Common/src/geometry/CGeometry.cpp +++ b/Common/src/geometry/CGeometry.cpp @@ -3969,5 +3969,81 @@ void CGeometry::ComputeWallDistance(const CConfig* const* config_container, CGeo } } } + + + su2vector>> WallNormal_container; + WallNormal_container.resize(nZone) = su2vector>(); + for (int iZone = 0; iZone < nZone; iZone++){ + const CConfig* config = config_container[iZone]; + const CGeometry* geometry = geometry_container[iZone][iInst][MESH_0]; + WallNormal_container[iZone].resize(geometry->GetnMarker()); + for (auto iMarker = 0; iMarker < geometry->GetnMarker(); iMarker++){ + if (config->GetViscous_Wall(iMarker)) { + WallNormal_container[iZone][iMarker].resize(geometry->GetnElem_Bound(iMarker), 3); + + // cout << "geometry->nVertex[iMarker] = " << geometry->nVertex[iMarker] << endl; + for (auto iElem = 0u; iElem < geometry->GetnElem_Bound(iMarker); iElem++) { + su2vector NormalHere; + NormalHere.resize(3) = su2double(0.0); + + for (unsigned short iNode = 0; iNode < geometry->bound[iMarker][iElem]->GetnNodes(); iNode++) { + // Extract global coordinate of the node + unsigned long iPointHere = geometry->bound[iMarker][iElem]->GetNode(iNode); + long iVertexHere = geometry->nodes->GetVertex(iPointHere, iMarker); + for (auto iDim = 0u; iDim < 3; iDim++) + NormalHere[iDim] += geometry->vertex[iMarker][iVertexHere]->GetNormal(iDim); + } + + for (auto iDim = 0u; iDim < 3; iDim++) + NormalHere[iDim] /= geometry->bound[iMarker][iElem]->GetnNodes(); + + su2double NormalMag = 0.0; + for (auto iDim = 0u; iDim < 3; iDim++) + NormalMag += NormalHere[iDim]*NormalHere[iDim]; + NormalMag = sqrt(NormalMag); + + for (auto iDim = 0u; iDim < 3; iDim++) + NormalHere[iDim] /= NormalMag; + + for (auto iDim = 0u; iDim < 3; iDim++) + WallNormal_container[iZone][iMarker](iElem, iDim) = NormalHere[iDim]; + + } + } else { + WallNormal_container[iZone][iMarker].resize(1, 3) = su2double(0.0); + } + } + } + + auto normal_i = + make_pair(nZone, [config_container,geometry_container,iInst,WallNormal_container](unsigned long iZone){ + const CConfig* config = config_container[iZone]; + const CGeometry* geometry = geometry_container[iZone][iInst][MESH_0]; + const auto nMarker = geometry->GetnMarker(); + const auto WallNormal = WallNormal_container[iZone]; + + return make_pair( nMarker, [config,geometry,WallNormal](unsigned long iMarker){ + auto nElem_Bou = geometry->GetnElem_Bound(iMarker); + if (!config->GetViscous_Wall(iMarker)) nElem_Bou = 1; + + return make_pair(nElem_Bou, [WallNormal,iMarker](unsigned long iElem){ + const auto dimensions = 3; + + return make_pair(dimensions, [WallNormal,iMarker,iElem](unsigned short iDim){ + + return WallNormal[iMarker](iElem, iDim); + }); + }); + }); + }); + + NdFlattener<4>Normals_Local(normal_i); + NdFlattener<5> Normals_global(Nd_MPI_Environment(), Normals_Local); + + + // use it to update roughnesses + for(int jZone=0; jZonenodes->SetWallNormals(Normals_global); + } } } diff --git a/Common/src/geometry/dual_grid/CPoint.cpp b/Common/src/geometry/dual_grid/CPoint.cpp index f11b7dd8945..ed8715e0105 100644 --- a/Common/src/geometry/dual_grid/CPoint.cpp +++ b/Common/src/geometry/dual_grid/CPoint.cpp @@ -132,6 +132,9 @@ void CPoint::FullAllocation(unsigned short imesh, const CConfig* config) { RoughnessHeight.resize(npoint) = su2double(0.0); SharpEdge_Distance.resize(npoint) = su2double(0.0); + + Normals.resize(npoint, 3) = su2double(0.0); + } void CPoint::SetElems(const vector >& elemsMatrix) { Elem = CCompressedSparsePatternL(elemsMatrix); } diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index 450112df8f7..9c80ec71ff3 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -729,6 +729,41 @@ class CNumerics { */ su2double GetIntermittencyEff() const { return intermittency_eff_i; } + /*! + * \brief Get the value of the Transition Momentum Thickness Reynolds number from correlations. + * \param[out] Corr_Rec_i - Value of the Transition Momentum Thickness Reynolds number at point i. + */ + inline virtual su2double GetCorr_Rec() {return 0.0;} + + /*! + * \brief Get the value of the Momentum Thickness Reynolds number. + * \param[out] re_t - Value of the Momentum Thickness Reynolds number at point i. + */ + inline virtual su2double GetRe_t() {return 0.0;} + inline virtual su2double GetTu() {return 0.0;} + inline virtual su2double GetLambda_theta() {return 0.0;} + inline virtual su2double Getduds() {return 0.0;} + inline virtual su2double GetRe_v() {return 0.0;} + inline virtual su2double GetProd() {return 0.0;} + inline virtual su2double GetDestr() {return 0.0;} + inline virtual su2double GetF_onset1() {return 0.0;} + inline virtual su2double GetF_onset2() {return 0.0;} + inline virtual su2double GetF_onset3() {return 0.0;} + inline virtual su2double GetF_onset() {return 0.0;} + + /*! + * \brief Set the value of the F2 blending function into SLM transition model. + * \param[in] val_F2 - F2 blending function. + */ + inline virtual void SetF2(su2double val_F2) {} + + /*! + * \brief Set the gradient of the auxiliary variables. + * \param[in] val_auxvar_grad_i - Gradient of the auxiliary variable at point i. + * \param[in] val_auxvar_grad_j - Gradient of the auxiliary variable at point j. + */ + inline virtual void SetAuxVar(su2double val_AuxVar) {} + /*! * \brief Set the gradient of the auxiliary variables. * \param[in] val_auxvar_grad_i - Gradient of the auxiliary variable at point i. diff --git a/SU2_CFD/include/numerics/turbulent/transition/trans_convection.hpp b/SU2_CFD/include/numerics/turbulent/transition/trans_convection.hpp index 669c3a3dfa4..0b5c70701d6 100644 --- a/SU2_CFD/include/numerics/turbulent/transition/trans_convection.hpp +++ b/SU2_CFD/include/numerics/turbulent/transition/trans_convection.hpp @@ -38,3 +38,10 @@ template using CUpwSca_TransLM = CUpwSca_TurbSST; +/*! + * \class CUpwSca_TransSLM + * \brief Re-use the SA convective fluxes for the scalar upwind discretization of Simplified LM transition model equations. + * \ingroup ConvDiscr + */ +template +using CUpwSca_TransSLM = CUpwSca_TurbSA; diff --git a/SU2_CFD/include/numerics/turbulent/transition/trans_correlations.hpp b/SU2_CFD/include/numerics/turbulent/transition/trans_correlations.hpp index b8102120701..817d832f245 100644 --- a/SU2_CFD/include/numerics/turbulent/transition/trans_correlations.hpp +++ b/SU2_CFD/include/numerics/turbulent/transition/trans_correlations.hpp @@ -25,6 +25,7 @@ */ #pragma once +//#include /*! * \class TransLMCorrelations @@ -189,4 +190,109 @@ class TransLMCorrelations { return F_length1; } + + /*! + * \brief Compute Re_theta_c from correlations for the Simplified LM model. + * \param[in] Tu - Turbulence intensity. + * \param[in] du_ds - Streamwise velocity gradient. + * \param[out] rethetac - Corrected value for Re_theta. + */ + su2double ReThetaC_Correlations_SLM(const su2double Tu_L, const su2double lambda_theta, const su2double wall_dist, const su2double VorticityMag, const su2double VelocityMag) const { + + su2double rethetac = 0.0; + + switch (options.Correlation_SLM) { + case TURB_TRANS_CORRELATION_SLM::MENTER_SLM: { + + /*-- Thwaites parameter ---*/ + su2double lambda_theta_local = lambda_theta; + + /*-- Function to sensitize the transition onset to the streamwise pressure gradient ---*/ + su2double FPG = 0.0; + const su2double C_PG1 = 14.68; + const su2double C_PG1_lim = 1.5; + const su2double C_PG2 = -7.34; + const su2double C_PG2_lim = 3.0; + const su2double C_PG3 = 0.0; + if (lambda_theta_local >= 0.0) { + FPG = min(1+ C_PG1 * lambda_theta_local, C_PG1_lim); + } else { + const su2double FirstTerm = C_PG2 * lambda_theta_local; + const su2double SecondTerm = C_PG3 * min(lambda_theta_local + 0.0681, 0.0); + FPG = min(1 + FirstTerm + SecondTerm, C_PG2_lim); + } + + FPG = max(FPG, 0.0); + + const su2double C_TU1 = 100.0; + const su2double C_TU2 = 1000.0; + const su2double C_TU3 = 1.0; + rethetac = C_TU1 + C_TU2 * exp(-C_TU3 * Tu_L * FPG); + + break; + } case TURB_TRANS_CORRELATION_SLM::CODER_SLM: { + + /*-- Local pressure gradient parameter ---*/ + const su2double H_c = max(min(wall_dist * VorticityMag / VelocityMag, 1.1542), 0.3823); + + /*-- Thwaites parameter ---*/ + su2double lambda_theta_local = 0.0; + const su2double H_c_delta = 0.587743 - H_c; + if ( H_c >= 0.587743 ) { + const su2double FirstTerm = 0.1919 * pow(H_c_delta, 3.0); + const su2double SecondTerm = 0.4182 * pow(H_c_delta, 2.0); + const su2double ThirdTerm = 0.2959 * H_c_delta; + lambda_theta_local = FirstTerm + SecondTerm + ThirdTerm; + } else { + const su2double FirstTerm = 4.7596 * pow(H_c_delta, 3.0); + const su2double SecondTerm = -0.3837 * pow(H_c_delta, 2.0); + const su2double ThirdTerm = 0.3575 * H_c_delta; + lambda_theta_local = FirstTerm + SecondTerm + ThirdTerm; + } + + /*-- Function to sensitize the transition onset to the streamwise pressure gradient ---*/ + su2double FPG = 0.0; + if (lambda_theta_local <= 0.0) { + const su2double FirstTerm = -12.986 * lambda_theta_local; + const su2double SecondTerm = -123.66 * pow(lambda_theta_local, 2.0); + const su2double ThirdTerm = -405.689 * pow(lambda_theta_local, 3.0); + FPG = 1 - (FirstTerm + SecondTerm + ThirdTerm) * exp(-pow(Tu_L/1.5,1.5)); + } else { + FPG = 1 + 0.275 * (1 - exp(-35.0 * lambda_theta_local)) * exp(-Tu_L/0.5); + } + + // This is not reported in the paper + //FPG = max(FPG, 0.0); + + if (Tu_L <= 1.3) { + const su2double FirstTerm = -589.428 * Tu_L; + const su2double SecondTerm = 0.2196 / max(pow(Tu_L, 2.0), 1e-12); + rethetac = 1173.51 + FirstTerm + SecondTerm; + } else { + rethetac = 331.50 * pow(Tu_L-0.5658, -0.671); + } + rethetac = rethetac * FPG; + + break; + } case TURB_TRANS_CORRELATION_SLM::MOD_EPPLER_SLM: { + + /*-- Local pressure gradient parameter ---*/ + const su2double H_c = max(min(wall_dist * VorticityMag / VelocityMag, 1.1542), 0.3823); + + /*-- H_32 Shape factor --*/ + const su2double H_32 = 1.515095 + 0.2041 * pow((1.1542 - H_c), 2.0956); + + rethetac = exp(127.94 * pow((H_32-1.515095), 2.0) + 6.774224); + + break; + } + case TURB_TRANS_CORRELATION_SLM::DEFAULT: + SU2_MPI::Error("Transition correlation for Simplified LM model is set to DEFAULT but no default value has ben set in the code.", + CURRENT_FUNCTION); + break; + } + + return rethetac; + } + }; diff --git a/SU2_CFD/include/numerics/turbulent/transition/trans_diffusion.hpp b/SU2_CFD/include/numerics/turbulent/transition/trans_diffusion.hpp index 50153606860..4beb7d374eb 100644 --- a/SU2_CFD/include/numerics/turbulent/transition/trans_diffusion.hpp +++ b/SU2_CFD/include/numerics/turbulent/transition/trans_diffusion.hpp @@ -103,3 +103,72 @@ class CAvgGrad_TransLM final : public CAvgGrad_Scalar { } }; + + +/*! + * \class CAvgGrad_TransSLM + * \brief Class for computing viscous term using average of gradient with correction (Simplified LM transition model). + * \ingroup ViscDiscr + * \author S. Kang. + */ +template +class CAvgGrad_TransSLM final : public CAvgGrad_Scalar { +private: + using Base = CAvgGrad_Scalar; + using Base::Laminar_Viscosity_i; + using Base::Laminar_Viscosity_j; + using Base::Eddy_Viscosity_i; + using Base::Eddy_Viscosity_j; + using Base::Density_i; + using Base::Density_j; + using Base::ScalarVar_i; + using Base::ScalarVar_j; + using Base::Proj_Mean_GradScalarVar; + using Base::proj_vector_ij; + using Base::Flux; + using Base::Jacobian_i; + using Base::Jacobian_j; + + /*! + * \brief Adds any extra variables to AD + */ + void ExtraADPreaccIn() override {} + + /*! + * \brief LM transition model specific steps in the ComputeResidual method + * \param[in] config - Definition of the particular problem. + */ + void FinishResidualCalc(const CConfig* config) override { + const bool implicit = config->GetKind_TimeIntScheme() == EULER_IMPLICIT; + + /*--- Compute mean effective dynamic viscosity ---*/ + const su2double diff_i_gamma = Laminar_Viscosity_i + Eddy_Viscosity_i; + const su2double diff_j_gamma = Laminar_Viscosity_j + Eddy_Viscosity_j; + + const su2double diff_gamma = 0.5*(diff_i_gamma + diff_j_gamma); + + Flux[0] = diff_gamma*Proj_Mean_GradScalarVar[0]; + + /*--- For Jacobians -> Use of TSL (Thin Shear Layer) approx. to compute derivatives of the gradients ---*/ + if (implicit) { + const su2double proj_on_rho_i = proj_vector_ij/Density_i; + Jacobian_i[0][0] = -diff_gamma*proj_on_rho_i; + + const su2double proj_on_rho_j = proj_vector_ij/Density_j; + Jacobian_j[0][0] = diff_gamma*proj_on_rho_j; + } + } + +public: + /*! + * \brief Constructor of the class. + * \param[in] val_nDim - Number of dimensions of the problem. + * \param[in] val_nVar - Number of variables of the problem. + * \param[in] correct_grad - Whether to correct gradient for skewness. + * \param[in] config - Definition of the particular problem. + */ + CAvgGrad_TransSLM(unsigned short val_nDim, unsigned short val_nVar, bool correct_grad, const CConfig* config) + : CAvgGrad_Scalar(val_nDim, val_nVar, correct_grad, config){ + } + +}; \ No newline at end of file diff --git a/SU2_CFD/include/numerics/turbulent/transition/trans_sources.hpp b/SU2_CFD/include/numerics/turbulent/transition/trans_sources.hpp index c1388c2df9c..d00ea1bb691 100644 --- a/SU2_CFD/include/numerics/turbulent/transition/trans_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/transition/trans_sources.hpp @@ -56,6 +56,17 @@ class CSourcePieceWise_TransLM final : public CNumerics { TURB_FAMILY TurbFamily; su2double hRoughness; + su2double Re_v_Here; + su2double Corr_Rec_Here; + su2double Prod_Here = 0.0; + su2double Destr_Here = 0.0; + su2double F_onset1_Here = 0.0; + su2double F_onset2_Here = 0.0; + su2double F_onset3_Here = 0.0; + su2double F_onset_Here = 0.0; + su2double lambda_theta_Here = 0.0; + su2double duds_Here = 0.0; + su2double IntermittencySep = 1.0; su2double IntermittencyEff = 1.0; @@ -92,7 +103,7 @@ class CSourcePieceWise_TransLM final : public CNumerics { * \return A lightweight const-view (read-only) of the residual/flux and Jacobians. */ ResidualType<> ComputeResidual(const CConfig* config) override { - /*--- ScalarVar[0] = k, ScalarVar[0] = w, TransVar[0] = gamma, and TransVar[0] = ReThetaT ---*/ + /*--- ScalarVar[0] = k, ScalarVar[0] = w, TransVar[0] = gamma, and TransVar[1] = ReThetaT ---*/ /*--- dU/dx = PrimVar_Grad[1][0] ---*/ AD::StartPreacc(); AD::SetPreaccIn(StrainMag_i); @@ -113,7 +124,7 @@ class CSourcePieceWise_TransLM final : public CNumerics { const su2double vel_v = V_i[1 + idx.Velocity()]; const su2double vel_w = (nDim == 3) ? V_i[2 + idx.Velocity()] : 0.0; - const su2double Velocity_Mag = sqrt(vel_u * vel_u + vel_v * vel_v + vel_w * vel_w); + const su2double Velocity_Mag = max(sqrt(vel_u * vel_u + vel_v * vel_v + vel_w * vel_w), 1e-20); AD::SetPreaccIn(V_i[idx.Density()], V_i[idx.LaminarViscosity()], V_i[idx.EddyViscosity()]); @@ -135,6 +146,8 @@ class CSourcePieceWise_TransLM final : public CNumerics { /*--- Corr_RetC correlation*/ const su2double Corr_Rec = TransCorrelations.ReThetaC_Correlations(Tu, TransVar_i[1]); + // AGGIUNTO PER DEBUG + Corr_Rec_Here = Corr_Rec; /*--- F_length correlation*/ const su2double Corr_F_length = TransCorrelations.FLength_Correlations(Tu, TransVar_i[1]); @@ -154,6 +167,9 @@ class CSourcePieceWise_TransLM final : public CNumerics { if (TurbFamily == TURB_FAMILY::SA) R_t = Eddy_Viscosity_i / Laminar_Viscosity_i; const su2double Re_v = Density_i * dist_i * dist_i * StrainMag_i / Laminar_Viscosity_i; + // AGGIUNTO PER DEBUG + Re_v_Here = Re_v; + const su2double F_onset1 = Re_v / (2.193 * Corr_Rec); su2double F_onset2 = 1.0; su2double F_onset3 = 1.0; @@ -166,6 +182,11 @@ class CSourcePieceWise_TransLM final : public CNumerics { F_onset3 = max(2.0 - pow(R_t / 2.5, 3.0), 0.0); } const su2double F_onset = max(F_onset2 - F_onset3, 0.0); + // AGGIUNTO PER DEBUG + F_onset1_Here = F_onset1; + F_onset2_Here = F_onset2; + F_onset3_Here = F_onset3; + F_onset_Here = F_onset; /*-- Gradient of velocity magnitude ---*/ @@ -186,7 +207,7 @@ class CSourcePieceWise_TransLM final : public CNumerics { /*-- Calculate blending function f_theta --*/ su2double time_scale = 500.0 * Laminar_Viscosity_i / Density_i / Velocity_Mag / Velocity_Mag; - if (options.LM2015) + if (options.CrossFlow) time_scale = min(time_scale, Density_i * LocalGridLength_i * LocalGridLength_i / (Laminar_Viscosity_i + Eddy_Viscosity_i)); const su2double theta_bl = TransVar_i[1] * Laminar_Viscosity_i / Density_i / Velocity_Mag; @@ -206,14 +227,14 @@ class CSourcePieceWise_TransLM final : public CNumerics { const su2double f_turb = exp(-pow(R_t / 4, 4)); su2double f_theta_2 = 0.0; - if (options.LM2015) + if (options.CrossFlow) f_theta_2 = min(f_wake * exp(-pow(dist_i / delta, 4.0)), 1.0); /*--- Corr_Ret correlation*/ const su2double Corr_Ret_lim = 20.0; su2double f_lambda = 1.0; - su2double Retheta_Error = 200.0, Retheta_old = 0.0; + su2double Retheta_Error = 200.0, Retheta_old = 1.0; su2double lambda = 0.0; su2double Corr_Ret = 20.0; @@ -245,9 +266,13 @@ class CSourcePieceWise_TransLM final : public CNumerics { Retheta_old = Corr_Ret; } + // DEBUG + lambda_theta_Here = lambda; + duds_Here = du_ds; + /*-- Corr_RetT_SCF Correlations--*/ su2double ReThetat_SCF = 0.0; - if (options.LM2015) { + if (options.CrossFlow) { su2double VelocityNormalized[3]; VelocityNormalized[0] = vel_u / Velocity_Mag; VelocityNormalized[1] = vel_v / Velocity_Mag; @@ -292,13 +317,17 @@ class CSourcePieceWise_TransLM final : public CNumerics { /*-- destruction term of Intermeittency(Gamma) --*/ const su2double Dg = c_a2 * Density_i * VorticityMag * TransVar_i[0] * f_turb * (c_e2 * TransVar_i[0] - 1.0); + // DEBUG + Prod_Here = Pg; + Destr_Here = Dg; + /*-- production term of ReThetaT --*/ const su2double PRethetat = c_theta * Density_i / time_scale * (Corr_Ret - TransVar_i[1]) * (1.0 - f_theta); /*-- destruction term of ReThetaT --*/ // It should not be with the minus sign but I put for consistency su2double DRethetat = 0.0; - if (options.LM2015) + if (options.CrossFlow) DRethetat = -c_theta * (Density_i / time_scale) * c_CF * min(ReThetat_SCF - TransVar_i[1], 0.0) * f_theta_2; /*--- Source ---*/ @@ -313,7 +342,7 @@ class CSourcePieceWise_TransLM final : public CNumerics { Jacobian_i[0][1] = 0.0; Jacobian_i[1][0] = 0.0; Jacobian_i[1][1] = -c_theta / time_scale * (1.0 - f_theta) * Volume; - if (options.LM2015 && ReThetat_SCF - TransVar_i[1] < 0) + if (options.CrossFlow && ReThetat_SCF - TransVar_i[1] < 0) Jacobian_i[1][1] += (c_theta / time_scale) * c_CF * f_theta_2 * Volume; } @@ -322,4 +351,293 @@ class CSourcePieceWise_TransLM final : public CNumerics { return ResidualType<>(Residual, Jacobian_i, nullptr); } + + inline su2double GetRe_v() override {return Re_v_Here;} + inline su2double GetCorr_Rec() override {return Corr_Rec_Here;} + inline su2double GetProd() override {return Prod_Here;} + inline su2double GetDestr() override {return Destr_Here;} + inline su2double GetF_onset1() override {return F_onset1_Here;} + inline su2double GetF_onset2() override {return F_onset2_Here;} + inline su2double GetF_onset3() override {return F_onset3_Here;} + inline su2double GetF_onset() override {return F_onset_Here;} + inline su2double GetLambda_theta() override {return lambda_theta_Here;} + inline su2double Getduds() override {return duds_Here;} + }; + +/*! + * \class CSourcePieceWise_TranSLM + * \brief Class for integrating the source terms of the Simplified LM transition model equations. + * \ingroup SourceDiscr + * \author S. Kang. + */ +template +class CSourcePieceWise_TransSLM final : public CNumerics { + private: + const FlowIndices idx; /*!< \brief Object to manage the access to the flow primitives. */ + + const LM_ParsedOptions options; + + /*--- LM Closure constants ---*/ + const su2double c_e1 = 1.0; + const su2double c_a1 = 2.0; + const su2double c_e2 = 50.0; + const su2double c_a2 = 0.06; + const su2double sigmaf = 1.0; + const su2double s1 = 2.0; + const su2double c_theta = 0.03; + const su2double c_CF = 0.6; + const su2double sigmat = 2.0; + + TURB_FAMILY TurbFamily; + su2double hRoughness; + + su2double IntermittencySep = 1.0; + su2double IntermittencyEff = 1.0; + + su2double Re_t; + su2double Corr_Rec = 1.0; + su2double AuxVar; + su2double F2; + su2double Tu_Here = 0.0; + su2double duds_Here = 0.0; + su2double lambda_theta_Here = 0.0; + su2double Re_v_Here = 0.0; + su2double Prod_Here = 0.0; + su2double Destr_Here = 0.0; + su2double F_onset1_Here = 0.0; + su2double F_onset2_Here = 0.0; + su2double F_onset3_Here = 0.0; + su2double F_onset_Here = 0.0; + + su2double Residual; + su2double* Jacobian_i; + su2double Jacobian_Buffer; // Static storage for the Jacobian (which needs to be pointer for return type). + + TransLMCorrelations TransCorrelations; + + public: + /*! + * \brief Constructor of the class. + * \param[in] val_nDim - Number of dimensions of the problem. + * \param[in] val_nVar - Number of variables of the problem. + * \param[in] config - Definition of the particular problem. + */ + CSourcePieceWise_TransSLM(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config) + : CNumerics(val_nDim, 1, config), idx(val_nDim, config->GetnSpecies()), options(config->GetLMParsedOptions()){ + /*--- "Allocate" the Jacobian using the static buffer. ---*/ + Jacobian_i = &Jacobian_Buffer; + + TurbFamily = TurbModelFamily(config->GetKind_Turb_Model()); + + hRoughness = config->GethRoughness(); + + TransCorrelations.SetOptions(options); + + } + + /*! + * \brief Residual for source term integration. + * \param[in] config - Definition of the particular problem. + * \return A lightweight const-view (read-only) of the residual/flux and Jacobians. + */ + ResidualType<> ComputeResidual(const CConfig* config) override { + /*--- ScalarVar[0] = k, ScalarVar[0] = w, TransVar[0] = gamma ---*/ + /*--- dU/dx = PrimVar_Grad[1][0] ---*/ + AD::StartPreacc(); + AD::SetPreaccIn(StrainMag_i); + AD::SetPreaccIn(ScalarVar_i, nVar); + AD::SetPreaccIn(ScalarVar_Grad_i, nVar, nDim); + AD::SetPreaccIn(TransVar_i, nVar); + AD::SetPreaccIn(TransVar_Grad_i, nVar, nDim); + AD::SetPreaccIn(Volume); + AD::SetPreaccIn(dist_i); + AD::SetPreaccIn(&V_i[idx.Velocity()], nDim); + AD::SetPreaccIn(PrimVar_Grad_i, nDim + idx.Velocity(), nDim); + AD::SetPreaccIn(Vorticity_i, 3); + + su2double VorticityMag = + sqrt(Vorticity_i[0] * Vorticity_i[0] + Vorticity_i[1] * Vorticity_i[1] + Vorticity_i[2] * Vorticity_i[2]); + + const su2double vel_u = V_i[idx.Velocity()]; + const su2double vel_v = V_i[1 + idx.Velocity()]; + const su2double vel_w = (nDim == 3) ? V_i[2 + idx.Velocity()] : 0.0; + + su2double Velocity[nDim]; + Velocity[0] = vel_u; + Velocity[1] = vel_v; + if(nDim == 3) Velocity[2] = vel_w; + + + const su2double Velocity_Mag = max(sqrt(vel_u * vel_u + vel_v * vel_v + vel_w * vel_w), 1e-20); + + AD::SetPreaccIn(V_i[idx.Density()], V_i[idx.LaminarViscosity()], V_i[idx.EddyViscosity()]); + + Density_i = V_i[idx.Density()]; + Laminar_Viscosity_i = V_i[idx.LaminarViscosity()]; + Eddy_Viscosity_i = V_i[idx.EddyViscosity()]; + + Residual = 0.0; + Jacobian_i[0] = 0.0; + + if (dist_i > 1e-10) { + su2double Tu_L = 1.0; + if (TurbFamily == TURB_FAMILY::KW) Tu_L = min(100.0 * sqrt(2.0 * ScalarVar_i[0] / 3.0) / (ScalarVar_i[1]*dist_i), 100.0); + // if (TurbFamily == TURB_FAMILY::KW) Tu_L = min(100.0 * sqrt(2.0 * ScalarVar_i[0] / 3.0) / (Velocity_Mag), 100.0); + if (TurbFamily == TURB_FAMILY::SA) Tu_L = config->GetTurbulenceIntensity_FreeStream() * 100; + + Tu_Here = Tu_L; + + /*--- F_length ---*/ + su2double F_length = 100.0; + + /*--- F_onset ---*/ + su2double R_t = 1.0; + if (TurbFamily == TURB_FAMILY::KW) R_t = Density_i * ScalarVar_i[0] / (Laminar_Viscosity_i * ScalarVar_i[1]); + if (TurbFamily == TURB_FAMILY::SA) R_t = Eddy_Viscosity_i / Laminar_Viscosity_i; + + /*-- Gradient of velocity magnitude ---*/ + + // su2double du_ds = 0.0; + // for (int i = 0; i < nDim; i++) + // for (int j = 0; j < nDim; j++) + // du_ds = du_ds + Velocity[i]*Velocity[j]*PrimVar_Grad_i[i][j]; + + // du_ds = du_ds/(Velocity_Mag*Velocity_Mag); + + // const su2double lambda_theta = -7.57e-3 * du_ds * dist_i * dist_i * Density_i / Laminar_Viscosity_i + 0.0128; + const su2double lambda_theta = max(min(-7.57e-3 * AuxVar * dist_i * dist_i * Density_i / Laminar_Viscosity_i + 0.0128, 1.0), -1.0); + // const su2double lambda_theta = du_ds * dist_i * dist_i * Density_i / Laminar_Viscosity_i; + // duds_Here = du_ds; + duds_Here = AuxVar; + lambda_theta_Here = lambda_theta; + + // const su2double lambda_theta = 7.57e-3 * AuxVar * dist_i * dist_i * Density_i / Laminar_Viscosity_i + 0.0128; + + /*--- Corr_RetC correlation*/ + Re_t = TransCorrelations.ReThetaC_Correlations_SLM(Tu_L, lambda_theta, dist_i, VorticityMag, Velocity_Mag); + Corr_Rec = Re_t; // If the MENTER_SLM correlation is used then they are the same thing + + if (options.Correlation_SLM == TURB_TRANS_CORRELATION_SLM::CODER_SLM || options.Correlation_SLM == TURB_TRANS_CORRELATION_SLM::MOD_EPPLER_SLM) { + // If these correlations are used, then the value of Corr_Rec has to be used instead of the TransVar[1] of the original LM model + F_length = TransCorrelations.FLength_Correlations(Tu_L, Re_t); + Corr_Rec = TransCorrelations.ReThetaC_Correlations(Tu_L, Re_t); + } + + + const su2double Re_v = Density_i * dist_i * dist_i * StrainMag_i / Laminar_Viscosity_i; + Re_v_Here = Re_v; + const su2double F_onset1 = Re_v / (2.2 * Corr_Rec); + su2double F_onset2 = 1.0; + su2double F_onset3 = 1.0; + if (TurbFamily == TURB_FAMILY::KW) { + F_onset2 = min(F_onset1, 2.0); + F_onset3 = max(1.0 - pow(R_t / 3.5, 3.0), 0.0); + } + if (TurbFamily == TURB_FAMILY::SA) { + F_onset2 = min(max(F_onset1, pow(F_onset1, 4.0)), 4.0); + F_onset3 = max(2.0 - pow(R_t / 2.5, 3.0), 0.0); + } + su2double F_onset = max(F_onset2 - F_onset3, 0.0); + + F_onset1_Here = F_onset1; + F_onset2_Here = F_onset2; + F_onset3_Here = F_onset3; + F_onset_Here = F_onset; + + if (options.CrossFlow) { + + // Computation of shape factor + const su2double k = 0.25 - lambda_theta; + const su2double FirstTerm = 4.14 * k; + const su2double SecondTerm = 83.5 * pow(k, 2.0); + const su2double ThirdTerm = 854.0 * pow(k, 3.0); + const su2double ForthTerm = 3337.0 * pow(k, 4.0); + const su2double FifthTerm = 4576.0 * pow(k, 5.0); + const su2double H = min(2.0 + FirstTerm - SecondTerm + ThirdTerm - ForthTerm + FifthTerm, 2.7); + + // Computation of critical cross flow Reynolds number + su2double Re_Crit_CF = 0.0; + if(H < 2.3) { + Re_Crit_CF = 150.0; + } else { + Re_Crit_CF = -(300.0/PI_NUMBER) * atan(0.106/(pow(H-2.3, 2.05))); + } + + // Helicity computation + su2double VelocityNormalized[3]; + VelocityNormalized[0] = vel_u / Velocity_Mag; + VelocityNormalized[1] = vel_v / Velocity_Mag; + if (nDim == 3) VelocityNormalized[2] = vel_w / Velocity_Mag; + + su2double StreamwiseVort = 0.0; + for (auto iDim = 0u; iDim < nDim; iDim++) { + StreamwiseVort += VelocityNormalized[iDim] * Vorticity_i[iDim]; + } + StreamwiseVort = abs(StreamwiseVort); + + const su2double H_CF = StreamwiseVort * dist_i / Velocity_Mag; + + // Computation of Delta_H_CF. Here I have included directly R_t as the ration between turb and lam viscosity + const su2double Delta_H_CF = H_CF * (1.0 + min(R_t, 0.4)); + + // Take into account for roughness + const su2double h_0 = 0.25e-6; + const su2double C_r = 2.0 - pow(0.5, config->GethRoughness()/h_0); + + // Construct Cross flow activation function + const su2double C_CF = 1.0; + const su2double f_CF = (C_CF * C_r * Delta_H_CF * Corr_Rec) / Re_Crit_CF; + const su2double F_onset_CF = min(max(0.0, f_CF - 1.0), 1.0); + + // Adjust onset function for intermittency + F_onset = max(F_onset, F_onset_CF); + + } + + const su2double f_turb = exp(-pow(R_t / 2, 4)); + + /*-- production term of Intermeittency(Gamma) --*/ + const su2double Pg = + F_length * Density_i * StrainMag_i * F_onset * TransVar_i[0] * (1.0 - TransVar_i[0]); + + /*-- destruction term of Intermeittency(Gamma) --*/ + const su2double Dg = c_a2 * Density_i * VorticityMag * TransVar_i[0] * f_turb * (c_e2 * TransVar_i[0] - 1.0); + + Prod_Here = Pg; + Destr_Here = Dg; + + /*--- Source ---*/ + Residual += (Pg - Dg) * Volume; + + /*--- Implicit part ---*/ + Jacobian_i[0] = (F_length * StrainMag_i * F_onset * (1 - 2*TransVar_i[0]) - + c_a2 * VorticityMag * f_turb * (2.0 * c_e2 * TransVar_i[0] - 1.0)) * + Volume; + + } + + AD::SetPreaccOut(Residual, nVar); + AD::EndPreacc(); + + return ResidualType<>(&Residual, &Jacobian_i, nullptr); + + } + + inline su2double GetRe_t() override {return Re_t;} + inline su2double GetCorr_Rec() override {return Corr_Rec;} + inline su2double GetTu() override {return Tu_Here;} + inline su2double GetLambda_theta() override {return lambda_theta_Here;} + inline su2double Getduds() override {return duds_Here;} + inline su2double GetRe_v() override {return Re_v_Here;} + inline su2double GetProd() override {return Prod_Here;} + inline su2double GetDestr() override {return Destr_Here;} + inline su2double GetF_onset1() override {return F_onset1_Here;} + inline su2double GetF_onset2() override {return F_onset2_Here;} + inline su2double GetF_onset3() override {return F_onset3_Here;} + inline su2double GetF_onset() override {return F_onset_Here;} + inline void SetAuxVar(su2double val_AuxVar) override { AuxVar = val_AuxVar;} + // non serve più + inline void SetF2(su2double val_F2) override { F2 = val_F2;} + +}; \ No newline at end of file diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp index ce75b38f115..b30eaab2ec4 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp @@ -910,6 +910,16 @@ class CSourcePieceWise_TurbSST final : public CNumerics { /*--- LM model coupling with production and dissipation term for k transport equation---*/ if (config->GetKind_Trans_Model() == TURB_TRANS_MODEL::LM) { pk = pk * eff_intermittency; + // Check if the Prod_lim_k has to be introduced based on input options + if ((config->GetLMParsedOptions()).SLM && (config->GetLMParsedOptions()).Correlation_SLM == TURB_TRANS_CORRELATION_SLM::MENTER_SLM) { + const su2double Re_theta_c_lim = 1100.0; + const su2double C_k = 1.0; + const su2double C_SEP = 1.0; + const su2double Re_v = Density_i * dist_i * dist_i * StrainMag_i / Laminar_Viscosity_i; + const su2double F_on_lim = min(max(Re_v/(2.2*Re_theta_c_lim)-1.0, 0.0), 3.0); + const su2double IntermittencyRelated = max(eff_intermittency-0.2, 0.0) * (1.0 - eff_intermittency); + pk = pk + 5*C_k * IntermittencyRelated * F_on_lim * max(3.0*C_SEP*Laminar_Viscosity_i - Eddy_Viscosity_i, 0.0) * StrainMag_i * VorticityMag; + } dk = min(max(eff_intermittency, 0.1), 1.0) * dk; } diff --git a/SU2_CFD/include/solvers/CTransLMSolver.hpp b/SU2_CFD/include/solvers/CTransLMSolver.hpp index 98595abb266..40515ed5054 100644 --- a/SU2_CFD/include/solvers/CTransLMSolver.hpp +++ b/SU2_CFD/include/solvers/CTransLMSolver.hpp @@ -42,6 +42,7 @@ class CTransLMSolver final : public CTurbSolver { LM_ParsedOptions options; TURB_FAMILY TurbFamily; + bool isSepNeeded; TransLMCorrelations TransCorrelations; diff --git a/SU2_CFD/include/variables/CTransLMVariable.hpp b/SU2_CFD/include/variables/CTransLMVariable.hpp index 96b556f5450..f1c3f280264 100644 --- a/SU2_CFD/include/variables/CTransLMVariable.hpp +++ b/SU2_CFD/include/variables/CTransLMVariable.hpp @@ -40,6 +40,23 @@ class CTransLMVariable final : public CTurbVariable { protected: VectorType Intermittency_Eff; VectorType Intermittency_Sep; + + VectorType Corr_Rec; + VectorType Re_t; + VectorType Tu; + VectorType Lambda_theta; + VectorType duds; + VectorType Re_v; + VectorType Prod; + VectorType Destr; + VectorType F_onset1; + VectorType F_onset2; + VectorType F_onset3; + VectorType F_onset; + + VectorType normal_x; + VectorType normal_y; + VectorType normal_z; public: /*! @@ -70,6 +87,28 @@ class CTransLMVariable final : public CTurbVariable { */ void SetIntermittencyEff(unsigned long iPoint, su2double val_Intermittency_sep) override; + /*! + * \brief Set Value of Transition Momentum Thickness Reynolds number from correlations. + */ + void SetCorr_Rec(unsigned long iPoint, su2double val_Corr_Rec) override; + + /*! + * \brief Set Value of Momentum Thickness Reynolds number from correlations (substitute to the second equation of original LM model). + */ + void SetRe_t(unsigned long iPoint, su2double val_Re_t) override; + void SetTu(unsigned long iPoint, su2double val_Tu) override; + void SetLambda_theta(unsigned long iPoint, su2double val_Lambda_theta) override; + void Setduds(unsigned long iPoint, su2double val_duds) override; + void SetRe_v(unsigned long iPoint, su2double val_Re_v) override; + void SetProd(unsigned long iPoint, su2double val_Prod) override; + void SetDestr(unsigned long iPoint, su2double val_Destr) override; + void SetF_onset1(unsigned long iPoint, su2double val_F_onset1) override; + void SetF_onset2(unsigned long iPoint, su2double val_F_onset2) override; + void SetF_onset3(unsigned long iPoint, su2double val_F_onset3) override; + void SetF_onset(unsigned long iPoint, su2double val_F_onset) override; + void SetNormal(unsigned long iPoint, su2double val_normal_x, su2double val_normal_y, su2double val_normal_z) override; + + /*! * \brief Calculate effective intermittency. */ @@ -80,4 +119,26 @@ class CTransLMVariable final : public CTurbVariable { */ inline su2double GetIntermittencySep(unsigned long iPoint) const override { return Intermittency_Sep(iPoint); } + /*! + * \brief Get Value of Transition Momentum Thickness Reynolds number from correlations. + */ + inline su2double GetCorr_Rec(unsigned long iPoint) const override { return Corr_Rec(iPoint); } + + /*! + * \brief Get Value of Momentum Thickness Reynolds number from correlations (substitute to the second equation of original LM model). + */ + inline su2double GetRe_t(unsigned long iPoint) const override { return Re_t(iPoint); } + inline su2double GetTu(unsigned long iPoint) const override { return Tu(iPoint); } + inline su2double GetLambda_theta(unsigned long iPoint) const override { return Lambda_theta(iPoint); } + inline su2double Getduds(unsigned long iPoint) const override { return duds(iPoint); } + inline su2double GetRe_v(unsigned long iPoint) const override { return Re_v(iPoint); } + inline su2double GetProd(unsigned long iPoint) const override { return Prod(iPoint); } + inline su2double GetDestr(unsigned long iPoint) const override { return Destr(iPoint); } + inline su2double GetF_onset1(unsigned long iPoint) const override { return F_onset1(iPoint); } + inline su2double GetF_onset2(unsigned long iPoint) const override { return F_onset2(iPoint); } + inline su2double GetF_onset3(unsigned long iPoint) const override { return F_onset3(iPoint); } + inline su2double GetF_onset(unsigned long iPoint) const override { return F_onset(iPoint); } + inline su2double GetNormal_x(unsigned long iPoint) const override {return normal_x(iPoint);}; + inline su2double GetNormal_y(unsigned long iPoint) const override {return normal_y(iPoint);}; + inline su2double GetNormal_z(unsigned long iPoint) const override {return normal_z(iPoint);}; }; diff --git a/SU2_CFD/include/variables/CTurbSSTVariable.hpp b/SU2_CFD/include/variables/CTurbSSTVariable.hpp index 8735ddc6de6..d745ca82d71 100644 --- a/SU2_CFD/include/variables/CTurbSSTVariable.hpp +++ b/SU2_CFD/include/variables/CTurbSSTVariable.hpp @@ -45,6 +45,7 @@ class CTurbSSTVariable final : public CTurbVariable { VectorType F2; /*!< \brief Menter blending function for blending of k-w and k-eps. */ VectorType CDkw; /*!< \brief Cross-diffusion. */ SST_ParsedOptions sstParsedOptions; + LM_ParsedOptions lmParsedOptions; public: /*! * \brief Constructor of the class. diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index 920f5680915..9c51ffafb39 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -1699,6 +1699,50 @@ class CVariable { */ inline virtual void SetIntermittencyEff(unsigned long iPoint, su2double val_Intermittency_eff) {} + /*! + * \brief Set Value of Transition Momentum Thickness Reynolds number from correlations. + */ + inline virtual void SetCorr_Rec(unsigned long iPoint, su2double val_Corr_Rec) {}; + inline virtual void SetTu(unsigned long iPoint, su2double val_Tu) {}; + inline virtual void SetLambda_theta(unsigned long iPoint, su2double val_Lambda_theta) {}; + inline virtual void Setduds(unsigned long iPoint, su2double val_duds) {}; + inline virtual void SetRe_v(unsigned long iPoint, su2double val_Re_v) {}; + inline virtual void SetProd(unsigned long iPoint, su2double val_Prod) {}; + inline virtual void SetDestr(unsigned long iPoint, su2double val_Destr) {}; + inline virtual void SetF_onset1(unsigned long iPoint, su2double val_F_onset1) {}; + inline virtual void SetF_onset2(unsigned long iPoint, su2double val_F_onset2) {}; + inline virtual void SetF_onset3(unsigned long iPoint, su2double val_F_onset3) {}; + inline virtual void SetF_onset(unsigned long iPoint, su2double val_F_onset1) {}; + inline virtual void SetNormal(unsigned long iPoint, su2double val_normal_x, su2double val_normal_y, su2double val_normal_z) {}; + + /*! + * \brief Get Value of Transition Momentum Thickness Reynolds number from correlations. + */ + inline virtual su2double GetCorr_Rec(unsigned long iPoint) const { return 0.0; } + inline virtual su2double GetTu(unsigned long iPoint) const { return 0.0; } + inline virtual su2double GetLambda_theta(unsigned long iPoint) const { return 0.0; } + inline virtual su2double Getduds(unsigned long iPoint) const { return 0.0; } + inline virtual su2double GetRe_v(unsigned long iPoint) const { return 0.0; } + inline virtual su2double GetProd(unsigned long iPoint) const { return 0.0; } + inline virtual su2double GetDestr(unsigned long iPoint) const { return 0.0; } + inline virtual su2double GetF_onset1(unsigned long iPoint) const { return 0.0; } + inline virtual su2double GetF_onset2(unsigned long iPoint) const { return 0.0; } + inline virtual su2double GetF_onset3(unsigned long iPoint) const { return 0.0; } + inline virtual su2double GetF_onset(unsigned long iPoint) const { return 0.0; } + inline virtual su2double GetNormal_x(unsigned long iPoint) const { return 0.0; } + inline virtual su2double GetNormal_y(unsigned long iPoint) const { return 0.0; } + inline virtual su2double GetNormal_z(unsigned long iPoint) const { return 0.0; } + + /*! + * \brief Set Value of Momentum Thickness Reynolds number from correlations (substitute to the second equation of original LM model). + */ + inline virtual void SetRe_t(unsigned long iPoint, su2double val_Re_t) {}; + + /*! + * \brief Get Value of Momentum Thickness Reynolds number from correlations (substitute to the second equation of original LM model). + */ + inline virtual su2double GetRe_t(unsigned long iPoint) const { return 0.0; } + /*! * \brief Set the value of the eddy viscosity. * \param[in] val_muT diff --git a/SU2_CFD/src/SU2_CFD.cpp b/SU2_CFD/src/SU2_CFD.cpp index 353a744c891..62c98ee28f7 100644 --- a/SU2_CFD/src/SU2_CFD.cpp +++ b/SU2_CFD/src/SU2_CFD.cpp @@ -33,7 +33,7 @@ #endif /* Include file, needed for the runtime NaN catching. You also have to include feenableexcept(...) below. */ -//#include +#include using namespace std; @@ -76,6 +76,7 @@ int main(int argc, char *argv[]) { /*--- Uncomment the following line if runtime NaN catching is desired. ---*/ // feenableexcept(FE_INVALID | FE_OVERFLOW | FE_DIVBYZERO ); + feenableexcept(FE_OVERFLOW | FE_DIVBYZERO ); /*--- Initialize libxsmm, if supported. ---*/ #ifdef HAVE_LIBXSMM diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index 65a25b6c9a2..6c4b4fb4b2d 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -1319,6 +1319,8 @@ void CDriver::InstantiateTransitionNumerics(unsigned short nVar_Trans, int offse const int visc_bound_term = VISC_BOUND_TERM + offset; const bool LM = config->GetKind_Trans_Model() == TURB_TRANS_MODEL::LM; + LM_ParsedOptions options; + if(LM) options = config->GetLMParsedOptions(); /*--- Definition of the convective scheme for each equation and mesh level ---*/ @@ -1328,7 +1330,10 @@ void CDriver::InstantiateTransitionNumerics(unsigned short nVar_Trans, int offse break; case SPACE_UPWIND : for (auto iMGlevel = 0u; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { - if (LM) numerics[iMGlevel][TRANS_SOL][conv_term] = new CUpwSca_TransLM(nDim, nVar_Trans, config); + if (LM){ + if (!options.SLM) numerics[iMGlevel][TRANS_SOL][conv_term] = new CUpwSca_TransLM(nDim, nVar_Trans, config); + if (options.SLM) numerics[iMGlevel][TRANS_SOL][conv_term] = new CUpwSca_TransSLM(nDim, nVar_Trans, config); + } } break; default: @@ -1339,7 +1344,10 @@ void CDriver::InstantiateTransitionNumerics(unsigned short nVar_Trans, int offse /*--- Definition of the viscous scheme for each equation and mesh level ---*/ for (auto iMGlevel = 0u; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { - if (LM) numerics[iMGlevel][TRANS_SOL][visc_term] = new CAvgGrad_TransLM(nDim, nVar_Trans, true, config); + if (LM){ + if (!options.SLM) numerics[iMGlevel][TRANS_SOL][visc_term] = new CAvgGrad_TransLM(nDim, nVar_Trans, true, config); + if (options.SLM) numerics[iMGlevel][TRANS_SOL][visc_term] = new CAvgGrad_TransSLM(nDim, nVar_Trans, true, config); + } } /*--- Definition of the source term integration scheme for each equation and mesh level ---*/ @@ -1347,7 +1355,10 @@ void CDriver::InstantiateTransitionNumerics(unsigned short nVar_Trans, int offse for (auto iMGlevel = 0u; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { auto& trans_source_first_term = numerics[iMGlevel][TRANS_SOL][source_first_term]; - if (LM) trans_source_first_term = new CSourcePieceWise_TransLM(nDim, nVar_Trans, config); + if (LM){ + if (!options.SLM) trans_source_first_term = new CSourcePieceWise_TransLM(nDim, nVar_Trans, config); + if (options.SLM) trans_source_first_term = new CSourcePieceWise_TransSLM(nDim, nVar_Trans, config); + } numerics[iMGlevel][TRANS_SOL][source_second_term] = new CSourceNothing(nDim, nVar_Trans, config); } @@ -1355,9 +1366,12 @@ void CDriver::InstantiateTransitionNumerics(unsigned short nVar_Trans, int offse /*--- Definition of the boundary condition method ---*/ for (auto iMGlevel = 0u; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { - if (LM) { + if (!options.SLM) { numerics[iMGlevel][TRANS_SOL][conv_bound_term] = new CUpwSca_TransLM(nDim, nVar_Trans, config); numerics[iMGlevel][TRANS_SOL][visc_bound_term] = new CAvgGrad_TransLM(nDim, nVar_Trans, false, config); + } else { + numerics[iMGlevel][TRANS_SOL][conv_bound_term] = new CUpwSca_TransSLM(nDim, nVar_Trans, config); + numerics[iMGlevel][TRANS_SOL][visc_bound_term] = new CAvgGrad_TransSLM(nDim, nVar_Trans, false, config); } } } diff --git a/SU2_CFD/src/output/CFlowIncOutput.cpp b/SU2_CFD/src/output/CFlowIncOutput.cpp index cd275b8f51e..c70551c14d4 100644 --- a/SU2_CFD/src/output/CFlowIncOutput.cpp +++ b/SU2_CFD/src/output/CFlowIncOutput.cpp @@ -211,7 +211,7 @@ void CFlowIncOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CSolv SetHistoryOutputValue("MAX_PRESSURE", log10(flow_solver->GetRes_Max(0))); SetHistoryOutputValue("MAX_VELOCITY-X", log10(flow_solver->GetRes_Max(1))); SetHistoryOutputValue("MAX_VELOCITY-Y", log10(flow_solver->GetRes_Max(2))); - if (nDim == 3) SetHistoryOutputValue("RMS_VELOCITY-Z", log10(flow_solver->GetRes_Max(3))); + if (nDim == 3) SetHistoryOutputValue("MAX_VELOCITY-Z", log10(flow_solver->GetRes_Max(3))); if (multiZone){ SetHistoryOutputValue("BGS_PRESSURE", log10(flow_solver->GetRes_BGS(0))); diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp index 5ffee0d23f5..baafc717594 100644 --- a/SU2_CFD/src/output/CFlowOutput.cpp +++ b/SU2_CFD/src/output/CFlowOutput.cpp @@ -981,7 +981,10 @@ void CFlowOutput::AddHistoryOutputFields_ScalarRMS_RES(const CConfig* config) { /// DESCRIPTION: Root-mean square residual of the intermittency (LM model). AddHistoryOutput("RMS_INTERMITTENCY", "rms[LM_1]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of intermittency (LM model).", HistoryFieldType::RESIDUAL); /// DESCRIPTION: Root-mean square residual of the momentum thickness Reynolds number (LM model). - AddHistoryOutput("RMS_RE_THETA_T", "rms[LM_2]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of momentum thickness Reynolds number (LM model).", HistoryFieldType::RESIDUAL); + if (!(config->GetLMParsedOptions()).SLM) { + /// DESCRIPTION: Root-mean square residual of the momentum thickness Reynolds number (LM model). + AddHistoryOutput("RMS_RE_THETA_T", "rms[LM_2]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of momentum thickness Reynolds number (LM model).", HistoryFieldType::RESIDUAL); + } break; case TURB_TRANS_MODEL::NONE: break; @@ -1037,7 +1040,10 @@ void CFlowOutput::AddHistoryOutputFields_ScalarMAX_RES(const CConfig* config) { /// DESCRIPTION: Maximum residual of the intermittency (LM model). AddHistoryOutput("MAX_INTERMITTENCY", "max[LM_1]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of the intermittency (LM model).", HistoryFieldType::RESIDUAL); /// DESCRIPTION: Maximum residual of the momentum thickness Reynolds number (LM model). - AddHistoryOutput("MAX_RE_THETA_T", "max[LM_2]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of the momentum thickness Reynolds number (LM model).", HistoryFieldType::RESIDUAL); + if (!(config->GetLMParsedOptions()).SLM) { + /// DESCRIPTION: Maximum residual of the momentum thickness Reynolds number (LM model). + AddHistoryOutput("MAX_RE_THETA_T", "max[LM_2]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of the momentum thickness Reynolds number (LM model).", HistoryFieldType::RESIDUAL); + } break; case TURB_TRANS_MODEL::NONE: @@ -1093,7 +1099,10 @@ void CFlowOutput::AddHistoryOutputFields_ScalarBGS_RES(const CConfig* config) { /// DESCRIPTION: Maximum residual of the intermittency (LM model). AddHistoryOutput("BGS_INTERMITTENCY", "bgs[LM_1]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the intermittency (LM model).", HistoryFieldType::RESIDUAL); /// DESCRIPTION: Maximum residual of the momentum thickness Reynolds number (LM model). - AddHistoryOutput("BGS_RE_THETA_T", "bgs[LM_2]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the momentum thickness Reynolds number (LM model).", HistoryFieldType::RESIDUAL); + if (!(config->GetLMParsedOptions()).SLM) { + /// DESCRIPTION: Maximum residual of the momentum thickness Reynolds number (LM model). + AddHistoryOutput("BGS_RE_THETA_T", "bgs[LM_2]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the momentum thickness Reynolds number (LM model).", HistoryFieldType::RESIDUAL); + } break; case TURB_TRANS_MODEL::NONE: break; @@ -1184,12 +1193,16 @@ void CFlowOutput::LoadHistoryDataScalar(const CConfig* config, const CSolver* co switch (config->GetKind_Trans_Model()) { case TURB_TRANS_MODEL::LM: SetHistoryOutputValue("RMS_INTERMITTENCY", log10(solver[TRANS_SOL]->GetRes_RMS(0))); - SetHistoryOutputValue("RMS_RE_THETA_T",log10(solver[TRANS_SOL]->GetRes_RMS(1))); SetHistoryOutputValue("MAX_INTERMITTENCY", log10(solver[TRANS_SOL]->GetRes_Max(0))); - SetHistoryOutputValue("MAX_RE_THETA_T", log10(solver[TRANS_SOL]->GetRes_Max(1))); + if (!(config->GetLMParsedOptions()).SLM) { + SetHistoryOutputValue("RMS_RE_THETA_T",log10(solver[TRANS_SOL]->GetRes_RMS(1))); + SetHistoryOutputValue("MAX_RE_THETA_T", log10(solver[TRANS_SOL]->GetRes_Max(1))); + } if (multiZone) { SetHistoryOutputValue("BGS_INTERMITTENCY", log10(solver[TRANS_SOL]->GetRes_BGS(0))); - SetHistoryOutputValue("BGS_RE_THETA_T", log10(solver[TRANS_SOL]->GetRes_BGS(1))); + if (!(config->GetLMParsedOptions()).SLM) { + SetHistoryOutputValue("BGS_RE_THETA_T", log10(solver[TRANS_SOL]->GetRes_BGS(1))); + } } SetHistoryOutputValue("LINSOL_ITER_TRANS", solver[TRANS_SOL]->GetIterLinSolver()); SetHistoryOutputValue("LINSOL_RESIDUAL_TRANS", log10(solver[TRANS_SOL]->GetResLinSolver())); @@ -1262,6 +1275,33 @@ void CFlowOutput::SetVolumeOutputFieldsScalarSolution(const CConfig* config){ case TURB_TRANS_MODEL::LM: AddVolumeOutput("INTERMITTENCY", "LM_gamma", "SOLUTION", "LM intermittency"); AddVolumeOutput("RE_THETA_T", "LM_Re_t", "SOLUTION", "LM RE_THETA_T"); + AddVolumeOutput("RE_V", "Re_v", "DEBUG", "LM Re_v"); + AddVolumeOutput("RE_THETA_CORR", "LM_Corr_Rec", "DEBUG", "LM RE_THETA_CORR"); + AddVolumeOutput("PROD", "LM_Prod", "DEBUG", "LM PROD"); + AddVolumeOutput("DESTR", "LM_Destr", "DEBUG", "LM DESTR"); + AddVolumeOutput("F_ONSET1", "LM_F_onset1", "DEBUG", "LM F_ONSET1"); + AddVolumeOutput("F_ONSET2", "LM_F_onset2", "DEBUG", "LM F_ONSET2"); + AddVolumeOutput("F_ONSET3", "LM_F_onset3", "DEBUG", "LM F_ONSET3"); + AddVolumeOutput("F_ONSET", "LM_F_onset", "DEBUG", "LM F_ONSET"); + AddVolumeOutput("STRAINMAG", "StrainMag", "DEBUG", "LM STRAINMAG"); + AddVolumeOutput("LAMBDA_THETA", "Lambda_theta", "DEBUG", "LM Lambda_theta"); + AddVolumeOutput("DU_DS", "du_ds", "DEBUG", "LM du_ds"); + if ((config->GetLMParsedOptions()).SLM) { + AddVolumeOutput("TU", "Tu", "SOLUTION", "LM Tu"); + AddVolumeOutput("NORMAL_X", "Normal_x", "DEBUG", "LM Normal_x"); + AddVolumeOutput("NORMAL_Y", "Normal_y", "DEBUG", "LM Normal_y"); + AddVolumeOutput("NORMAL_Z", "Normal_z", "DEBUG", "LM Normal_z"); + if (!((config->GetLMParsedOptions()).Correlation_SLM == TURB_TRANS_CORRELATION_SLM::MENTER_SLM)) { + AddVolumeOutput("INTERMITTENCY_SEP", "LM_gamma_sep", "PRIMITIVE", "LM intermittency"); + AddVolumeOutput("INTERMITTENCY_EFF", "LM_gamma_eff", "PRIMITIVE", "LM RE_THETA_T"); + } + } + + if (!(config->GetLMParsedOptions()).SLM) { + AddVolumeOutput("INTERMITTENCY_SEP", "LM_gamma_sep", "PRIMITIVE", "LM intermittency"); + AddVolumeOutput("INTERMITTENCY_EFF", "LM_gamma_eff", "PRIMITIVE", "LM RE_THETA_T"); + } + AddVolumeOutput("TURB_INDEX", "Turb_index", "PRIMITIVE", "Turbulence index"); break; case TURB_TRANS_MODEL::NONE: @@ -1334,7 +1374,9 @@ void CFlowOutput::SetVolumeOutputFieldsScalarResidual(const CConfig* config) { switch (config->GetKind_Trans_Model()) { case TURB_TRANS_MODEL::LM: AddVolumeOutput("RES_INTERMITTENCY", "Residual_LM_intermittency", "RESIDUAL", "Residual of LM intermittency"); - AddVolumeOutput("RES_RE_THETA_T", "Residual_LM_RE_THETA_T", "RESIDUAL", "Residual of LM RE_THETA_T"); + if (!(config->GetLMParsedOptions()).SLM) { + AddVolumeOutput("RES_RE_THETA_T", "Residual_LM_RE_THETA_T", "RESIDUAL", "Residual of LM RE_THETA_T"); + } break; case TURB_TRANS_MODEL::NONE: @@ -1553,12 +1595,39 @@ void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* con switch (config->GetKind_Trans_Model()) { case TURB_TRANS_MODEL::LM: SetVolumeOutputValue("INTERMITTENCY", iPoint, Node_Trans->GetSolution(iPoint, 0)); - SetVolumeOutputValue("RE_THETA_T", iPoint, Node_Trans->GetSolution(iPoint, 1)); - SetVolumeOutputValue("INTERMITTENCY_SEP", iPoint, Node_Trans->GetIntermittencySep(iPoint)); - SetVolumeOutputValue("INTERMITTENCY_EFF", iPoint, Node_Trans->GetIntermittencyEff(iPoint)); + SetVolumeOutputValue("RE_V", iPoint, Node_Trans->GetRe_v(iPoint)); + SetVolumeOutputValue("RE_THETA_CORR", iPoint, Node_Trans->GetCorr_Rec(iPoint)); + SetVolumeOutputValue("PROD", iPoint, Node_Trans->GetProd(iPoint)); + SetVolumeOutputValue("DESTR", iPoint, Node_Trans->GetDestr(iPoint)); + SetVolumeOutputValue("F_ONSET1", iPoint, Node_Trans->GetF_onset1(iPoint)); + SetVolumeOutputValue("F_ONSET2", iPoint, Node_Trans->GetF_onset2(iPoint)); + SetVolumeOutputValue("F_ONSET3", iPoint, Node_Trans->GetF_onset3(iPoint)); + SetVolumeOutputValue("F_ONSET", iPoint, Node_Trans->GetF_onset(iPoint)); + SetVolumeOutputValue("STRAINMAG", iPoint, Node_Flow->GetStrainMag(iPoint)); + SetVolumeOutputValue("LAMBDA_THETA", iPoint, Node_Trans->GetLambda_theta(iPoint)); + SetVolumeOutputValue("DU_DS", iPoint, Node_Trans->Getduds(iPoint)); + if (!(config->GetLMParsedOptions()).SLM) { + SetVolumeOutputValue("RE_THETA_T", iPoint, Node_Trans->GetSolution(iPoint, 1)); + } else { + SetVolumeOutputValue("RE_THETA_T", iPoint, Node_Trans->GetRe_t(iPoint)); + SetVolumeOutputValue("TU", iPoint, Node_Trans->GetTu(iPoint)); + SetVolumeOutputValue("NORMAL_X", iPoint, Node_Trans->GetNormal_x(iPoint)); + SetVolumeOutputValue("NORMAL_Y", iPoint, Node_Trans->GetNormal_y(iPoint)); + SetVolumeOutputValue("NORMAL_Z", iPoint, Node_Trans->GetNormal_z(iPoint)); + if (!((config->GetLMParsedOptions()).Correlation_SLM == TURB_TRANS_CORRELATION_SLM::MENTER_SLM)) { + SetVolumeOutputValue("INTERMITTENCY_SEP", iPoint, Node_Trans->GetIntermittencySep(iPoint)); + SetVolumeOutputValue("INTERMITTENCY_EFF", iPoint, Node_Trans->GetIntermittencyEff(iPoint)); + } + } + if (!(config->GetLMParsedOptions()).SLM) { + SetVolumeOutputValue("INTERMITTENCY_SEP", iPoint, Node_Trans->GetIntermittencySep(iPoint)); + SetVolumeOutputValue("INTERMITTENCY_EFF", iPoint, Node_Trans->GetIntermittencyEff(iPoint)); + } SetVolumeOutputValue("TURB_INDEX", iPoint, Node_Turb->GetTurbIndex(iPoint)); SetVolumeOutputValue("RES_INTERMITTENCY", iPoint, trans_solver->LinSysRes(iPoint, 0)); - SetVolumeOutputValue("RES_RE_THETA_T", iPoint, trans_solver->LinSysRes(iPoint, 1)); + if (!(config->GetLMParsedOptions()).SLM) { + SetVolumeOutputValue("RES_RE_THETA_T", iPoint, trans_solver->LinSysRes(iPoint, 1)); + } break; case TURB_TRANS_MODEL::NONE: break; @@ -2667,7 +2736,7 @@ void CFlowOutput::WriteForcesBreakdown(const CConfig* config, const CSolver* flo case TURB_TRANS_MODEL::NONE: break; case TURB_TRANS_MODEL::LM: file << "Langtry and Menter's transition"; - if (config->GetLMParsedOptions().LM2015) { + if (config->GetLMParsedOptions().CrossFlow) { file << " w/ cross-flow corrections (2015)\n"; } else { file << " (2009)\n"; diff --git a/SU2_CFD/src/solvers/CTransLMSolver.cpp b/SU2_CFD/src/solvers/CTransLMSolver.cpp index 902a67925c4..3d774401642 100644 --- a/SU2_CFD/src/solvers/CTransLMSolver.cpp +++ b/SU2_CFD/src/solvers/CTransLMSolver.cpp @@ -49,6 +49,14 @@ CTransLMSolver::CTransLMSolver(CGeometry *geometry, CConfig *config, unsigned sh /*--- Dimension of the problem --> 2 Transport equations (intermittency, Reth) ---*/ nVar = 2; nPrimVar = 2; + + /*--- Check if Simplified version is used ---*/ + options = config->GetLMParsedOptions(); + if (options.SLM) { + nVar = 1; + nPrimVar = 1; + } + nPoint = geometry->GetnPoint(); nPointDomain = geometry->GetnPointDomain(); @@ -60,11 +68,15 @@ CTransLMSolver::CTransLMSolver(CGeometry *geometry, CConfig *config, unsigned sh nDim = geometry->GetnDim(); - /*--- Define variables needed for transition from config file */ - options = config->GetLMParsedOptions(); + /*--- Define variables needed for transition from config file ---*/ TransCorrelations.SetOptions(options); TurbFamily = TurbModelFamily(config->GetKind_Turb_Model()); + isSepNeeded = true; + // If the Simplified model is used coupled with SST then we do not need the separation induced intermittency + // due to the added production term to k equation + if (options.SLM && options.Correlation_SLM == TURB_TRANS_CORRELATION_SLM::MENTER_SLM) isSepNeeded = false; + /*--- Single grid simulation ---*/ if (iMesh == MESH_0) { @@ -101,13 +113,17 @@ CTransLMSolver::CTransLMSolver(CGeometry *geometry, CConfig *config, unsigned sh lowerlimit[0] = 1.0e-4; upperlimit[0] = 5.0; - lowerlimit[1] = 1.0e-4; - upperlimit[1] = 1.0e15; + if (!options.SLM) { + lowerlimit[1] = 1.0e-4; + upperlimit[1] = 1.0e15; + } /*--- Far-field flow state quantities and initialization. ---*/ const su2double Intensity = config->GetTurbulenceIntensity_FreeStream()*100.0; const su2double Intermittency_Inf = 1.0; + Solution_Inf[0] = Intermittency_Inf; + su2double ReThetaT_Inf = 100.0; /*--- Momentum thickness Reynolds number, initialized from freestream turbulent intensity*/ @@ -124,10 +140,13 @@ CTransLMSolver::CTransLMSolver(CGeometry *geometry, CConfig *config, unsigned sh } Solution_Inf[0] = Intermittency_Inf; - Solution_Inf[1] = ReThetaT_Inf; + if (!options.SLM) { + Solution_Inf[1] = ReThetaT_Inf; + } /*--- Initialize the solution to the far-field state everywhere. ---*/ nodes = new CTransLMVariable(Intermittency_Inf, ReThetaT_Inf, 1.0, 1.0, nPoint, nDim, nVar, config); + SetBaseClassPointerToNodes(); /*--- MPI solution ---*/ @@ -155,7 +174,7 @@ CTransLMSolver::CTransLMSolver(CGeometry *geometry, CConfig *config, unsigned sh Inlet_TurbVars[iMarker].resize(nVertex[iMarker],nVar); for (unsigned long iVertex = 0; iVertex < nVertex[iMarker]; ++iVertex) { Inlet_TurbVars[iMarker](iVertex,0) = Intermittency_Inf; - Inlet_TurbVars[iMarker](iVertex,1) = ReThetaT_Inf; + if (!options.SLM) Inlet_TurbVars[iMarker](iVertex,1) = ReThetaT_Inf; } } @@ -178,6 +197,38 @@ void CTransLMSolver::Preprocessing(CGeometry *geometry, CSolver **solver_contain /*--- Upwind second order reconstruction and gradients ---*/ CommonPreprocessing(geometry, config, Output); + + if (options.SLM && options.Correlation_SLM == TURB_TRANS_CORRELATION_SLM::MENTER_SLM) { + + auto* flowNodes = su2staticcast_p(solver_container[FLOW_SOL]->GetNodes()); + SU2_OMP_FOR_STAT(omp_chunk_size) + for (unsigned long iPoint = 0; iPoint < nPoint; iPoint ++) { + auto Normal = geometry->nodes->GetNormal(iPoint); + nodes->SetAuxVar(iPoint, 0, flowNodes->GetProjVel(iPoint, Normal)); + nodes->SetNormal(iPoint, Normal[0], Normal[1], Normal[2]); + } + END_SU2_OMP_FOR + + if (config->GetKind_Gradient_Method() == GREEN_GAUSS) { + SetAuxVar_Gradient_GG(geometry, config); + } + if (config->GetKind_Gradient_Method() == WEIGHTED_LEAST_SQUARES) { + SetAuxVar_Gradient_LS(geometry, config); + } + + + SU2_OMP_FOR_STAT(omp_chunk_size) + for (unsigned long iPoint = 0; iPoint < nPoint; iPoint ++) { + su2double AuxVarHere = 0.0; + auto Normal = geometry->nodes->GetNormal(iPoint); + for (auto iDim = 0u; iDim < nDim; iDim++) + AuxVarHere += Normal[iDim] * nodes->GetAuxVarGradient(iPoint, 0, iDim); + nodes->SetAuxVar(iPoint, 0, AuxVarHere); + } + END_SU2_OMP_FOR + + } + } void CTransLMSolver::Postprocessing(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iMesh) { @@ -191,73 +242,103 @@ void CTransLMSolver::Postprocessing(CGeometry *geometry, CSolver **solver_contai SetSolution_Gradient_LS(geometry, config, -1); } + AD::StartNoSharedReading(); - auto* flowNodes = su2staticcast_p(solver_container[FLOW_SOL]->GetNodes()); - auto* turbNodes = su2staticcast_p(solver_container[TURB_SOL]->GetNodes()); - - SU2_OMP_FOR_STAT(omp_chunk_size) - for (unsigned long iPoint = 0; iPoint < nPoint; iPoint ++) { - - // Here the nodes already have the new solution, thus I have to compute everything from scratch - - const su2double rho = flowNodes->GetDensity(iPoint); - const su2double mu = flowNodes->GetLaminarViscosity(iPoint); - const su2double muT = turbNodes->GetmuT(iPoint); - const su2double dist = geometry->nodes->GetWall_Distance(iPoint); - su2double VorticityMag = GeometryToolbox::Norm(3, flowNodes->GetVorticity(iPoint)); - su2double StrainMag =flowNodes->GetStrainMag(iPoint); - VorticityMag = max(VorticityMag, 1e-12); - StrainMag = max(StrainMag, 1e-12); // safety against division by zero - const su2double Intermittency = nodes->GetSolution(iPoint,0); - const su2double Re_t = nodes->GetSolution(iPoint,1); - const su2double Re_v = rho*dist*dist*StrainMag/mu; - const su2double vel_u = flowNodes->GetVelocity(iPoint, 0); - const su2double vel_v = flowNodes->GetVelocity(iPoint, 1); - const su2double vel_w = (nDim ==3) ? flowNodes->GetVelocity(iPoint, 2) : 0.0; - const su2double VelocityMag = sqrt(vel_u*vel_u + vel_v*vel_v + vel_w*vel_w); - su2double omega = 0.0; - su2double k = 0.0; - if(TurbFamily == TURB_FAMILY::KW){ - omega = turbNodes->GetSolution(iPoint,1); - k = turbNodes->GetSolution(iPoint,0); - } - su2double Tu = 1.0; - if(TurbFamily == TURB_FAMILY::KW) - Tu = max(100.0*sqrt( 2.0 * k / 3.0 ) / VelocityMag,0.027); - if(TurbFamily == TURB_FAMILY::SA) - Tu = config->GetTurbulenceIntensity_FreeStream()*100; - - const su2double Corr_Rec = TransCorrelations.ReThetaC_Correlations(Tu, Re_t); - - su2double R_t = 1.0; - if(TurbFamily == TURB_FAMILY::KW) - R_t = rho*k/ mu/ omega; - if(TurbFamily == TURB_FAMILY::SA) - R_t = muT/ mu; - - const su2double f_reattach = exp(-pow(R_t/20,4)); - su2double f_wake = 0.0; - if(TurbFamily == TURB_FAMILY::KW){ - const su2double re_omega = rho*omega*dist*dist/mu; - f_wake = exp(-pow(re_omega/(1.0e+05),2)); + + if(isSepNeeded){ + + AD::StartNoSharedReading(); + auto* flowNodes = su2staticcast_p(solver_container[FLOW_SOL]->GetNodes()); + auto* turbNodes = su2staticcast_p(solver_container[TURB_SOL]->GetNodes()); + + SU2_OMP_FOR_STAT(omp_chunk_size) + for (unsigned long iPoint = 0; iPoint < nPoint; iPoint ++) { + + // Here the nodes already have the new solution, thus I have to compute everything from scratch + + const su2double rho = flowNodes->GetDensity(iPoint); + const su2double mu = flowNodes->GetLaminarViscosity(iPoint); + const su2double muT = turbNodes->GetmuT(iPoint); + const su2double dist = geometry->nodes->GetWall_Distance(iPoint); + su2double VorticityMag = GeometryToolbox::Norm(3, flowNodes->GetVorticity(iPoint)); + su2double StrainMag =flowNodes->GetStrainMag(iPoint); + VorticityMag = max(VorticityMag, 1e-12); + StrainMag = max(StrainMag, 1e-12); // safety against division by zero + const su2double Intermittency = nodes->GetSolution(iPoint,0); + const su2double Re_v = rho*dist*dist*StrainMag/mu; + const su2double VelocityMag = max(sqrt(flowNodes->GetVelocity2(iPoint)), 1e-20); + su2double omega = 0.0; + su2double k = 0.0; + if(TurbFamily == TURB_FAMILY::KW){ + omega = turbNodes->GetSolution(iPoint,1); + k = turbNodes->GetSolution(iPoint,0); + } + + su2double Re_t = 0.0; + su2double Corr_Rec = 0.0; + + if (options.SLM) { + Re_t = nodes->GetRe_t(iPoint); + Corr_Rec = nodes->GetCorr_Rec(iPoint); + } else { + Re_t = nodes->GetSolution(iPoint,1); + su2double Tu = 1.0; + if(TurbFamily == TURB_FAMILY::KW) + Tu = max(100.0*sqrt( 2.0 * k / 3.0 ) / VelocityMag,0.027); + if(TurbFamily == TURB_FAMILY::SA) + Tu = config->GetTurbulenceIntensity_FreeStream()*100; + + Corr_Rec = TransCorrelations.ReThetaC_Correlations(Tu, Re_t); + } + + + // cout << Re_t << " " << Corr_Rec << endl; + // if (geometry->nodes->GetDomain(iPoint)) cout << "Point is on boundary" << endl; + + su2double R_t = 1.0; + if(TurbFamily == TURB_FAMILY::KW) + R_t = rho*k/ mu/ omega; + if(TurbFamily == TURB_FAMILY::SA) + R_t = muT/ mu; + + const su2double f_reattach = exp(-pow(R_t/20,4)); + + su2double f_wake = 0.0; + if(TurbFamily == TURB_FAMILY::KW){ + const su2double re_omega = rho*omega*dist*dist/mu; + f_wake = exp(-pow(re_omega/(1.0e+05),2)); + } + if(TurbFamily == TURB_FAMILY::SA) + f_wake = 1.0; + + //cout << "StrainMag = " << StrainMag << " rho = " << rho << " dist = " << dist << " Re_v = " << Re_v << " Corr_Rec = " << Corr_Rec << endl; + + const su2double theta_bl = Re_t*mu / rho /VelocityMag; + const su2double delta_bl = 7.5*theta_bl; + const su2double delta = 50.0*VorticityMag*dist/VelocityMag*delta_bl + 1e-20; + const su2double var1 = (Intermittency-1.0/50.0)/(1.0-1.0/50.0); + const su2double var2 = 1.0 - pow(var1,2.0); + const su2double f_theta = min(max(f_wake*exp(-pow(dist/delta, 4)), var2), 1.0); + su2double Intermittency_Sep = 2.0*max(0.0, Re_v/(3.235*Corr_Rec)-1.0)*f_reattach; + //if (Intermittency_Sep>1.0) cout << "StrainMag = " << StrainMag << " rho = " << rho << " dist = " << dist << " Re_v = " << Re_v << " Corr_Rec = " << Corr_Rec << " Intermittency: " << Intermittency_Sep << " f_reattach = " << f_reattach << endl; + Intermittency_Sep = min(Intermittency_Sep,2.0)*f_theta; + Intermittency_Sep = min(max(0.0, Intermittency_Sep), 2.0); + nodes -> SetIntermittencySep(iPoint, Intermittency_Sep); + nodes -> SetIntermittencyEff(iPoint, Intermittency_Sep); + } - if(TurbFamily == TURB_FAMILY::SA) - f_wake = 1.0; - - const su2double theta_bl = Re_t*mu / rho /VelocityMag; - const su2double delta_bl = 7.5*theta_bl; - const su2double delta = 50.0*VorticityMag*dist/VelocityMag*delta_bl + 1e-20; - const su2double var1 = (Intermittency-1.0/50.0)/(1.0-1.0/50.0); - const su2double var2 = 1.0 - pow(var1,2.0); - const su2double f_theta = min(max(f_wake*exp(-pow(dist/delta, 4)), var2), 1.0); - su2double Intermittency_Sep = 2.0*max(0.0, Re_v/(3.235*Corr_Rec)-1.0)*f_reattach; - Intermittency_Sep = min(Intermittency_Sep,2.0)*f_theta; - Intermittency_Sep = min(max(0.0, Intermittency_Sep), 2.0); - nodes -> SetIntermittencySep(iPoint, Intermittency_Sep); - nodes -> SetIntermittencyEff(iPoint, Intermittency_Sep); + END_SU2_OMP_FOR } - END_SU2_OMP_FOR + else { + + // Effective intermittency and separation intermittency are not taken into account here! There is the added production term in the k-equation + SU2_OMP_FOR_STAT(omp_chunk_size) + for (unsigned long iPoint = 0; iPoint < nPoint; iPoint ++) { + nodes -> SetIntermittencyEff(iPoint, nodes->GetSolution(iPoint,0)); + } + END_SU2_OMP_FOR + } AD::EndNoSharedReading(); } @@ -333,15 +414,36 @@ void CTransLMSolver::Source_Residual(CGeometry *geometry, CSolver **solver_conta /*--- Set coordinate (for debugging) ---*/ numerics->SetCoord(geometry->nodes->GetCoord(iPoint), nullptr); - if (options.LM2015) { - /*--- Set local grid length (for LM2015)*/ + if (options.CrossFlow && !options.SLM) { + /*--- Set local grid length (for LM2015 cross-flow corrections)*/ numerics->SetLocalGridLength(geometry->nodes->GetMaxLength(iPoint)); } + if(options.SLM) { + if (options.Correlation_SLM == TURB_TRANS_CORRELATION_SLM::MENTER_SLM) numerics->SetAuxVar(nodes->GetAuxVar(iPoint, 0)); + numerics->SetF2(turbNodes->GetF2blending(iPoint)); + } + /*--- Compute the source term ---*/ auto residual = numerics->ComputeResidual(config); + if(options.SLM) { + nodes->SetRe_t(iPoint, numerics->GetRe_t()); + nodes->SetTu(iPoint, numerics->GetTu()); + } + + nodes->SetRe_v(iPoint, numerics->GetRe_v()); + nodes->SetCorr_Rec(iPoint, numerics->GetCorr_Rec()); + nodes->SetProd(iPoint, numerics->GetProd()); + nodes->SetDestr(iPoint, numerics->GetDestr()); + nodes->SetF_onset1(iPoint, numerics->GetF_onset1()); + nodes->SetF_onset2(iPoint, numerics->GetF_onset2()); + nodes->SetF_onset3(iPoint, numerics->GetF_onset3()); + nodes->SetF_onset(iPoint, numerics->GetF_onset()); + nodes->SetLambda_theta(iPoint, numerics->GetLambda_theta()); + nodes->Setduds(iPoint, numerics->Getduds()); + /*--- Subtract residual and the Jacobian ---*/ LinSysRes.SubtractBlock(iPoint, residual); @@ -540,6 +642,7 @@ void CTransLMSolver::LoadRestart(CGeometry** geometry, CSolver*** solver, CConfi const auto index = counter * Restart_Vars[1] + skipVars; for (auto iVar = 0u; iVar < nVar; iVar++) nodes->SetSolution(iPoint_Local, iVar, Restart_Data[index + iVar]); + // I am really not sure about these, I have to check if they are actually saved nodes ->SetIntermittencySep(iPoint_Local, Restart_Data[index + 2]); nodes ->SetIntermittencyEff(iPoint_Local, Restart_Data[index + 3]); diff --git a/SU2_CFD/src/solvers/CTurbSASolver.cpp b/SU2_CFD/src/solvers/CTurbSASolver.cpp index 40e19bd024c..a85fd084ef5 100644 --- a/SU2_CFD/src/solvers/CTurbSASolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSASolver.cpp @@ -269,7 +269,7 @@ void CTurbSASolver::Postprocessing(CGeometry *geometry, CSolver **solver_contain } shearStress = sqrt(shearStress); - FrictionVelocity = sqrt(shearStress/flowNodes->GetDensity(iPoint)); + FrictionVelocity = sqrt(shearStress/max(flowNodes->GetDensity(iPoint), 1e-20)); } else { su2double VorticityMag = max(GeometryToolbox::Norm(3, flowNodes->GetVorticity(iPoint)), 1e-12); FrictionVelocity = sqrt(flowNodes->GetLaminarViscosity(iPoint)*VorticityMag); @@ -277,7 +277,7 @@ void CTurbSASolver::Postprocessing(CGeometry *geometry, CSolver **solver_contain const su2double wall_dist = geometry->nodes->GetWall_Distance(jPoint); const su2double Derivative = nodes->GetSolution(jPoint, 0) / wall_dist; - const su2double turbulence_index = Derivative / (FrictionVelocity * 0.41); + const su2double turbulence_index = Derivative / max((FrictionVelocity * 0.41), 1e-20); nodes->SetTurbIndex(iPoint, turbulence_index); diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp index 63c29102175..eeb1b2f75e3 100644 --- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp @@ -268,10 +268,10 @@ void CTurbSSTSolver::Postprocessing(CGeometry *geometry, CSolver **solver_contai } shearStress = sqrt(shearStress); - const su2double FrictionVelocity = sqrt(shearStress/flowNodes->GetDensity(iPoint)); + const su2double FrictionVelocity = sqrt(shearStress/max(flowNodes->GetDensity(iPoint),1e-20)); const su2double wall_dist = geometry->nodes->GetWall_Distance(jPoint); const su2double Derivative = flowNodes->GetLaminarViscosity(jPoint) * pow(nodes->GetSolution(jPoint, 0), 0.673) / wall_dist; - const su2double turbulence_index = 6.1 * Derivative / pow(FrictionVelocity, 2.346); + const su2double turbulence_index = 6.1 * Derivative / max(pow(FrictionVelocity, 2.346),1e-20); nodes->SetTurbIndex(iPoint, turbulence_index); diff --git a/SU2_CFD/src/variables/CIncNSVariable.cpp b/SU2_CFD/src/variables/CIncNSVariable.cpp index 7aee6086d05..53300da5521 100644 --- a/SU2_CFD/src/variables/CIncNSVariable.cpp +++ b/SU2_CFD/src/variables/CIncNSVariable.cpp @@ -106,7 +106,6 @@ bool CIncNSVariable::SetPrimVar(unsigned long iPoint, su2double eddy_visc, su2do SetVelocity(iPoint); /*--- Set laminar viscosity ---*/ - SetLaminarViscosity(iPoint, FluidModel->GetLaminarViscosity()); /*--- Set eddy viscosity locally and in the fluid model. ---*/ diff --git a/SU2_CFD/src/variables/CTransLMVariable.cpp b/SU2_CFD/src/variables/CTransLMVariable.cpp index 551968b1c35..dca6f89113b 100644 --- a/SU2_CFD/src/variables/CTransLMVariable.cpp +++ b/SU2_CFD/src/variables/CTransLMVariable.cpp @@ -31,10 +31,19 @@ CTransLMVariable::CTransLMVariable(su2double Intermittency, su2double ReThetaT, su2double gammaSep, su2double gammaEff, unsigned long npoint, unsigned long ndim, unsigned long nvar, CConfig *config) : CTurbVariable(npoint, ndim, nvar, config) { - for(unsigned long iPoint=0; iPointGetLMParsedOptions(); + + if (!options.SLM) { + for(unsigned long iPoint=0; iPointGetSSTParsedOptions(); + lmParsedOptions = config->GetLMParsedOptions(); for(unsigned long iPoint=0; iPoint