From e652a53383686d2ff91cd47e41575dba9126c8a2 Mon Sep 17 00:00:00 2001 From: "oosun93@pusan.ac.kr" Date: Mon, 6 Jan 2025 23:50:48 +0900 Subject: [PATCH 1/6] Primary SA_AFT_2019b push Primary back up --- Common/include/CConfig.hpp | 27 +- Common/include/option_structure.hpp | 60 ++ Common/src/CConfig.cpp | 23 + .../turbulent/transition/trans_convection.hpp | 7 + .../transition/trans_correlations.hpp | 180 +++++ .../turbulent/transition/trans_diffusion.hpp | 75 +++ .../turbulent/transition/trans_sources.hpp | 149 +++++ .../numerics/turbulent/turb_sources.hpp | 23 +- SU2_CFD/include/solvers/CSolver.hpp | 18 + SU2_CFD/include/solvers/CTransAFTSolver.hpp | 218 +++++++ .../include/variables/CTransAFTVariable.hpp | 146 +++++ SU2_CFD/include/variables/CVariable.hpp | 114 ++++ SU2_CFD/src/drivers/CDriver.cpp | 10 +- SU2_CFD/src/iteration/CFluidIteration.cpp | 2 +- SU2_CFD/src/meson.build | 2 + SU2_CFD/src/output/CFlowOutput.cpp | 88 +++ SU2_CFD/src/solvers/CEulerSolver.cpp | 20 +- SU2_CFD/src/solvers/CSolverFactory.cpp | 7 + SU2_CFD/src/solvers/CTransAFTSolver.cpp | 615 ++++++++++++++++++ SU2_CFD/src/solvers/CTurbSASolver.cpp | 4 +- SU2_CFD/src/variables/CTransAFTVariable.cpp | 85 +++ config_template.cfg | 7 +- 22 files changed, 1869 insertions(+), 11 deletions(-) create mode 100644 SU2_CFD/include/solvers/CTransAFTSolver.hpp create mode 100644 SU2_CFD/include/variables/CTransAFTVariable.hpp create mode 100644 SU2_CFD/src/solvers/CTransAFTSolver.cpp create mode 100644 SU2_CFD/src/variables/CTransAFTVariable.cpp diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index 0c32f044288..65b9fe09bad 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -611,6 +611,7 @@ class CConfig { TURB_SGS_MODEL Kind_SGS_Model; /*!< \brief LES SGS model definition. */ TURB_TRANS_MODEL Kind_Trans_Model; /*!< \brief Transition model definition. */ TURB_TRANS_CORRELATION Kind_Trans_Correlation; /*!< \brief Transition correlation model definition. */ + AFT_CORRELATION Kind_AFT_Correlation; /*!< \brief Transition correlation of AFT model definition. */ su2double hRoughness; /*!< \brief RMS roughness for Transition model. */ unsigned short Kind_ActDisk, Kind_Engine_Inflow, *Kind_Data_Riemann, @@ -734,10 +735,12 @@ class CConfig { string *Config_Filenames; /*!< \brief List of names for configuration files. */ SST_OPTIONS *SST_Options; /*!< \brief List of modifications/corrections/versions of SST turbulence model.*/ SA_OPTIONS *SA_Options; /*!< \brief List of modifications/corrections/versions of SA turbulence model.*/ - LM_OPTIONS *LM_Options; /*!< \brief List of modifications/corrections/versions of SA turbulence model.*/ + LM_OPTIONS *LM_Options; /*!< \brief List of modifications/corrections/versions of LM transition model.*/ + AFT_OPTIONS *AFT_Options; /*!< \brief List of modifications/corrections/versions of AFT transition model.*/ unsigned short nSST_Options; /*!< \brief Number of SST options specified. */ unsigned short nSA_Options; /*!< \brief Number of SA options specified. */ - unsigned short nLM_Options; /*!< \brief Number of SA options specified. */ + unsigned short nLM_Options; /*!< \brief Number of LM options specified. */ + unsigned short nAFT_Options; /*!< \brief Number of AFT options specified. */ WALL_FUNCTIONS *Kind_WallFunctions; /*!< \brief The kind of wall function to use for the corresponding markers. */ unsigned short **IntInfo_WallFunctions; /*!< \brief Additional integer information for the wall function markers. */ su2double **DoubleInfo_WallFunctions; /*!< \brief Additional double information for the wall function markers. */ @@ -885,6 +888,7 @@ class CConfig { Tke_FreeStream, /*!< \brief Total turbulent kinetic energy of the fluid. */ Intermittency_FreeStream, /*!< \brief Freestream intermittency (for sagt transition model) of the fluid. */ ReThetaT_FreeStream, /*!< \brief Freestream Transition Momentum Thickness Reynolds Number (for LM transition model) of the fluid. */ + N_Critcal, /*!< \brief Critical N-factor (for AFT model). */ NuFactor_FreeStream, /*!< \brief Ratio of turbulent to laminar viscosity. */ NuFactor_Engine, /*!< \brief Ratio of turbulent to laminar viscosity at the engine. */ KFactor_LowerLimit, /*!< \Non dimensional coefficient for lower limit of K in SST model. */ @@ -1178,6 +1182,7 @@ class CConfig { SST_ParsedOptions sstParsedOptions; /*!< \brief Additional parameters for the SST turbulence model. */ SA_ParsedOptions saParsedOptions; /*!< \brief Additional parameters for the SA turbulence model. */ LM_ParsedOptions lmParsedOptions; /*!< \brief Additional parameters for the LM transition model. */ + AFT_ParsedOptions aftParsedOptions; /*!< \brief Additional parameters for the AFT transition model. */ su2double uq_delta_b; /*!< \brief Parameter used to perturb eigenvalues of Reynolds Stress Matrix */ unsigned short eig_val_comp; /*!< \brief Parameter used to determine type of eigenvalue perturbation */ su2double uq_urlx; /*!< \brief Under-relaxation factor */ @@ -2007,6 +2012,12 @@ class CConfig { */ su2double GetReThetaT_FreeStream() const { return ReThetaT_FreeStream; } + /*! + * \brief Get the value of the critical N-factor. + * \return The critical N-factor. + */ + su2double GetN_Critical(void) const { return N_Critcal; } + /*! * \brief Get the value of the non-dimensionalized freestream turbulence intensity. * \return Non-dimensionalized freestream intensity. @@ -2757,6 +2768,12 @@ class CConfig { */ void SetReThetaT_FreeStream(su2double val_ReThetaT_freestream) { ReThetaT_FreeStream = val_ReThetaT_freestream; } + /*! + * \brief Set the freestream momentum thickness Reynolds number. + * \param[in] val_ReThetaT_freestream - Value of the freestream momentum thickness Reynolds number. + */ + void SetN_Crtical(su2double val_N_critcal) { N_Critcal = val_N_critcal; } + /*! * \brief Set the non-dimensional freestream energy. * \param[in] val_energy_freestreamnd - Value of the non-dimensional freestream energy. @@ -9920,4 +9937,10 @@ class CConfig { */ LM_ParsedOptions GetLMParsedOptions() const { return lmParsedOptions; } + /*! + * \brief Get parsed AFT option data structure. + * \return AFT option data structure. + */ + AFT_ParsedOptions GetAFTParsedOptions() const { return aftParsedOptions; } + }; diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index ddf570ab697..70592980ecc 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -1205,10 +1205,12 @@ inline SA_ParsedOptions ParseSAOptions(const SA_OPTIONS *SA_Options, unsigned sh enum class TURB_TRANS_MODEL { NONE, /*!< \brief No transition model. */ LM, /*!< \brief Kind of transition model (Langtry-Menter (LM) for SST and Spalart-Allmaras). */ + AFT, /*!< \brief Kind of transition model (Amplification Factor Transport model for Spalart-Allmaras). */ }; static const MapType Trans_Model_Map = { MakePair("NONE", TURB_TRANS_MODEL::NONE) MakePair("LM", TURB_TRANS_MODEL::LM) + MakePair("AFT", TURB_TRANS_MODEL::AFT) }; /*! @@ -1324,6 +1326,64 @@ inline LM_ParsedOptions ParseLMOptions(const LM_OPTIONS *LM_Options, unsigned sh return LMParsedOptions; } +/*! + * \brief AFT Options + */ +enum class AFT_OPTIONS { + NONE, /*!< \brief No option / default. */ + AFT2019b /*!< \brief using AFT2019b model. */ +}; + +static const MapType AFT_Options_Map = { + MakePair("NONE", AFT_OPTIONS::NONE) + MakePair("AFT2019b", AFT_OPTIONS::AFT2019b) +}; + +/*! + * \brief Types of transition correlations + */ +enum class AFT_CORRELATION { + NONE, /*!< \brief No option / default. */ + AFT2019b /*!< \brief Kind of transition correlation model (AFT2019b). */ +}; + +/*! + * \brief Structure containing parsed AFT options. + */ +struct AFT_ParsedOptions { + AFT_OPTIONS version = AFT_OPTIONS::NONE; /*!< \brief AFT base model. */ + AFT_CORRELATION Correlation = AFT_CORRELATION::NONE; +}; + +/*! + * \brief Function to parse AFT options. + * \param[in] AFT_Options - Selected AFT option from config. + * \param[in] nAFT_Options - Number of options selected. + * \param[in] rank - MPI rank. + * \return Struct with AFT options. + */ +inline AFT_ParsedOptions ParseAFTOptions(const AFT_OPTIONS *AFT_Options, unsigned short nAFT_Options, int rank) { + AFT_ParsedOptions AFTParsedOptions; + + auto IsPresent = [&](AFT_OPTIONS option) { + const auto aft_options_end = AFT_Options + nAFT_Options; + return std::find(AFT_Options, aft_options_end, option) != aft_options_end; + }; + + int NFoundCorrelations = 0; + if (IsPresent(AFT_OPTIONS::AFT2019b)) { + AFTParsedOptions.Correlation = AFT_CORRELATION::AFT2019b; + AFTParsedOptions.version = AFT_OPTIONS::AFT2019b; + NFoundCorrelations++; + } + + if (NFoundCorrelations > 1) { + SU2_MPI::Error("Two correlations selected for AFT_OPTIONS. Please choose only one.", CURRENT_FUNCTION); + } + + return AFTParsedOptions; +} + /*! * \brief types of species transport models */ diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index c5bd653b04a..37c00b0ab35 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -1123,6 +1123,8 @@ void CConfig::SetConfig_Options() { addEnumOption("KIND_TRANS_MODEL", Kind_Trans_Model, Trans_Model_Map, TURB_TRANS_MODEL::NONE); /*!\brief SST_OPTIONS \n DESCRIPTION: Specify LM transition model options/correlations. \n Options: see \link LM_Options_Map \endlink \n DEFAULT: NONE \ingroup Config*/ addEnumListOption("LM_OPTIONS", nLM_Options, LM_Options, LM_Options_Map); + /*!\brief AFT_OPTIONS \n DESCRIPTION: Specify AFT transition model options/correlations. \n Options: see \link AFT_Options_Map \endlink \n DEFAULT: NONE \ingroup Config*/ + addEnumListOption("AFT_OPTIONS", nAFT_Options, AFT_Options, AFT_Options_Map); /*!\brief HROUGHNESS \n DESCRIPTION: Value of RMS roughness for transition model \n DEFAULT: 1E-6 \ingroup Config*/ addDoubleOption("HROUGHNESS", hRoughness, 1e-6); @@ -1419,6 +1421,8 @@ void CConfig::SetConfig_Options() { /* DESCRIPTION: */ addDoubleOption("FREESTREAM_TURBULENCEINTENSITY", TurbIntensityAndViscRatioFreeStream[0], 0.05); /* DESCRIPTION: */ + addDoubleOption("N_CRITICAL", N_Critcal, 0.0); + /* DESCRIPTION: */ addDoubleOption("FREESTREAM_NU_FACTOR", NuFactor_FreeStream, 3.0); /* DESCRIPTION: */ addDoubleOption("LOWER_LIMIT_K_FACTOR", KFactor_LowerLimit, 1.0e-15); @@ -3521,6 +3525,8 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i if (lmParsedOptions.LM2015 && val_nDim == 2) { SU2_MPI::Error("LM2015 is available only for 3D problems", CURRENT_FUNCTION); } + } else if (Kind_Trans_Model == TURB_TRANS_MODEL::AFT) { + aftParsedOptions = ParseAFTOptions(AFT_Options, nAFT_Options, rank); } /*--- Set the boolean Wall_Functions equal to true if there is a @@ -6304,6 +6310,10 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { } break; } + case TURB_TRANS_MODEL::AFT:{ + cout << "Transition model: Amplification Factor Transport model"; + break; + } } if (Kind_Trans_Model == TURB_TRANS_MODEL::LM) { @@ -6325,6 +6335,19 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { break; } } + if (Kind_Trans_Model == TURB_TRANS_MODEL::AFT) { + + switch (aftParsedOptions.Correlation) { + case AFT_CORRELATION::AFT2019b: + switch (Kind_Turb_Model) { + case TURB_MODEL::NONE: SU2_MPI::Error("No turbulence model has been selected but AFT transition model is active.", CURRENT_FUNCTION); break; + case TURB_MODEL::SST: SU2_MPI::Error("k-w SST turbulence model has been selected but AFT transition model is active.", CURRENT_FUNCTION); break; + } + cout << "-2019b" << endl; break; + if(!saParsedOptions.ft2) SU2_MPI::Error("ft2 option of SA model has been not selected.", CURRENT_FUNCTION); + case AFT_CORRELATION::NONE: SU2_MPI::Error("NONE has been selected.", CURRENT_FUNCTION); break; + } + } cout << "Hybrid RANS/LES: "; switch (Kind_HybridRANSLES) { case NO_HYBRIDRANSLES: cout << "No Hybrid RANS/LES" << endl; break; diff --git a/SU2_CFD/include/numerics/turbulent/transition/trans_convection.hpp b/SU2_CFD/include/numerics/turbulent/transition/trans_convection.hpp index 82912d00c74..27ab21469b2 100644 --- a/SU2_CFD/include/numerics/turbulent/transition/trans_convection.hpp +++ b/SU2_CFD/include/numerics/turbulent/transition/trans_convection.hpp @@ -38,3 +38,10 @@ template using CUpwSca_TransLM = CUpwSca_TurbSST; +/*! + * \class CUpwSca_TransAFT + * \brief Re-use the SST convective fluxes for the scalar upwind discretization of AFT transition model equations. + * \ingroup ConvDiscr + */ +template +using CUpwSca_TransAFT = CUpwSca_TurbSST; \ No newline at end of file diff --git a/SU2_CFD/include/numerics/turbulent/transition/trans_correlations.hpp b/SU2_CFD/include/numerics/turbulent/transition/trans_correlations.hpp index 617258174c9..af0a434c90d 100644 --- a/SU2_CFD/include/numerics/turbulent/transition/trans_correlations.hpp +++ b/SU2_CFD/include/numerics/turbulent/transition/trans_correlations.hpp @@ -190,3 +190,183 @@ class TransLMCorrelations { return F_length1; } }; + +/*! + * \class TransAFTCorrelations + * \brief Class for AFT model's correlation functions. + * \ingroup SourceDiscr + * \author S. Kang. + */ +class TransAFTCorrelations { + private: + + AFT_ParsedOptions options; + + public: + + /*! + * \brief Set AFT options. + * \param[in] val_options - AFT options structure. + */ + void SetOptions(const AFT_ParsedOptions val_options){ + options = val_options; + } + + /*! + * \brief Compute H12 from correlations. + * \param[in] HL - Local shape factor. + * \param[out] H12 - Integrated shape factor. + */ + su2double H12_Correlations(const su2double HL) const { + su2double H12 = 0.0; + + switch (options.Correlation) { + case AFT_CORRELATION::AFT2019b: { + H12 = 0.26 * HL + 2.4; + break; + } + + case AFT_CORRELATION::NONE: + SU2_MPI::Error("Transition correlation is set to DEFAULT but no default value has ben set in the code.", + CURRENT_FUNCTION); + break; + } + H12 = min(max(H12, 2.2), 20.0); + return H12; + } + + /*! + * \brief Compute dN/dRet from correlations. + * \param[in] H12 - Integreated shape factor. + * \param[out] dNdRet - N-factor gradient for 1st mode. + */ + su2double dNdRet_Correlations(const su2double H12) const { + su2double dNdRet = 0.0; + + switch (options.Correlation) { + case AFT_CORRELATION::AFT2019b: { + dNdRet = 0.028*(H12 - 1.0) - 0.0345 * exp(-pow(3.87/(H12 -1.0) - 2.52, 2.0)); + break; + } + + case AFT_CORRELATION::NONE: + SU2_MPI::Error("Transition correlation is set to DEFAULT but no default value has ben set in the code.", + CURRENT_FUNCTION); + break; + } + return dNdRet; + } + + /*! + * \brief Compute Ret0 from correlations. + * \param[in] H12 - Integreated Shape Factor. + * \param[out] Ret0 - N-factor growth onset momentumthickness Reynolds number. + */ + su2double Ret0_Correlations(const su2double H12) const { + su2double Ret0 = 0.0; + + switch (options.Correlation) { + case AFT_CORRELATION::AFT2019b: { + Ret0 = 0.7 * tanh( 14.0 / (H12 - 1.0) - 9.24) + 2.492 / pow(H12 - 1.0, 0.43) + 0.62; + Ret0 = pow(10, Ret0); + break; + } + + case AFT_CORRELATION::NONE: + SU2_MPI::Error("Transition correlation is set to DEFAULT but no default value has ben set in the code.", + CURRENT_FUNCTION); + break; + } + return Ret0; + } + + /*! + * \brief Compute D from correlations. + * \param[in] H12 - Integreated shape factor. + * \param[out] D - Model correlation. + */ + su2double D_Correlations(const su2double H12) const { + su2double D_H12 = 0.0; + + switch (options.Correlation) { + case AFT_CORRELATION::AFT2019b: { + D_H12 = 2.4 * H12 / (H12 - 1.0); + break; + } + + case AFT_CORRELATION::NONE: + SU2_MPI::Error("Transition correlation is set to DEFAULT but no default value has ben set in the code.", + CURRENT_FUNCTION); + break; + } + return D_H12; + } + + /*! + * \brief Compute l from correlations. + * \param[in] H12 - Integreated Shape Factor. + * \param[out] l(H12) - l function. + */ + su2double l_Correlations(const su2double H12) const { + su2double l_H12 = 0.0; + + switch (options.Correlation) { + case AFT_CORRELATION::AFT2019b: { + l_H12 = (6.54 * H12 - 14.07) / pow(H12, 2.0); + break; + } + + case AFT_CORRELATION::NONE: + SU2_MPI::Error("Transition correlation is set to DEFAULT but no default value has ben set in the code.", + CURRENT_FUNCTION); + break; + } + return l_H12; + } + + /*! + * \brief Compute m from correlations. + * \param[in] H12 - Integreated Shape Factor. + * \param[out] m(H12) - m function. + */ + su2double m_Correlations(const su2double H12, const su2double l_Correlation) const { + su2double m_H12 = 0.0; + + switch (options.Correlation) { + case AFT_CORRELATION::AFT2019b: { + m_H12 = (0.058 * pow(H12 - 4, 2.0) / (H12 - 1.0) - 0.068) / l_Correlation; + break; + } + + case AFT_CORRELATION::NONE: + SU2_MPI::Error("Transition correlation is set to DEFAULT but no default value has ben set in the code.", + CURRENT_FUNCTION); + break; + } + return m_H12; + } + + /*! + * \brief Compute kv from correlations. + * \param[in] H12 - Integreated Shape Factor. + * \param[out] kv - Kinetic Shape Factor. + */ + su2double kv_Correlations(const su2double H12) const { + su2double kv = 0.0; + + switch (options.Correlation) { + case AFT_CORRELATION::AFT2019b: { + kv = 1.0 / (0.4036 * pow(H12, 2.0) - 2.5394 * H12 + 4.3273); + break; + } + + case AFT_CORRELATION::NONE: + SU2_MPI::Error("Transition correlation is set to DEFAULT but no default value has ben set in the code.", + CURRENT_FUNCTION); + break; + } + return kv; + } + + +}; diff --git a/SU2_CFD/include/numerics/turbulent/transition/trans_diffusion.hpp b/SU2_CFD/include/numerics/turbulent/transition/trans_diffusion.hpp index 372d89b253c..31e42d1e828 100644 --- a/SU2_CFD/include/numerics/turbulent/transition/trans_diffusion.hpp +++ b/SU2_CFD/include/numerics/turbulent/transition/trans_diffusion.hpp @@ -103,3 +103,78 @@ class CAvgGrad_TransLM final : public CAvgGrad_Scalar { } }; + + +/*! + * \class CAvgGrad_TransAFT + * \brief Class for computing viscous term using average of gradient with correction (AFT transition model). + * \ingroup ViscDiscr + * \author S. Kang. + */ +template +class CAvgGrad_TransAFT final : public CAvgGrad_Scalar { +private: + using Base = CAvgGrad_Scalar; + using Base::Laminar_Viscosity_i; + using Base::Laminar_Viscosity_j; + using Base::Eddy_Viscosity_i; + using Base::Eddy_Viscosity_j; + using Base::Density_i; + using Base::Density_j; + using Base::ScalarVar_i; + using Base::ScalarVar_j; + using Base::Proj_Mean_GradScalarVar; + using Base::proj_vector_ij; + using Base::Flux; + using Base::Jacobian_i; + using Base::Jacobian_j; + + /*! + * \brief Adds any extra variables to AD + */ + void ExtraADPreaccIn() override {} + + /*! + * \brief LM transition model specific steps in the ComputeResidual method + * \param[in] config - Definition of the particular problem. + */ + void FinishResidualCalc(const CConfig* config) override { + const bool implicit = config->GetKind_TimeIntScheme() == EULER_IMPLICIT; + + /*--- Compute mean effective dynamic viscosity ---*/ + const su2double diff_i_AF = Laminar_Viscosity_i + Eddy_Viscosity_i; + const su2double diff_j_AF = Laminar_Viscosity_j + Eddy_Viscosity_j; + const su2double diff_i_Gamma = Laminar_Viscosity_i + Eddy_Viscosity_i; + const su2double diff_j_Gamma = Laminar_Viscosity_j + Eddy_Viscosity_j; + + const su2double diff_AF = 0.5*(diff_i_AF + diff_j_AF); + const su2double diff_Gamma = 0.5*(diff_i_Gamma + diff_j_Gamma); + + Flux[0] = diff_AF*Proj_Mean_GradScalarVar[0]; + Flux[1] = diff_Gamma*Proj_Mean_GradScalarVar[1]; + + /*--- For Jacobians -> Use of TSL (Thin Shear Layer) approx. to compute derivatives of the gradients ---*/ + if (implicit) { + const su2double proj_on_rho_i = proj_vector_ij/Density_i; + Jacobian_i[0][0] = -diff_AF*proj_on_rho_i; Jacobian_i[0][1] = 0.0; + Jacobian_i[1][0] = 0.0; Jacobian_i[1][1] = -diff_Gamma*proj_on_rho_i; + + const su2double proj_on_rho_j = proj_vector_ij/Density_j; + Jacobian_j[0][0] = diff_AF*proj_on_rho_j; Jacobian_j[0][1] = 0.0; + Jacobian_j[1][0] = 0.0; Jacobian_j[1][1] = diff_Gamma*proj_on_rho_j; + } + } + +public: + /*! + * \brief Constructor of the class. + * \param[in] val_nDim - Number of dimensions of the problem. + * \param[in] val_nVar - Number of variables of the problem. + * \param[in] correct_grad - Whether to correct gradient for skewness. + * \param[in] config - Definition of the particular problem. + */ + CAvgGrad_TransAFT(unsigned short val_nDim, unsigned short val_nVar, bool correct_grad, const CConfig* config) + : CAvgGrad_Scalar(val_nDim, val_nVar, correct_grad, config){ + } + +}; \ No newline at end of file diff --git a/SU2_CFD/include/numerics/turbulent/transition/trans_sources.hpp b/SU2_CFD/include/numerics/turbulent/transition/trans_sources.hpp index 4030e35bab9..0e652501f8d 100644 --- a/SU2_CFD/include/numerics/turbulent/transition/trans_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/transition/trans_sources.hpp @@ -323,3 +323,152 @@ class CSourcePieceWise_TransLM final : public CNumerics { return ResidualType<>(Residual, Jacobian_i, nullptr); } }; + +/*! + * \class CSourcePieceWise_TranAFT + * \brief Class for integrating the source terms of the AFT transition model equations. + * \ingroup SourceDiscr + * \author S. Kang. + */ +template +class CSourcePieceWise_TransAFT final : public CNumerics { + private: + const FlowIndices idx; /*!< \brief Object to manage the access to the flow primitives. */ + + const AFT_ParsedOptions options; + + /*--- AFT Closure constants ---*/ + const su2double c_1 = 100.0; + const su2double c_2 = 0.06; + const su2double c_3 = 50.0; + + TURB_FAMILY TurbFamily; + su2double Residual[2]; + su2double* Jacobian_i[2]; + su2double Jacobian_Buffer[4]; // Static storage for the Jacobian (which needs to be pointer for return type). + + TransAFTCorrelations TransCorrelations; + + public: + /*! + * \brief Constructor of the class. + * \param[in] val_nDim - Number of dimensions of the problem. + * \param[in] val_nVar - Number of variables of the problem. + * \param[in] config - Definition of the particular problem. + */ + CSourcePieceWise_TransAFT(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config) + : CNumerics(val_nDim, 2, config), idx(val_nDim, config->GetnSpecies()), options(config->GetAFTParsedOptions()){ + /*--- "Allocate" the Jacobian using the static buffer. ---*/ + Jacobian_i[0] = Jacobian_Buffer; + Jacobian_i[1] = Jacobian_Buffer + 2; + + 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[0] = w, TransVar[0] = Amplification Factor, and TransVar[1] = Gamma ---*/ + /*--- dU/dx = PrimVar_Grad[1][0] ---*/ + AD::StartPreacc(); + AD::SetPreaccIn(StrainMag_i); + AD::SetPreaccIn(ScalarVar_i, nVar); + AD::SetPreaccIn(ScalarVar_Grad_i, nVar, nDim); + AD::SetPreaccIn(TransVar_i, nVar); + AD::SetPreaccIn(TransVar_Grad_i, nVar, nDim); + AD::SetPreaccIn(Volume); + AD::SetPreaccIn(dist_i); + AD::SetPreaccIn(&V_i[idx.Velocity()], nDim); + AD::SetPreaccIn(PrimVar_Grad_i, nDim + idx.Velocity(), nDim); + AD::SetPreaccIn(Vorticity_i, 3); + AD::SetPreaccIn(AuxVar_Grad_i, 2); + + su2double VorticityMag = + sqrt(Vorticity_i[0] * Vorticity_i[0] + Vorticity_i[1] * Vorticity_i[1] + Vorticity_i[2] * Vorticity_i[2]); + + const su2double vel_u = V_i[idx.Velocity()]; + const su2double vel_v = V_i[1 + idx.Velocity()]; + const su2double vel_w = (nDim == 3) ? V_i[2 + idx.Velocity()] : 0.0; + + const su2double Velocity_Mag = sqrt(vel_u * vel_u + vel_v * vel_v + vel_w * vel_w); + const su2double Critical_N_Factor = config->GetN_Critical(); + + 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.0; + Residual[1] = 0.0; + Jacobian_i[0][0] = 0.0; + Jacobian_i[0][1] = 0.0; + Jacobian_i[1][0] = 0.0; + Jacobian_i[1][1] = 0.0; + + if (dist_i > 1e-10) { + su2double HLGradTerm = 0.0; + HLGradTerm = AuxVar_Grad_i[1][0] * AuxVar_Grad_i[0][0] + AuxVar_Grad_i[1][1] * AuxVar_Grad_i[0][1]; + if(nDim == 3) { + HLGradTerm += AuxVar_Grad_i[1][2] * AuxVar_Grad_i[0][2]; + } + if(HLGradTerm > 0) { + su2double tt = 0; + } + + const su2double HL = dist_i * dist_i / Laminar_Viscosity_i * HLGradTerm; + /*--- Cal H12, dNdRet, Ret0, D_H12, l_H12, m_H12, kv ---*/ + const su2double H12 = TransCorrelations.H12_Correlations(HL); + const su2double dNdRet = TransCorrelations.dNdRet_Correlations(H12); + const su2double Ret0 = TransCorrelations.Ret0_Correlations(H12); + const su2double D_H12 = TransCorrelations.D_Correlations(H12); + const su2double l_H12 = TransCorrelations.l_Correlations(H12); + const su2double m_H12 = TransCorrelations.m_Correlations(H12, l_H12); + const su2double kv = TransCorrelations.kv_Correlations(H12); + const su2double Rev = Density_i * dist_i * dist_i * StrainMag_i / (Laminar_Viscosity_i + Eddy_Viscosity_i); + const su2double Rev0 = kv * Ret0; + + if(H12 != 2.2 && H12 != 20.0 && dist_i < 0.01){ + su2double tt2 = 0; + } + + /*--- production term of the amplification factor ---*/ + const su2double F_growth = max(D_H12 * (1.0 + m_H12) / 2.0 * l_H12, 0.0); + const su2double F_crit = (Rev < Rev0) ? 0.0 : 1.0; + const su2double PAF = Density_i * VorticityMag * F_crit * F_growth * dNdRet; + + /*--- production term of the intermittency(Gamma) ---*/ + const su2double RT = Eddy_Viscosity_i / Laminar_Viscosity_i; + const su2double F_onset1 = TransVar_i[0] / Critical_N_Factor; + const su2double F_onset2 = min(F_onset1, 2.0); + const su2double F_onset3 = max(1.0 - pow(RT / 3.5, 3), 0.0); + const su2double F_onset = max(F_onset2 - F_onset3, 0.0); + const su2double F_turb = exp(-pow(RT / 2.0,4)); + const su2double Pg = c_1 * Density_i * StrainMag_i * F_onset * (1.0 - exp(TransVar_i[1])); + + /*-- destruction term of Intermeittency(Gamma) --*/ + const su2double Dg = c_2 * Density_i * VorticityMag * F_turb * (c_3 * exp(TransVar_i[1]) - 1.0); + + /*--- Source ---*/ + Residual[0] += (PAF) * Volume; + Residual[1] += (Pg - Dg) * Volume; + + /*--- Implicit part ---*/ + Jacobian_i[0][0] = 0.0; + Jacobian_i[0][1] = 0.0; + Jacobian_i[1][0] = 0.0; + Jacobian_i[1][1] = ( -c_1 * StrainMag_i * F_onset * exp(TransVar_i[1]) -c_2 * VorticityMag * F_turb * c_3 * exp(TransVar_i[1]) ) * Volume; + } + + AD::SetPreaccOut(Residual, nVar); + AD::EndPreacc(); + + return ResidualType<>(Residual, Jacobian_i, nullptr); + } +}; diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp index c9464b1dfc3..2bfbe830d87 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp @@ -58,6 +58,7 @@ struct CSAVariables { su2double Omega, dist_i_2, inv_k2_d2, inv_Shat, g_6, norm2_Grad; su2double intermittency, interDestrFactor; + su2double lnintermittecncy; }; /*! @@ -80,6 +81,7 @@ class CSourceBase_TurbSA : public CNumerics { const bool axisymmetric = false; bool transition_LM; + bool transition_AFT; /*! * \brief Add contribution from diffusion due to axisymmetric formulation to 2D residual @@ -123,7 +125,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_AFT(config->GetKind_Trans_Model() == TURB_TRANS_MODEL::AFT) { /*--- 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. ---*/ Jacobian_i = &Jacobian_Buffer; @@ -144,6 +147,7 @@ class CSourceBase_TurbSA : public CNumerics { AD::SetPreaccIn(Vorticity_i, 3); AD::SetPreaccIn(PrimVar_Grad_i + idx.Velocity(), nDim, nDim); AD::SetPreaccIn(ScalarVar_Grad_i[0], nDim); + AD::SetPreaccIn(intermittency_i); /*--- Common auxiliary variables and constants of the model. ---*/ CSAVariables var; @@ -192,7 +196,13 @@ class CSourceBase_TurbSA : public CNumerics { } /*--- Compute ft2 term ---*/ - ft2::get(var); + if(transition_AFT) { + var.lnintermittecncy = exp(intermittency_i); + } + else { + var.lnintermittecncy = 1.0; + } + ft2::get(var, transition_AFT); /*--- Compute auxiliary function r ---*/ rFunc::get(ScalarVar_i[0], var); @@ -307,12 +317,13 @@ struct Edw { * \brief SA classes to set the ft2 term and its derivative. * \ingroup SourceDiscr * \param[in,out] var: Common SA variables struct. + * \param[in] transition_AFT: Common AFT model struct. */ struct ft2 { /*! \brief No-ft2 term. */ struct Zero { - static void get(CSAVariables& var) { + static void get(CSAVariables& var, bool transition_AFT) { var.ft2 = 0.0; var.d_ft2 = 0.0; } @@ -320,10 +331,14 @@ struct Zero { /*! \brief Non-zero ft2 term according to the literature. */ struct Nonzero { - static void get(CSAVariables& var) { + static void get(CSAVariables& var, bool transition_AFT) { const su2double xsi2 = pow(var.Ji, 2); var.ft2 = var.ct3 * exp(-var.ct4 * xsi2); var.d_ft2 = -2.0 * var.ct4 * var.Ji * var.ft2 * var.d_Ji; + if(transition_AFT) { + var.ft2 = var.ct3 * exp(1.0 - var.lnintermittecncy); + var.d_ft2 = 0.0; + } } }; }; diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index b6c0d45d219..d968d05c6e3 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -3170,6 +3170,24 @@ class CSolver { */ inline virtual su2double GetReThetaT_Inf() const { return 0; } + /*! + * \brief Get value of the Amplification Factor. + * \return Value of the Amplification Factor. + */ + inline virtual su2double GetAF_Inf() const { return 0; } + + /*! + * \brief Get value of the natural logarithnm Intermittency. + * \return Value of the natural logarithnm Intermittency. + */ + inline virtual su2double GetLnIntermittency_Inf() const { return 0; } + + /*! + * \brief Get value of the critical N-factor. + * \return Value of the critical N-factor. + */ + inline virtual su2double N_Critical() const { return 0; } + /*! * \brief A virtual member. * \return Value of the sensitivity coefficient for the Young Modulus E diff --git a/SU2_CFD/include/solvers/CTransAFTSolver.hpp b/SU2_CFD/include/solvers/CTransAFTSolver.hpp new file mode 100644 index 00000000000..e9140a7abb0 --- /dev/null +++ b/SU2_CFD/include/solvers/CTransAFTSolver.hpp @@ -0,0 +1,218 @@ +/*! + * \file CTransLMSolver.hpp + * \brief Headers of the CTransLMSolver class + * \author A. Aranake + * \version 8.1.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2024, 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 "CTurbSolver.hpp" +#include "../numerics/turbulent/transition/trans_correlations.hpp" + +/*! + * \class CTransAFTSolver + * \brief Main class for defining the transition model solver. + * \ingroup Turbulence_Model + * \author S. Kang. + */ + +class CTransAFTSolver final : public CTurbSolver { +private: + + AFT_ParsedOptions options; + TURB_FAMILY TurbFamily; + + TransAFTCorrelations TransCorrelations; + +public: + /*! + * \overload + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + * \param[in] iMesh - Index of the mesh in multigrid computations. + */ + CTransAFTSolver(CGeometry *geometry, CConfig *config, unsigned short iMesh); + + /*! + * \brief Destructor of the class. + */ + ~CTransAFTSolver() = default; + + /*! + * \brief Restart residual and compute gradients. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] config - Definition of the particular problem. + * \param[in] iMesh - Index of the mesh in multigrid computations. + * \param[in] iRKStep - Current step of the Runge-Kutta iteration. + * \param[in] RunTime_EqSystem - System of equations which is going to be solved. + * \param[in] Output - boolean to determine whether to print output. + */ + void Preprocessing(CGeometry *geometry, + CSolver **solver_container, + CConfig *config, + unsigned short iMesh, + unsigned short iRKStep, + unsigned short RunTime_EqSystem, + bool Output) override; + + /*! + * \brief Computes the effective intermtittency. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] config - Definition of the particular problem. + * \param[in] iMesh - Index of the mesh in multigrid computations. + */ + void Postprocessing(CGeometry *geometry, + CSolver **solver_container, + CConfig *config, + unsigned short iMesh) override; + + /*! + * \brief Compute the viscous flux for the LM equation at a particular edge. + * \param[in] iEdge - Edge for which we want to compute the flux + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] numerics - Description of the numerical method. + * \param[in] config - Definition of the particular problem. + * \note Calls a generic implementation after defining a SolverSpecificNumerics object. + */ + void Viscous_Residual(const unsigned long iEdge, const CGeometry* geometry, CSolver** solver_container, + CNumerics* numerics, const CConfig* config) override; + + /*! + * \brief Source term computation. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] numerics_container - Description of the numerical method. + * \param[in] config - Definition of the particular problem. + * \param[in] iMesh - Index of the mesh in multigrid computations. + */ + void Source_Residual(CGeometry *geometry, + CSolver **solver_container, + CNumerics **numerics_container, + CConfig *config, + unsigned short iMesh) override; + + /*! + * \brief Source term computation. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] numerics - Description of the numerical method. + * \param[in] config - Definition of the particular problem. + * \param[in] iMesh - Index of the mesh in multigrid computations. + */ + void Source_Template(CGeometry *geometry, + CSolver **solver_container, + CNumerics *numerics, + CConfig *config, + unsigned short iMesh) override; + + /*! + * \brief Impose the Langtry Menter transition wall boundary condition. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] conv_numerics - Description of the numerical method. + * \param[in] visc_numerics - Description of the numerical method. + * \param[in] config - Definition of the particular problem. + * \param[in] val_marker - Surface marker where the boundary condition is applied. + */ + void BC_HeatFlux_Wall(CGeometry *geometry, + CSolver **solver_container, + CNumerics *conv_numerics, + CNumerics *visc_numerics, + CConfig *config, + unsigned short val_marker) override; + + /*! + * \brief Impose the Navier-Stokes wall boundary condition. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] conv_numerics - Description of the numerical method. + * \param[in] visc_numerics - Description of the numerical method. + * \param[in] config - Definition of the particular problem. + * \param[in] val_marker - Surface marker where the boundary condition is applied. + */ + void BC_Isothermal_Wall(CGeometry *geometry, + CSolver **solver_container, + CNumerics *conv_numerics, + CNumerics *visc_numerics, + CConfig *config, + unsigned short val_marker) override; + + /*! + * \brief Impose the inlet boundary condition. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] conv_numerics - Description of the numerical method. + * \param[in] visc_numerics - Description of the numerical method. + * \param[in] config - Definition of the particular problem. + * \param[in] val_marker - Surface marker where the boundary condition is applied. + */ + void BC_Inlet(CGeometry *geometry, + CSolver **solver_container, + CNumerics *conv_numerics, + CNumerics *visc_numerics, + CConfig *config, + unsigned short val_marker) override; + + /*! + * \brief Impose the outlet boundary condition. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] conv_numerics - Description of the numerical method. + * \param[in] visc_numerics - Description of the numerical method. + * \param[in] config - Definition of the particular problem. + * \param[in] val_marker - Surface marker where the boundary condition is applied. + */ + void BC_Outlet(CGeometry *geometry, + CSolver **solver_container, + CNumerics *conv_numerics, + CNumerics *visc_numerics, + CConfig *config, + unsigned short val_marker) override; + + /*! + * \brief Get the value of the intermittency. + * \return Value of the turbulent kinetic energy. + */ + inline su2double GetAF_Inf(void) const override { return Solution_Inf[0]; } + + /*! + * \brief Get the value of the intermittency. + * \return Value of the turbulent kinetic energy. + */ + inline su2double GetLnIntermittency_Inf(void) const override { return Solution_Inf[1]; } + + /*! + * \brief Load a solution from a restart file. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver - Container vector with all of the solvers. + * \param[in] config - Definition of the particular problem. + * \param[in] val_iter - Current external iteration number. + * \param[in] val_update_geo - Flag for updating coords and grid velocity. + */ + void LoadRestart(CGeometry** geometry, CSolver*** solver, CConfig* config, int val_iter, bool val_update_geo) final; + +}; diff --git a/SU2_CFD/include/variables/CTransAFTVariable.hpp b/SU2_CFD/include/variables/CTransAFTVariable.hpp new file mode 100644 index 00000000000..98daf2c3c32 --- /dev/null +++ b/SU2_CFD/include/variables/CTransAFTVariable.hpp @@ -0,0 +1,146 @@ +/*! + * \file CTransAFTVariable.hpp + * \brief Declaration of the variables of the transition model. + * \author S. Kang. + * \version 8.1.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2024, 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 "CTurbVariable.hpp" + +/*! + * \class CTransAFTVariable + * \brief Transition model variables. + * \ingroup Turbulence_Model + * \author S. Kang. + */ + +class CTransAFTVariable final : public CTurbVariable { +protected: + VectorType TempVar1, TempVar2, TempVar3, TempVar4, TempVar5, TempVar6; + VectorType TempVar7, TempVar8, TempVar9, TempVar10, TempVar11; + VectorType TempVar12, TempVar13, TempVar14, TempVar15; + +public: + /*! + * \brief Constructor of the class. + * \param[in] AF - Amplification Factor (initialization value). + * \param[in] LnIntermittency - Natural logarithm intermittency (initialization value). + * \param[in] npoint - Number of points/nodes/vertices in the domain. + * \param[in] ndim - Number of dimensions of the problem. + * \param[in] nvar - Number of variables of the problem. + * \param[in] config - Definition of the particular problem. + */ + CTransAFTVariable(su2double AF, su2double LnIntermittency, unsigned long npoint, unsigned long ndim, unsigned long nvar, CConfig *config); + + /*! + * \brief Destructor of the class. + */ + ~CTransAFTVariable() override = default; + + void SetAFT_Wonder_Func(unsigned long iPoint, su2double tempVar1, su2double tempVar2, su2double tempVar3, su2double tempVar4, su2double tempVar5, + su2double tempVar6, su2double tempVar7, su2double tempVar8, su2double tempVar9, su2double tempVar10, + su2double tempVar11, su2double tempVar12, su2double tempVar13, su2double tempVar14, su2double tempVar15) override; + + /*! + * \brief Value of Wonder Variable. + */ + inline su2double GetAFT_Wonder_Func_var1(unsigned long iPoint) const override { return TempVar1(iPoint); } + + /*! + * \brief Value of Wonder Variable. + */ + inline su2double GetAFT_Wonder_Func_var2(unsigned long iPoint) const override { return TempVar2(iPoint); } + + /*! + * \brief Value of Wonder Variable. + */ + inline su2double GetAFT_Wonder_Func_var3(unsigned long iPoint) const override { return TempVar3(iPoint); } + + /*! + * \brief Value of Wonder Variable. + */ + inline su2double GetAFT_Wonder_Func_var4(unsigned long iPoint) const override { return TempVar4(iPoint); } + + /*! + * \brief Value of Wonder Variable. + */ + inline su2double GetAFT_Wonder_Func_var5(unsigned long iPoint) const override { return TempVar5(iPoint); } + + /*! + * \brief Value of Wonder Variable. + */ + inline su2double GetAFT_Wonder_Func_var6(unsigned long iPoint) const override { return TempVar6(iPoint); } + + /*! + * \brief Value of Wonder Variable. + */ + inline su2double GetAFT_Wonder_Func_var7(unsigned long iPoint) const override { return TempVar7(iPoint); } + + /*! + * \brief Value of Wonder Variable. + */ + inline su2double GetAFT_Wonder_Func_var8(unsigned long iPoint) const override { return TempVar8(iPoint); } + + /*! + * \brief Value of Wonder Variable. + */ + inline su2double GetAFT_Wonder_Func_var9(unsigned long iPoint) const override { return TempVar9(iPoint); } + + /*! + * \brief Value of Wonder Variable. + */ + inline su2double GetAFT_Wonder_Func_var10(unsigned long iPoint) const override { return TempVar10(iPoint); } + + /*! + * \brief Value of Wonder Variable. + */ + inline su2double GetAFT_Wonder_Func_var11(unsigned long iPoint) const override { return TempVar11(iPoint); } + + + /*! + * \brief Value of Wonder Variable. + */ + inline su2double GetAFT_Wonder_Func_var12(unsigned long iPoint) const override { return TempVar12(iPoint); } + + + /*! + * \brief Value of Wonder Variable. + */ + inline su2double GetAFT_Wonder_Func_var13(unsigned long iPoint) const override { return TempVar13(iPoint); } + + + /*! + * \brief Value of Wonder Variable. + */ + inline su2double GetAFT_Wonder_Func_var14(unsigned long iPoint) const override { return TempVar14(iPoint); } + + + /*! + * \brief Value of Wonder Variable. + */ + inline su2double GetAFT_Wonder_Func_var15(unsigned long iPoint) const override { return TempVar15(iPoint); } + + +}; diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index bc855bc5dca..a64075a4687 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -1675,6 +1675,120 @@ class CVariable { */ inline virtual void SetIntermittency(unsigned long iPoint, su2double val_Intermittency) {} + /*! + * \brief Set the variables for the AFT model. + * \param[in] val_iPoint - Value of the iPoint. + * \param[in] tempvar1 - tempvariable. + * \param[in] tempvar2 - tempvariable. + * \param[in] tempvar3 - tempvariable. + * \param[in] tempvar4 - tempvariable. + * \param[in] tempvar5 - tempvariable. + * \param[in] tempvar6 - tempvariable. + * \param[in] tempvar7 - tempvariable. + * \param[in] tempvar8 - tempvariable. + * \param[in] tempvar9 - tempvariable. + * \param[in] tempvar10 - tempvariable. + * \param[in] tempvar11 - tempvariable. + */ + inline virtual void SetAFT_Wonder_Func(unsigned long iPoint, su2double tempVar1, su2double tempVar2, su2double tempVar3, + su2double tempVar4, su2double tempVar5, su2double tempVar6, su2double tempVar7, + su2double tempVar8, su2double tempVar9, su2double tempVar10, su2double tempVar11, + su2double tempVar12, su2double tempVar13, su2double tempVar14, su2double tempVar15) {} + + /*! + * \brief Get the value of the Length scale. + * \return the value of the Length scale. + */ + inline virtual su2double GetAFT_Wonder_Func_var1(unsigned long iPoint) const {return 0.0;} + + /*! + * \brief Get the value of the Length scale. + * \return the value of the Length scale. + */ + inline virtual su2double GetAFT_Wonder_Func_var2(unsigned long iPoint) const {return 0.0;} + + /*! + * \brief Get the value of the Length scale. + * \return the value of the Length scale. + */ + inline virtual su2double GetAFT_Wonder_Func_var3(unsigned long iPoint) const {return 0.0;} + + /*! + * \brief Get the value of the Length scale. + * \return the value of the Length scale. + */ + inline virtual su2double GetAFT_Wonder_Func_var4(unsigned long iPoint) const {return 0.0;} + + /*! + * \brief Get the value of the Length scale. + * \return the value of the Length scale. + */ + inline virtual su2double GetAFT_Wonder_Func_var5(unsigned long iPoint) const {return 0.0;} + + /*! + * \brief Get the value of the Length scale. + * \return the value of the Length scale. + */ + inline virtual su2double GetAFT_Wonder_Func_var6(unsigned long iPoint) const {return 0.0;} + + /*! + * \brief Get the value of the Length scale. + * \return the value of the Length scale. + */ + inline virtual su2double GetAFT_Wonder_Func_var7(unsigned long iPoint) const {return 0.0;} + + /*! + * \brief Get the value of the Length scale. + * \return the value of the Length scale. + */ + inline virtual su2double GetAFT_Wonder_Func_var8(unsigned long iPoint) const {return 0.0;} + + /*! + * \brief Get the value of the Length scale. + * \return the value of the Length scale. + */ + inline virtual su2double GetAFT_Wonder_Func_var9(unsigned long iPoint) const {return 0.0;} + + /*! + * \brief Get the value of the Length scale. + * \return the value of the Length scale. + */ + inline virtual su2double GetAFT_Wonder_Func_var10(unsigned long iPoint) const {return 0.0;} + + /*! + * \brief Get the value of the Length scale. + * \return the value of the Length scale. + */ + inline virtual su2double GetAFT_Wonder_Func_var11(unsigned long iPoint) const {return 0.0;} + + + /*! + * \brief Get the value of the Length scale. + * \return the value of the Length scale. + */ + inline virtual su2double GetAFT_Wonder_Func_var12(unsigned long iPoint) const {return 0.0;} + + + /*! + * \brief Get the value of the Length scale. + * \return the value of the Length scale. + */ + inline virtual su2double GetAFT_Wonder_Func_var13(unsigned long iPoint) const {return 0.0;} + + + /*! + * \brief Get the value of the Length scale. + * \return the value of the Length scale. + */ + inline virtual su2double GetAFT_Wonder_Func_var14(unsigned long iPoint) const {return 0.0;} + + + /*! + * \brief Get the value of the Length scale. + * \return the value of the Length scale. + */ + inline virtual su2double GetAFT_Wonder_Func_var15(unsigned long iPoint) const {return 0.0;} + /*! * \brief Get the value of the separation intermittency. * \return the value of the separation intermittency. diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index 015e808b626..6126b74a6f7 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -1319,6 +1319,7 @@ void CDriver::InstantiateTransitionNumerics(unsigned short nVar_Trans, int offse const int visc_bound_term = VISC_BOUND_TERM + offset; const bool LM = config->GetKind_Trans_Model() == TURB_TRANS_MODEL::LM; + const bool AFT = config->GetKind_Trans_Model() == TURB_TRANS_MODEL::AFT; /*--- Definition of the convective scheme for each equation and mesh level ---*/ @@ -1329,6 +1330,7 @@ void CDriver::InstantiateTransitionNumerics(unsigned short nVar_Trans, int offse case SPACE_UPWIND : for (auto iMGlevel = 0u; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { if (LM) numerics[iMGlevel][TRANS_SOL][conv_term] = new CUpwSca_TransLM(nDim, nVar_Trans, config); + if (AFT) numerics[iMGlevel][TRANS_SOL][conv_term] = new CUpwSca_TransAFT(nDim, nVar_Trans, config); } break; default: @@ -1340,6 +1342,7 @@ void CDriver::InstantiateTransitionNumerics(unsigned short nVar_Trans, int offse for (auto iMGlevel = 0u; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { if (LM) numerics[iMGlevel][TRANS_SOL][visc_term] = new CAvgGrad_TransLM(nDim, nVar_Trans, true, config); + if (AFT) numerics[iMGlevel][TRANS_SOL][visc_term] = new CAvgGrad_TransAFT(nDim, nVar_Trans, true, config); } /*--- Definition of the source term integration scheme for each equation and mesh level ---*/ @@ -1348,6 +1351,7 @@ void CDriver::InstantiateTransitionNumerics(unsigned short nVar_Trans, int offse 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 (AFT) trans_source_first_term = new CSourcePieceWise_TransAFT(nDim, nVar_Trans, config); numerics[iMGlevel][TRANS_SOL][source_second_term] = new CSourceNothing(nDim, nVar_Trans, config); } @@ -1359,6 +1363,10 @@ void CDriver::InstantiateTransitionNumerics(unsigned short nVar_Trans, int offse numerics[iMGlevel][TRANS_SOL][conv_bound_term] = new CUpwSca_TransLM(nDim, nVar_Trans, config); numerics[iMGlevel][TRANS_SOL][visc_bound_term] = new CAvgGrad_TransLM(nDim, nVar_Trans, false, config); } + if (AFT) { + numerics[iMGlevel][TRANS_SOL][conv_bound_term] = new CUpwSca_TransAFT(nDim, nVar_Trans, config); + numerics[iMGlevel][TRANS_SOL][visc_bound_term] = new CAvgGrad_TransAFT(nDim, nVar_Trans, false, config); + } } } /*--- Explicit instantiation of the template above, needed because it is defined in a cpp file, instead of hpp. ---*/ @@ -1486,7 +1494,7 @@ void CDriver::InitializeNumerics(CConfig *config, CGeometry **geometry, CSolver case MAIN_SOLVER::RANS: case MAIN_SOLVER::DISC_ADJ_RANS: ns = compressible = turbulent = true; - transition = (config->GetKind_Trans_Model() == TURB_TRANS_MODEL::LM); + transition = (config->GetKind_Trans_Model() != TURB_TRANS_MODEL::NONE); species = config->GetKind_Species_Model() != SPECIES_MODEL::NONE; break; case MAIN_SOLVER::INC_EULER: diff --git a/SU2_CFD/src/iteration/CFluidIteration.cpp b/SU2_CFD/src/iteration/CFluidIteration.cpp index 4ca6c45a172..8525c5e7164 100644 --- a/SU2_CFD/src/iteration/CFluidIteration.cpp +++ b/SU2_CFD/src/iteration/CFluidIteration.cpp @@ -85,7 +85,7 @@ void CFluidIteration::Iterate(COutput* output, CIntegration**** integration, CGe /*--- Solve transition model ---*/ - if (config[val_iZone]->GetKind_Trans_Model() == TURB_TRANS_MODEL::LM) { + if (config[val_iZone]->GetKind_Trans_Model() != TURB_TRANS_MODEL::NONE) { config[val_iZone]->SetGlobalParam(main_solver, RUNTIME_TRANS_SYS); integration[val_iZone][val_iInst][TRANS_SOL]->SingleGrid_Iteration(geometry, solver, numerics, config, RUNTIME_TRANS_SYS, val_iZone, val_iInst); diff --git a/SU2_CFD/src/meson.build b/SU2_CFD/src/meson.build index 3b40822a34f..5175a6c09a2 100644 --- a/SU2_CFD/src/meson.build +++ b/SU2_CFD/src/meson.build @@ -56,6 +56,7 @@ su2_cfd_src += files(['output/COutputFactory.cpp', su2_cfd_src += files(['variables/CIncNSVariable.cpp', 'variables/CTransLMVariable.cpp', + 'variables/CTransAFTVariable.cpp', 'variables/CAdjEulerVariable.cpp', 'variables/CHeatVariable.cpp', 'variables/CTurbVariable.cpp', @@ -112,6 +113,7 @@ su2_cfd_src += files(['solvers/CSolverFactory.cpp', 'solvers/CTemplateSolver.cpp', 'solvers/CSpeciesSolver.cpp', 'solvers/CSpeciesFlameletSolver.cpp', + 'solvers/CTransAFTSolver.cpp', 'solvers/CTransLMSolver.cpp', 'solvers/CTurbSolver.cpp', 'solvers/CTurbSASolver.cpp', diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp index 3c651269b02..c9fc690284d 100644 --- a/SU2_CFD/src/output/CFlowOutput.cpp +++ b/SU2_CFD/src/output/CFlowOutput.cpp @@ -984,6 +984,13 @@ void CFlowOutput::AddHistoryOutputFields_ScalarRMS_RES(const CConfig* config) { 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::AFT: + /// DESCRIPTION: Root-mean square residual of the intermittency (LM model). + AddHistoryOutput("RMS_AF", "rms[AFT_1]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of amplification factor (AFT model).", HistoryFieldType::RESIDUAL); + /// DESCRIPTION: Root-mean square residual of natural logarithm intermittency (AFT model). + AddHistoryOutput("RMS_LNINTERMITTENCY", "rms[AFT_2]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of natural logarithm intermittency (AFT model).", HistoryFieldType::RESIDUAL); + break; + case TURB_TRANS_MODEL::NONE: break; } @@ -1040,6 +1047,13 @@ void CFlowOutput::AddHistoryOutputFields_ScalarMAX_RES(const CConfig* config) { 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::AFT: + /// DESCRIPTION: Maximum residual of the amplification factor (AFT model). + AddHistoryOutput("MAX_AF", "max[AFT_1]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of the amplification factor (AFT model).", HistoryFieldType::RESIDUAL); + /// DESCRIPTION: Maximum residual of the natural logarithm intermittency (AFT model). + AddHistoryOutput("MAX_LNINTERMITTENCY", "max[AFT_2]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of the natural logarithm intermittency (AFT model).", HistoryFieldType::RESIDUAL); + break; + case TURB_TRANS_MODEL::NONE: break; } @@ -1096,6 +1110,13 @@ void CFlowOutput::AddHistoryOutputFields_ScalarBGS_RES(const CConfig* config) { 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::AFT: + /// DESCRIPTION: Maximum residual of amplification factor (AFT model). + AddHistoryOutput("BGS_AF", "bgs[AFT_1]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the amplification factor (AFT model).", HistoryFieldType::RESIDUAL); + /// DESCRIPTION: Maximum residual of the natural logarithm intermittency (AFT model). + AddHistoryOutput("BGS_LNINTERMITTENCY", "bgs[AFT_2]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the the natural logarithm intermittency (AFT model).", HistoryFieldType::RESIDUAL); + break; + case TURB_TRANS_MODEL::NONE: break; } @@ -1195,6 +1216,19 @@ void CFlowOutput::LoadHistoryDataScalar(const CConfig* config, const CSolver* co SetHistoryOutputValue("LINSOL_RESIDUAL_TRANS", log10(solver[TRANS_SOL]->GetResLinSolver())); break; + case TURB_TRANS_MODEL::AFT: + SetHistoryOutputValue("RMS_AF", log10(solver[TRANS_SOL]->GetRes_RMS(0))); + SetHistoryOutputValue("RMS_LNINTERMITTENCY",log10(solver[TRANS_SOL]->GetRes_RMS(1))); + SetHistoryOutputValue("MAX_AF", log10(solver[TRANS_SOL]->GetRes_Max(0))); + SetHistoryOutputValue("MAX_LNINTERMITTENCY", log10(solver[TRANS_SOL]->GetRes_Max(1))); + if (multiZone) { + SetHistoryOutputValue("BGS_AF", log10(solver[TRANS_SOL]->GetRes_BGS(0))); + SetHistoryOutputValue("BGS_LNINTERMITTENCY", log10(solver[TRANS_SOL]->GetRes_BGS(1))); + } + SetHistoryOutputValue("LINSOL_ITER_TRANS", solver[TRANS_SOL]->GetIterLinSolver()); + SetHistoryOutputValue("LINSOL_RESIDUAL_TRANS", log10(solver[TRANS_SOL]->GetResLinSolver())); + break; + case TURB_TRANS_MODEL::NONE: break; } @@ -1264,6 +1298,11 @@ void CFlowOutput::SetVolumeOutputFieldsScalarSolution(const CConfig* config){ AddVolumeOutput("RE_THETA_T", "LM_Re_t", "SOLUTION", "LM RE_THETA_T"); break; + case TURB_TRANS_MODEL::AFT: + AddVolumeOutput("AF", "AFT_AF", "SOLUTION", "AFT_AF"); + AddVolumeOutput("LNINTERMITTENCY", "AFT_LNINTERMITTENCY", "SOLUTION", "AFT_LNINTERMITTENCY"); + break; + case TURB_TRANS_MODEL::NONE: break; } @@ -1337,6 +1376,11 @@ void CFlowOutput::SetVolumeOutputFieldsScalarResidual(const CConfig* config) { AddVolumeOutput("RES_RE_THETA_T", "Residual_LM_RE_THETA_T", "RESIDUAL", "Residual of LM RE_THETA_T"); break; + case TURB_TRANS_MODEL::AFT: + AddVolumeOutput("RES_AF", "Residual_AFT_AF", "RESIDUAL", "Residual of AFT amplification factor"); + AddVolumeOutput("RES_LNINTERMITTENCY", "Residual_AFT_LNINTERMITTENCY", "RESIDUAL", "Residual of AFT LNINTERMITTENCY"); + break; + case TURB_TRANS_MODEL::NONE: break; } @@ -1407,6 +1451,25 @@ void CFlowOutput::SetVolumeOutputFieldsScalarPrimitive(const CConfig* config) { AddVolumeOutput("TURB_INDEX", "Turb_index", "PRIMITIVE", "Turbulence index"); break; + case TURB_TRANS_MODEL::AFT: + //nodes -> SetAFT_Wonder_Func(iPoint, HL, H12, dNdRet, Ret0, D_H12, l_H12, m_H12, kv, Rev, Rev0, F_growth, F_crit, PAF, Pg, Dg); + AddVolumeOutput("HL", "HL", "PRIMITIVE", "HL"); + AddVolumeOutput("H12", "H12", "PRIMITIVE", "H12"); + AddVolumeOutput("dNdRet", "dNdRet", "PRIMITIVE", "dNdRet"); + AddVolumeOutput("Ret0", "Ret0", "PRIMITIVE", "Ret0"); + AddVolumeOutput("D_H12", "D_H12", "PRIMITIVE", "D_H12"); + AddVolumeOutput("l_H12", "l_H12", "PRIMITIVE", "l_H12"); + AddVolumeOutput("m_H12", "m_H12", "PRIMITIVE", "m_H12"); + AddVolumeOutput("kv", "kv", "PRIMITIVE", "kv"); + AddVolumeOutput("Rev", "Rev", "PRIMITIVE", "Rev"); + AddVolumeOutput("Rev0", "Rev0", "PRIMITIVE", "Rev0"); + AddVolumeOutput("F_growth", "F_growth", "PRIMITIVE", "F_growth"); + AddVolumeOutput("F_crit", "F_crit", "PRIMITIVE", "F_crit"); + AddVolumeOutput("PAF", "PAF", "PRIMITIVE", "PAF"); + AddVolumeOutput("Pg", "Pg", "PRIMITIVE", "Pg"); + AddVolumeOutput("Dg", "Dg", "PRIMITIVE", "Dg"); + break; + case TURB_TRANS_MODEL::NONE: break; } @@ -1561,6 +1624,25 @@ void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* con SetVolumeOutputValue("RES_RE_THETA_T", iPoint, trans_solver->LinSysRes(iPoint, 1)); break; + case TURB_TRANS_MODEL::AFT: + SetVolumeOutputValue("AF", iPoint, Node_Trans->GetSolution(iPoint, 0)); + SetVolumeOutputValue("LNINTERMITTENCY", iPoint, Node_Trans->GetSolution(iPoint, 1)); + SetVolumeOutputValue("HL", iPoint, Node_Trans->GetAFT_Wonder_Func_var1(iPoint)); + SetVolumeOutputValue("H12", iPoint, Node_Trans->GetAFT_Wonder_Func_var2(iPoint)); + SetVolumeOutputValue("dNdRet", iPoint, Node_Trans->GetAFT_Wonder_Func_var3(iPoint)); + SetVolumeOutputValue("Ret0", iPoint, Node_Trans->GetAFT_Wonder_Func_var4(iPoint)); + SetVolumeOutputValue("D_H12", iPoint, Node_Trans->GetAFT_Wonder_Func_var5(iPoint)); + SetVolumeOutputValue("l_H12", iPoint, Node_Trans->GetAFT_Wonder_Func_var6(iPoint)); + SetVolumeOutputValue("m_H12", iPoint, Node_Trans->GetAFT_Wonder_Func_var7(iPoint)); + SetVolumeOutputValue("kv", iPoint, Node_Trans->GetAFT_Wonder_Func_var8(iPoint)); + SetVolumeOutputValue("Rev", iPoint, Node_Trans->GetAFT_Wonder_Func_var9(iPoint)); + SetVolumeOutputValue("Rev0", iPoint, Node_Trans->GetAFT_Wonder_Func_var10(iPoint)); + SetVolumeOutputValue("F_growth", iPoint, Node_Trans->GetAFT_Wonder_Func_var11(iPoint)); + SetVolumeOutputValue("F_crit", iPoint, Node_Trans->GetAFT_Wonder_Func_var12(iPoint)); + SetVolumeOutputValue("PAF", iPoint, Node_Trans->GetAFT_Wonder_Func_var13(iPoint)); + SetVolumeOutputValue("Pg", iPoint, Node_Trans->GetAFT_Wonder_Func_var14(iPoint)); + SetVolumeOutputValue("Dg", iPoint, Node_Trans->GetAFT_Wonder_Func_var15(iPoint)); + case TURB_TRANS_MODEL::NONE: break; } @@ -2673,6 +2755,12 @@ void CFlowOutput::WriteForcesBreakdown(const CConfig* config, const CSolver* flo file << " (2009)\n"; } break; + case TURB_TRANS_MODEL::AFT: + file << "Coder's Amplification Factor Tranport model"; + if (config->GetAFTParsedOptions().version == AFT_OPTIONS::AFT2019b ) { + file << " : SA-AFT2019b\n"; + } + break; } } break; diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 01626563d40..268dec6eadb 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -788,7 +788,7 @@ void CEulerSolver::SetNondimensionalization(CConfig *config, unsigned short iMes su2double Temperature_FreeStream = 0.0, Mach2Vel_FreeStream = 0.0, ModVel_FreeStream = 0.0, Energy_FreeStream = 0.0, ModVel_FreeStreamND = 0.0, Velocity_Reynolds = 0.0, Omega_FreeStream = 0.0, Omega_FreeStreamND = 0.0, Viscosity_FreeStream = 0.0, - Density_FreeStream = 0.0, Pressure_FreeStream = 0.0, Tke_FreeStream = 0.0, Re_ThetaT_FreeStream = 0.0, + Density_FreeStream = 0.0, Pressure_FreeStream = 0.0, Tke_FreeStream = 0.0, Re_ThetaT_FreeStream = 0.0, N_crit = 0.0, Length_Ref = 0.0, Density_Ref = 0.0, Pressure_Ref = 0.0, Velocity_Ref = 0.0, Temperature_Ref = 0.0, Time_Ref = 0.0, Omega_Ref = 0.0, Force_Ref = 0.0, Gas_Constant_Ref = 0.0, Viscosity_Ref = 0.0, Conductivity_Ref = 0.0, Energy_Ref= 0.0, @@ -811,6 +811,7 @@ void CEulerSolver::SetNondimensionalization(CConfig *config, unsigned short iMes bool gravity = config->GetGravityForce(); bool turbulent = (config->GetKind_Turb_Model() != TURB_MODEL::NONE); bool tkeNeeded = (turbulent && config->GetKind_Turb_Model() == TURB_MODEL::SST); + bool NcritNeeded = (config->GetN_Critical() == 0.0); bool free_stream_temp = (config->GetKind_FreeStreamOption() == FREESTREAM_OPTION::TEMPERATURE_FS); bool reynolds_init = (config->GetKind_InitOption() == REYNOLDS); bool aeroelastic = config->GetAeroelastic_Simulation(); @@ -1073,6 +1074,19 @@ void CEulerSolver::SetNondimensionalization(CConfig *config, unsigned short iMes } config->SetReThetaT_FreeStream(Re_ThetaT_FreeStream); + if (NcritNeeded) { + if (tkeNeeded) { + N_crit = -8.43 - 2.4 * log( (2.5 * tanh( config->GetTurbulenceIntensity_FreeStream() *100 / 2.5)) / 100.0); + } + else { + N_crit = 9.0; + } + } + else { + N_crit = config->GetN_Critical(); + } + config->SetN_Crtical(N_crit); + const su2double MassDiffusivityND = config->GetDiffusivity_Constant() / (Velocity_Ref * Length_Ref); config->SetDiffusivity_ConstantND(MassDiffusivityND); @@ -1373,6 +1387,10 @@ void CEulerSolver::SetNondimensionalization(CConfig *config, unsigned short iMes NonDimTable << "Moment. Thick. Re" << "-" << "-" << "-" << config->GetReThetaT_FreeStream(); Unit.str(""); } + if (config-> GetKind_Trans_Model() == TURB_TRANS_MODEL::AFT) { + NonDimTable << "Critical N-Factor" << "-" << "-" << "-" << config->GetN_Critical(); + Unit.str(""); + } } if (config->GetKind_Species_Model() != SPECIES_MODEL::NONE) { if (config->GetSystemMeasurements() == SI) Unit << "m^2/s"; diff --git a/SU2_CFD/src/solvers/CSolverFactory.cpp b/SU2_CFD/src/solvers/CSolverFactory.cpp index 743775ad139..3c20a8c0750 100644 --- a/SU2_CFD/src/solvers/CSolverFactory.cpp +++ b/SU2_CFD/src/solvers/CSolverFactory.cpp @@ -36,6 +36,7 @@ #include "../../include/solvers/CTurbSASolver.hpp" #include "../../include/solvers/CTurbSSTSolver.hpp" #include "../../include/solvers/CTransLMSolver.hpp" +#include "../../include/solvers/CTransAFTSolver.hpp" #include "../../include/solvers/CAdjEulerSolver.hpp" #include "../../include/solvers/CAdjNSSolver.hpp" #include "../../include/solvers/CAdjTurbSolver.hpp" @@ -383,6 +384,12 @@ CSolver* CSolverFactory::CreateTransSolver(TURB_TRANS_MODEL kindTransModel, CSol transSolver->Postprocessing(geometry, solver, config, iMGLevel); solver[FLOW_SOL]->Preprocessing(geometry, solver, config, iMGLevel, NO_RK_ITER, RUNTIME_FLOW_SYS, false); break; + case TURB_TRANS_MODEL::AFT : + transSolver = new CTransAFTSolver(geometry, config, iMGLevel); + solver[FLOW_SOL]->Preprocessing(geometry, solver, config, iMGLevel, NO_RK_ITER, RUNTIME_FLOW_SYS, false); + transSolver->Postprocessing(geometry, solver, config, iMGLevel); + solver[FLOW_SOL]->Preprocessing(geometry, solver, config, iMGLevel, NO_RK_ITER, RUNTIME_FLOW_SYS, false); + break; case TURB_TRANS_MODEL::NONE: break; } diff --git a/SU2_CFD/src/solvers/CTransAFTSolver.cpp b/SU2_CFD/src/solvers/CTransAFTSolver.cpp new file mode 100644 index 00000000000..f5c4f402d13 --- /dev/null +++ b/SU2_CFD/src/solvers/CTransAFTSolver.cpp @@ -0,0 +1,615 @@ +/*! + * \file CTransAFTSolver.cpp + * \brief Main subroutines for Langtry-Menter Transition model solver. + * \author S. Kang. + * \version 8.1.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2024, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#include "../../include/solvers/CTransAFTSolver.hpp" +#include "../../include/variables/CTransAFTVariable.hpp" +#include "../../include/variables/CFlowVariable.hpp" +#include "../../include/variables/CTurbSAVariable.hpp" +#include "../../../Common/include/parallelization/omp_structure.hpp" +#include "../../../Common/include/toolboxes/geometry_toolbox.hpp" + +/*--- This is the implementation of the Amplification Factor Transport(SA-AFT2019b) transition model. + The main reference for this model is:Coder, In AIAA scitech 2019 forum. + DOI: https://doi.org/10.2514/6.2019-0039 ---*/ + +// Note: TransAFT seems to use rho*AF, rho*ln(Intermittency) as Solution variables, thus Conservative=true + +CTransAFTSolver::CTransAFTSolver(CGeometry *geometry, CConfig *config, unsigned short iMesh) + : CTurbSolver(geometry, config, true) { + unsigned long iPoint; + ifstream restart_file; + string text_line; + + bool multizone = config->GetMultizone_Problem(); + + /*--- Dimension of the problem --> 2 Transport equations (AF, ln(intermittency)) ---*/ + nVar = 2; + nPrimVar = 2; + nPoint = geometry->GetnPoint(); + nPointDomain = geometry->GetnPointDomain(); + + /*--- Initialize nVarGrad for deallocation ---*/ + + nVarGrad = nVar; + + /*--- Define geometry constants in the solver structure ---*/ + + nDim = geometry->GetnDim(); + + /*--- Define variables needed for transition from config file */ + options = config->GetAFTParsedOptions(); + TransCorrelations.SetOptions(options); + + /*--- Single grid simulation ---*/ + + if (iMesh == MESH_0) { + + /*--- Define some auxiliary vector related with the residual ---*/ + + Residual_RMS.resize(nVar,0.0); + Residual_Max.resize(nVar,0.0); + Point_Max.resize(nVar,0); + Point_Max_Coord.resize(nVar,nDim) = su2double(0.0); + + /*--- Initialization of the structure of the whole Jacobian ---*/ + + if (rank == MASTER_NODE) cout << "Initialize Jacobian structure (AFT transition model)." << endl; + Jacobian.Initialize(nPoint, nPointDomain, nVar, nVar, true, geometry, config, ReducerStrategy); + LinSysSol.Initialize(nPoint, nPointDomain, nVar, 0.0); + LinSysRes.Initialize(nPoint, nPointDomain, nVar, 0.0); + System.SetxIsZero(true); + + if (ReducerStrategy) + EdgeFluxes.Initialize(geometry->GetnEdge(), geometry->GetnEdge(), nVar, nullptr); + + /*--- Initialize the BGS residuals in multizone problems. ---*/ + if (multizone){ + Residual_BGS.resize(nVar,0.0); + Residual_Max_BGS.resize(nVar,0.0); + Point_Max_BGS.resize(nVar,0); + Point_Max_Coord_BGS.resize(nVar,nDim) = su2double(0.0); + } + + } + + /*--- Initialize lower and upper limits---*/ + lowerlimit[0] = 0.0; + upperlimit[0] = 100.0; + + lowerlimit[1] = -20.0; + upperlimit[1] = 1.0e-12; + + /*--- Far-field flow state quantities and initialization. ---*/ + const su2double AF_Inf = 0.0; + const su2double lnIntermittency_Inf = 0.0; + + Solution_Inf[0] = AF_Inf; + Solution_Inf[1] = lnIntermittency_Inf; + + /*--- Initialize the solution to the far-field state everywhere. ---*/ + nodes = new CTransAFTVariable(AF_Inf, lnIntermittency_Inf, nPoint, nDim, nVar, config); + SetBaseClassPointerToNodes(); + + /*--- MPI solution ---*/ + + InitiateComms(geometry, config, MPI_QUANTITIES::SOLUTION); + CompleteComms(geometry, config, MPI_QUANTITIES::SOLUTION); + + /*--- Initializate quantities for SlidingMesh Interface ---*/ + + SlidingState.resize(nMarker); + SlidingStateNodes.resize(nMarker); + + for (unsigned long iMarker = 0; iMarker < nMarker; iMarker++) { + if (config->GetMarker_All_KindBC(iMarker) == FLUID_INTERFACE) { + SlidingState[iMarker].resize(nVertex[iMarker], nPrimVar+1) = nullptr; + SlidingStateNodes[iMarker].resize(nVertex[iMarker],0); + } + } + + /*-- Allocation of inlets has to happen in derived classes (not CTurbSolver), + due to arbitrary number of turbulence variables ---*/ + + Inlet_TurbVars.resize(nMarker); + for (unsigned long iMarker = 0; iMarker < nMarker; iMarker++) { + Inlet_TurbVars[iMarker].resize(nVertex[iMarker],nVar); + for (unsigned long iVertex = 0; iVertex < nVertex[iMarker]; ++iVertex) { + Inlet_TurbVars[iMarker](iVertex,0) = AF_Inf; + Inlet_TurbVars[iMarker](iVertex,1) = lnIntermittency_Inf; + } + } + + const su2double CFL = config->GetCFL(MGLevel)*config->GetCFLRedCoeff_Turb(); + for (iPoint = 0; iPoint < nPoint; iPoint++) { + nodes->SetLocalCFL(iPoint, CFL); + nodes->SetAuxVar(iPoint, 0, geometry->nodes->GetWall_Distance(iPoint)); + } + Min_CFL_Local = CFL; + Max_CFL_Local = CFL; + Avg_CFL_Local = CFL; + + /*--- Add the solver name. ---*/ + SolverName = "AFT model"; + +} + +void CTransAFTSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container, CConfig *config, + unsigned short iMesh, unsigned short iRKStep, unsigned short RunTime_EqSystem, bool Output) { + SU2_OMP_SAFE_GLOBAL_ACCESS(config->SetGlobalParam(config->GetKind_Solver(), RunTime_EqSystem);) + + /*--- Upwind second order reconstruction and gradients ---*/ + CommonPreprocessing(geometry, config, Output); +} + +void CTransAFTSolver::Postprocessing(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iMesh) { + + /*--- Compute LM model gradients. ---*/ + + if (config->GetKind_Gradient_Method() == GREEN_GAUSS) { + SetSolution_Gradient_GG(geometry, config, -1); + } + if (config->GetKind_Gradient_Method() == WEIGHTED_LEAST_SQUARES) { + SetSolution_Gradient_LS(geometry, config, -1); + } + + AD::StartNoSharedReading(); + auto* flowNodes = su2staticcast_p(solver_container[FLOW_SOL]->GetNodes()); + auto* turbNodes = su2staticcast_p(solver_container[TURB_SOL]->GetNodes()); + + SU2_OMP_FOR_STAT(omp_chunk_size) + for (unsigned long iPoint = 0; iPoint < nPoint; iPoint ++) { + + // Here the nodes already have the new solution, thus I have to compute everything from scratch + const su2double AF = nodes->GetSolution(iPoint,0); + const su2double lnIntermittency = nodes->GetSolution(iPoint,1); + const su2double Critical_N_Factor = config->GetN_Critical(); + /* Ampification Factor term */ + const su2double vel_u = flowNodes->GetVelocity(iPoint, 0); + const su2double vel_v = flowNodes->GetVelocity(iPoint, 1); + const su2double vel_w = (nDim == 3) ? flowNodes->GetVelocity(iPoint, 2) : 0.0; + const su2double Velocity_Mag = sqrt(vel_u * vel_u + vel_v * vel_v + vel_w * vel_w); + const su2double dist_i = geometry->nodes->GetWall_Distance(iPoint); + const su2double StrainMag_i = flowNodes->GetStrainMag(iPoint); + const su2double Eddy_Viscosity_i = turbNodes->GetmuT(iPoint); + const su2double Laminar_Viscosity_i = flowNodes->GetLaminarViscosity(iPoint); + const su2double Density_i = flowNodes->GetDensity(iPoint); + const su2double c_1 = 100.0; + const su2double c_2 = 0.06; + const su2double c_3 = 50.0; + su2double Temp1 = nodes->GetAuxVarGradient(iPoint, 0, 0); + su2double Temp2 = nodes->GetAuxVarGradient(iPoint, 0, 1); + su2double Temp3 = flowNodes->GetDensity(iPoint) * flowNodes->GetVelocity(iPoint, 0) * nodes->GetAuxVarGradient(iPoint, 0, 0); + su2double Temp4 = flowNodes->GetDensity(iPoint) * flowNodes->GetVelocity(iPoint, 1) * nodes->GetAuxVarGradient(iPoint, 0, 1); + + if(nDim == 2) { + nodes->SetAuxVar(iPoint, 1, Temp3 + Temp4); + } + else { + nodes->SetAuxVar(iPoint, 1, Temp3 + Temp4 + flowNodes->GetDensity(iPoint) * flowNodes->GetVelocity(iPoint, 2) * nodes->GetAuxVarGradient(iPoint, 0, 2)); + } + + su2double VorticityMag = 0.0; + su2double HLGradTerm = 0.0; + HLGradTerm = nodes->GetAuxVarGradient(iPoint, 1, 0) * nodes->GetAuxVarGradient(iPoint, 0, 0) + nodes->GetAuxVarGradient(iPoint, 1, 1) * nodes->GetAuxVarGradient(iPoint, 0, 1); + + if(nDim == 2) { + VorticityMag = sqrt(flowNodes->GetVorticity(iPoint)[0] * flowNodes->GetVorticity(iPoint)[0] + flowNodes->GetVorticity(iPoint)[1] * flowNodes->GetVorticity(iPoint)[1] ); + + } + else { + VorticityMag = sqrt(flowNodes->GetVorticity(iPoint)[0] * flowNodes->GetVorticity(iPoint)[0] + flowNodes->GetVorticity(iPoint)[1] * flowNodes->GetVorticity(iPoint)[1] + + flowNodes->GetVorticity(iPoint)[2] * flowNodes->GetVorticity(iPoint)[2] ); + HLGradTerm += nodes->GetAuxVarGradient(iPoint, 1, 2) * nodes->GetAuxVarGradient(iPoint, 0, 2); + } + + /*--- Cal H12, Hk, dNdRet, Ret0 ---*/ + const su2double HL = dist_i * dist_i / Laminar_Viscosity_i * HLGradTerm; + const su2double H12 = TransCorrelations.H12_Correlations(HL); + const su2double dNdRet = TransCorrelations.dNdRet_Correlations(H12); + const su2double Ret0 = TransCorrelations.Ret0_Correlations(H12); + const su2double D_H12 = TransCorrelations.D_Correlations(H12); + const su2double l_H12 = TransCorrelations.l_Correlations(H12); + const su2double m_H12 = TransCorrelations.m_Correlations(H12, l_H12); + const su2double kv = TransCorrelations.kv_Correlations(H12); + const su2double Rev = Density_i * dist_i * dist_i * StrainMag_i / (Laminar_Viscosity_i + Eddy_Viscosity_i); + const su2double Rev0 = kv * Ret0; + + /*--- production term of the amplification factor ---*/ + const su2double F_growth = max(D_H12 * (1.0 + m_H12) / 2.0 * l_H12, 0.0); + const su2double F_crit = (Rev < Rev0) ? 0.0 : 1.0; + const su2double PAF = Density_i * VorticityMag * F_crit * F_growth * dNdRet; + + /*--- production term of the intermittency(Gamma) ---*/ + const su2double RT = Eddy_Viscosity_i / Laminar_Viscosity_i; + const su2double F_onset1 = AF / Critical_N_Factor; + const su2double F_onset2 = min(F_onset1, 2.0); + const su2double F_onset3 = max(1.0 - pow(RT / 3.5, 3), 0.0); + const su2double F_onset = max(F_onset2 - F_onset3, 0.0); + const su2double F_turb = exp(-pow(RT / 2.0,4)); + const su2double Pg = c_1 * Density_i * StrainMag_i * F_onset * (1.0 - exp(lnIntermittency)); + + /*-- destruction term of Intermeittency(Gamma) --*/ + const su2double Dg = c_2 * Density_i * VorticityMag * F_turb * (c_3 * exp(lnIntermittency) - 1.0); + nodes -> SetAFT_Wonder_Func(iPoint, HL, H12, dNdRet, Ret0, D_H12, l_H12, m_H12, kv, Rev, Rev0, F_growth, F_crit, PAF, Pg, Dg); + } + END_SU2_OMP_FOR + + AD::EndNoSharedReading(); + + switch (config->GetKind_Gradient_Method()) { + case GREEN_GAUSS: + SetAuxVar_Gradient_GG(geometry, config); + break; + case WEIGHTED_LEAST_SQUARES: + SetAuxVar_Gradient_LS(geometry, config); + default: + break; + } +} + + +void CTransAFTSolver::Viscous_Residual(const unsigned long iEdge, const CGeometry* geometry, CSolver** solver_container, + CNumerics* numerics, const CConfig* config) { + + /*--- Define an object to set solver specific numerics contribution. ---*/ + + auto SolverSpecificNumerics = [&](unsigned long iPoint, unsigned long jPoint) {}; + + /*--- Now instantiate the generic implementation with the functor above. ---*/ + + Viscous_Residual_impl(SolverSpecificNumerics, iEdge, geometry, solver_container, numerics, config); +} + + +void CTransAFTSolver::Source_Residual(CGeometry *geometry, CSolver **solver_container, + CNumerics **numerics_container, CConfig *config, unsigned short iMesh) { + + const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); + + auto* flowNodes = su2staticcast_p(solver_container[FLOW_SOL]->GetNodes()); + //auto* turbNodes = su2staticcast_p(solver_container[TURB_SOL]->GetNodes()); + CVariable* turbNodes = solver_container[TURB_SOL]->GetNodes(); + + /*--- Pick one numerics object per thread. ---*/ + auto* numerics = numerics_container[SOURCE_FIRST_TERM + omp_get_thread_num()*MAX_TERMS]; + + /*--- Loop over all points. ---*/ + + AD::StartNoSharedReading(); + + SU2_OMP_FOR_DYN(omp_chunk_size) + for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) { + + + + /*--- Conservative variables w/o reconstruction ---*/ + + numerics->SetPrimitive(flowNodes->GetPrimitive(iPoint), nullptr); + + /*--- Gradient of the primitive and conservative variables ---*/ + + numerics->SetPrimVarGradient(flowNodes->GetGradient_Primitive(iPoint), nullptr); + + /*--- Turbulent variables w/o reconstruction, and its gradient ---*/ + /*--- ScalarVar & ScalarVarGradient : Turbulence model solution(k&w) ---*/ + + numerics->SetScalarVar(turbNodes->GetSolution(iPoint), nullptr); + numerics->SetScalarVarGradient(turbNodes->GetGradient(iPoint), nullptr); + + /*--- Transition variables w/o reconstruction, and its gradient ---*/ + + numerics->SetTransVar(nodes->GetSolution(iPoint), nullptr); + numerics->SetTransVarGradient(nodes->GetGradient(iPoint), nullptr); + + /*--- Set volume ---*/ + + numerics->SetVolume(geometry->nodes->GetVolume(iPoint)); + + /*--- Set distance to the surface ---*/ + + numerics->SetDistance(geometry->nodes->GetWall_Distance(iPoint), 0.0); + + /*--- Set vorticity and strain rate magnitude ---*/ + + numerics->SetVorticity(flowNodes->GetVorticity(iPoint), nullptr); + + numerics->SetStrainMag(flowNodes->GetStrainMag(iPoint), 0.0); + + /*--- Set coordinate (for debugging) ---*/ + numerics->SetCoord(geometry->nodes->GetCoord(iPoint), nullptr); + numerics->SetAuxVarGrad(nodes->GetAuxVarGradient(iPoint), nullptr); + + su2double Temp1 = flowNodes->GetDensity(iPoint) * flowNodes->GetVelocity(iPoint, 0) * nodes->GetAuxVarGradient(iPoint, 0, 0); + su2double Temp2 = flowNodes->GetDensity(iPoint) * flowNodes->GetVelocity(iPoint, 1) * nodes->GetAuxVarGradient(iPoint, 0, 1); + su2double Temp3 = nodes->GetAuxVarGradient(iPoint, 0, 0); + su2double Temp4 = nodes->GetAuxVarGradient(iPoint, 0, 1); + su2double Temp5 = nodes->GetAuxVarGradient(iPoint, 1, 0); + su2double Temp6 = nodes->GetAuxVarGradient(iPoint, 1, 1); + + + /*--- Compute the source term ---*/ + auto residual = numerics->ComputeResidual(config); + + /*--- Subtract residual and the Jacobian ---*/ + + LinSysRes.SubtractBlock(iPoint, residual); + if (implicit) Jacobian.SubtractBlock2Diag(iPoint, residual.jacobian_i); + + } + END_SU2_OMP_FOR + + AD::EndNoSharedReading(); + +} + +void CTransAFTSolver::Source_Template(CGeometry *geometry, CSolver **solver_container, CNumerics *numerics, + CConfig *config, unsigned short iMesh) { +} + +void CTransAFTSolver::BC_HeatFlux_Wall(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, + CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { + + const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); + + SU2_OMP_FOR_STAT(OMP_MIN_SIZE) + for (auto iVertex = 0u; iVertex < geometry->nVertex[val_marker]; iVertex++) { + + const auto iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); + + /*--- Check if the node belongs to the domain (i.e, not a halo node) ---*/ + + if (geometry->nodes->GetDomain(iPoint)) { + + /*--- Allocate the value at the infinity ---*/ + + auto V_infty = solver_container[FLOW_SOL]->GetCharacPrimVar(val_marker, iVertex); + + /*--- Retrieve solution at the farfield boundary node ---*/ + + auto V_domain = solver_container[FLOW_SOL]->GetNodes()->GetPrimitive(iPoint); + + conv_numerics->SetPrimitive(V_domain, V_infty); + + /*--- Set turbulent variable at the wall, and at infinity ---*/ + + conv_numerics->SetScalarVar(nodes->GetSolution(iPoint), Solution_Inf); + + /*--- Set Normal (it is necessary to change the sign) ---*/ + /*--- It's mean wall normal zero flux. */ + + su2double Normal[MAXNDIM] = {0.0}; + for (auto iDim = 0u; iDim < nDim; iDim++) + Normal[iDim] = -geometry->vertex[val_marker][iVertex]->GetNormal(iDim); + conv_numerics->SetNormal(Normal); + + /*--- Grid Movement ---*/ + + if (dynamic_grid) + conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), + geometry->nodes->GetGridVel(iPoint)); + + /*--- Compute residuals and Jacobians ---*/ + + auto residual = conv_numerics->ComputeResidual(config); + + /*--- Add residuals and Jacobians ---*/ + + LinSysRes.AddBlock(iPoint, residual); + if (implicit) Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i); + } + } + END_SU2_OMP_FOR + +} + +void CTransAFTSolver::BC_Isothermal_Wall(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, + CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { + + BC_HeatFlux_Wall(geometry, solver_container, conv_numerics, visc_numerics, config, val_marker); + +} + +void CTransAFTSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, + unsigned short val_marker) { + const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); + + /*--- Loop over all the vertices on this boundary marker ---*/ + + SU2_OMP_FOR_STAT(OMP_MIN_SIZE) + for (auto iVertex = 0u; iVertex < geometry->nVertex[val_marker]; iVertex++) { + + const auto iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); + + /*--- Check if the node belongs to the domain (i.e., not a halo node) ---*/ + + if (geometry->nodes->GetDomain(iPoint)) { + + /*--- Normal vector for this vertex (negate for outward convention) ---*/ + + su2double Normal[MAXNDIM] = {0.0}; + for (auto iDim = 0u; iDim < nDim; iDim++) + Normal[iDim] = -geometry->vertex[val_marker][iVertex]->GetNormal(iDim); + conv_numerics->SetNormal(Normal); + + /*--- Allocate the value at the inlet ---*/ + + auto V_inlet = solver_container[FLOW_SOL]->GetCharacPrimVar(val_marker, iVertex); + + /*--- Retrieve solution at the farfield boundary node ---*/ + + auto V_domain = solver_container[FLOW_SOL]->GetNodes()->GetPrimitive(iPoint); + + /*--- Set various quantities in the solver class ---*/ + + conv_numerics->SetPrimitive(V_domain, V_inlet); + + /*--- 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)); + } + + /*--- Set the LM variable states. ---*/ + /*--- Load the inlet transition LM model variables (uniform by default). ---*/ + + conv_numerics->SetScalarVar(nodes->GetSolution(iPoint), Inlet_Vars); + + /*--- Set various other quantities in the solver class ---*/ + + if (dynamic_grid) + conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), + geometry->nodes->GetGridVel(iPoint)); + + /*--- Compute the residual using an upwind scheme ---*/ + + auto residual = conv_numerics->ComputeResidual(config); + LinSysRes.AddBlock(iPoint, residual); + + /*--- Jacobian contribution for implicit integration ---*/ + + if (implicit) Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i); + } + } + END_SU2_OMP_FOR + +} + +void CTransAFTSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, + CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { + BC_Far_Field(geometry, solver_container, conv_numerics, visc_numerics, config, val_marker); +} + +void CTransAFTSolver::LoadRestart(CGeometry** geometry, CSolver*** solver, CConfig* config, int val_iter, + bool val_update_geo) { + + const string restart_filename = config->GetFilename(config->GetSolution_FileName(), "", val_iter); + + /*--- To make this routine safe to call in parallel most of it can only be executed by one thread. ---*/ + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { + /*--- Read the restart data from either an ASCII or binary SU2 file. ---*/ + + if (config->GetRead_Binary_Restart()) { + Read_SU2_Restart_Binary(geometry[MESH_0], config, restart_filename); + } else { + Read_SU2_Restart_ASCII(geometry[MESH_0], config, restart_filename); + } + + /*--- Skip flow variables ---*/ + + unsigned short skipVars = nDim + solver[MESH_0][FLOW_SOL]->GetnVar() + solver[MESH_0][TURB_SOL] ->GetnVar(); + + /*--- Adjust the number of solution variables in the incompressible + restart. We always carry a space in nVar for the energy equation in the + mean flow solver, but we only write it to the restart if it is active. + Therefore, we must reduce skipVars here if energy is inactive so that + the turbulent variables are read correctly. ---*/ + + const bool incompressible = (config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE); + const bool energy = config->GetEnergy_Equation(); + const bool weakly_coupled_heat = config->GetWeakly_Coupled_Heat(); + + if (incompressible && ((!energy) && (!weakly_coupled_heat))) skipVars--; + + /*--- Load data from the restart into correct containers. ---*/ + + unsigned long counter = 0; + for (auto iPoint_Global = 0ul; iPoint_Global < geometry[MESH_0]->GetGlobal_nPointDomain(); iPoint_Global++) { + /*--- Retrieve local index. If this node from the restart file lives + on the current processor, we will load and instantiate the vars. ---*/ + + const auto iPoint_Local = geometry[MESH_0]->GetGlobal_to_Local_Point(iPoint_Global); + + if (iPoint_Local > -1) { + /*--- We need to store this point's data, so jump to the correct + offset in the buffer of data from the restart file and load it. ---*/ + + 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 ->SetAFT_Wonder_Func(iPoint_Local, Restart_Data[index + 2], Restart_Data[index + 3] + , Restart_Data[index + 4], Restart_Data[index + 5], Restart_Data[index + 6], Restart_Data[index + 7] + , Restart_Data[index + 8], Restart_Data[index + 9], Restart_Data[index + 10], Restart_Data[index + 11] + , Restart_Data[index + 12], Restart_Data[index + 13], Restart_Data[index + 14], Restart_Data[index + 15] + , Restart_Data[index + 16]); + + /*--- Increment the overall counter for how many points have been loaded. ---*/ + counter++; + } + } + + /*--- Detect a wrong solution file ---*/ + + if (counter != nPointDomain) { + SU2_MPI::Error(string("The solution file ") + restart_filename + string(" does not match with the mesh file!\n") + + string("This can be caused by empty lines at the end of the file."), + CURRENT_FUNCTION); + } + + } // end SU2_OMP_MASTER, pre and postprocessing are thread-safe. + END_SU2_OMP_SAFE_GLOBAL_ACCESS + + /*--- MPI solution and compute the eddy viscosity ---*/ + + solver[MESH_0][TRANS_SOL]->InitiateComms(geometry[MESH_0], config, MPI_QUANTITIES::SOLUTION); + solver[MESH_0][TRANS_SOL]->CompleteComms(geometry[MESH_0], config, MPI_QUANTITIES::SOLUTION); + + /*--- For turbulent+species simulations the solver Pre-/Postprocessing is done by the species solver. ---*/ + if (config->GetKind_Species_Model() == SPECIES_MODEL::NONE) { + solver[MESH_0][FLOW_SOL]->Preprocessing(geometry[MESH_0], solver[MESH_0], config, MESH_0, NO_RK_ITER, + RUNTIME_FLOW_SYS, false); + solver[MESH_0][TURB_SOL]->Postprocessing(geometry[MESH_0], solver[MESH_0], config, MESH_0); + solver[MESH_0][TRANS_SOL]->Postprocessing(geometry[MESH_0], solver[MESH_0], config, MESH_0); + } + + /*--- Interpolate the solution down to the coarse multigrid levels ---*/ + + for (auto iMesh = 1u; iMesh <= config->GetnMGLevels(); iMesh++) { + + MultigridRestriction(*geometry[iMesh - 1], solver[iMesh - 1][TRANS_SOL]->GetNodes()->GetSolution(), + *geometry[iMesh], solver[iMesh][TRANS_SOL]->GetNodes()->GetSolution()); + solver[iMesh][TRANS_SOL]->InitiateComms(geometry[iMesh], config, MPI_QUANTITIES::SOLUTION); + solver[iMesh][TRANS_SOL]->CompleteComms(geometry[iMesh], config, MPI_QUANTITIES::SOLUTION); + + if (config->GetKind_Species_Model() == SPECIES_MODEL::NONE) { + solver[iMesh][FLOW_SOL]->Preprocessing(geometry[iMesh], solver[iMesh], config, iMesh, NO_RK_ITER, RUNTIME_FLOW_SYS, + false); + solver[iMesh][TRANS_SOL]->Postprocessing(geometry[iMesh], solver[iMesh], config, iMesh); + } + } + + /*--- Go back to single threaded execution. ---*/ + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { + /*--- Delete the class memory that is used to load the restart. ---*/ + + Restart_Vars.clear(); + Restart_Data.clear(); + } + 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 6b465e802fd..ed4b56d9eed 100644 --- a/SU2_CFD/src/solvers/CTurbSASolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSASolver.cpp @@ -380,9 +380,11 @@ void CTurbSASolver::Source_Residual(CGeometry *geometry, CSolver **solver_contai /*--- Effective Intermittency ---*/ - if (config->GetKind_Trans_Model() != TURB_TRANS_MODEL::NONE) { + if (config->GetKind_Trans_Model() == TURB_TRANS_MODEL::LM) { numerics->SetIntermittencyEff(solver_container[TRANS_SOL]->GetNodes()->GetIntermittencyEff(iPoint)); numerics->SetIntermittency(solver_container[TRANS_SOL]->GetNodes()->GetSolution(iPoint, 0)); + } else if (config->GetKind_Trans_Model() == TURB_TRANS_MODEL::AFT) { + numerics->SetIntermittency(solver_container[TRANS_SOL]->GetNodes()->GetSolution(iPoint, 1)); } if (axisymmetric) { diff --git a/SU2_CFD/src/variables/CTransAFTVariable.cpp b/SU2_CFD/src/variables/CTransAFTVariable.cpp new file mode 100644 index 00000000000..2d0a7b9fa2b --- /dev/null +++ b/SU2_CFD/src/variables/CTransAFTVariable.cpp @@ -0,0 +1,85 @@ +/*! + * \file CTransAFTVariable.cpp + * \brief Definition of the solution fields. + * \author S. Kang + * \version 8.1.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2024, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ +#include "../../include/variables/CTransAFTVariable.hpp" + +CTransAFTVariable::CTransAFTVariable(su2double AF, su2double LnIntermittency, unsigned long npoint, unsigned long ndim, unsigned long nvar, CConfig *config) + : CTurbVariable(npoint, ndim, nvar, config) { + + for(unsigned long iPoint=0; iPoint Date: Tue, 7 Jan 2025 10:37:17 +0900 Subject: [PATCH 2/6] Test_SA_AFT2019b Test for SAAFT2019b --- SU2_CFD/include/numerics/turbulent/turb_sources.hpp | 2 +- SU2_CFD/src/solvers/CTransAFTSolver.cpp | 10 ++-------- 2 files changed, 3 insertions(+), 9 deletions(-) diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp index 2bfbe830d87..6ad99b0f2fe 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp @@ -336,7 +336,7 @@ struct Nonzero { var.ft2 = var.ct3 * exp(-var.ct4 * xsi2); var.d_ft2 = -2.0 * var.ct4 * var.Ji * var.ft2 * var.d_Ji; if(transition_AFT) { - var.ft2 = var.ct3 * exp(1.0 - var.lnintermittecncy); + var.ft2 = var.ct3 * (1.0 - var.lnintermittecncy); var.d_ft2 = 0.0; } } diff --git a/SU2_CFD/src/solvers/CTransAFTSolver.cpp b/SU2_CFD/src/solvers/CTransAFTSolver.cpp index f5c4f402d13..d91d959b8bd 100644 --- a/SU2_CFD/src/solvers/CTransAFTSolver.cpp +++ b/SU2_CFD/src/solvers/CTransAFTSolver.cpp @@ -212,17 +212,11 @@ void CTransAFTSolver::Postprocessing(CGeometry *geometry, CSolver **solver_conta nodes->SetAuxVar(iPoint, 1, Temp3 + Temp4 + flowNodes->GetDensity(iPoint) * flowNodes->GetVelocity(iPoint, 2) * nodes->GetAuxVarGradient(iPoint, 0, 2)); } - su2double VorticityMag = 0.0; + su2double VorticityMag = GeometryToolbox::Norm(3, flowNodes->GetVorticity(iPoint)); su2double HLGradTerm = 0.0; HLGradTerm = nodes->GetAuxVarGradient(iPoint, 1, 0) * nodes->GetAuxVarGradient(iPoint, 0, 0) + nodes->GetAuxVarGradient(iPoint, 1, 1) * nodes->GetAuxVarGradient(iPoint, 0, 1); - if(nDim == 2) { - VorticityMag = sqrt(flowNodes->GetVorticity(iPoint)[0] * flowNodes->GetVorticity(iPoint)[0] + flowNodes->GetVorticity(iPoint)[1] * flowNodes->GetVorticity(iPoint)[1] ); - - } - else { - VorticityMag = sqrt(flowNodes->GetVorticity(iPoint)[0] * flowNodes->GetVorticity(iPoint)[0] + flowNodes->GetVorticity(iPoint)[1] * flowNodes->GetVorticity(iPoint)[1] - + flowNodes->GetVorticity(iPoint)[2] * flowNodes->GetVorticity(iPoint)[2] ); + if(nDim == 3) { HLGradTerm += nodes->GetAuxVarGradient(iPoint, 1, 2) * nodes->GetAuxVarGradient(iPoint, 0, 2); } From 448f5b91a2ec961113bebb929b1d86a6af842be6 Mon Sep 17 00:00:00 2001 From: "oosun93@pusan.ac.kr" Date: Tue, 7 Jan 2025 19:55:41 +0900 Subject: [PATCH 3/6] Renew_Min_Max Renew_Min_Max --- SU2_CFD/src/solvers/CTransAFTSolver.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/SU2_CFD/src/solvers/CTransAFTSolver.cpp b/SU2_CFD/src/solvers/CTransAFTSolver.cpp index d91d959b8bd..0e93a4051c9 100644 --- a/SU2_CFD/src/solvers/CTransAFTSolver.cpp +++ b/SU2_CFD/src/solvers/CTransAFTSolver.cpp @@ -98,10 +98,10 @@ CTransAFTSolver::CTransAFTSolver(CGeometry *geometry, CConfig *config, unsigned /*--- Initialize lower and upper limits---*/ lowerlimit[0] = 0.0; - upperlimit[0] = 100.0; + upperlimit[0] = 20.0; - lowerlimit[1] = -20.0; - upperlimit[1] = 1.0e-12; + lowerlimit[1] = -4.0; + upperlimit[1] = 1.0; /*--- Far-field flow state quantities and initialization. ---*/ const su2double AF_Inf = 0.0; From a583d7c240cf7586ff81f59d95d93bb07b08e88d Mon Sep 17 00:00:00 2001 From: "oosun93@pusan.ac.kr" Date: Tue, 14 Jan 2025 22:02:35 +0900 Subject: [PATCH 4/6] HL update HL --- .../numerics/turbulent/transition/trans_sources.hpp | 4 ++-- SU2_CFD/src/solvers/CTransAFTSolver.cpp | 12 ++++++------ 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/SU2_CFD/include/numerics/turbulent/transition/trans_sources.hpp b/SU2_CFD/include/numerics/turbulent/transition/trans_sources.hpp index 0e652501f8d..a6c845314b6 100644 --- a/SU2_CFD/include/numerics/turbulent/transition/trans_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/transition/trans_sources.hpp @@ -422,7 +422,7 @@ class CSourcePieceWise_TransAFT final : public CNumerics { su2double tt = 0; } - const su2double HL = dist_i * dist_i / Laminar_Viscosity_i * HLGradTerm; + const su2double HL = dist_i * dist_i * Density_i / Laminar_Viscosity_i * HLGradTerm; /*--- Cal H12, dNdRet, Ret0, D_H12, l_H12, m_H12, kv ---*/ const su2double H12 = TransCorrelations.H12_Correlations(HL); const su2double dNdRet = TransCorrelations.dNdRet_Correlations(H12); @@ -463,7 +463,7 @@ class CSourcePieceWise_TransAFT final : public CNumerics { Jacobian_i[0][0] = 0.0; Jacobian_i[0][1] = 0.0; Jacobian_i[1][0] = 0.0; - Jacobian_i[1][1] = ( -c_1 * StrainMag_i * F_onset * exp(TransVar_i[1]) -c_2 * VorticityMag * F_turb * c_3 * exp(TransVar_i[1]) ) * Volume; + Jacobian_i[1][1] = (-c_1 * StrainMag_i * F_onset * exp(TransVar_i[1]) - c_2 * VorticityMag * F_turb * c_3 * exp(TransVar_i[1])) * Volume; } AD::SetPreaccOut(Residual, nVar); diff --git a/SU2_CFD/src/solvers/CTransAFTSolver.cpp b/SU2_CFD/src/solvers/CTransAFTSolver.cpp index 0e93a4051c9..a9494fae717 100644 --- a/SU2_CFD/src/solvers/CTransAFTSolver.cpp +++ b/SU2_CFD/src/solvers/CTransAFTSolver.cpp @@ -202,14 +202,14 @@ void CTransAFTSolver::Postprocessing(CGeometry *geometry, CSolver **solver_conta const su2double c_3 = 50.0; su2double Temp1 = nodes->GetAuxVarGradient(iPoint, 0, 0); su2double Temp2 = nodes->GetAuxVarGradient(iPoint, 0, 1); - su2double Temp3 = flowNodes->GetDensity(iPoint) * flowNodes->GetVelocity(iPoint, 0) * nodes->GetAuxVarGradient(iPoint, 0, 0); - su2double Temp4 = flowNodes->GetDensity(iPoint) * flowNodes->GetVelocity(iPoint, 1) * nodes->GetAuxVarGradient(iPoint, 0, 1); + su2double Temp3 = flowNodes->GetVelocity(iPoint, 0) * nodes->GetAuxVarGradient(iPoint, 0, 0); + su2double Temp4 = flowNodes->GetVelocity(iPoint, 1) * nodes->GetAuxVarGradient(iPoint, 0, 1); if(nDim == 2) { nodes->SetAuxVar(iPoint, 1, Temp3 + Temp4); } else { - nodes->SetAuxVar(iPoint, 1, Temp3 + Temp4 + flowNodes->GetDensity(iPoint) * flowNodes->GetVelocity(iPoint, 2) * nodes->GetAuxVarGradient(iPoint, 0, 2)); + nodes->SetAuxVar(iPoint, 1, Temp3 + Temp4 + flowNodes->GetVelocity(iPoint, 2) * nodes->GetAuxVarGradient(iPoint, 0, 2)); } su2double VorticityMag = GeometryToolbox::Norm(3, flowNodes->GetVorticity(iPoint)); @@ -221,7 +221,7 @@ void CTransAFTSolver::Postprocessing(CGeometry *geometry, CSolver **solver_conta } /*--- Cal H12, Hk, dNdRet, Ret0 ---*/ - const su2double HL = dist_i * dist_i / Laminar_Viscosity_i * HLGradTerm; + const su2double HL = dist_i * dist_i * Density_i / Laminar_Viscosity_i * HLGradTerm; const su2double H12 = TransCorrelations.H12_Correlations(HL); const su2double dNdRet = TransCorrelations.dNdRet_Correlations(H12); const su2double Ret0 = TransCorrelations.Ret0_Correlations(H12); @@ -337,8 +337,8 @@ void CTransAFTSolver::Source_Residual(CGeometry *geometry, CSolver **solver_cont numerics->SetCoord(geometry->nodes->GetCoord(iPoint), nullptr); numerics->SetAuxVarGrad(nodes->GetAuxVarGradient(iPoint), nullptr); - su2double Temp1 = flowNodes->GetDensity(iPoint) * flowNodes->GetVelocity(iPoint, 0) * nodes->GetAuxVarGradient(iPoint, 0, 0); - su2double Temp2 = flowNodes->GetDensity(iPoint) * flowNodes->GetVelocity(iPoint, 1) * nodes->GetAuxVarGradient(iPoint, 0, 1); + su2double Temp1 = flowNodes->GetVelocity(iPoint, 0) * nodes->GetAuxVarGradient(iPoint, 0, 0); + su2double Temp2 = flowNodes->GetVelocity(iPoint, 1) * nodes->GetAuxVarGradient(iPoint, 0, 1); su2double Temp3 = nodes->GetAuxVarGradient(iPoint, 0, 0); su2double Temp4 = nodes->GetAuxVarGradient(iPoint, 0, 1); su2double Temp5 = nodes->GetAuxVarGradient(iPoint, 1, 0); From 5349d973b1830d429dd431b30008e4ba5528c956 Mon Sep 17 00:00:00 2001 From: "oosun93@pusan.ac.kr" Date: Thu, 16 Jan 2025 20:48:12 +0900 Subject: [PATCH 5/6] Updat_AFT2017b AFT2017b --- Common/include/option_structure.hpp | 9 ++++ Common/src/CConfig.cpp | 7 ++++ .../transition/trans_correlations.hpp | 41 ++++++++++++++++++- .../turbulent/transition/trans_sources.hpp | 2 +- SU2_CFD/src/solvers/CTransAFTSolver.cpp | 2 +- 5 files changed, 57 insertions(+), 4 deletions(-) diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 70592980ecc..a2d29a48287 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -1331,11 +1331,13 @@ inline LM_ParsedOptions ParseLMOptions(const LM_OPTIONS *LM_Options, unsigned sh */ enum class AFT_OPTIONS { NONE, /*!< \brief No option / default. */ + AFT2017b, /*!< \brief using AFT2017b model. */ AFT2019b /*!< \brief using AFT2019b model. */ }; static const MapType AFT_Options_Map = { MakePair("NONE", AFT_OPTIONS::NONE) + MakePair("AFT2017b", AFT_OPTIONS::AFT2017b) MakePair("AFT2019b", AFT_OPTIONS::AFT2019b) }; @@ -1344,6 +1346,7 @@ static const MapType AFT_Options_Map = { */ enum class AFT_CORRELATION { NONE, /*!< \brief No option / default. */ + AFT2017b, /*!< \brief Kind of transition correlation model (AFT2017b). */ AFT2019b /*!< \brief Kind of transition correlation model (AFT2019b). */ }; @@ -1371,6 +1374,12 @@ inline AFT_ParsedOptions ParseAFTOptions(const AFT_OPTIONS *AFT_Options, unsigne }; int NFoundCorrelations = 0; + if (IsPresent(AFT_OPTIONS::AFT2017b)) { + AFTParsedOptions.Correlation = AFT_CORRELATION::AFT2017b; + AFTParsedOptions.version = AFT_OPTIONS::AFT2017b; + NFoundCorrelations++; + } + if (IsPresent(AFT_OPTIONS::AFT2019b)) { AFTParsedOptions.Correlation = AFT_CORRELATION::AFT2019b; AFTParsedOptions.version = AFT_OPTIONS::AFT2019b; diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 37c00b0ab35..4aa7e3bc0c0 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -6338,6 +6338,13 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { if (Kind_Trans_Model == TURB_TRANS_MODEL::AFT) { switch (aftParsedOptions.Correlation) { + case AFT_CORRELATION::AFT2017b: + switch (Kind_Turb_Model) { + case TURB_MODEL::NONE: SU2_MPI::Error("No turbulence model has been selected but AFT transition model is active.", CURRENT_FUNCTION); break; + case TURB_MODEL::SST: SU2_MPI::Error("k-w SST turbulence model has been selected but AFT transition model is active.", CURRENT_FUNCTION); break; + } + cout << "-2017b" << endl; break; + if(!saParsedOptions.ft2) SU2_MPI::Error("ft2 option of SA model has been not selected.", CURRENT_FUNCTION); case AFT_CORRELATION::AFT2019b: switch (Kind_Turb_Model) { case TURB_MODEL::NONE: SU2_MPI::Error("No turbulence model has been selected but AFT transition model is active.", CURRENT_FUNCTION); break; diff --git a/SU2_CFD/include/numerics/turbulent/transition/trans_correlations.hpp b/SU2_CFD/include/numerics/turbulent/transition/trans_correlations.hpp index af0a434c90d..5934af3ece9 100644 --- a/SU2_CFD/include/numerics/turbulent/transition/trans_correlations.hpp +++ b/SU2_CFD/include/numerics/turbulent/transition/trans_correlations.hpp @@ -221,8 +221,15 @@ class TransAFTCorrelations { su2double H12 = 0.0; switch (options.Correlation) { + case AFT_CORRELATION::AFT2017b: { + H12 = 0.376960 + sqrt( (HL + 2.453432) / 0.653181 ); + H12 = min(max(H12, 2.2), 20.0); + break; + } + case AFT_CORRELATION::AFT2019b: { H12 = 0.26 * HL + 2.4; + H12 = min(max(H12, 2.2), 20.0); break; } @@ -231,7 +238,6 @@ class TransAFTCorrelations { CURRENT_FUNCTION); break; } - H12 = min(max(H12, 2.2), 20.0); return H12; } @@ -244,6 +250,11 @@ class TransAFTCorrelations { su2double dNdRet = 0.0; switch (options.Correlation) { + case AFT_CORRELATION::AFT2017b: { + dNdRet = 0.028*(H12 - 1.0) - 0.0345 * exp(-pow(3.87/(H12 -1.0) - 2.52, 2.0)); + break; + } + case AFT_CORRELATION::AFT2019b: { dNdRet = 0.028*(H12 - 1.0) - 0.0345 * exp(-pow(3.87/(H12 -1.0) - 2.52, 2.0)); break; @@ -266,6 +277,12 @@ class TransAFTCorrelations { su2double Ret0 = 0.0; switch (options.Correlation) { + case AFT_CORRELATION::AFT2017b: { + Ret0 = 0.7 * tanh( 14.0 / (H12 - 1.0) - 9.24) + 2.492 / pow(H12 - 1.0, 0.43) + 0.62; + Ret0 = pow(10, Ret0); + break; + } + case AFT_CORRELATION::AFT2019b: { Ret0 = 0.7 * tanh( 14.0 / (H12 - 1.0) - 9.24) + 2.492 / pow(H12 - 1.0, 0.43) + 0.62; Ret0 = pow(10, Ret0); @@ -289,6 +306,11 @@ class TransAFTCorrelations { su2double D_H12 = 0.0; switch (options.Correlation) { + case AFT_CORRELATION::AFT2017b: { + D_H12 = H12 / ( 0.5482 * H12 - 0.5185); + break; + } + case AFT_CORRELATION::AFT2019b: { D_H12 = 2.4 * H12 / (H12 - 1.0); break; @@ -311,6 +333,11 @@ class TransAFTCorrelations { su2double l_H12 = 0.0; switch (options.Correlation) { + case AFT_CORRELATION::AFT2017b: { + l_H12 = (6.54 * H12 - 14.07) / pow(H12, 2.0); + break; + } + case AFT_CORRELATION::AFT2019b: { l_H12 = (6.54 * H12 - 14.07) / pow(H12, 2.0); break; @@ -333,11 +360,16 @@ class TransAFTCorrelations { su2double m_H12 = 0.0; switch (options.Correlation) { - case AFT_CORRELATION::AFT2019b: { + case AFT_CORRELATION::AFT2017b: { m_H12 = (0.058 * pow(H12 - 4, 2.0) / (H12 - 1.0) - 0.068) / l_Correlation; break; } + case AFT_CORRELATION::AFT2019b: { + m_H12 = (0.058 * pow(H12 - 4.0, 2.0) / (H12 - 1.0) - 0.068) / l_Correlation; + break; + } + case AFT_CORRELATION::NONE: SU2_MPI::Error("Transition correlation is set to DEFAULT but no default value has ben set in the code.", CURRENT_FUNCTION); @@ -355,6 +387,11 @@ class TransAFTCorrelations { su2double kv = 0.0; switch (options.Correlation) { + case AFT_CORRELATION::AFT2017b: { + kv = 0.246175 * pow(H12, 2.0) - 0.141831 * H12 + 0.008886; + break; + } + case AFT_CORRELATION::AFT2019b: { kv = 1.0 / (0.4036 * pow(H12, 2.0) - 2.5394 * H12 + 4.3273); break; diff --git a/SU2_CFD/include/numerics/turbulent/transition/trans_sources.hpp b/SU2_CFD/include/numerics/turbulent/transition/trans_sources.hpp index a6c845314b6..9dfdf570727 100644 --- a/SU2_CFD/include/numerics/turbulent/transition/trans_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/transition/trans_sources.hpp @@ -422,7 +422,7 @@ class CSourcePieceWise_TransAFT final : public CNumerics { su2double tt = 0; } - const su2double HL = dist_i * dist_i * Density_i / Laminar_Viscosity_i * HLGradTerm; + su2double HL = min( max(dist_i * dist_i * Density_i / Laminar_Viscosity_i * HLGradTerm, -0.25), 200.0); /*--- Cal H12, dNdRet, Ret0, D_H12, l_H12, m_H12, kv ---*/ const su2double H12 = TransCorrelations.H12_Correlations(HL); const su2double dNdRet = TransCorrelations.dNdRet_Correlations(H12); diff --git a/SU2_CFD/src/solvers/CTransAFTSolver.cpp b/SU2_CFD/src/solvers/CTransAFTSolver.cpp index a9494fae717..b4ac57bec82 100644 --- a/SU2_CFD/src/solvers/CTransAFTSolver.cpp +++ b/SU2_CFD/src/solvers/CTransAFTSolver.cpp @@ -221,7 +221,7 @@ void CTransAFTSolver::Postprocessing(CGeometry *geometry, CSolver **solver_conta } /*--- Cal H12, Hk, dNdRet, Ret0 ---*/ - const su2double HL = dist_i * dist_i * Density_i / Laminar_Viscosity_i * HLGradTerm; + su2double HL = min( max(dist_i * dist_i * Density_i / Laminar_Viscosity_i * HLGradTerm, -0.25), 200.0); const su2double H12 = TransCorrelations.H12_Correlations(HL); const su2double dNdRet = TransCorrelations.dNdRet_Correlations(H12); const su2double Ret0 = TransCorrelations.Ret0_Correlations(H12); From 0de4e65684a1fea091770442ad9ef74b2be9b731 Mon Sep 17 00:00:00 2001 From: "oosun93@pusan.ac.kr" Date: Thu, 16 Jan 2025 23:02:38 +0900 Subject: [PATCH 6/6] SA_AFT2019b Test_Strain --- .../turbulent/transition/trans_sources.hpp | 24 ++++++++++++------- SU2_CFD/src/solvers/CTransAFTSolver.cpp | 11 --------- 2 files changed, 16 insertions(+), 19 deletions(-) diff --git a/SU2_CFD/include/numerics/turbulent/transition/trans_sources.hpp b/SU2_CFD/include/numerics/turbulent/transition/trans_sources.hpp index 9dfdf570727..a6241afb1dc 100644 --- a/SU2_CFD/include/numerics/turbulent/transition/trans_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/transition/trans_sources.hpp @@ -418,11 +418,22 @@ class CSourcePieceWise_TransAFT final : public CNumerics { if(nDim == 3) { HLGradTerm += AuxVar_Grad_i[1][2] * AuxVar_Grad_i[0][2]; } - if(HLGradTerm > 0) { - su2double tt = 0; + su2double HL = 0.0; + switch (options.Correlation) + { + case AFT_CORRELATION::AFT2017b : + HL = min( max(dist_i * dist_i * Density_i / Laminar_Viscosity_i * HLGradTerm, -0.25), 200.0); + break; + + case AFT_CORRELATION::AFT2019b : + HL = dist_i * dist_i * Density_i / Laminar_Viscosity_i * HLGradTerm; + break; + + default: HL = 0.0; + break; } - su2double HL = min( max(dist_i * dist_i * Density_i / Laminar_Viscosity_i * HLGradTerm, -0.25), 200.0); + /*--- Cal H12, dNdRet, Ret0, D_H12, l_H12, m_H12, kv ---*/ const su2double H12 = TransCorrelations.H12_Correlations(HL); const su2double dNdRet = TransCorrelations.dNdRet_Correlations(H12); @@ -431,15 +442,12 @@ class CSourcePieceWise_TransAFT final : public CNumerics { const su2double l_H12 = TransCorrelations.l_Correlations(H12); const su2double m_H12 = TransCorrelations.m_Correlations(H12, l_H12); const su2double kv = TransCorrelations.kv_Correlations(H12); - const su2double Rev = Density_i * dist_i * dist_i * StrainMag_i / (Laminar_Viscosity_i + Eddy_Viscosity_i); + const su2double Rev = Density_i * dist_i * dist_i * StrainMag_i / (Laminar_Viscosity_i + Eddy_Viscosity_i) * sqrt(2.0); const su2double Rev0 = kv * Ret0; - if(H12 != 2.2 && H12 != 20.0 && dist_i < 0.01){ - su2double tt2 = 0; - } /*--- production term of the amplification factor ---*/ - const su2double F_growth = max(D_H12 * (1.0 + m_H12) / 2.0 * l_H12, 0.0); + const su2double F_growth = D_H12 * (1.0 + m_H12) / 2.0 * l_H12; const su2double F_crit = (Rev < Rev0) ? 0.0 : 1.0; const su2double PAF = Density_i * VorticityMag * F_crit * F_growth * dNdRet; diff --git a/SU2_CFD/src/solvers/CTransAFTSolver.cpp b/SU2_CFD/src/solvers/CTransAFTSolver.cpp index b4ac57bec82..8449b1986c1 100644 --- a/SU2_CFD/src/solvers/CTransAFTSolver.cpp +++ b/SU2_CFD/src/solvers/CTransAFTSolver.cpp @@ -200,8 +200,6 @@ void CTransAFTSolver::Postprocessing(CGeometry *geometry, CSolver **solver_conta const su2double c_1 = 100.0; const su2double c_2 = 0.06; const su2double c_3 = 50.0; - su2double Temp1 = nodes->GetAuxVarGradient(iPoint, 0, 0); - su2double Temp2 = nodes->GetAuxVarGradient(iPoint, 0, 1); su2double Temp3 = flowNodes->GetVelocity(iPoint, 0) * nodes->GetAuxVarGradient(iPoint, 0, 0); su2double Temp4 = flowNodes->GetVelocity(iPoint, 1) * nodes->GetAuxVarGradient(iPoint, 0, 1); @@ -285,7 +283,6 @@ void CTransAFTSolver::Source_Residual(CGeometry *geometry, CSolver **solver_cont const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); auto* flowNodes = su2staticcast_p(solver_container[FLOW_SOL]->GetNodes()); - //auto* turbNodes = su2staticcast_p(solver_container[TURB_SOL]->GetNodes()); CVariable* turbNodes = solver_container[TURB_SOL]->GetNodes(); /*--- Pick one numerics object per thread. ---*/ @@ -337,14 +334,6 @@ void CTransAFTSolver::Source_Residual(CGeometry *geometry, CSolver **solver_cont numerics->SetCoord(geometry->nodes->GetCoord(iPoint), nullptr); numerics->SetAuxVarGrad(nodes->GetAuxVarGradient(iPoint), nullptr); - su2double Temp1 = flowNodes->GetVelocity(iPoint, 0) * nodes->GetAuxVarGradient(iPoint, 0, 0); - su2double Temp2 = flowNodes->GetVelocity(iPoint, 1) * nodes->GetAuxVarGradient(iPoint, 0, 1); - su2double Temp3 = nodes->GetAuxVarGradient(iPoint, 0, 0); - su2double Temp4 = nodes->GetAuxVarGradient(iPoint, 0, 1); - su2double Temp5 = nodes->GetAuxVarGradient(iPoint, 1, 0); - su2double Temp6 = nodes->GetAuxVarGradient(iPoint, 1, 1); - - /*--- Compute the source term ---*/ auto residual = numerics->ComputeResidual(config);