diff --git a/.github/workflows/regression.yml b/.github/workflows/regression.yml index e1985a4911f2..fa6381b48e6d 100644 --- a/.github/workflows/regression.yml +++ b/.github/workflows/regression.yml @@ -233,7 +233,7 @@ jobs: uses: docker://ghcr.io/su2code/su2/test-su2:260405-0054 with: # -t -c - args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} + args: -b ${{github.ref}} -t develop -c feature_Trans_SLM -s ${{matrix.testscript}} - name: Cleanup uses: docker://ghcr.io/su2code/su2/test-su2:260405-0054 with: @@ -282,7 +282,7 @@ jobs: uses: docker://ghcr.io/su2code/su2/test-su2:260405-0054 with: # -t -c - args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} -a "--tapetests" + args: -b ${{github.ref}} -t develop -c feature_Trans_SLM -s ${{matrix.testscript}} -a "--tapetests" - name: Cleanup uses: docker://ghcr.io/su2code/su2/test-su2:260405-0054 with: @@ -330,7 +330,7 @@ jobs: PMIX_MCA_gds: hash with: # -t -c - args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} -a "--tsan" + args: -b ${{github.ref}} -t develop -c feature_Trans_SLM -s ${{matrix.testscript}} -a "--tsan" - name: Cleanup uses: docker://ghcr.io/su2code/su2/test-su2-tsan:260405-0054 with: @@ -375,7 +375,7 @@ jobs: uses: docker://ghcr.io/su2code/su2/test-su2-asan:260405-0054 with: # -t -c - args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} -a "--asan" + args: -b ${{github.ref}} -t develop -c feature_Trans_SLM -s ${{matrix.testscript}} -a "--asan" - name: Cleanup uses: docker://ghcr.io/su2code/su2/test-su2-asan:260405-0054 with: diff --git a/Common/include/geometry/dual_grid/CPoint.hpp b/Common/include/geometry/dual_grid/CPoint.hpp index bac83ccbfd2d..b23c30b39a43 100644 --- a/Common/include/geometry/dual_grid/CPoint.hpp +++ b/Common/include/geometry/dual_grid/CPoint.hpp @@ -111,6 +111,7 @@ class CPoint { 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 before solver iteration. */ @@ -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,23 @@ 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 (auto iDim = 0u; 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. @@ -902,4 +941,19 @@ 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 < GlobalIndex.size(); ++iPoint) { + auto rankID = ClosestWall_Rank[iPoint]; + auto zoneID = ClosestWall_Zone[iPoint]; + auto markerID = ClosestWall_Marker[iPoint]; + auto elementID = ClosestWall_Elem[iPoint]; + if (rankID >= 0) SetNormal(iPoint, normals[rankID][zoneID][markerID][elementID]); + } + } }; diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 351e534671ae..7dbe6bdabcae 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -1373,27 +1373,36 @@ 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. */ + 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("LM2015", LM_OPTIONS::CROSSFLOW) // name of the option before the one-equation model was added + MakePair("SLM", LM_OPTIONS::SLM) 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) }; @@ -1411,13 +1420,34 @@ 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 Lower limit of the roughness height (HROUGHNESS, same units) in the cross-flow corrections of the LM models, + * where it enters a logarithm (Langtry et al., log(h/theta_t)) or a ratio (Vallinayagam Pillai and Lardeau, + * h/h0). The papers give no calibration limit for small heights, so a small value keeps these terms finite + * for smooth surfaces (h = 0). h and theta_t are both in mesh length units (the reference length of the + * non-dimensionalization is 1), so log(h/theta_t) is consistent for any REF_DIMENSIONALIZATION. + */ +constexpr double LM_CROSSFLOW_MIN_ROUGHNESS = 1e-8; + /*! * \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 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; }; /*! @@ -1435,7 +1465,9 @@ 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.CrossFlow = IsPresent(LM_OPTIONS::CROSSFLOW); int NFoundCorrelations = 0; if (IsPresent(LM_OPTIONS::MALAN)) { @@ -1467,10 +1499,28 @@ 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); } + if (NFoundCorrelations_SLM > 1) { + SU2_MPI::Error("Two SLM correlations selected for LM_OPTIONS. Please choose only one.", CURRENT_FUNCTION); + } + if (LMParsedOptions.Correlation == TURB_TRANS_CORRELATION::DEFAULT){ if (Kind_Turb_Model == TURB_MODEL::SST) { LMParsedOptions.Correlation = TURB_TRANS_CORRELATION::MENTER_LANGTRY; @@ -1479,6 +1529,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 1ab7ced2d354..e0d741451a59 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -3844,13 +3844,40 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i SU2_MPI::Error("COMPRESSIBILITY-WILCOX only supported for SOLVER= RANS", CURRENT_FUNCTION); } + /*--- The transition solver has no discrete adjoint: the adjoint solvers do not create it, and the turbulence + * sources would read a missing transition solution. ---*/ + if (DiscreteAdjoint && Kind_Trans_Model != TURB_TRANS_MODEL::NONE) { + SU2_MPI::Error("KIND_TRANS_MODEL= LM is not available for discrete adjoint simulations (MATH_PROBLEM= DISCRETE_ADJOINT).", + CURRENT_FUNCTION); + } + /*--- Postprocess LM_OPTIONS into structure. ---*/ 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); + } + + /*--- The simplified (one-equation) model is available only for the combinations found in the literature: + * SST + MENTER_SLM (+ CROSSFLOW): Menter et al., Flow Turbul. Combust. 95, 2015, with the cross-flow extension of + * Vallinayagam Pillai and Lardeau, AIAA 2017-3159. + * SA + MENTER_SLM (+ CROSSFLOW): Lee and Baeder, AIAA 2021-1532. + * The CODER_SLM and MOD_EPPLER_SLM correlations (Coder and Maughmer, AIAA 2012-672) are defined for the + * Langtry-Menter intermittency equation, not for the one-equation model implemented here. ---*/ + if (lmParsedOptions.SLM) { + if (lmParsedOptions.Correlation_SLM != TURB_TRANS_CORRELATION_SLM::MENTER_SLM) { + SU2_MPI::Error("The CODER_SLM and MOD_EPPLER_SLM correlations of LM_OPTIONS are not available with the " + "one-equation (SLM) transition model:\nCoder and Maughmer (AIAA 2012-672) use them with the " + "Langtry-Menter intermittency equation.\nSupported combinations: SLM with MENTER_SLM and " + "KIND_TURB_MODEL= SST (Menter et al. 2015) or SA (Lee and Baeder, AIAA 2021-1532), " + "with or without CROSSFLOW.", CURRENT_FUNCTION); + } + if (Kind_Turb_Model != TURB_MODEL::SST && Kind_Turb_Model != TURB_MODEL::SA) { + SU2_MPI::Error("The one-equation (SLM) transition model is available only with KIND_TURB_MODEL= SST or SA.", + CURRENT_FUNCTION); + } } } @@ -6706,33 +6733,79 @@ 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"; + if (!lmParsedOptions.SLM) { + cout << " (2015)"; + } + cout << endl; + cout << "Roughness height of the cross-flow model (HROUGHNESS, in mesh length units): limited to at least " + << LM_CROSSFLOW_MIN_ROUGHNESS << " to keep log(h/theta_t) and h/h0 finite (the papers give no\n" + << "calibration limit for small heights); "; + if (hRoughness < LM_CROSSFLOW_MIN_ROUGHNESS) + cout << "the given value " << hRoughness << " is below it, " << LM_CROSSFLOW_MIN_ROUGHNESS << " is used." << endl; + else + cout << "the given value " << hRoughness << " is used." << endl; + } else { + if (!lmParsedOptions.SLM) { + cout << " (2009)"; + } + cout << endl; } break; } } - if (Kind_Trans_Model == TURB_TRANS_MODEL::LM) { - cout << "Correlation Functions: "; - switch (lmParsedOptions.Correlation) { - case TURB_TRANS_CORRELATION::MALAN: cout << "Malan et al. (2009)" << endl; break; - case TURB_TRANS_CORRELATION::SULUKSNA: cout << "Suluksna et al. (2009)" << endl; break; - case TURB_TRANS_CORRELATION::KRAUSE: cout << "Krause et al. (2008)" << endl; break; - case TURB_TRANS_CORRELATION::KRAUSE_HYPER: cout << "Krause et al. (2008, paper)" << endl; break; - case TURB_TRANS_CORRELATION::MEDIDA_BAEDER: cout << "Medida and Baeder (2011)" << endl; break; - case TURB_TRANS_CORRELATION::MEDIDA: cout << "Medida PhD (2014)" << endl; break; - case TURB_TRANS_CORRELATION::MENTER_LANGTRY: cout << "Menter and Langtry (2009)" << endl; break; - case TURB_TRANS_CORRELATION::DEFAULT: - switch (Kind_Turb_Model) { - case TURB_MODEL::SA: cout << "Malan et al. (2009)" << endl; break; - case TURB_MODEL::SST: cout << "Menter and Langtry (2009)" << endl; break; - case TURB_MODEL::NONE: SU2_MPI::Error("No turbulence model has been selected but LM transition model is active.", CURRENT_FUNCTION); break; - } - break; + if (Kind_Trans_Model == TURB_TRANS_MODEL::LM) { + if (!lmParsedOptions.SLM){ + cout << "Correlation Functions: "; + switch (lmParsedOptions.Correlation) { + case TURB_TRANS_CORRELATION::MALAN: cout << "Malan et al. (2009)" << endl; break; + case TURB_TRANS_CORRELATION::SULUKSNA: cout << "Suluksna et al. (2009)" << endl; break; + case TURB_TRANS_CORRELATION::KRAUSE: cout << "Krause et al. (2008)" << endl; break; + case TURB_TRANS_CORRELATION::KRAUSE_HYPER: cout << "Krause et al. (2008, paper)" << endl; break; + case TURB_TRANS_CORRELATION::MEDIDA_BAEDER: cout << "Medida and Baeder (2011)" << endl; break; + case TURB_TRANS_CORRELATION::MEDIDA: cout << "Medida PhD (2014)" << endl; break; + case TURB_TRANS_CORRELATION::MENTER_LANGTRY: cout << "Menter and Langtry (2009)" << endl; break; + case TURB_TRANS_CORRELATION::DEFAULT: + switch (Kind_Turb_Model) { + case TURB_MODEL::SA: cout << "Malan et al. (2009)" << endl; break; + case TURB_MODEL::SST: cout << "Menter and Langtry (2009)" << endl; break; + case TURB_MODEL::NONE: SU2_MPI::Error("No turbulence model has been selected but LM transition model is active.", CURRENT_FUNCTION); break; + } + break; + } + } + else { + 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; + } + const bool menterSLM = lmParsedOptions.Correlation_SLM == TURB_TRANS_CORRELATION_SLM::MENTER_SLM || + lmParsedOptions.Correlation_SLM == TURB_TRANS_CORRELATION_SLM::DEFAULT; + if (menterSLM && Kind_Turb_Model == TURB_MODEL::SST) { + cout << "WARNING: Menter et al. (2015, Eq. 24) compute the k production of the SST model with the\n" + " Kato-Launder form P_k = mu_t*S*Omega and without the SST production limiter.\n"; + if (sstParsedOptions.production != SST_OPTIONS::KL) + cout << " The selected SST options do not use Kato-Launder, add KATO-LAUNDER to SST_OPTIONS.\n"; + cout << " The SST production limiter is active." << endl; + } } } cout << "Hybrid RANS/LES: "; diff --git a/Common/src/geometry/CGeometry.cpp b/Common/src/geometry/CGeometry.cpp index 360da1aaaa1e..1d223cb2cb0b 100644 --- a/Common/src/geometry/CGeometry.cpp +++ b/Common/src/geometry/CGeometry.cpp @@ -4538,6 +4538,52 @@ su2double NearestNeighborDistance(CGeometry* geometry, const CConfig* config, co const su2double Vol = geometry->nodes->GetVolume(iPoint) + geometry->nodes->GetPeriodicVolume(iPoint); return 2 * Vol / GeometryToolbox::Norm(3, Normal); } + +/*! \brief Average the vertex normals of a wall element and normalize the result. */ +std::array WallElementNormal(const CGeometry* geometry, unsigned short iMarker, unsigned long iElem) { + std::array normal = {}; + const auto* element = geometry->bound[iMarker][iElem]; + for (auto iNode = 0u; iNode < element->GetnNodes(); ++iNode) { + const auto iPoint = element->GetNode(iNode); + const auto iVertex = geometry->nodes->GetVertex(iPoint, iMarker); + for (auto iDim = 0u; iDim < geometry->GetnDim(); ++iDim) + normal[iDim] += geometry->vertex[iMarker][iVertex]->GetNormal(iDim); + } + for (auto& component : normal) component /= element->GetnNodes(); + const auto magnitude = GeometryToolbox::Norm(3, normal.data()); + for (auto& component : normal) component /= magnitude; + return normal; +} + +/*! \brief Collect unit wall-element normals for each marker of one zone. */ +su2vector> CollectWallNormals(const CGeometry* geometry, const CConfig* config) { + su2vector> normals; + normals.resize(geometry->GetnMarker()); + for (auto iMarker = 0u; iMarker < geometry->GetnMarker(); ++iMarker) { + if (!config->GetViscous_Wall(iMarker)) { + normals[iMarker].resize(1, 3) = su2double(0.0); + continue; + } + normals[iMarker].resize(geometry->GetnElem_Bound(iMarker), 3); + for (auto iElem = 0ul; iElem < geometry->GetnElem_Bound(iMarker); ++iElem) { + const auto normal = WallElementNormal(geometry, iMarker, iElem); + for (auto iDim = 0u; iDim < 3; ++iDim) normals[iMarker](iElem, iDim) = normal[iDim]; + } + } + return normals; +} + +/*! \brief Adapt one marker's normal matrix to the recursive NdFlattener interface. */ +auto WallNormalMatrix(const su2matrix& normals) { + return make_pair(normals.rows(), [&normals](auto iElem) { + return make_pair(normals.cols(), [&normals, iElem](auto iDim) { return normals(iElem, iDim); }); + }); +} + +/*! \brief Adapt a zone's marker normals without copying the matrices. */ +auto WallNormalMarkers(const su2vector>& normals) { + return make_pair(normals.size(), [&normals](auto iMarker) { return WallNormalMatrix(normals[iMarker]); }); +} } // namespace void CGeometry::ComputeWallDistance(const CConfig* const* config_container, CGeometry**** geometry_container, @@ -4639,5 +4685,25 @@ void CGeometry::ComputeWallDistance(const CConfig* const* config_container, CGeo END_SU2_OMP_FOR } } + + MAIN_SOLVER kindSolver = config_container[ZONE_0]->GetKind_Solver(); + if (!wallDistanceNeeded[ZONE_0] || kindSolver == MAIN_SOLVER::FEM_LES || kindSolver == MAIN_SOLVER::FEM_RANS) { + continue; + } else { + su2vector>> wallNormals; + wallNormals.resize(nZone); + for (auto iZone = 0; iZone < nZone; ++iZone) + wallNormals[iZone] = CollectWallNormals(geometry_container[iZone][iInst][MESH_0], config_container[iZone]); + + auto normal_i = make_pair(nZone, [&wallNormals](auto iZone) { return WallNormalMarkers(wallNormals[iZone]); }); + + NdFlattener<4> Normals_Local(normal_i); + NdFlattener<5> Normals_global(Nd_MPI_Environment(), Normals_Local); + + // Use the gathered normals of the nearest wall elements. + for (int jZone = 0; jZone < nZone; jZone++) { + geometry_container[jZone][iInst][MESH_0]->nodes->SetWallNormals(Normals_global); + } + } } } diff --git a/Common/src/geometry/dual_grid/CPoint.cpp b/Common/src/geometry/dual_grid/CPoint.cpp index 24616df65cc6..0f935659b729 100644 --- a/Common/src/geometry/dual_grid/CPoint.cpp +++ b/Common/src/geometry/dual_grid/CPoint.cpp @@ -127,6 +127,8 @@ 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 ec511a3fa1db..593061f334af 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -28,6 +28,8 @@ #pragma once +#include "../transition_data.hpp" + #include #include #include @@ -808,6 +810,21 @@ class CNumerics { */ su2double GetIntermittencyEff() const { return intermittency_eff_i; } + /*! \brief Diagnostic values produced by the transition-model source term. */ + virtual const TransitionLMData* GetTransitionData() const { return nullptr; } + + /*! + * \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 cross-flow strength Psi = |n . grad(e_omega)| d_w of the one-equation transition model. + */ + inline virtual void SetCrossFlowStrength(su2double val_Psi) {} + /*! * \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_correlations.hpp b/SU2_CFD/include/numerics/turbulent/transition/trans_correlations.hpp index 4a6be1c6c2a6..eff4fe6c75f9 100644 --- a/SU2_CFD/include/numerics/turbulent/transition/trans_correlations.hpp +++ b/SU2_CFD/include/numerics/turbulent/transition/trans_correlations.hpp @@ -37,7 +37,126 @@ class TransLMCorrelations { LM_ParsedOptions options; + su2double MenterSLM(const su2double Tu_L, const su2double lambda_theta, const bool spalartAllmaras) const { + su2double rethetac = 0.0; + + /*-- 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); + + /*--- Menter et al. (2015), Eq. 14. ---*/ + su2double C_TU1 = 100.0; + su2double C_TU2 = 1000.0; + const su2double C_TU3 = 1.0; + + if (spalartAllmaras) { + /*--- Lee and Baeder, AIAA 2021-1532, Eqs. 11-13: the constants of Colonia et al. (Eq. 12) below Tu = 0.51%, + * the original ones (Eq. 11) above Tu = 2%, linearly blended in between. Tu_L is the freestream Tu. ---*/ + const su2double Tu_blend = min(max(Tu_L, 0.51), 2.0); + C_TU1 = (100.0 - 163.0) / (2.0 - 0.51) * (Tu_blend - 2.0) + 100.0; + C_TU2 = (1000.0 - 1002.25) / (2.0 - 0.51) * (Tu_blend - 2.0) + 1000.0; + } + rethetac = C_TU1 + C_TU2 * exp(-C_TU3 * Tu_L * FPG); + + return rethetac; + + } + + su2double CoderSLM(const su2double Tu_L, const su2double wall_dist, const su2double VorticityMag, const su2double VelocityMag) const { + su2double rethetac = 0.0; + + /*-- 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); + } + + 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; + + return rethetac; + + } + + su2double ModifiedEpplerSLM(const su2double wall_dist, const su2double VorticityMag, const su2double VelocityMag) const { + su2double rethetac = 0.0; + + /*-- 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); + + return rethetac; + + } + + public: + /*! \brief Langtry-Menter (2009) farfield correlation, with Tu in percent. + * See https://tmbwg.github.io/turbmodels/langtrymenter_4eqn.html (farfield boundary condition). + * The lower Tu limit avoids the inverse-square singularity. */ + static su2double FreestreamReThetaT(su2double intensity) { + constexpr passivedouble minIntensity = 0.027; + constexpr passivedouble intensitySwitch = 1.3; + constexpr passivedouble lowIntensityOffset = 1173.51; + constexpr passivedouble lowIntensitySlope = 589.428; + constexpr passivedouble inverseSquareCoefficient = 0.2196; + constexpr passivedouble highIntensityCoefficient = 331.5; + constexpr passivedouble highIntensityOffset = 0.5658; + constexpr passivedouble highIntensityExponent = -0.671; + intensity = max(intensity, minIntensity); + if (intensity <= intensitySwitch) + return lowIntensityOffset - lowIntensitySlope * intensity + inverseSquareCoefficient / (intensity * intensity); + return highIntensityCoefficient * pow(intensity - highIntensityOffset, highIntensityExponent); + } /*! * \brief Set LM options. @@ -189,4 +308,36 @@ class TransLMCorrelations { return F_length1; } + + /*! + * \brief Compute Re_theta_c from correlations for the Simplified LM model. + * \param[in] Tu_L - Turbulence intensity in percent. + * \param[in] lambda_theta - Local pressure-gradient parameter. + * \param[in] wall_dist - Wall distance. + * \param[in] VorticityMag - Vorticity magnitude. + * \param[in] VelocityMag - Velocity magnitude. + * \param[in] spalartAllmaras - Use the SA coefficients of Lee and Baeder. + * \return Critical transition Reynolds number. + */ + su2double ReThetaC_Correlations_SLM(const su2double Tu_L, const su2double lambda_theta, const su2double wall_dist, + const su2double VorticityMag, const su2double VelocityMag, + const bool spalartAllmaras = false) const { + + switch (options.Correlation_SLM) { + case TURB_TRANS_CORRELATION_SLM::MENTER_SLM: + return MenterSLM(Tu_L, lambda_theta, spalartAllmaras); + case TURB_TRANS_CORRELATION_SLM::CODER_SLM: + return CoderSLM(Tu_L, wall_dist, VorticityMag, VelocityMag); + case TURB_TRANS_CORRELATION_SLM::MOD_EPPLER_SLM: + return ModifiedEpplerSLM(wall_dist, VorticityMag, VelocityMag); + + 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 0.0; + } + }; diff --git a/SU2_CFD/include/numerics/turbulent/transition/trans_edge_flux.hpp b/SU2_CFD/include/numerics/turbulent/transition/trans_edge_flux.hpp index 927fcb6f5b02..daca1d982629 100644 --- a/SU2_CFD/include/numerics/turbulent/transition/trans_edge_flux.hpp +++ b/SU2_CFD/include/numerics/turbulent/transition/trans_edge_flux.hpp @@ -32,7 +32,8 @@ /*! * \class CScalarFlux_TransLM * \ingroup ViscDiscr - * \brief Convection and diffusion of the Langtry-Menter transition model, conservative with a + * \brief Convection and diffusion of the Langtry-Menter transition model (nVar = 2) and of its simplified, + * one-equation version (nVar = 1, intermittency only), conservative with a * diagonal, i/j-symmetric diffusion matrix. The coefficients depend only on the flow's * mu/mu_t, not on the transported gamma/Re_theta, so no coefficientJacobians override * is needed. @@ -67,12 +68,16 @@ class CScalarFlux_TransLM const Double diff_i_gamma = mu_i + muT_i; const Double diff_j_gamma = mu_j + muT_j; - const Double diff_i_ReThetaT = 2.0 * (mu_i + muT_i); - const Double diff_j_ReThetaT = 2.0 * (mu_j + muT_j); Vector D; D(0) = 0.5 * (diff_i_gamma + diff_j_gamma); - D(1) = 0.5 * (diff_i_ReThetaT + diff_j_ReThetaT); + + /*--- The simplified (one-equation) model only transports the intermittency. ---*/ + if constexpr (nVar > 1) { + const Double diff_i_ReThetaT = 2.0 * (mu_i + muT_i); + const Double diff_j_ReThetaT = 2.0 * (mu_j + muT_j); + D(1) = 0.5 * (diff_i_ReThetaT + diff_j_ReThetaT); + } return {D, D}; } }; diff --git a/SU2_CFD/include/numerics/turbulent/transition/trans_sources.hpp b/SU2_CFD/include/numerics/turbulent/transition/trans_sources.hpp index 73f37cc8defe..d9de034322b8 100644 --- a/SU2_CFD/include/numerics/turbulent/transition/trans_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/transition/trans_sources.hpp @@ -47,17 +47,14 @@ class CSourcePieceWise_TransLM final : public CNumerics { 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; + TransitionLMData TransitionData; + su2double Residual[2]; su2double* Jacobian_i[2]; @@ -65,6 +62,165 @@ class CSourcePieceWise_TransLM final : public CNumerics { TransLMCorrelations TransCorrelations; + /*! \brief Streamwise derivative of the velocity magnitude. */ + su2double StreamwiseVelocityGradient(const su2double Velocity_Mag) const { + 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; + /*-- Gradient of velocity magnitude ---*/ + + su2double dU_dx = 0.5 / Velocity_Mag * (2. * vel_u * PrimVar_Grad_i[1][0] + 2. * vel_v * PrimVar_Grad_i[2][0]); + if (nDim == 3) dU_dx += 0.5 / Velocity_Mag * (2. * vel_w * PrimVar_Grad_i[3][0]); + + su2double dU_dy = 0.5 / Velocity_Mag * (2. * vel_u * PrimVar_Grad_i[1][1] + 2. * vel_v * PrimVar_Grad_i[2][1]); + if (nDim == 3) dU_dy += 0.5 / Velocity_Mag * (2. * vel_w * PrimVar_Grad_i[3][1]); + + su2double dU_dz = 0.0; + if (nDim == 3) + dU_dz = + 0.5 / Velocity_Mag * + (2. * vel_u * PrimVar_Grad_i[1][2] + 2. * vel_v * PrimVar_Grad_i[2][2] + 2. * vel_w * PrimVar_Grad_i[3][2]); + + su2double du_ds = vel_u / Velocity_Mag * dU_dx + vel_v / Velocity_Mag * dU_dy; + if (nDim == 3) du_ds += vel_w / Velocity_Mag * dU_dz; + return du_ds; + } + + /*! \brief Fixed-point pressure-gradient correction of the transition Reynolds number. */ + su2double TransitionReynolds(const su2double Tu, const su2double du_ds, const su2double Velocity_Mag) { + /*--- Corr_Ret correlation*/ + const su2double Corr_Ret_lim = 20.0; + su2double f_lambda = 1.0; + + su2double Retheta_Error = 200.0, Retheta_old = 1.0; + su2double lambda = 0.0; + su2double Corr_Ret = 20.0; + + for (int iter = 0; iter < 100; iter++) { + su2double theta = Corr_Ret * Laminar_Viscosity_i / Density_i / Velocity_Mag; + lambda = Density_i * theta * theta / Laminar_Viscosity_i * du_ds; + lambda = min(max(-0.1, lambda), 0.1); + + if (lambda <= 0.0) { + f_lambda = 1. - (-12.986 * lambda - 123.66 * lambda * lambda - 405.689 * lambda * lambda * lambda) * + exp(-pow(Tu / 1.5, 1.5)); + } else { + f_lambda = 1. + 0.275 * (1. - exp(-35. * lambda)) * exp(-Tu / 0.5); + } + + if (Tu <= 1.3) { + Corr_Ret = f_lambda * (1173.51 - 589.428 * Tu + 0.2196 / Tu / Tu); + } else { + Corr_Ret = 331.5 * f_lambda * pow(Tu - 0.5658, -0.671); + } + Corr_Ret = max(Corr_Ret, Corr_Ret_lim); + + Retheta_Error = fabs(Retheta_old - Corr_Ret) / Retheta_old; + + if (Retheta_Error < 0.0000001) { + break; + } + + Retheta_old = Corr_Ret; + } + + // Stored for the output + TransitionData.pressureGradient = lambda; + TransitionData.streamwiseVelocityGradient = du_ds; + return Corr_Ret; + } + + /*! \brief Fixed-point stationary cross-flow correlation of LM2015. */ + su2double CrossFlowReynolds(const su2double Velocity_Mag) const { + 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 ReThetat_SCF = 0.0; + 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; + const su2double DeltaH_CF = H_CF * (1.0 + min(Eddy_Viscosity_i / Laminar_Viscosity_i, 0.4)); + const su2double DeltaH_CF_Minus = max(-1.0 * (0.1066 - DeltaH_CF), 0.0); + const su2double DeltaH_CF_Plus = max(0.1066 - DeltaH_CF, 0.0); + const su2double fDeltaH_CF_Minus = 75.0 * tanh(DeltaH_CF_Minus / 0.0125); + const su2double fDeltaH_CF_Plus = 6200 * DeltaH_CF_Plus + 50000 * DeltaH_CF_Plus * DeltaH_CF_Plus; + + const su2double toll = 1e-5; + su2double error = toll + 1.0; + su2double thetat_SCF = 0.0; + su2double rethetat_SCF_old = 20.0; + const int nMax = 100; + + int iter; + for (iter = 0; iter < nMax && error > toll; iter++) { + thetat_SCF = rethetat_SCF_old * Laminar_Viscosity_i / (Density_i * (Velocity_Mag / 0.82)); + thetat_SCF = max(1e-20, thetat_SCF); + + ReThetat_SCF = -35.088 * log(max(hRoughness, su2double(LM_CROSSFLOW_MIN_ROUGHNESS)) / thetat_SCF) + 319.51 + + fDeltaH_CF_Plus - fDeltaH_CF_Minus; + + error = abs(ReThetat_SCF - rethetat_SCF_old) / rethetat_SCF_old; + + rethetat_SCF_old = ReThetat_SCF; + } + return ReThetat_SCF; + } + + /*! \brief Transition onset and length correlations, including near-wall blending. */ + std::pair IntermittencyOnset(const su2double Tu, const su2double R_t) { + /*--- Corr_RetC correlation*/ + const su2double Corr_Rec = TransCorrelations.ReThetaC_Correlations(Tu, TransVar_i[1]); + // Stored for the output + TransitionData.criticalReynolds = Corr_Rec; + + /*--- F_length correlation*/ + const su2double Corr_F_length = TransCorrelations.FLength_Correlations(Tu, TransVar_i[1]); + + /*--- F_length ---*/ + su2double F_length = 0.0; + if (TurbFamily == TURB_FAMILY::KW) { + const su2double r_omega = Density_i * dist_i * dist_i * ScalarVar_i[1] / Laminar_Viscosity_i; + const su2double f_sub = exp(-pow(r_omega / 200.0, 2)); + F_length = Corr_F_length * (1. - f_sub) + 40.0 * f_sub; + } + if (TurbFamily == TURB_FAMILY::SA) F_length = Corr_F_length; + + /*--- F_onset ---*/ + + const su2double Re_v = Density_i * dist_i * dist_i * StrainMag_i / Laminar_Viscosity_i; + // Stored for the output + TransitionData.vorticityReynolds = Re_v; + + const su2double F_onset1 = Re_v / (2.193 * Corr_Rec); + su2double F_onset2 = 1.0; + su2double F_onset3 = 1.0; + if (TurbFamily == TURB_FAMILY::KW) { + F_onset2 = min(max(F_onset1, pow(F_onset1, 4.0)), 2.0); + F_onset3 = max(1.0 - pow(R_t / 2.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); + } + const su2double F_onset = max(F_onset2 - F_onset3, 0.0); + // Stored for the output + TransitionData.onset1 = F_onset1; + TransitionData.onset2 = F_onset2; + TransitionData.onset3 = F_onset3; + TransitionData.onset = F_onset; + return {F_length, F_onset}; + } + + public: /*! * \brief Constructor of the class. @@ -92,12 +248,14 @@ 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[1] = omega, TransVar[0] = gamma, and TransVar[1] = ReThetaT ---*/ /*--- dU/dx = PrimVar_Grad[1][0] ---*/ AD::StartPreacc(); AD::SetPreaccIn(StrainMag_i); - AD::SetPreaccIn(ScalarVar_i, nVar); - AD::SetPreaccIn(ScalarVar_Grad_i, nVar, nDim); + /*--- The turbulence model has one variable (SA) or two (k-omega). ---*/ + const unsigned short nTurbVar = (TurbFamily == TURB_FAMILY::SA) ? 1 : 2; + AD::SetPreaccIn(ScalarVar_i, nTurbVar); + AD::SetPreaccIn(ScalarVar_Grad_i, nTurbVar, nDim); AD::SetPreaccIn(TransVar_i, nVar); AD::SetPreaccIn(TransVar_Grad_i, nVar, nDim); AD::SetPreaccIn(Volume); @@ -106,14 +264,9 @@ class CSourcePieceWise_TransLM final : public CNumerics { 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 VorticityMag = max(GeometryToolbox::Norm(3, Vorticity_i), 1e-20); - 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; - - const su2double Velocity_Mag = sqrt(vel_u * vel_u + vel_v * vel_v + vel_w * vel_w); + const su2double Velocity_Mag = max(GeometryToolbox::Norm(nDim, &V_i[idx.Velocity()]), 1e-20); AD::SetPreaccIn(V_i[idx.Density()], V_i[idx.LaminarViscosity()], V_i[idx.EddyViscosity()]); @@ -131,62 +284,19 @@ class CSourcePieceWise_TransLM final : public CNumerics { if (dist_i > 1e-10) { su2double Tu = 1.0; if (TurbFamily == TURB_FAMILY::KW) Tu = max(100.0 * sqrt(2.0 * ScalarVar_i[0] / 3.0) / Velocity_Mag, 0.027); - if (TurbFamily == TURB_FAMILY::SA) Tu = config->GetTurbulenceIntensity_FreeStream() * 100; - - /*--- Corr_RetC correlation*/ - const su2double Corr_Rec = TransCorrelations.ReThetaC_Correlations(Tu, TransVar_i[1]); - - /*--- F_length correlation*/ - const su2double Corr_F_length = TransCorrelations.FLength_Correlations(Tu, TransVar_i[1]); + /*--- The same lower limit as for k-omega, the correlations divide by Tu. ---*/ + if (TurbFamily == TURB_FAMILY::SA) Tu = max(config->GetTurbulenceIntensity_FreeStream() * 100, 0.027); - /*--- F_length ---*/ - su2double F_length = 0.0; - if (TurbFamily == TURB_FAMILY::KW) { - const su2double r_omega = Density_i * dist_i * dist_i * ScalarVar_i[1] / Laminar_Viscosity_i; - const su2double f_sub = exp(-pow(r_omega / 200.0, 2)); - F_length = Corr_F_length * (1. - f_sub) + 40.0 * f_sub; - } - if (TurbFamily == TURB_FAMILY::SA) F_length = Corr_F_length; - - /*--- 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; + const auto [F_length, F_onset] = IntermittencyOnset(Tu, R_t); - const su2double Re_v = Density_i * dist_i * dist_i * StrainMag_i / Laminar_Viscosity_i; - const su2double F_onset1 = Re_v / (2.193 * Corr_Rec); - su2double F_onset2 = 1.0; - su2double F_onset3 = 1.0; - if (TurbFamily == TURB_FAMILY::KW) { - F_onset2 = min(max(F_onset1, pow(F_onset1, 4.0)), 2.0); - F_onset3 = max(1.0 - pow(R_t / 2.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); - } - const su2double F_onset = max(F_onset2 - F_onset3, 0.0); - - /*-- Gradient of velocity magnitude ---*/ - - su2double dU_dx = 0.5 / Velocity_Mag * (2. * vel_u * PrimVar_Grad_i[1][0] + 2. * vel_v * PrimVar_Grad_i[2][0]); - if (nDim == 3) dU_dx += 0.5 / Velocity_Mag * (2. * vel_w * PrimVar_Grad_i[3][0]); - - su2double dU_dy = 0.5 / Velocity_Mag * (2. * vel_u * PrimVar_Grad_i[1][1] + 2. * vel_v * PrimVar_Grad_i[2][1]); - if (nDim == 3) dU_dy += 0.5 / Velocity_Mag * (2. * vel_w * PrimVar_Grad_i[3][1]); - - su2double dU_dz = 0.0; - if (nDim == 3) - dU_dz = - 0.5 / Velocity_Mag * - (2. * vel_u * PrimVar_Grad_i[1][2] + 2. * vel_v * PrimVar_Grad_i[2][2] + 2. * vel_w * PrimVar_Grad_i[3][2]); - - su2double du_ds = vel_u / Velocity_Mag * dU_dx + vel_v / Velocity_Mag * dU_dy; - if (nDim == 3) du_ds += vel_w / Velocity_Mag * dU_dz; + const su2double du_ds = StreamwiseVelocityGradient(Velocity_Mag); /*-- 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,99 +316,30 @@ 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 lambda = 0.0; - su2double Corr_Ret = 20.0; + const su2double Corr_Ret = TransitionReynolds(Tu, du_ds, Velocity_Mag); - for (int iter = 0; iter < 100; iter++) { - su2double theta = Corr_Ret * Laminar_Viscosity_i / Density_i / Velocity_Mag; - lambda = Density_i * theta * theta / Laminar_Viscosity_i * du_ds; - lambda = min(max(-0.1, lambda), 0.1); + const su2double ReThetat_SCF = options.CrossFlow ? CrossFlowReynolds(Velocity_Mag) : su2double(0.0); - if (lambda <= 0.0) { - f_lambda = 1. - (-12.986 * lambda - 123.66 * lambda * lambda - 405.689 * lambda * lambda * lambda) * - exp(-pow(Tu / 1.5, 1.5)); - } else { - f_lambda = 1. + 0.275 * (1. - exp(-35. * lambda)) * exp(-Tu / 0.5); - } - - if (Tu <= 1.3) { - Corr_Ret = f_lambda * (1173.51 - 589.428 * Tu + 0.2196 / Tu / Tu); - } else { - Corr_Ret = 331.5 * f_lambda * pow(Tu - 0.5658, -0.671); - } - Corr_Ret = max(Corr_Ret, Corr_Ret_lim); - - Retheta_Error = fabs(Retheta_old - Corr_Ret) / Retheta_old; - - if (Retheta_Error < 0.0000001) { - break; - } - - Retheta_old = Corr_Ret; - } - - /*-- Corr_RetT_SCF Correlations--*/ - su2double ReThetat_SCF = 0.0; - if (options.LM2015) { - 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; - const su2double DeltaH_CF = H_CF * (1.0 + min(Eddy_Viscosity_i / Laminar_Viscosity_i, 0.4)); - const su2double DeltaH_CF_Minus = max(-1.0 * (0.1066 - DeltaH_CF), 0.0); - const su2double DeltaH_CF_Plus = max(0.1066 - DeltaH_CF, 0.0); - const su2double fDeltaH_CF_Minus = 75.0 * tanh(DeltaH_CF_Minus / 0.0125); - const su2double fDeltaH_CF_Plus = 6200 * DeltaH_CF_Plus + 50000 * DeltaH_CF_Plus * DeltaH_CF_Plus; - - const su2double toll = 1e-5; - su2double error = toll + 1.0; - su2double thetat_SCF = 0.0; - su2double rethetat_SCF_old = 20.0; - const int nMax = 100; - - int iter; - for (iter = 0; iter < nMax && error > toll; iter++) { - thetat_SCF = rethetat_SCF_old * Laminar_Viscosity_i / (Density_i * (Velocity_Mag / 0.82)); - thetat_SCF = max(1e-20, thetat_SCF); - - ReThetat_SCF = -35.088 * log(hRoughness / thetat_SCF) + 319.51 + fDeltaH_CF_Plus - fDeltaH_CF_Minus; - - error = abs(ReThetat_SCF - rethetat_SCF_old) / rethetat_SCF_old; - - rethetat_SCF_old = ReThetat_SCF; - } - } - - /*-- production term of Intermeittency(Gamma) --*/ + /*-- production term of intermittency (gamma) --*/ const su2double Pg = F_length * c_a1 * Density_i * StrainMag_i * sqrt(F_onset * TransVar_i[0]) * (1.0 - c_e1 * TransVar_i[0]); - /*-- destruction term of Intermeittency(Gamma) --*/ + /*-- destruction term of intermittency (gamma) --*/ const su2double Dg = c_a2 * Density_i * VorticityMag * TransVar_i[0] * f_turb * (c_e2 * TransVar_i[0] - 1.0); + // Stored for the output + TransitionData.production = Pg; + TransitionData.destruction = 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 +354,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 +363,292 @@ class CSourcePieceWise_TransLM final : public CNumerics { return ResidualType<>(Residual, Jacobian_i, nullptr); } + + const TransitionLMData* GetTransitionData() const override { return &TransitionData; } + +}; + +/*! + * \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_e2 = 50.0; + const su2double c_a2 = 0.06; + + TURB_FAMILY TurbFamily; + + + TransitionLMData TransitionData; + su2double AuxVar = 0.0; /*!< \brief Wall-normal derivative of the wall-normal velocity (Menter correlation). */ + su2double CrossFlowPsi = 0.0; /*!< \brief Wall-normal change of the vorticity direction times the wall distance. */ + su2double Residual; + su2double* Jacobian_i; + su2double Jacobian_Buffer; // Static storage for the Jacobian (which needs to be pointer for return type). + + TransLMCorrelations TransCorrelations; + + /*! + * \brief Magnitude of the streamwise vorticity, |U/|U| . Omega|. + */ + su2double StreamwiseVorticity(su2double vel_u, su2double vel_v, su2double vel_w, su2double Velocity_Mag) const { + const su2double velocity[3] = {vel_u, vel_v, vel_w}; + su2double streamwiseVort = 0.0; + for (auto iDim = 0u; iDim < nDim; iDim++) streamwiseVort += velocity[iDim] / Velocity_Mag * Vorticity_i[iDim]; + return fabs(streamwiseVort); + } + + /*! + * \brief Stationary cross-flow Reynolds number of Langtry et al., Lee and Baeder (AIAA 2021-1532), Eqs. 36-42. + * Eq. 36 is implicit in theta_t and is solved by fixed-point iteration, as in the LM2015 option of + * CSourcePieceWise_TransLM. + */ + su2double StationaryCrossFlowReynolds(su2double vel_u, su2double vel_v, su2double vel_w, su2double Velocity_Mag, + su2double R_t, const CConfig* config) const { + const su2double H_CF = StreamwiseVorticity(vel_u, vel_v, vel_w, Velocity_Mag) * dist_i / Velocity_Mag; // Eq. 37 + const su2double DeltaH_CF = H_CF * (1.0 + min(R_t, 0.4)); // Eq. 38 + const su2double DeltaH_CF_Plus = max(0.1066 - DeltaH_CF, 0.0); // Eq. 39 + const su2double fDeltaH_CF_Plus = 6200 * DeltaH_CF_Plus + 50000 * DeltaH_CF_Plus * DeltaH_CF_Plus; // Eq. 40 + const su2double DeltaH_CF_Minus = max(-(0.1066 - DeltaH_CF), 0.0); // Eq. 41 + const su2double fDeltaH_CF_Minus = 75.0 * tanh(DeltaH_CF_Minus / 0.0125); // Eq. 42 + + const su2double h = max(config->GethRoughness(), su2double(LM_CROSSFLOW_MIN_ROUGHNESS)); + su2double reScf = 20.0, error = 1.0; + for (int iter = 0; iter < 100 && error > 1e-5; iter++) { + const su2double theta_t = max(1e-20, reScf * Laminar_Viscosity_i / (Density_i * Velocity_Mag / 0.82)); + const su2double reScfNew = -35.088 * log(h / theta_t) + 319.51 + fDeltaH_CF_Plus - fDeltaH_CF_Minus; + error = fabs(reScfNew - reScf) / fabs(reScf); + reScf = reScfNew; + } + return max(reScf, 1e-20); + } + + /*! \brief Menter-Smirnov cross-flow onset for the SA-based simplified model. */ + su2double SACrossFlowOnset(const su2double Re_v) const { + /*--- Menter and Smirnov C1-based cross-flow criterion, Lee and Baeder (2021), Eqs. 25-35. ---*/ + const su2double lambda_CF = min(max(-7.57e-3 * AuxVar * dist_i * dist_i * Density_i / Laminar_Viscosity_i + 0.0174, 0.0), 0.0477); + const su2double g_CF = min(max(27864.0 * pow(lambda_CF, 3) - 1962.0 * pow(lambda_CF, 2) + 54.3 * lambda_CF + 1.0, 1.0), 2.3); + const su2double G_CF = 0.684 / g_CF; + const su2double C_RSF = 1.35; // Calibrated for SA by Lee and Baeder (1.0 in the original model). + const su2double T_C1 = C_RSF / 150.0 * G_CF * CrossFlowPsi * Re_v; + const su2double F_onset_CF = min(max(100.0 * (T_C1 - 1.0), 0.0), 1.0); + return F_onset_CF; + + } + + /*! \brief Cross-flow onset for the SST-based simplified model. */ + su2double SSTCrossFlowOnset(const su2double lambda_theta, const su2double vel_u, const su2double vel_v, const su2double vel_w, const su2double Velocity_Mag, const CConfig* config) const { + + /*--- Vallinayagam Pillai and Lardeau, "Accounting crossflow effects in one-equation local correlation-based + * transition model", AIAA 2017-3159 (equation numbers below). ---*/ + + // Shape factor, Eq. 3 with k = 0.25 - lambda (lambda_theta_L of the gamma model), limited to 2.7 above which + // the crossflow criterion does not apply (text after Eq. 3) + 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); + + // Critical crossflow Reynolds number, Eq. 2 + su2double Re_Crit_CF = 0.0; + if(H < 2.3) { + Re_Crit_CF = 150.0; + } else { + // Eq. 2 is printed with a minus sign, which gives -150 at H = 2.3 and negative critical Reynolds numbers; + // the positive sign makes it continuous with the value 150 below H = 2.3 + Re_Crit_CF = (300.0/PI_NUMBER) * atan(0.106/(pow(H-2.3, 2.05))); + } + + const su2double H_CF = StreamwiseVorticity(vel_u, vel_v, vel_w, Velocity_Mag) * dist_i / Velocity_Mag; + + // Crossflow strength H_cf, Eqs. 4-6, and Delta_H_cf, Eq. 8 + const su2double Delta_H_CF = H_CF * (1.0 + min(Eddy_Viscosity_i / Laminar_Viscosity_i, 0.4)); + + // Roughness, Eq. 9, h0 = 0.25 micrometers: HROUGHNESS and the mesh must be in meters + const su2double h_0 = 0.25e-6; + const su2double C_r = 2.0 - pow(0.5, max(config->GethRoughness(), LM_CROSSFLOW_MIN_ROUGHNESS)/h_0); + + // f_cf, Eqs. 7 and 10, with C_cf = 1 + const su2double C_CF = 1.0; + const su2double f_CF = (C_CF * C_r * Delta_H_CF * TransitionData.criticalReynolds) / Re_Crit_CF; + const su2double F_onset_CF = min(max(0.0, f_CF - 1.0), 1.0); + + // Eqs. 17-18 + return F_onset_CF; + + + } + + 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()); + + + 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[1] = omega, TransVar[0] = gamma ---*/ + /*--- dU/dx = PrimVar_Grad[1][0] ---*/ + AD::StartPreacc(); + AD::SetPreaccIn(StrainMag_i); + /*--- ScalarVar_i holds the turbulence variables (k and omega for SST), nVar is the single transition variable. ---*/ + const unsigned short nVarTurb = (TurbFamily == TURB_FAMILY::KW) ? 2 : 1; + AD::SetPreaccIn(ScalarVar_i, nVarTurb); + AD::SetPreaccIn(ScalarVar_Grad_i, nVarTurb, nDim); + AD::SetPreaccIn(AuxVar); + AD::SetPreaccIn(CrossFlowPsi); + 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); + + const su2double VorticityMag = GeometryToolbox::Norm(3, Vorticity_i); + + 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; + + 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; + // Local value of the Turbulence intensity that makes it galileian invariant. Look at Eq. 7 in https://doi.org/10.1007/s10494-015-9622-4 + 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::SA) Tu_L = config->GetTurbulenceIntensity_FreeStream() * 100; + + TransitionData.turbulenceIntensity = Tu_L; + + /*--- F_length ---*/ + const su2double F_length = 100.0; // Menter et al. (2015), Eq. 6, and Lee and Baeder (2021), Eq. 7. + + /*--- 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]); + /*--- Lee and Baeder (AIAA 2021-1532), Eq. 6: k and omega are replaced by the viscosity ratio for SA. ---*/ + if (TurbFamily == TURB_FAMILY::SA) R_t = Eddy_Viscosity_i / Laminar_Viscosity_i; + + /*--- Menter et al. (2015), Eqs. 12-13, and Lee and Baeder (2021), Eq. 9. ---*/ + 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); + + TransitionData.streamwiseVelocityGradient = AuxVar; + TransitionData.pressureGradient = lambda_theta; + + /*--- Critical Reynolds number, Menter et al. (2015), Eq. 14, and Lee and Baeder (2021), Eqs. 10-13. ---*/ + TransitionData.momentumThicknessReynolds = TransCorrelations.ReThetaC_Correlations_SLM(Tu_L, lambda_theta, dist_i, VorticityMag, Velocity_Mag, + TurbFamily == TURB_FAMILY::SA); + TransitionData.criticalReynolds = TransitionData.momentumThicknessReynolds; + + const su2double Re_v = Density_i * dist_i * dist_i * StrainMag_i / Laminar_Viscosity_i; + TransitionData.vorticityReynolds = Re_v; + + /*--- Menter et al. (2015), Eqs. 4-5, and Lee and Baeder (2021), Eqs. 4-5 (the same for SST and SA). ---*/ + su2double F_onset1 = Re_v / (2.2 * TransitionData.criticalReynolds); + + if (options.CrossFlow && TurbFamily == TURB_FAMILY::SA) { + /*--- Langtry et al. stationary cross-flow criterion, Lee and Baeder (2021), Eqs. 36-43. ---*/ + const su2double Re_scf = StationaryCrossFlowReynolds(vel_u, vel_v, vel_w, Velocity_Mag, R_t, config); + const su2double c_scf = 0.94; + F_onset1 = max(F_onset1, c_scf * Re_v / (2.2 * Re_scf)); + } + + const su2double F_onset2 = min(F_onset1, 2.0); + const su2double F_onset3 = max(1.0 - pow(R_t / 3.5, 3.0), 0.0); + su2double F_onset = max(F_onset2 - F_onset3, 0.0); + + TransitionData.onset1 = F_onset1; + TransitionData.onset2 = F_onset2; + TransitionData.onset3 = F_onset3; + + if (options.CrossFlow && TurbFamily == TURB_FAMILY::SA) { + F_onset = max(F_onset, SACrossFlowOnset(Re_v)); + } + + if (options.CrossFlow && TurbFamily == TURB_FAMILY::KW) { + F_onset = max(F_onset, SSTCrossFlowOnset(lambda_theta, vel_u, vel_v, vel_w, Velocity_Mag, config)); + } + + /*--- Output value, including the cross-flow corrections. ---*/ + TransitionData.onset = F_onset; + + /*--- Menter et al. (2015), Eq. 5, and Lee and Baeder (2021), Eq. 6. ---*/ + const su2double f_turb = exp(-pow(R_t / 2, 4)); + + /*--- Production and destruction of the intermittency, Menter et al. (2015), Eqs. 2-3, and Lee and Baeder + * (2021), Eqs. 2-3, written as gamma*(P - D) with (Eqs. 22-23) ---*/ + const su2double P = F_length * StrainMag_i * (1.0 - TransVar_i[0]) * F_onset; + const su2double D = c_a2 * VorticityMag * f_turb * (c_e2 * TransVar_i[0] - 1.0); + const su2double Pg = Density_i * TransVar_i[0] * P; + const su2double Dg = Density_i * TransVar_i[0] * D; + + TransitionData.production = Pg; + TransitionData.destruction = Dg; + + /*--- Source ---*/ + Residual += (Pg - Dg) * Volume; + + /*--- Implicit part, per unit rho*gamma. ---*/ + if (TurbFamily == TURB_FAMILY::SA) { + /*--- Positivity of the implicit operator, Lee and Baeder (2021), Eq. 24. ---*/ + const su2double dP = -F_length * StrainMag_i * F_onset; + const su2double dD = c_a2 * VorticityMag * f_turb * c_e2; + Jacobian_i[0] = -(max(D - P, 0.0) + max(dD - dP, 0.0) * TransVar_i[0]) * Volume; + } else { + 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); + AD::EndPreacc(); + + return ResidualType<>(&Residual, &Jacobian_i, nullptr); + + } + + const TransitionLMData* GetTransitionData() const override { return &TransitionData; } + inline void SetAuxVar(su2double val_AuxVar) override { AuxVar = val_AuxVar;} + inline void SetCrossFlowStrength(su2double val_Psi) override { CrossFlowPsi = val_Psi; } + }; diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp index 9f72b6e37cc4..5cc6c2440e44 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp @@ -80,6 +80,7 @@ class CSourceBase_TurbSA : public CNumerics { const bool axisymmetric = false; bool transition_LM; + bool transition_SLM; /*!< \brief One-equation (simplified) LM model. */ /*! * \brief Add contribution from diffusion due to axisymmetric formulation to 2D residual @@ -195,7 +196,8 @@ class CSourceBase_TurbSA : public CNumerics { idx(nDim, config->GetnSpecies()), options(config->GetSAParsedOptions()), axisymmetric(config->GetAxisymmetric()), - transition_LM(config->GetKind_Trans_Model() == TURB_TRANS_MODEL::LM) { + transition_LM(config->GetKind_Trans_Model() == TURB_TRANS_MODEL::LM), + transition_SLM(transition_LM && config->GetLMParsedOptions().SLM) { /*--- Setup the Jacobian pointer, we need to return su2double** but we know * the Jacobian is 1x1 so we use this trick to avoid heap allocation. ---*/ /*--- Setup the Jacobian pointer (size increased for Stochastic Backscatter Model). ---*/ @@ -219,6 +221,7 @@ class CSourceBase_TurbSA : public CNumerics { AD::SetPreaccIn(PrimVar_Grad_i + idx.Velocity(), nDim, nDim); AD::SetPreaccIn(ScalarVar_Grad_i, nVar, nDim); AD::SetPreaccIn(stochSource, 3); + if (transition_LM) AD::SetPreaccIn(intermittency_i, intermittency_eff_i); /*--- Common auxiliary variables and constants of the model. ---*/ CSAVariables var; @@ -323,6 +326,18 @@ class CSourceBase_TurbSA : public CNumerics { var.intermittency = intermittency_eff_i; var.interDestrFactor = 1; + } else if (transition_SLM) { + + /*--- Lee and Baeder (AIAA 2021-1532), Eqs. 17-18: the scaled intermittency, which is zero in the laminar + * boundary layer, multiplies the production, and max(gamma_s, 0.1) the destruction. ---*/ + const su2double c_e2 = 50.0; + const su2double limitedIntermittency = min(intermittency_i, 1.0); + const su2double scaledIntermittency = (limitedIntermittency - 1.0 / c_e2) / (1.0 - 1.0 / c_e2); + const su2double cappedIntermittency = min(scaledIntermittency, 1.0); + const su2double gamma_s = max(cappedIntermittency, 0.0); + var.intermittency = gamma_s; + var.interDestrFactor = max(gamma_s, 0.1); + } else if (transition_LM){ var.intermittency = intermittency_eff_i; @@ -1001,6 +1016,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 7338679c1125..a2b5fffc0111 100644 --- a/SU2_CFD/include/solvers/CTransLMSolver.hpp +++ b/SU2_CFD/include/solvers/CTransLMSolver.hpp @@ -42,9 +42,13 @@ class CTransLMSolver final : public CTurbSolver { LM_ParsedOptions options; TURB_FAMILY TurbFamily; + bool isSepNeeded; TransLMCorrelations TransCorrelations; + /*! \brief Compute separation-induced and effective intermittency after the solution update. */ + void SetSeparationIntermittency(CGeometry* geometry, CSolver** solver_container, const CConfig* config); + /*! * \brief Resolve the compile-time parameters of CScalarFlux_TransLM and run one of this solver's * boundaries through the shared boundary flux pass. diff --git a/SU2_CFD/include/transition_data.hpp b/SU2_CFD/include/transition_data.hpp new file mode 100644 index 000000000000..bef7dd9466d6 --- /dev/null +++ b/SU2_CFD/include/transition_data.hpp @@ -0,0 +1,45 @@ +/*! + * \file transition_data.hpp + * \brief Diagnostic values shared by the LM transition numerics and solver. + * \version 8.5.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#pragma once + +#include "../../Common/include/basic_types/datatype_structure.hpp" + +/*! \brief Diagnostic values of the two-equation and simplified LM transition models. */ +struct TransitionLMData { + su2double criticalReynolds = 1.0; + su2double momentumThicknessReynolds = 0.0; + su2double turbulenceIntensity = 0.0; + su2double pressureGradient = 0.0; + su2double streamwiseVelocityGradient = 0.0; + su2double vorticityReynolds = 0.0; + su2double production = 0.0; + su2double destruction = 0.0; + su2double onset1 = 0.0; + su2double onset2 = 0.0; + su2double onset3 = 0.0; + su2double onset = 0.0; +}; diff --git a/SU2_CFD/include/variables/CTransLMVariable.hpp b/SU2_CFD/include/variables/CTransLMVariable.hpp index f2877c1cb155..e8f5dba9a05c 100644 --- a/SU2_CFD/include/variables/CTransLMVariable.hpp +++ b/SU2_CFD/include/variables/CTransLMVariable.hpp @@ -40,6 +40,9 @@ class CTransLMVariable final : public CTurbVariable { protected: VectorType Intermittency_Eff; VectorType Intermittency_Sep; + + std::vector TransitionData; + MatrixType WallNormal; public: /*! @@ -80,4 +83,8 @@ class CTransLMVariable final : public CTurbVariable { */ inline su2double GetIntermittencySep(unsigned long iPoint) const override { return Intermittency_Sep(iPoint); } + TransitionLMData* GetTransitionData(unsigned long iPoint) override { return &TransitionData[iPoint]; } + const TransitionLMData* GetTransitionData(unsigned long iPoint) const override { return &TransitionData[iPoint]; } + su2double* GetTransitionWallNormal(unsigned long iPoint) override { return WallNormal[iPoint]; } + const su2double* GetTransitionWallNormal(unsigned long iPoint) const override { return WallNormal[iPoint]; } }; diff --git a/SU2_CFD/include/variables/CTurbSSTVariable.hpp b/SU2_CFD/include/variables/CTurbSSTVariable.hpp index 0ec4abe466e9..9ecd804ba932 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 739a1af75fa5..292cccd555b5 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -36,6 +36,7 @@ #include #include "../../../Common/include/CConfig.hpp" +#include "../transition_data.hpp" #include "../../../Common/include/containers/container_decorators.hpp" class CFluidModel; @@ -1767,6 +1768,14 @@ class CVariable { */ inline virtual void SetIntermittencyEff(unsigned long iPoint, su2double val_Intermittency_eff) {} + /*! \brief Access the transition-model diagnostics at a point. */ + virtual TransitionLMData* GetTransitionData(unsigned long iPoint) { return nullptr; } + virtual const TransitionLMData* GetTransitionData(unsigned long iPoint) const { return nullptr; } + + /*! \brief Normal of the nearest wall element, stored by the simplified LM model. */ + virtual su2double* GetTransitionWallNormal(unsigned long iPoint) { return nullptr; } + virtual const su2double* GetTransitionWallNormal(unsigned long iPoint) const { return nullptr; } + /*! * \brief Set the value of the eddy viscosity. * \param[in] val_muT diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index 29fb325fa233..0d7b9f248668 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -1324,6 +1324,8 @@ void CDriver::InstantiateTransitionNumerics(unsigned short nVar_Trans, int offse const int source_second_term = SOURCE_SECOND_TERM + offset; const bool LM = config->GetKind_Trans_Model() == TURB_TRANS_MODEL::LM; + LM_ParsedOptions options; + if(LM) options = config->GetLMParsedOptions(); /*--- LM drives its interior loop and boundaries through its own CScalarFlux_TransLM edge kernel * (see CTransLMSolver), so conv_term/visc_term/conv_bound_term/visc_bound_term are never set @@ -1345,7 +1347,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); } diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp index 2ac2693b9952..e805e80ec5c9 100644 --- a/SU2_CFD/src/output/CFlowOutput.cpp +++ b/SU2_CFD/src/output/CFlowOutput.cpp @@ -1055,7 +1055,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; @@ -1122,7 +1125,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: @@ -1187,7 +1193,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; @@ -1292,12 +1301,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())); @@ -1378,7 +1391,29 @@ void CFlowOutput::SetVolumeOutputFieldsScalarSolution(const CConfig* config){ switch (config->GetKind_Trans_Model()) { case TURB_TRANS_MODEL::LM: AddVolumeOutput("INTERMITTENCY", "LM_gamma", "SOLUTION", "LM intermittency"); - AddVolumeOutput("RE_THETA_T", "LM_Re_t", "SOLUTION", "LM RE_THETA_T"); + /*--- Re_theta_t is a solution variable only of the two-equation model, the simplified model computes it + * algebraically, so it must not enter the restart files. ---*/ + if (!(config->GetLMParsedOptions()).SLM) { + AddVolumeOutput("RE_THETA_T", "LM_Re_t", "SOLUTION", "LM RE_THETA_T"); + } else { + AddVolumeOutput("RE_THETA_T", "LM_Re_t", "PRIMITIVE", "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("LAMBDA_THETA", "Lambda_theta", "DEBUG", "LM Lambda_theta"); + AddVolumeOutput("DU_DS", "du_ds", "DEBUG", "LM du_ds"); + if ((config->GetLMParsedOptions()).SLM) { + AddVolumeOutput("TU", "Tu", "PRIMITIVE", "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"); + } break; case TURB_TRANS_MODEL::NONE: @@ -1463,7 +1498,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: @@ -1695,15 +1732,37 @@ void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* con } switch (config->GetKind_Trans_Model()) { - case TURB_TRANS_MODEL::LM: + case TURB_TRANS_MODEL::LM: { + const auto& data = *Node_Trans->GetTransitionData(iPoint); SetVolumeOutputValue("INTERMITTENCY", iPoint, Node_Trans->GetSolution(iPoint, 0)); - SetVolumeOutputValue("RE_THETA_T", iPoint, Node_Trans->GetSolution(iPoint, 1)); + SetVolumeOutputValue("RE_V", iPoint, data.vorticityReynolds); + SetVolumeOutputValue("RE_THETA_CORR", iPoint, data.criticalReynolds); + SetVolumeOutputValue("PROD", iPoint, data.production); + SetVolumeOutputValue("DESTR", iPoint, data.destruction); + SetVolumeOutputValue("F_ONSET1", iPoint, data.onset1); + SetVolumeOutputValue("F_ONSET2", iPoint, data.onset2); + SetVolumeOutputValue("F_ONSET3", iPoint, data.onset3); + SetVolumeOutputValue("F_ONSET", iPoint, data.onset); + SetVolumeOutputValue("LAMBDA_THETA", iPoint, data.pressureGradient); + SetVolumeOutputValue("DU_DS", iPoint, data.streamwiseVelocityGradient); + if (!(config->GetLMParsedOptions()).SLM) { + SetVolumeOutputValue("RE_THETA_T", iPoint, Node_Trans->GetSolution(iPoint, 1)); + } else { + SetVolumeOutputValue("RE_THETA_T", iPoint, data.momentumThicknessReynolds); + SetVolumeOutputValue("TU", iPoint, data.turbulenceIntensity); + SetVolumeOutputValue("NORMAL_X", iPoint, Node_Trans->GetTransitionWallNormal(iPoint)[0]); + SetVolumeOutputValue("NORMAL_Y", iPoint, Node_Trans->GetTransitionWallNormal(iPoint)[1]); + SetVolumeOutputValue("NORMAL_Z", iPoint, Node_Trans->GetTransitionWallNormal(iPoint)[2]); + } 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; } @@ -2835,7 +2894,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/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 32c265e55656..e3312e5d9c54 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -26,6 +26,7 @@ */ #include "../../include/solvers/CEulerSolver.hpp" +#include "../../include/numerics/turbulent/transition/trans_correlations.hpp" #include "../../include/variables/CNSVariable.hpp" #include "../../../Common/include/toolboxes/geometry_toolbox.hpp" #include "../../../Common/include/toolboxes/printing_toolbox.hpp" @@ -1097,18 +1098,7 @@ void CEulerSolver::SetNondimensionalization(CConfig *config, unsigned short iMes Omega_FreeStreamND = Density_FreeStreamND*Tke_FreeStreamND/max(Viscosity_FreeStreamND*config->GetTurb2LamViscRatio_FreeStream(), EPS); config->SetOmega_FreeStreamND(Omega_FreeStreamND); - if (config->GetTurbulenceIntensity_FreeStream() *100 <= 1.3) { - if (config->GetTurbulenceIntensity_FreeStream() *100 >=0.027) { - Re_ThetaT_FreeStream = (1173.51-589.428*config->GetTurbulenceIntensity_FreeStream() *100+0.2196/ - (config->GetTurbulenceIntensity_FreeStream() *100*config->GetTurbulenceIntensity_FreeStream() *100)); - } - else { - Re_ThetaT_FreeStream = (1173.51-589.428*config->GetTurbulenceIntensity_FreeStream() *100+0.2196/(0.27*0.27)); - } - } - else { - Re_ThetaT_FreeStream = 331.5*pow(config->GetTurbulenceIntensity_FreeStream() *100-0.5658,-0.671); - } + Re_ThetaT_FreeStream = TransLMCorrelations::FreestreamReThetaT(config->GetTurbulenceIntensity_FreeStream() * 100.0); config->SetReThetaT_FreeStream(Re_ThetaT_FreeStream); const su2double MassDiffusivityND = config->GetDiffusivity_Constant() / (Velocity_Ref * Length_Ref); @@ -9180,7 +9170,7 @@ void CEulerSolver::PreprocessAverage(CSolver **solver, CGeometry *geometry, CCon Allreduce_inplace(nDim, TotalAreaVelocity); delete [] buffer; - + #endif /*--- initialize spanwise average quantities ---*/ @@ -9929,4 +9919,4 @@ void CEulerSolver::ComputeTurboBladePerformance(CGeometry* geometry, CConfig* co } TurbomachineryPerformance->ComputeTurbomachineryPerformance(bladePrimitives, iBlade); } -} \ No newline at end of file +} diff --git a/SU2_CFD/src/solvers/CTransLMSolver.cpp b/SU2_CFD/src/solvers/CTransLMSolver.cpp index d7404faa95e5..141dc538162e 100644 --- a/SU2_CFD/src/solvers/CTransLMSolver.cpp +++ b/SU2_CFD/src/solvers/CTransLMSolver.cpp @@ -53,6 +53,14 @@ CTransLMSolver::CTransLMSolver(CGeometry *geometry, CConfig *config, const CSolv /*--- 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,19 @@ CTransLMSolver::CTransLMSolver(CGeometry *geometry, CConfig *config, const CSolv nVarGrad = nVar; - /*--- Define variables needed for transition from config file */ - options = config->GetLMParsedOptions(); + /*--- Define geometry constants in the solver structure ---*/ + + nDim = geometry->GetnDim(); + + /*--- 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, and every level of a Full-MG startup, which solves the transition * equations on whichever grid the flow solver is currently on. ---*/ @@ -104,33 +120,27 @@ CTransLMSolver::CTransLMSolver(CGeometry *geometry, CConfig *config, const CSolv 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; - su2double ReThetaT_Inf = 100.0; + Solution_Inf[0] = Intermittency_Inf; - /*--- Momentum thickness Reynolds number, initialized from freestream turbulent intensity*/ - if (Intensity <= 1.3) { - if(Intensity >=0.027) { - ReThetaT_Inf = (1173.51-589.428*Intensity+0.2196/(Intensity*Intensity)); - } - else { - ReThetaT_Inf = (1173.51-589.428*Intensity+0.2196/(0.27*0.27)); - } - } - else if(Intensity>1.3) { - ReThetaT_Inf = 331.5*pow(Intensity-0.5658,-0.671); - } + const su2double ReThetaT_Inf = + TransLMCorrelations::FreestreamReThetaT(config->GetTurbulenceIntensity_FreeStream() * 100.0); 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(); /*--- Ghost states for boundary conditions, sized to the largest marker (see BoundaryFluxResidual). ---*/ @@ -165,7 +175,7 @@ CTransLMSolver::CTransLMSolver(CGeometry *geometry, CConfig *config, const CSolv 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; } } @@ -189,6 +199,57 @@ 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) { + + /*--- With cross-flow and SA, the auxiliary variables 1-3 are the direction of the vorticity, e_omega, of which + * the gradient gives the cross-flow strength of Lee and Baeder (AIAA 2021-1532), Eqs. 31-33. ---*/ + const bool crossFlowSA = options.CrossFlow && TurbFamily == TURB_FAMILY::SA; + + 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)); + for (auto iDim = 0u; iDim < nDim; ++iDim) nodes->GetTransitionWallNormal(iPoint)[iDim] = Normal[iDim]; + if (crossFlowSA) { + const auto Vorticity = flowNodes->GetVorticity(iPoint); + const su2double VorticityMag = GeometryToolbox::Norm(3, Vorticity); + for (auto iDim = 0u; iDim < 3; iDim++) + nodes->SetAuxVar(iPoint, 1 + iDim, VorticityMag > 1e-12 ? su2double(Vorticity[iDim] / VorticityMag) : su2double(0.0)); + } + } + 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); + } + + /*--- After the gradients, auxiliary variable 0 becomes dV/dy = grad(n . V) . n and, with cross-flow and SA, + * auxiliary variable 1 becomes Psi = |n . grad(e_omega)| d_w (Lee and Baeder, Eqs. 32-33). ---*/ + 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); + + if (crossFlowSA) { + su2double phi[3] = {0.0, 0.0, 0.0}; + for (auto iVar = 0u; iVar < 3; iVar++) + for (auto iDim = 0u; iDim < nDim; iDim++) + phi[iVar] += Normal[iDim] * nodes->GetAuxVarGradient(iPoint, 1 + iVar, iDim); + nodes->SetAuxVar(iPoint, 1, GeometryToolbox::Norm(3, phi) * geometry->nodes->GetWall_Distance(iPoint)); + } + } + END_SU2_OMP_FOR + + } + } void CTransLMSolver::Postprocessing(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iMesh) { @@ -204,6 +265,25 @@ void CTransLMSolver::Postprocessing(CGeometry *geometry, CSolver **solver_contai } AD::StartNoSharedReading(); + + if (isSepNeeded) { + SetSeparationIntermittency(geometry, solver_container, config); + } + 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(); +} + +void CTransLMSolver::SetSeparationIntermittency(CGeometry* geometry, CSolver** solver_container, const CConfig* config) { + auto* flowNodes = su2staticcast_p(solver_container[FLOW_SOL]->GetNodes()); auto* turbNodes = su2staticcast_p(solver_container[TURB_SOL]->GetNodes()); @@ -221,7 +301,6 @@ void CTransLMSolver::Postprocessing(CGeometry *geometry, CSolver **solver_contai 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); @@ -233,13 +312,23 @@ void CTransLMSolver::Postprocessing(CGeometry *geometry, CSolver **solver_contai 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 Re_t = 0.0; + su2double Corr_Rec = 0.0; + + if (options.SLM) { + Re_t = nodes->GetTransitionData(iPoint)->momentumThicknessReynolds; + Corr_Rec = nodes->GetTransitionData(iPoint)->criticalReynolds; + } 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); + } su2double R_t = 1.0; if(TurbFamily == TURB_FAMILY::KW) @@ -248,6 +337,7 @@ void CTransLMSolver::Postprocessing(CGeometry *geometry, CSolver **solver_contai 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; @@ -271,10 +361,8 @@ void CTransLMSolver::Postprocessing(CGeometry *geometry, CSolver **solver_contai } END_SU2_OMP_FOR - AD::EndNoSharedReading(); } - void CTransLMSolver::Upwind_Residual(CGeometry* geometry, CSolver** solver_container, CNumerics**, CConfig* config, unsigned short iMesh) { SU2_ZONE_SCOPED @@ -283,14 +371,14 @@ void CTransLMSolver::Upwind_Residual(CGeometry* geometry, CSolver** solver_conta * so there is no accurate-Jacobian correction to apply. ---*/ const auto opt = ScalarFluxOptions::Interior(*config, config->GetBounded_Turb()); - DispatchScheme(config, [&](auto tag) { + DispatchScheme(config, [&](auto tag) { EdgeFluxResidual(geometry, solver_container, config, opt); }); } void CTransLMSolver::BoundaryFlux(CGeometry* geometry, CSolver** solver_container, CConfig* config, const ScalarFluxOptions& opt, unsigned short val_marker) { - DispatchScheme(config, [&](auto tag) { + DispatchScheme(config, [&](auto tag) { BoundaryFluxResidual(geometry, solver_container, config, opt, val_marker); }); } @@ -315,7 +403,6 @@ void CTransLMSolver::Source_Residual(CGeometry *geometry, CSolver **solver_conta for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) { - /*--- Conservative variables w/o reconstruction ---*/ numerics->SetPrimitive(flowNodes->GetPrimitive(iPoint), nullptr); @@ -352,15 +439,22 @@ 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)); + if (options.CrossFlow && TurbFamily == TURB_FAMILY::SA) numerics->SetCrossFlowStrength(nodes->GetAuxVar(iPoint, 1)); + } + /*--- Compute the source term ---*/ auto residual = numerics->ComputeResidual(config); + *nodes->GetTransitionData(iPoint) = *numerics->GetTransitionData(); + /*--- Subtract residual and the Jacobian ---*/ LinSysRes.SubtractBlock(iPoint, residual); @@ -424,16 +518,10 @@ void CTransLMSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, C for (auto iVertex = 0u; iVertex < geometry->nVertex[val_marker]; iVertex++) { const auto* V_inlet = flowSolver->GetCharacPrimVar(val_marker, iVertex); - /*--- Non-dimensionalize Inlet_TurbVars if Inlet-Files are used. ---*/ - su2double Inlet_Vars[MAXNVAR]; - Inlet_Vars[0] = Inlet_TurbVars[val_marker][iVertex][0]; - Inlet_Vars[1] = Inlet_TurbVars[val_marker][iVertex][1]; - if (config->GetInlet_Profile_From_File()) { - Inlet_Vars[0] /= pow(config->GetVelocity_Ref(), 2); - Inlet_Vars[1] *= config->GetViscosity_Ref() / (config->GetDensity_Ref() * pow(config->GetVelocity_Ref(), 2)); - } - - for (auto iVar = 0u; iVar < nVar; iVar++) ghostNodes->SetSolution(iVertex, iVar, Inlet_Vars[iVar]); + /*--- The transition variables are dimensionless and are not read from inlet profile files: the free-stream + values set in the constructor are used (nVar of them, one for the simplified model). ---*/ + for (auto iVar = 0u; iVar < nVar; iVar++) + ghostNodes->SetSolution(iVertex, iVar, Inlet_TurbVars[val_marker][iVertex][iVar]); SetGhostPrimitives(iVertex, V_inlet); @@ -487,13 +575,12 @@ void CTransLMSolver::BC_Fluid_Interface(CGeometry *geometry, CSolver **solver_co /*--- LM's diffusion coefficients read no auxiliary ghost field. ---*/ const auto fillGhostExtras = [](unsigned long, unsigned long) {}; - DispatchScheme(config, [&](auto tag) { + DispatchScheme(config, [&](auto tag) { FluidInterfaceFluxResidual(geometry, solver_container, config, optConv, optVisc, fillGhostExtras); }); } - void CTransLMSolver::LoadRestart(CGeometry** geometry, CSolver*** solver, CConfig* config, int val_iter, bool val_update_geo) { SU2_ZONE_SCOPED @@ -543,8 +630,8 @@ 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]); - nodes->SetIntermittencySep(iPoint_Local, Restart_Data[index + 2]); - nodes->SetIntermittencyEff(iPoint_Local, Restart_Data[index + 3]); + /*--- The separation and effective intermittencies are not in the restart (compact restarts only have the + * solution variables), they are recomputed by Postprocessing below. ---*/ /*--- Increment the overall counter for how many points have been loaded. ---*/ counter++; @@ -600,4 +687,4 @@ void CTransLMSolver::LoadRestart(CGeometry** geometry, CSolver*** solver, CConfi } END_SU2_OMP_SAFE_GLOBAL_ACCESS -} \ No newline at end of file +} diff --git a/SU2_CFD/src/solvers/CTurbSASolver.cpp b/SU2_CFD/src/solvers/CTurbSASolver.cpp index ff580cf56b05..3124b974d3bb 100644 --- a/SU2_CFD/src/solvers/CTurbSASolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSASolver.cpp @@ -310,26 +310,16 @@ void CTurbSASolver::Postprocessing(CGeometry *geometry, CSolver **solver_contain if (geometry->nodes->GetDomain(iPoint)) { const auto jPoint = geometry->vertex[iMarker][iVertex]->GetNormal_Neighbor(); - su2double FrictionVelocity = 0.0; - /*--- Formulation varies for 2D and 3D problems: in 3D the friction velocity is assumed to be sqrt(mu * |Omega|) - (provided by the reference paper https://doi.org/10.2514/6.1992-439), whereas in 2D we have to use the - standard definition sqrt(c_f / rho) since Omega = 0. ---*/ - if(nDim == 2){ - su2double shearStress = 0.0; - for(auto iDim = 0u; iDim < nDim; iDim++) { - shearStress += pow(solver_container[FLOW_SOL]->GetCSkinFriction(iMarker, iVertex, iDim), 2.0); - } - shearStress = sqrt(shearStress); - - FrictionVelocity = sqrt(shearStress/flowNodes->GetDensity(iPoint)); - } else { - su2double VorticityMag = max(GeometryToolbox::Norm(3, flowNodes->GetVorticity(iPoint)), 1e-12); - FrictionVelocity = sqrt(flowNodes->GetLaminarViscosity(iPoint)*VorticityMag); - } + /*--- Friction velocity u_tau = sqrt(tau_w / rho) = sqrt(nu |Omega|) at the wall, where the wall shear stress + is mu |Omega| (https://doi.org/10.2514/6.1992-439). The vorticity is also defined in 2D (z component). ---*/ + const su2double VorticityMag = max(GeometryToolbox::Norm(3, flowNodes->GetVorticity(iPoint)), 1e-12); + const su2double kinematicViscosity = + flowNodes->GetLaminarViscosity(iPoint) / max(flowNodes->GetDensity(iPoint), 1e-20); + const su2double FrictionVelocity = sqrt(kinematicViscosity * VorticityMag); 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/variables/CTransLMVariable.cpp b/SU2_CFD/src/variables/CTransLMVariable.cpp index 9ab1e10bec80..dfa0969a6c07 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; iPointGetKind_Turb_Model()) == TURB_FAMILY::SA; + nAuxVar = crossFlowSA ? 4 : 1; + Grad_AuxVar.resize(nPoint, nAuxVar, nDim, su2double(0.0)); + AuxVar.resize(nPoint, nAuxVar) = su2double(0.0); + WallNormal.resize(nPoint, 3) = su2double(0.0); + } + } void CTransLMVariable::SetIntermittencyEff(unsigned long iPoint, su2double val_Intermittency_sep) { @@ -54,4 +75,4 @@ void CTransLMVariable::SetIntermittencyEff(unsigned long iPoint, su2double val_I void CTransLMVariable::SetIntermittencySep(unsigned long iPoint, su2double val_Intermittency_sep) { Intermittency_Sep(iPoint) = val_Intermittency_sep; -} \ No newline at end of file +} diff --git a/SU2_CFD/src/variables/CTurbSSTVariable.cpp b/SU2_CFD/src/variables/CTurbSSTVariable.cpp index a0a08669c766..bf04228e3a05 100644 --- a/SU2_CFD/src/variables/CTurbSSTVariable.cpp +++ b/SU2_CFD/src/variables/CTurbSSTVariable.cpp @@ -33,6 +33,7 @@ CTurbSSTVariable::CTurbSSTVariable(su2double kine, su2double omega, su2double mu : CTurbVariable(npoint, ndim, nvar, config) { sstParsedOptions = config->GetSSTParsedOptions(); + lmParsedOptions = config->GetLMParsedOptions(); for(unsigned long iPoint=0; iPoint #include #include "../../../SU2_CFD/include/numerics/CNumerics.hpp" +#include "../../../SU2_CFD/include/numerics/turbulent/transition/trans_correlations.hpp" #include "../../../SU2_CFD/include/numerics/NEMO/NEMO_diffusion.hpp" TEST_CASE("NTS blending has a minimum of 0.05", "[Upwind/central blending]") { @@ -243,3 +244,20 @@ TEST_CASE("NEMO corrected viscous residual returns distinct i and j Jacobians", } } } + +TEST_CASE("LM freestream and simplified correlations retain their limits", "[transition]") { + // Values from the published correlations, including both branches and the low-Tu limit. + REQUIRE(TransLMCorrelations::FreestreamReThetaT(0.5) == Approx(879.6744)); + REQUIRE(TransLMCorrelations::FreestreamReThetaT(2.0) == Approx(260.25459846968556)); + REQUIRE(TransLMCorrelations::FreestreamReThetaT(0.0) == TransLMCorrelations::FreestreamReThetaT(0.027)); + + LM_ParsedOptions options{}; + options.Correlation_SLM = TURB_TRANS_CORRELATION_SLM::MENTER_SLM; + TransLMCorrelations correlations; + correlations.SetOptions(options); + REQUIRE(correlations.ReThetaC_Correlations_SLM(1.0, 0.0, 0.1, 1.0, 1.0) == Approx(467.87944117144235)); + REQUIRE(correlations.ReThetaC_Correlations_SLM(0.1, 0.0, 0.1, 1.0, 1.0, true) == Approx(1069.8733022265405)); + // Above the SA blend, the original SST coefficients are recovered. + REQUIRE(correlations.ReThetaC_Correlations_SLM(2.0, 0.0, 0.1, 1.0, 1.0, true) == + correlations.ReThetaC_Correlations_SLM(2.0, 0.0, 0.1, 1.0, 1.0)); +} diff --git a/config_template.cfg b/config_template.cfg index 5e611f22ed48..5818927cd73d 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -33,7 +33,12 @@ KIND_TRANS_MODEL= NONE % Value of RMS roughness for transition model HROUGHNESS= 1.0e-6 % -% Specify versions/correlations of the LM model (LM2015, MALAN, SULUKSNA, KRAUSE, KRAUSE_HYPER, MEDIDA, MEDIDA_BAEDER, MENTER_LANGTRY) +% Specify versions/correlations of the LM model: +% CROSSFLOW: cross-flow corrections (LM2015 is accepted as its old name). +% SLM: one-equation (simplified) model of Menter et al. 2015 instead of the two-equation model. +% Correlations of the two-equation model (one of): MALAN, SULUKSNA, KRAUSE, KRAUSE_HYPER, MEDIDA, MEDIDA_BAEDER, +% MENTER_LANGTRY (default MENTER_LANGTRY with SST, MALAN with SA). +% Correlations of the one-equation model (one of): MENTER_SLM (default), CODER_SLM, MOD_EPPLER_SLM. LM_OPTIONS= NONE % % Specify subgrid scale model(NONE, IMPLICIT_LES, SMAGORINSKY, WALE, VREMAN) diff --git a/meson_scripts/init.py b/meson_scripts/init.py index 7a9ce0f0e887..81267ebb3132 100755 --- a/meson_scripts/init.py +++ b/meson_scripts/init.py @@ -26,7 +26,7 @@ # You should have received a copy of the GNU Lesser General Public # License along with SU2. If not, see . -import sys, os, subprocess, urllib.request, zipfile, time +import sys, os, subprocess, urllib.parse, urllib.request, zipfile, time def remove_file(path, retries=3, sleep=0.1): @@ -309,7 +309,10 @@ def _extract_member(self, member, targetpath, pwd): if not os.path.exists(filepath) and not os.path.exists(alt_filepath): try: - urllib.request.urlretrieve(url, filename) + if urllib.parse.urlsplit(url).scheme != "https": + raise ValueError("Dependency download URLs must use HTTPS") + # B310 audited: URLs use HTTPS; urllib rejects file/custom redirects. + urllib.request.urlretrieve(url, filename) # nosec B310 except Exception as e: print(e) print("Download of module " + name + " failed.")