diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index 4439410efff..7d8882b1759 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -902,6 +902,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. */ + A_r_FreeStream, /*!< \brief Freestream Roughness Amplification factor (for LMROUGH transition model) of the fluid. */ 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. */ @@ -2026,6 +2027,12 @@ class CConfig { * \return Freestream momentum thickness Reynolds number. */ su2double GetReThetaT_FreeStream() const { return ReThetaT_FreeStream; } + + /*! + * \brief Get the value of the freestream amplification factor for LMROUGH transition model. + * \return Freestream roughness amplification factor. + */ + su2double GetA_r_FreeStream(void) const { return A_r_FreeStream; } /*! * \brief Get the value of the non-dimensionalized freestream turbulence intensity. diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index bf99257b138..d9e838af592 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -1225,6 +1225,7 @@ enum class LM_OPTIONS { MEDIDA_BAEDER,/*!< \brief Kind of transition correlation model (Medida-Baeder). */ MEDIDA, /*!< \brief Kind of transition correlation model (Medida). */ MENTER_LANGTRY, /*!< \brief Kind of transition correlation model (Menter-Langtry). */ + LMROUGH, /*!< \brief Kind of roughness induced transition model. */ DEFAULT /*!< \brief Kind of transition correlation model (Menter-Langtry if SST, MALAN if SA). */ }; @@ -1237,6 +1238,7 @@ static const MapType LM_Options_Map = { MakePair("KRAUSE_HYPER", LM_OPTIONS::KRAUSE_HYPER) MakePair("MEDIDA_BAEDER", LM_OPTIONS::MEDIDA_BAEDER) MakePair("MENTER_LANGTRY", LM_OPTIONS::MENTER_LANGTRY) + MakePair("LMROUGH", LM_OPTIONS::LMROUGH) MakePair("DEFAULT", LM_OPTIONS::DEFAULT) }; @@ -1260,6 +1262,7 @@ enum class TURB_TRANS_CORRELATION { struct LM_ParsedOptions { LM_OPTIONS version = LM_OPTIONS::NONE; /*!< \brief LM base model. */ bool LM2015 = false; /*!< \brief Use cross-flow corrections. */ + bool LMROUGH = false; /*!< \brief Use roughness induced transition version. */ TURB_TRANS_CORRELATION Correlation = TURB_TRANS_CORRELATION::DEFAULT; }; @@ -1279,6 +1282,7 @@ inline LM_ParsedOptions ParseLMOptions(const LM_OPTIONS *LM_Options, unsigned sh }; LMParsedOptions.LM2015 = IsPresent(LM_OPTIONS::LM2015); + LMParsedOptions.LMROUGH = IsPresent(LM_OPTIONS::LMROUGH); int NFoundCorrelations = 0; if (IsPresent(LM_OPTIONS::MALAN)) { diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 3cadb5be5ef..1f46e5a8472 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -6414,6 +6414,8 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { cout << "Transition model: Langtry and Menter's 4 equation model"; if (lmParsedOptions.LM2015) { cout << " w/ cross-flow corrections (2015)" << endl; + } else if (lmParsedOptions.LMROUGH) { + cout << " w/ roughness-induced transition effects " << endl; } else { cout << " (2009)" << endl; } diff --git a/SU2_CFD/include/numerics/scalar/scalar_convection.hpp b/SU2_CFD/include/numerics/scalar/scalar_convection.hpp index 91b8429e0aa..a129fd49660 100644 --- a/SU2_CFD/include/numerics/scalar/scalar_convection.hpp +++ b/SU2_CFD/include/numerics/scalar/scalar_convection.hpp @@ -47,7 +47,7 @@ template class CUpwScalar : public CNumerics { protected: - enum : unsigned short {MAXNVAR = 8}; + enum : unsigned short {MAXNVAR = 16}; const FlowIndices idx; /*!< \brief Object to manage the access to the flow primitives. */ su2double a0 = 0.0; /*!< \brief The maximum of the face-normal velocity and 0. */ diff --git a/SU2_CFD/include/numerics/scalar/scalar_diffusion.hpp b/SU2_CFD/include/numerics/scalar/scalar_diffusion.hpp index 13ee01f7eda..bbae37e7a82 100644 --- a/SU2_CFD/include/numerics/scalar/scalar_diffusion.hpp +++ b/SU2_CFD/include/numerics/scalar/scalar_diffusion.hpp @@ -59,7 +59,7 @@ struct CNoFlowIndices { template class CAvgGrad_Scalar : public CNumerics { protected: - enum : unsigned short {MAXNVAR = 8}; + enum : unsigned short {MAXNVAR = 16}; const FlowIndices idx; /*!< \brief Object to manage the access to the flow primitives. */ su2double Proj_Mean_GradScalarVar[MAXNVAR]; /*!< \brief Mean_gradScalarVar DOT normal, corrected if required. */ diff --git a/SU2_CFD/include/numerics/turbulent/transition/trans_convection.hpp b/SU2_CFD/include/numerics/turbulent/transition/trans_convection.hpp index 2e0368e3d2e..b10a3fe0bfc 100644 --- a/SU2_CFD/include/numerics/turbulent/transition/trans_convection.hpp +++ b/SU2_CFD/include/numerics/turbulent/transition/trans_convection.hpp @@ -38,3 +38,11 @@ template using CUpwSca_TransLM = CUpwSca_TurbSST; +/*! + * \class CUpwSca_TransLMROUGH + * \brief Use modified convective fluxes for the scalar upwind discretization of LMROUGH transition model equations to handle 3 variables. + * \ingroup ConvDiscr + */ +template +using CUpwSca_TransLMROUGH = CUpwSca_TurbLMROUGH; + diff --git a/SU2_CFD/include/numerics/turbulent/transition/trans_diffusion.hpp b/SU2_CFD/include/numerics/turbulent/transition/trans_diffusion.hpp index 3c4ce9b6fae..d7bfd64535c 100644 --- a/SU2_CFD/include/numerics/turbulent/transition/trans_diffusion.hpp +++ b/SU2_CFD/include/numerics/turbulent/transition/trans_diffusion.hpp @@ -103,3 +103,86 @@ class CAvgGrad_TransLM final : public CAvgGrad_Scalar { } }; + +/*! + * \class CAvgGrad_TransLMROUGH + * \brief Class for computing viscous term using average of gradient with correction (LMROUGH transition model). + * \ingroup ViscDiscr + * \author M. Schifone. + */ +template +class CAvgGrad_TransLMROUGH 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 LMROUGH transition model specific steps in the ComputeResidual method + * \param[in] config - Definition of the particular problem. + */ + void FinishResidualCalc(const CConfig* config) override { + const bool implicit = config->GetKind_TimeIntScheme() == EULER_IMPLICIT; + + /*--- Compute mean effective dynamic viscosity ---*/ + const su2double diff_i_gamma = Laminar_Viscosity_i + Eddy_Viscosity_i; + const su2double diff_j_gamma = Laminar_Viscosity_j + Eddy_Viscosity_j; + const su2double diff_i_ReThetaT = 2.0*(Laminar_Viscosity_i + Eddy_Viscosity_i); + const su2double diff_j_ReThetaT = 2.0*(Laminar_Viscosity_j + Eddy_Viscosity_j); + // Add diffusion terms in the LMROUGH model for A_r transport equation + const su2double diff_i_A_r = 10.0*(Laminar_Viscosity_i + Eddy_Viscosity_i); + const su2double diff_j_A_r = 10.0*(Laminar_Viscosity_j + Eddy_Viscosity_j); + + const su2double diff_gamma = 0.5*(diff_i_gamma + diff_j_gamma); + const su2double diff_ReThetaT = 0.5*(diff_i_ReThetaT + diff_j_ReThetaT); + //Add computation of diffusion term for LMROUGH + const su2double diff_A_r = 0.5*(diff_i_A_r + diff_j_A_r); + + Flux[0] = diff_gamma*Proj_Mean_GradScalarVar[0]; + Flux[1] = diff_ReThetaT*Proj_Mean_GradScalarVar[1]; + Flux[2] = diff_A_r*Proj_Mean_GradScalarVar[2]; + + /*--- For Jacobians -> Use of TSL (Thin Shear Layer) approx. to compute derivatives of the gradients ---*/ + //modifications for the LMROUGH model + if (implicit) { + const su2double proj_on_rho_i = proj_vector_ij/Density_i; + Jacobian_i[0][0] = -diff_gamma*proj_on_rho_i; Jacobian_i[0][1] = 0.0; Jacobian_i[0][2] = 0.0; + Jacobian_i[1][0] = 0.0; Jacobian_i[1][1] = -diff_ReThetaT*proj_on_rho_i; Jacobian_i[1][2] = 0.0; + Jacobian_i[2][0] = 0.0; Jacobian_i[2][1] = 0.0; Jacobian_i[2][2] = -diff_A_r*proj_on_rho_i; + + const su2double proj_on_rho_j = proj_vector_ij/Density_j; + Jacobian_j[0][0] = diff_gamma*proj_on_rho_j; Jacobian_j[0][1] = 0.0; Jacobian_j[0][2] = 0.0; + Jacobian_j[1][0] = 0.0; Jacobian_j[1][1] = diff_ReThetaT*proj_on_rho_j; Jacobian_j[1][2] = 0.0; + Jacobian_j[2][0] = 0.0; Jacobian_j[2][1] = 0.0; Jacobian_j[2][2] = diff_A_r*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_TransLMROUGH(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 bb15656d739..bacde4c10d5 100644 --- a/SU2_CFD/include/numerics/turbulent/transition/trans_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/transition/trans_sources.hpp @@ -323,3 +323,341 @@ class CSourcePieceWise_TransLM final : public CNumerics { return ResidualType<>(Residual, Jacobian_i, nullptr); } }; + +/*! + * \class CSourcePieceWise_TranLMROUGH + * \brief Class for integrating the source terms of the LMROUGH transition model equations. + * \ingroup SourceDiscr + * \author M. Schifone. + */ +template +class CSourcePieceWise_TransLMROUGH final : public CNumerics { + private: + const FlowIndices idx; /*!< \brief Object to manage the access to the flow primitives. */ + + const LM_ParsedOptions options; + /*--- LM Closure constants ---*/ + const su2double c_e1 = 1.0; + const su2double c_a1 = 2.0; + const su2double c_e2 = 50.0; + const su2double c_a2 = 0.06; + const su2double sigmaf = 1.0; + const su2double s1 = 2.0; + const su2double c_theta = 0.03; + const su2double c_CF = 0.6; + const su2double sigmat = 2.0; + const su2double sigma_A_r = 10.0; + const su2double c_Ar1 = 8.0; + const su2double c_Ar2 = 0.0005; + const su2double c_Ar3 = 2.0; + + + TURB_FAMILY TurbFamily; + su2double hRoughness; + + su2double IntermittencySep = 1.0; + su2double IntermittencyEff = 1.0; + + su2double Residual[3]; + su2double* Jacobian_i[3]; + su2double Jacobian_Buffer[9]; // Static storage for the Jacobian (which needs to be pointer for return type). nb contains the partial derivatives of residual wrt trans variables + + TransLMCorrelations TransCorrelations; + + public: + /*! + * \brief Constructor of the class. + * \param[in] val_nDim - Number of dimensions of the problem. + * \param[in] val_nVar - Number of variables of the problem. + * \param[in] config - Definition of the particular problem. + */ + CSourcePieceWise_TransLMROUGH(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config) + : CNumerics(val_nDim, 3, config), idx(val_nDim, config->GetnSpecies()), options(config->GetLMParsedOptions()){ + /*--- "Allocate" the Jacobian using the static buffer. ---*/ + Jacobian_i[0] = Jacobian_Buffer; + Jacobian_i[1] = Jacobian_Buffer + 3; + Jacobian_i[2] = Jacobian_Buffer + 6; + + TurbFamily = TurbModelFamily(config->GetKind_Turb_Model()); + + hRoughness = config->GethRoughness(); + + TransCorrelations.SetOptions(options); + + } + + /*! + * \brief Residual for source term integration. + * \param[in] config - Definition of the particular problem. + * \return A lightweight const-view (read-only) of the residual/flux and Jacobians. + */ + ResidualType<> ComputeResidual(const CConfig* config) override { + /*--- ScalarVar[0] = k, ScalarVar[0] = w, TransVar[0] = gamma, TransVar[1] = ReThetaT and TransVar[2] = A_r ---*/ + /*--- dU/dx = PrimVar_Grad[1][0] ---*/ + AD::StartPreacc(); + AD::SetPreaccIn(StrainMag_i); + AD::SetPreaccIn(ScalarVar_i, nVar); + AD::SetPreaccIn(ScalarVar_Grad_i, nVar, nDim); + AD::SetPreaccIn(TransVar_i, nVar); + AD::SetPreaccIn(TransVar_Grad_i, nVar, nDim); + AD::SetPreaccIn(Volume); + AD::SetPreaccIn(dist_i); + AD::SetPreaccIn(&V_i[idx.Velocity()], nDim); + AD::SetPreaccIn(PrimVar_Grad_i, nDim + idx.Velocity(), nDim); + AD::SetPreaccIn(Vorticity_i, 3); + + su2double VorticityMag = + sqrt(Vorticity_i[0] * Vorticity_i[0] + Vorticity_i[1] * Vorticity_i[1] + Vorticity_i[2] * Vorticity_i[2]); + + const su2double vel_u = V_i[idx.Velocity()]; + const su2double vel_v = V_i[1 + idx.Velocity()]; + const su2double vel_w = (nDim == 3) ? V_i[2 + idx.Velocity()] : 0.0; + + const su2double Velocity_Mag = max(sqrt(vel_u * vel_u + vel_v * vel_v + vel_w * vel_w), 1e-20); + + AD::SetPreaccIn(V_i[idx.Density()], V_i[idx.LaminarViscosity()], V_i[idx.EddyViscosity()]); + + Density_i = V_i[idx.Density()]; + Laminar_Viscosity_i = V_i[idx.LaminarViscosity()]; + Eddy_Viscosity_i = V_i[idx.EddyViscosity()]; + + Residual[0] = 0.0; + Residual[1] = 0.0; + Residual[2] = 0.0; + Jacobian_i[0][0] = 0.0; + Jacobian_i[0][1] = 0.0; + Jacobian_i[0][2] = 0.0; + Jacobian_i[1][0] = 0.0; + Jacobian_i[1][1] = 0.0; + Jacobian_i[1][2] = 0.0; + Jacobian_i[2][0] = 0.0; + Jacobian_i[2][1] = 0.0; + Jacobian_i[2][2] = 0.0; + + if (dist_i > 1e-10) { + su2double Tu = 1.0; + if (TurbFamily == TURB_FAMILY::KW) Tu = max(100.0 * sqrt(2.0 * ScalarVar_i[0] / 3.0) / Velocity_Mag, 0.027); //computed as in LM class + if (TurbFamily == TURB_FAMILY::SA) Tu = config->GetTurbulenceIntensity_FreeStream() * 100; + + /*--- Roughness influence function for LMROUGH transition model ---*/ + su2double F_Ar1 = 0.0; + const su2double C_Ar = sqrt(c_Ar3 / (3.0 * c_Ar2)); + if (TransVar_i[2] < C_Ar) { + // Cubic growth for A_r < C_Ar + F_Ar1 = c_Ar2 * pow(TransVar_i[2], 3); + } else { + // Linear growth for A_r >= C_Ar + F_Ar1 = c_Ar3 * (TransVar_i[2] - C_Ar) + (c_Ar2 * pow(C_Ar, 3)); + } + + /*--- Corr_RetC correlation*/ + const su2double Corr_Rec = TransCorrelations.ReThetaC_Correlations(Tu, TransVar_i[1]); + + /*--- F_length correlation*/ + const su2double Corr_F_length = TransCorrelations.FLength_Correlations(Tu, TransVar_i[1]); + + /*--- F_length ---*/ + su2double F_length = 0.0; + if (TurbFamily == TURB_FAMILY::KW) { + const su2double r_omega = Density_i * dist_i * dist_i * ScalarVar_i[1] / Laminar_Viscosity_i; + const su2double f_sub = exp(-pow(r_omega / 200.0, 2)); + F_length = Corr_F_length * (1. - f_sub) + 40.0 * f_sub; + } + if (TurbFamily == TURB_FAMILY::SA) F_length = Corr_F_length; + + /*--- F_onset ---*/ + su2double R_t = 1.0; + if (TurbFamily == TURB_FAMILY::KW) R_t = Density_i * ScalarVar_i[0] / Laminar_Viscosity_i / ScalarVar_i[1]; + if (TurbFamily == TURB_FAMILY::SA) R_t = Eddy_Viscosity_i / Laminar_Viscosity_i; + + const su2double Re_v = Density_i * dist_i * dist_i * StrainMag_i / Laminar_Viscosity_i; + const su2double F_onset1 = Re_v / (2.193 * Corr_Rec); + su2double F_onset2 = 1.0; + su2double F_onset3 = 1.0; + if (TurbFamily == TURB_FAMILY::KW) { + F_onset2 = min(max(F_onset1, pow(F_onset1, 4.0)), 2.0); + F_onset3 = max(1.0 - pow(R_t / 2.5, 3.0), 0.0); + } + if (TurbFamily == TURB_FAMILY::SA) { + F_onset2 = min(max(F_onset1, pow(F_onset1, 4.0)), 4.0); + F_onset3 = max(2.0 - pow(R_t / 2.5, 3.0), 0.0); + } + const su2double F_onset = max(F_onset2 - F_onset3, 0.0); + + /*-- Gradient of velocity magnitude ---*/ + + su2double dU_dx = 0.5 / Velocity_Mag * (2. * vel_u * PrimVar_Grad_i[1][0] + 2. * vel_v * PrimVar_Grad_i[2][0]); + if (nDim == 3) dU_dx += 0.5 / Velocity_Mag * (2. * vel_w * PrimVar_Grad_i[3][0]); + + su2double dU_dy = 0.5 / Velocity_Mag * (2. * vel_u * PrimVar_Grad_i[1][1] + 2. * vel_v * PrimVar_Grad_i[2][1]); + if (nDim == 3) dU_dy += 0.5 / Velocity_Mag * (2. * vel_w * PrimVar_Grad_i[3][1]); + + su2double dU_dz = 0.0; + if (nDim == 3) + dU_dz = + 0.5 / Velocity_Mag * + (2. * vel_u * PrimVar_Grad_i[1][2] + 2. * vel_v * PrimVar_Grad_i[2][2] + 2. * vel_w * PrimVar_Grad_i[3][2]); + + su2double du_ds = vel_u / Velocity_Mag * dU_dx + vel_v / Velocity_Mag * dU_dy; + if (nDim == 3) du_ds += vel_w / Velocity_Mag * dU_dz; + + /*-- Calculate blending function f_theta --*/ + su2double time_scale = 500.0 * Laminar_Viscosity_i / Density_i / Velocity_Mag / Velocity_Mag; + if (options.LM2015) + time_scale = min(time_scale, + Density_i * LocalGridLength_i * LocalGridLength_i / (Laminar_Viscosity_i + Eddy_Viscosity_i)); + const su2double theta_bl = TransVar_i[1] * Laminar_Viscosity_i / Density_i / Velocity_Mag; + const su2double delta_bl = 7.5 * theta_bl; + const su2double delta = 50.0 * VorticityMag * dist_i / Velocity_Mag * delta_bl + 1e-20; + + su2double f_wake = 0.0; + if (TurbFamily == TURB_FAMILY::KW) { + const su2double re_omega = Density_i * ScalarVar_i[1] * dist_i * dist_i / Laminar_Viscosity_i; + f_wake = exp(-pow(re_omega / (1.0e+05), 2)); + } + if (TurbFamily == TURB_FAMILY::SA) f_wake = 1.0; + + const su2double var1 = (TransVar_i[0] - 1.0 / c_e2) / (1.0 - 1.0 / c_e2); + const su2double var2 = 1.0 - pow(var1, 2.0); + const su2double f_theta = min(max(f_wake * exp(-pow(dist_i / delta, 4)), var2), 1.0); + const su2double f_turb = exp(-pow(R_t / 4, 4)); + + su2double f_theta_2 = 0.0; + if (options.LM2015) + f_theta_2 = min(f_wake * exp(-pow(dist_i / delta, 4.0)), 1.0); + + /*--- Corr_Ret correlation*/ + const su2double Corr_Ret_lim = 20.0; + su2double f_lambda = 1.0; + + su2double Retheta_Error = 200.0, Retheta_old = 0.0; + su2double lambda = 0.0; + su2double Corr_Ret = 20.0; + + for (int iter = 0; iter < 100; iter++) { + su2double theta = Corr_Ret * Laminar_Viscosity_i / Density_i / Velocity_Mag; + lambda = Density_i * theta * theta / Laminar_Viscosity_i * du_ds; + lambda = min(max(-0.1, lambda), 0.1); + + if (lambda <= 0.0) { + f_lambda = 1. - (-12.986 * lambda - 123.66 * lambda * lambda - 405.689 * lambda * lambda * lambda) * + exp(-pow(Tu / 1.5, 1.5)); + } else { + f_lambda = 1. + 0.275 * (1. - exp(-35. * lambda)) * exp(-Tu / 0.5); + } + + if (Tu <= 1.3) { + Corr_Ret = f_lambda * (1173.51 - 589.428 * Tu + 0.2196 / Tu / Tu); + } else { + Corr_Ret = 331.5 * f_lambda * pow(Tu - 0.5658, -0.671); + } + Corr_Ret = max(Corr_Ret, Corr_Ret_lim); + + Retheta_Error = fabs(Retheta_old - Corr_Ret) / Retheta_old; + + if (Retheta_Error < 0.0000001) { + break; + } + + Retheta_old = Corr_Ret; + } + + /*-- Corr_RetT_SCF Correlations--*/ + su2double ReThetat_SCF = 0.0; + if (options.LM2015) { + su2double VelocityNormalized[3]; + VelocityNormalized[0] = vel_u / Velocity_Mag; + VelocityNormalized[1] = vel_v / Velocity_Mag; + if (nDim == 3) VelocityNormalized[2] = vel_w / Velocity_Mag; + + su2double StreamwiseVort = 0.0; + for (auto iDim = 0u; iDim < nDim; iDim++) { + StreamwiseVort += VelocityNormalized[iDim] * Vorticity_i[iDim]; + } + StreamwiseVort = abs(StreamwiseVort); + + const su2double H_CF = StreamwiseVort * dist_i / Velocity_Mag; + const su2double DeltaH_CF = H_CF * (1.0 + min(Eddy_Viscosity_i / Laminar_Viscosity_i, 0.4)); + const su2double DeltaH_CF_Minus = max(-1.0 * (0.1066 - DeltaH_CF), 0.0); + const su2double DeltaH_CF_Plus = max(0.1066 - DeltaH_CF, 0.0); + const su2double fDeltaH_CF_Minus = 75.0 * tanh(DeltaH_CF_Minus / 0.0125); + const su2double fDeltaH_CF_Plus = 6200 * DeltaH_CF_Plus + 50000 * DeltaH_CF_Plus * DeltaH_CF_Plus; + + const su2double toll = 1e-5; + su2double error = toll + 1.0; + su2double thetat_SCF = 0.0; + su2double rethetat_SCF_old = 20.0; + const int nMax = 100; + + int iter; + for (iter = 0; iter < nMax && error > toll; iter++) { + thetat_SCF = rethetat_SCF_old * Laminar_Viscosity_i / (Density_i * (Velocity_Mag / 0.82)); + thetat_SCF = max(1e-20, thetat_SCF); + + ReThetat_SCF = -35.088 * log(hRoughness / thetat_SCF) + 319.51 + fDeltaH_CF_Plus - fDeltaH_CF_Minus; + + error = abs(ReThetat_SCF - rethetat_SCF_old) / rethetat_SCF_old; + + rethetat_SCF_old = ReThetat_SCF; + } + } + + /*-- production term of Intermeittency(Gamma) --*/ + const su2double Pg = + F_length * c_a1 * Density_i * StrainMag_i * sqrt(F_onset * TransVar_i[0]) * (1.0 - c_e1 * TransVar_i[0]); + + /*-- destruction term of Intermeittency(Gamma) --*/ + const su2double Dg = c_a2 * Density_i * VorticityMag * TransVar_i[0] * f_turb * (c_e2 * TransVar_i[0] - 1.0); + + /*-- production term of ReThetaT modified for LMROUGHm model--*/ + su2double PRethetat = 0.0; + PRethetat = c_theta *( Density_i / time_scale) * ((Corr_Ret - TransVar_i[1]) * (1.0 - f_theta) - ((TransVar_i[1]-30.0)/TransVar_i[1])*F_Ar1); + + /*-- destruction term of ReThetaT --*/ + // It should not be with the minus sign but I put for consistency + su2double DRethetat = 0.0; + if (options.LM2015) { + DRethetat = -c_theta * (Density_i / time_scale) * c_CF * min(ReThetat_SCF - TransVar_i[1], 0.0) * f_theta_2; + } + /* To enable the production term defined in LM solver class*/ + su2double PAr = 0.0; + su2double DAr = 0.0; + + /*--- Source ---*/ + Residual[0] += (Pg - Dg) * Volume; + Residual[1] += (PRethetat - DRethetat) * Volume; + Residual[2] += (PAr - DAr)* Volume; + // No production term for A_r in this model + + /*--- Implicit part ---*/ + Jacobian_i[0][0] = (F_length * c_a1 * StrainMag_i * sqrt(F_onset) * + (0.5 * pow(TransVar_i[0], -0.5) - 1.5 * c_e1 * pow(TransVar_i[0], 0.5)) - + c_a2 * VorticityMag * f_turb * (2.0 * c_e2 * TransVar_i[0] - 1.0)) * + Volume; + Jacobian_i[0][1] = 0.0; + Jacobian_i[0][2] = 0.0; + Jacobian_i[1][0] = 0.0; + Jacobian_i[1][1] = -c_theta / time_scale * (1.0 - f_theta) * Volume; + if (options.LM2015 && ReThetat_SCF - TransVar_i[1] < 0) { + Jacobian_i[1][1] += (c_theta / time_scale) * c_CF * f_theta_2 * Volume; + } + + if (TransVar_i[2] < C_Ar) { + Jacobian_i[1][2] += -(c_theta / time_scale) * (((TransVar_i[1]-30.0)/TransVar_i[1])*3 * c_Ar2 * pow(TransVar_i[2], 2)) * Volume; + } else { + Jacobian_i[1][2] += -(c_theta / time_scale) * (((TransVar_i[1]-30.0)/TransVar_i[1])*c_Ar3) * Volume; + } + + /*--- Jacobian for A_r (Residual[2]) ---*/ + //Since the source term for A_r (Residual[2]) is zero in this model, derivatives are also zero. + Jacobian_i[2][0] = 0.0; + Jacobian_i[2][1] = 0.0; + Jacobian_i[2][2] = 0.0; + + } + AD::SetPreaccOut(Residual, nVar); + AD::EndPreacc(); + + return ResidualType<>(Residual, Jacobian_i, nullptr); + } +}; \ No newline at end of file diff --git a/SU2_CFD/include/numerics/turbulent/turb_convection.hpp b/SU2_CFD/include/numerics/turbulent/turb_convection.hpp index 03c33d38138..769c107867e 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_convection.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_convection.hpp @@ -128,3 +128,70 @@ class CUpwSca_TurbSST final : public CUpwScalar { CUpwSca_TurbSST(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config) : CUpwScalar(val_nDim, val_nVar, config) { bounded_scalar = config->GetBounded_Turb(); } }; +/*! + * \class CUpwSca_TurbLMROUGH + * \brief Class for doing a scalar upwind solver for the LMROUGH transition model equations. + * \ingroup ConvDiscr + * \author M. Schifone. + */ +template +class CUpwSca_TurbLMROUGH final : public CUpwScalar { +private: + using Base = CUpwScalar; + using Base::nDim; + using Base::V_i; + using Base::V_j; + using Base::a0; + using Base::a1; + using Base::Flux; + using Base::Jacobian_i; + using Base::Jacobian_j; + using Base::ScalarVar_i; + using Base::ScalarVar_j; + using Base::TransVar_i; + using Base::TransVar_j; + using Base::idx; + using Base::bounded_scalar; + + /*! + * \brief Adds any extra variables to AD + */ + void ExtraADPreaccIn() override {} + + /*! + * \brief SST specific steps in the ComputeResidual method + * \param[in] config - Definition of the particular problem. + */ + + void FinishResidualCalc(const CConfig* config) override { + + + + + + Flux[0] = a0*V_i[idx.Density()]*ScalarVar_i[0] + a1*V_j[idx.Density()]*ScalarVar_j[0]; + Flux[1] = a0*V_i[idx.Density()]*ScalarVar_i[1] + a1*V_j[idx.Density()]*ScalarVar_j[1]; + Flux[2] = a0*V_i[idx.Density()]*ScalarVar_i[2] + a1*V_j[idx.Density()]*ScalarVar_j[2]; + + + + Jacobian_i[0][0] = a0; Jacobian_i[0][1] = 0.0; Jacobian_i[0][2] = 0.0; + Jacobian_i[1][0] = 0.0; Jacobian_i[1][1] = a0; Jacobian_i[1][2] = 0.0; + Jacobian_i[2][0] = 0.0; Jacobian_i[2][1] = 0.0; Jacobian_i[2][2] = a0; + + Jacobian_j[0][0] = a1; Jacobian_j[0][1] = 0.0; Jacobian_j[0][2] = 0.0; + Jacobian_j[1][0] = 0.0; Jacobian_j[1][1] = a1; Jacobian_j[1][2] = 0.0; + Jacobian_j[2][0] = 0.0; Jacobian_j[2][1] = 0.0; Jacobian_j[2][2] = a1; + + } + +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. + */ + CUpwSca_TurbLMROUGH(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config) + : CUpwScalar(val_nDim, val_nVar, config) { bounded_scalar = config->GetBounded_Turb(); } +}; diff --git a/SU2_CFD/include/solvers/CScalarSolver.hpp b/SU2_CFD/include/solvers/CScalarSolver.hpp index e5daba900bb..315ef6a276b 100644 --- a/SU2_CFD/include/solvers/CScalarSolver.hpp +++ b/SU2_CFD/include/solvers/CScalarSolver.hpp @@ -44,9 +44,9 @@ template class CScalarSolver : public CSolver { protected: enum : size_t { MAXNDIM = 3 }; /*!< \brief Max number of space dimensions, used in some static arrays. */ - static constexpr size_t MAXNVAR = VariableType::MAXNVAR; /*!< \brief Max number of variables, for static arrays. */ + //static constexpr size_t MAXNVAR = VariableType::MAXNVAR; /*!< \brief Max number of variables, for static arrays. */ enum : size_t { MAXNVARFLOW = 12 }; /*!< \brief Max number of flow variables, used in some static arrays. */ - + enum : size_t { MAXNVAR = 10 }; /*!< \brief Max number of flow variables, used in some static arrays. */ enum : size_t { OMP_MAX_SIZE = 512 }; /*!< \brief Max chunk size for light point loops. */ enum : size_t { OMP_MIN_SIZE = 32 }; /*!< \brief Min chunk size for edge loops (max is color group size). */ diff --git a/SU2_CFD/include/solvers/CScalarSolver.inl b/SU2_CFD/include/solvers/CScalarSolver.inl index 1af1f828d86..7810bce777f 100644 --- a/SU2_CFD/include/solvers/CScalarSolver.inl +++ b/SU2_CFD/include/solvers/CScalarSolver.inl @@ -485,6 +485,7 @@ void CScalarSolver::PrepareImplicitIteration(CGeometry* geometry, /*--- "Add" residual at (iPoint,iVar) to local residual variables. ---*/ ResidualReductions_PerThread(iPoint, iVar, LinSysRes[total_index], resRMS, resMax, idxMax); + if ( idxMax[iVar] > nPointDomain ) idxMax[iVar] = nPointDomain-1; } } END_SU2_OMP_FOR diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index 7c38e054d4a..f75c32bb927 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -3169,6 +3169,12 @@ class CSolver { */ inline virtual su2double GetReThetaT_Inf() const { return 0; } + /*! + * \brief Get value of the roughness amplification factor in the freestream. + * \return Value of the roughness amplification factor in the freestream. + */ + inline virtual su2double GetA_r_Inf() const { return 0.0; } + /*! * \brief A virtual member. * \return Value of the sensitivity coefficient for the Young Modulus E diff --git a/SU2_CFD/include/solvers/CTransLMSolver.hpp b/SU2_CFD/include/solvers/CTransLMSolver.hpp index 7f52a072767..bb97693f265 100644 --- a/SU2_CFD/include/solvers/CTransLMSolver.hpp +++ b/SU2_CFD/include/solvers/CTransLMSolver.hpp @@ -205,6 +205,12 @@ class CTransLMSolver final : public CTurbSolver { */ inline su2double GetReThetaT_Inf(void) const override { return Solution_Inf[1]; } + /*! + * \brief Get the value of the roughness amplification factor in the freestream. + * \return Value of the freestream roughness amplification factor. + */ + inline su2double GetA_r_Inf(void) const override { return Solution_Inf[2]; } + /*! * \brief Load a solution from a restart file. * \param[in] geometry - Geometrical definition of the problem. diff --git a/SU2_CFD/include/variables/CTransLMVariable.hpp b/SU2_CFD/include/variables/CTransLMVariable.hpp index bcc3bf94543..f83198cff4e 100644 --- a/SU2_CFD/include/variables/CTransLMVariable.hpp +++ b/SU2_CFD/include/variables/CTransLMVariable.hpp @@ -53,7 +53,7 @@ class CTransLMVariable final : public CTurbVariable { * \param[in] nvar - Number of variables of the problem. * \param[in] config - Definition of the particular problem. */ - CTransLMVariable(su2double Intermittency, su2double ReThetaT, su2double gammaSep, su2double gammaEff, unsigned long npoint, unsigned long ndim, unsigned long nvar, CConfig *config); + CTransLMVariable(su2double Intermittency, su2double ReThetaT, su2double gammaSep, su2double gammaEff, su2double A_r, unsigned long npoint, unsigned long ndim, unsigned long nvar, CConfig *config); /*! * \brief Destructor of the class. diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index dc905d79fcb..6de1e2ef89a 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -1679,6 +1679,18 @@ class CVariable { */ inline virtual su2double GetmuT(unsigned long iPoint) const { return 0.0; } + /*! + * \brief Get the value of the roughness amplification factor for LMROUGH transition model. + * \return the value of the roughness amplification factor. + */ + inline virtual su2double GetA_r(unsigned long iPoint) const { return 0.0; } + + /*! + * \brief Set the roughness amplification factor, modification for the LMROUGH transition model. + * \param[in] val_dist - Value of the A_r. + */ + inline virtual void SetA_r(unsigned long iPoint, su2double val_A_r) {} + /*! * \brief Get the value of the intermittency. * \return the value of the intermittency. diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index 7e17ed9c601..cc1c9a7f781 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -1321,7 +1321,8 @@ void CDriver::InstantiateTransitionNumerics(unsigned short nVar_Trans, int offse const int visc_bound_term = VISC_BOUND_TERM + offset; const bool LM = config->GetKind_Trans_Model() == TURB_TRANS_MODEL::LM; - + LM_ParsedOptions options; + if (LM) options = config->GetLMParsedOptions(); /*--- Definition of the convective scheme for each equation and mesh level ---*/ switch (config->GetKind_ConvNumScheme_Turb()) { @@ -1330,7 +1331,10 @@ void CDriver::InstantiateTransitionNumerics(unsigned short nVar_Trans, int offse break; case SPACE_UPWIND : for (auto iMGlevel = 0u; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { - if (LM) numerics[iMGlevel][TRANS_SOL][conv_term] = new CUpwSca_TransLM(nDim, nVar_Trans, config); + if (LM) { + if (options.LMROUGH) numerics[iMGlevel][TRANS_SOL][conv_term] = new CUpwSca_TransLMROUGH(nDim,nVar_Trans, config); + else numerics[iMGlevel][TRANS_SOL][conv_term] = new CUpwSca_TransLM(nDim, nVar_Trans, config); + } } break; default: @@ -1341,7 +1345,10 @@ void CDriver::InstantiateTransitionNumerics(unsigned short nVar_Trans, int offse /*--- Definition of the viscous scheme for each equation and mesh level ---*/ for (auto iMGlevel = 0u; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { - if (LM) numerics[iMGlevel][TRANS_SOL][visc_term] = new CAvgGrad_TransLM(nDim, nVar_Trans, true, config); + if (LM) { + if (options.LMROUGH) numerics[iMGlevel][TRANS_SOL][visc_term] = new CAvgGrad_TransLMROUGH(nDim, nVar_Trans, true, config); + else numerics[iMGlevel][TRANS_SOL][visc_term] = new CAvgGrad_TransLM(nDim, nVar_Trans, true, config); + } } /*--- Definition of the source term integration scheme for each equation and mesh level ---*/ @@ -1349,15 +1356,21 @@ void CDriver::InstantiateTransitionNumerics(unsigned short nVar_Trans, int offse for (auto iMGlevel = 0u; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { auto& trans_source_first_term = numerics[iMGlevel][TRANS_SOL][source_first_term]; - if (LM) trans_source_first_term = new CSourcePieceWise_TransLM(nDim, nVar_Trans, config); - + if (LM) { + if (options.LMROUGH) trans_source_first_term = new CSourcePieceWise_TransLMROUGH(nDim, nVar_Trans, config); + else trans_source_first_term = new CSourcePieceWise_TransLM(nDim, nVar_Trans, config); + } numerics[iMGlevel][TRANS_SOL][source_second_term] = new CSourceNothing(nDim, nVar_Trans, config); } /*--- Definition of the boundary condition method ---*/ for (auto iMGlevel = 0u; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { - if (LM) { + if (options.LMROUGH) { + numerics[iMGlevel][TRANS_SOL][conv_bound_term] = new CUpwSca_TransLMROUGH(nDim, nVar_Trans, config); + numerics[iMGlevel][TRANS_SOL][visc_bound_term] = new CAvgGrad_TransLMROUGH(nDim, nVar_Trans, false, config); + } + else { 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); } diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp index b355c0698af..481d52ef1c4 100644 --- a/SU2_CFD/src/output/CFlowOutput.cpp +++ b/SU2_CFD/src/output/CFlowOutput.cpp @@ -980,8 +980,15 @@ void CFlowOutput::AddHistoryOutputFields_ScalarRMS_RES(const CConfig* config) { case TURB_TRANS_MODEL::LM: /// DESCRIPTION: Root-mean square residual of the intermittency (LM model). AddHistoryOutput("RMS_INTERMITTENCY", "rms[LM_1]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of intermittency (LM model).", HistoryFieldType::RESIDUAL); - /// DESCRIPTION: Root-mean square residual of the momentum thickness Reynolds number (LM model). - AddHistoryOutput("RMS_RE_THETA_T", "rms[LM_2]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of momentum thickness Reynolds number (LM model).", HistoryFieldType::RESIDUAL); + if ((config->GetLMParsedOptions()).LMROUGH) { + /// DESCRIPTION: Root-mean square residual of the momentum thickness Reynolds number (LM model). + AddHistoryOutput("RMS_RE_THETA_T", "rms[LM_2]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of momentum thickness Reynolds number (LM model).", HistoryFieldType::RESIDUAL); + /// DESCRIPTION: Root-mean square residual of the roughness amplification factor (LMROUGH model). + AddHistoryOutput("RMS_A_R", "rms[LM_3]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of momentum thickness Reynolds number (LM model).", HistoryFieldType::RESIDUAL); + } else { + /// DESCRIPTION: Root-mean square residual of the momentum thickness Reynolds number (LM model). + AddHistoryOutput("RMS_RE_THETA_T", "rms[LM_2]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of momentum thickness Reynolds number (LM model).", HistoryFieldType::RESIDUAL); + } break; case TURB_TRANS_MODEL::NONE: break; @@ -1037,8 +1044,15 @@ void CFlowOutput::AddHistoryOutputFields_ScalarMAX_RES(const CConfig* config) { case TURB_TRANS_MODEL::LM: /// DESCRIPTION: Maximum residual of the intermittency (LM model). AddHistoryOutput("MAX_INTERMITTENCY", "max[LM_1]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of the intermittency (LM model).", HistoryFieldType::RESIDUAL); - /// DESCRIPTION: Maximum residual of the momentum thickness Reynolds number (LM model). - AddHistoryOutput("MAX_RE_THETA_T", "max[LM_2]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of the momentum thickness Reynolds number (LM model).", HistoryFieldType::RESIDUAL); + if ((config->GetLMParsedOptions()).LMROUGH) { + /// DESCRIPTION: Maximum residual of the momentum thickness Reynolds number (LM model). + AddHistoryOutput("MAX_RE_THETA_T", "max[LM_2]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of the momentum thickness Reynolds number (LM model).", HistoryFieldType::RESIDUAL); + /// DESCRIPTION: Maximum residual of the roughness amplification factor (LMROUGH model). + AddHistoryOutput("MAX_A_R", "max[LM_3]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of the roughness amplification factor (LMROUGH model).", HistoryFieldType::RESIDUAL); + } else { + /// DESCRIPTION: Maximum residual of the momentum thickness Reynolds number (LM model). + AddHistoryOutput("MAX_RE_THETA_T", "max[LM_2]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of the momentum thickness Reynolds number (LM model).", HistoryFieldType::RESIDUAL); + } break; case TURB_TRANS_MODEL::NONE: @@ -1094,8 +1108,15 @@ void CFlowOutput::AddHistoryOutputFields_ScalarBGS_RES(const CConfig* config) { case TURB_TRANS_MODEL::LM: /// DESCRIPTION: Maximum residual of the intermittency (LM model). AddHistoryOutput("BGS_INTERMITTENCY", "bgs[LM_1]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the intermittency (LM model).", HistoryFieldType::RESIDUAL); - /// DESCRIPTION: Maximum residual of the momentum thickness Reynolds number (LM model). - AddHistoryOutput("BGS_RE_THETA_T", "bgs[LM_2]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the momentum thickness Reynolds number (LM model).", HistoryFieldType::RESIDUAL); + if ((config->GetLMParsedOptions()).LMROUGH) { + /// DESCRIPTION: Maximum residual of the momentum thickness Reynolds number (LM model). + AddHistoryOutput("BGS_RE_THETA_T", "bgs[LM_2]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the momentum thickness Reynolds number (LM model).", HistoryFieldType::RESIDUAL); + /// DESCRIPTION: Maximum residual of the roughness amplification factor (LMROUGH model). + AddHistoryOutput("BGS_A_R", "bgs[LM_3]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the roughness amplification factor (LMROUGH model).", HistoryFieldType::RESIDUAL); + } else { + /// DESCRIPTION: Maximum residual of the momentum thickness Reynolds number (LM model). + AddHistoryOutput("BGS_RE_THETA_T", "bgs[LM_2]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the momentum thickness Reynolds number (LM model).", HistoryFieldType::RESIDUAL); + } break; case TURB_TRANS_MODEL::NONE: break; @@ -1187,12 +1208,24 @@ void CFlowOutput::LoadHistoryDataScalar(const CConfig* config, const CSolver* co switch (config->GetKind_Trans_Model()) { case TURB_TRANS_MODEL::LM: SetHistoryOutputValue("RMS_INTERMITTENCY", log10(solver[TRANS_SOL]->GetRes_RMS(0))); - SetHistoryOutputValue("RMS_RE_THETA_T",log10(solver[TRANS_SOL]->GetRes_RMS(1))); SetHistoryOutputValue("MAX_INTERMITTENCY", log10(solver[TRANS_SOL]->GetRes_Max(0))); - SetHistoryOutputValue("MAX_RE_THETA_T", log10(solver[TRANS_SOL]->GetRes_Max(1))); + if ((config->GetLMParsedOptions()).LMROUGH) { + SetHistoryOutputValue("RMS_RE_THETA_T", log10(solver[TRANS_SOL]->GetRes_RMS(1))); + SetHistoryOutputValue("MAX_RE_THETA_T", log10(solver[TRANS_SOL]->GetRes_Max(1))); + SetHistoryOutputValue("RMS_A_R", log10(solver[TRANS_SOL]->GetRes_RMS(2))); + SetHistoryOutputValue("MAX_A_R", log10(solver[TRANS_SOL]->GetRes_Max(2))); + } else { + SetHistoryOutputValue("RMS_RE_THETA_T",log10(solver[TRANS_SOL]->GetRes_RMS(1))); + SetHistoryOutputValue("MAX_RE_THETA_T", log10(solver[TRANS_SOL]->GetRes_Max(1))); + } if (multiZone) { SetHistoryOutputValue("BGS_INTERMITTENCY", log10(solver[TRANS_SOL]->GetRes_BGS(0))); - SetHistoryOutputValue("BGS_RE_THETA_T", log10(solver[TRANS_SOL]->GetRes_BGS(1))); + if ((config->GetLMParsedOptions()).LMROUGH) { + SetHistoryOutputValue("BGS_RE_THETA_T", log10(solver[TRANS_SOL]->GetRes_BGS(1))); + SetHistoryOutputValue("BGS_A_R", log10(solver[TRANS_SOL]->GetRes_BGS(2))); + } else { + SetHistoryOutputValue("BGS_RE_THETA_T", log10(solver[TRANS_SOL]->GetRes_BGS(1))); + } } SetHistoryOutputValue("LINSOL_ITER_TRANS", solver[TRANS_SOL]->GetIterLinSolver()); SetHistoryOutputValue("LINSOL_RESIDUAL_TRANS", log10(solver[TRANS_SOL]->GetResLinSolver())); @@ -1266,6 +1299,9 @@ void CFlowOutput::SetVolumeOutputFieldsScalarSolution(const CConfig* config){ case TURB_TRANS_MODEL::LM: AddVolumeOutput("INTERMITTENCY", "LM_gamma", "SOLUTION", "LM intermittency"); AddVolumeOutput("RE_THETA_T", "LM_Re_t", "SOLUTION", "LM RE_THETA_T"); + if ((config->GetLMParsedOptions()).LMROUGH){ + AddVolumeOutput("A_R", "A_r", "SOLUTION", "LM A_r"); + } break; case TURB_TRANS_MODEL::NONE: @@ -1341,7 +1377,12 @@ void CFlowOutput::SetVolumeOutputFieldsScalarResidual(const CConfig* config) { switch (config->GetKind_Trans_Model()) { case TURB_TRANS_MODEL::LM: AddVolumeOutput("RES_INTERMITTENCY", "Residual_LM_intermittency", "RESIDUAL", "Residual of LM intermittency"); - AddVolumeOutput("RES_RE_THETA_T", "Residual_LM_RE_THETA_T", "RESIDUAL", "Residual of LM RE_THETA_T"); + if ((config->GetLMParsedOptions()).LMROUGH) { + AddVolumeOutput("RES_RE_THETA_T", "Residual_LM_RE_THETA_T", "RESIDUAL", "Residual of LM RE_THETA_T"); + AddVolumeOutput("RES_A_R", "Residual_LM_A_R", "RESIDUAL", "Residual of LM A_r"); + } else { + AddVolumeOutput("RES_RE_THETA_T", "Residual_LM_RE_THETA_T", "RESIDUAL", "Residual of LM RE_THETA_T"); + } break; case TURB_TRANS_MODEL::NONE: @@ -1568,12 +1609,22 @@ void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* con switch (config->GetKind_Trans_Model()) { case TURB_TRANS_MODEL::LM: SetVolumeOutputValue("INTERMITTENCY", iPoint, Node_Trans->GetSolution(iPoint, 0)); - SetVolumeOutputValue("RE_THETA_T", iPoint, Node_Trans->GetSolution(iPoint, 1)); + if ((config->GetLMParsedOptions()).LMROUGH) { + SetVolumeOutputValue("RE_THETA_T", iPoint, Node_Trans->GetSolution(iPoint, 1)); + SetVolumeOutputValue("A_R", iPoint, Node_Trans->GetSolution(iPoint, 2)); + } else { + SetVolumeOutputValue("RE_THETA_T", iPoint, Node_Trans->GetSolution(iPoint, 1)); + } SetVolumeOutputValue("INTERMITTENCY_SEP", iPoint, Node_Trans->GetIntermittencySep(iPoint)); SetVolumeOutputValue("INTERMITTENCY_EFF", iPoint, Node_Trans->GetIntermittencyEff(iPoint)); SetVolumeOutputValue("TURB_INDEX", iPoint, Node_Turb->GetTurbIndex(iPoint)); SetVolumeOutputValue("RES_INTERMITTENCY", iPoint, trans_solver->LinSysRes(iPoint, 0)); - SetVolumeOutputValue("RES_RE_THETA_T", iPoint, trans_solver->LinSysRes(iPoint, 1)); + if ((config->GetLMParsedOptions()).LMROUGH) { + SetVolumeOutputValue("RES_RE_THETA_T", iPoint, trans_solver->LinSysRes(iPoint, 1)); + SetVolumeOutputValue("RES_A_R", iPoint, trans_solver->LinSysRes(iPoint, 2)); + } else { + SetVolumeOutputValue("RES_RE_THETA_T", iPoint, trans_solver->LinSysRes(iPoint, 1)); + } break; case TURB_TRANS_MODEL::NONE: break; diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 80f905ec2d3..79be29e3928 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -1372,6 +1372,10 @@ void CEulerSolver::SetNondimensionalization(CConfig *config, unsigned short iMes Unit.str(""); NonDimTable << "Moment. Thick. Re" << "-" << "-" << "-" << config->GetReThetaT_FreeStream(); Unit.str(""); + if ((config->GetLMParsedOptions()).LMROUGH) { + NonDimTable << "Amplification Factor" << "-" << "-" << "-" << config->GetA_r_FreeStream(); + Unit.str(""); + } } } if (config->GetKind_Species_Model() != SPECIES_MODEL::NONE) { diff --git a/SU2_CFD/src/solvers/CTransLMSolver.cpp b/SU2_CFD/src/solvers/CTransLMSolver.cpp index e8a6dcf389d..68e94805e89 100644 --- a/SU2_CFD/src/solvers/CTransLMSolver.cpp +++ b/SU2_CFD/src/solvers/CTransLMSolver.cpp @@ -49,6 +49,13 @@ CTransLMSolver::CTransLMSolver(CGeometry *geometry, CConfig *config, unsigned sh /*--- Dimension of the problem --> 2 Transport equations (intermittency, Reth) ---*/ nVar = 2; nPrimVar = 2; + + options = config->GetLMParsedOptions(); + if (options.LMROUGH) { + nVar = 3; + nPrimVar = 3; + } + nPoint = geometry->GetnPoint(); nPointDomain = geometry->GetnPointDomain(); @@ -101,14 +108,21 @@ CTransLMSolver::CTransLMSolver(CGeometry *geometry, CConfig *config, unsigned sh lowerlimit[0] = 1.0e-4; upperlimit[0] = 5.0; - lowerlimit[1] = 1.0e-4; + lowerlimit[1] = 20.0; //Suggested lower limit to ensure robustness from Langtry Menter reference upperlimit[1] = 1.0e15; + if (options.LMROUGH) { + lowerlimit[1] = 20.0; + upperlimit[1] = 1.0e15; + lowerlimit[2] = 1.0e-4; + upperlimit[2] = 1.0e4; + } /*--- Far-field flow state quantities and initialization. ---*/ const su2double Intensity = config->GetTurbulenceIntensity_FreeStream()*100.0; const su2double Intermittency_Inf = 1.0; su2double ReThetaT_Inf = 100.0; + const su2double A_r_Inf = 1.0e-4; /*--- Momentum thickness Reynolds number, initialized from freestream turbulent intensity*/ if (Intensity <= 1.3) { @@ -125,9 +139,12 @@ CTransLMSolver::CTransLMSolver(CGeometry *geometry, CConfig *config, unsigned sh Solution_Inf[0] = Intermittency_Inf; Solution_Inf[1] = ReThetaT_Inf; + if (options.LMROUGH) { + Solution_Inf[2] = A_r_Inf; + } /*--- Initialize the solution to the far-field state everywhere. ---*/ - nodes = new CTransLMVariable(Intermittency_Inf, ReThetaT_Inf, 1.0, 1.0, nPoint, nDim, nVar, config); + nodes = new CTransLMVariable(Intermittency_Inf, ReThetaT_Inf, 1.0, 1.0, A_r_Inf, nPoint, nDim, nVar, config); SetBaseClassPointerToNodes(); /*--- MPI solution ---*/ @@ -156,6 +173,9 @@ CTransLMSolver::CTransLMSolver(CGeometry *geometry, CConfig *config, unsigned sh for (unsigned long iVertex = 0; iVertex < nVertex[iMarker]; ++iVertex) { Inlet_TurbVars[iMarker](iVertex,0) = Intermittency_Inf; Inlet_TurbVars[iMarker](iVertex,1) = ReThetaT_Inf; + if (options.LMROUGH) { + Inlet_TurbVars[iMarker](iVertex, 2) = A_r_Inf; + } } } @@ -368,6 +388,12 @@ void CTransLMSolver::BC_HeatFlux_Wall(CGeometry *geometry, CSolver **solver_cont const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); + bool rough_wall = false; + string Marker_Tag = config->GetMarker_All_TagBound(val_marker); + WALL_TYPE WallType; su2double Roughness_Height; + tie(WallType, Roughness_Height) = config->GetWallRoughnessProperties(Marker_Tag); + if (WallType == WALL_TYPE::ROUGH) rough_wall = true; + SU2_OMP_FOR_STAT(OMP_MIN_SIZE) for (auto iVertex = 0u; iVertex < geometry->nVertex[val_marker]; iVertex++) { @@ -412,7 +438,50 @@ void CTransLMSolver::BC_HeatFlux_Wall(CGeometry *geometry, CSolver **solver_cont /*--- Add residuals and Jacobians ---*/ LinSysRes.AddBlock(iPoint, residual); - if (implicit) Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i); + //if (implicit) Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i); + + options = config->GetLMParsedOptions(); + if (rough_wall) { + + /*--- Set wall values ---*/ + su2double density = solver_container[FLOW_SOL]->GetNodes()->GetDensity(iPoint); + su2double laminar_viscosity = solver_container[FLOW_SOL]->GetNodes()->GetLaminarViscosity(iPoint); + su2double WallShearStress = solver_container[FLOW_SOL]->GetWallShearStress(val_marker, iVertex); + + /*--- Compute non-dimensional velocity ---*/ + su2double FrictionVel = sqrt(fabs(WallShearStress)/density); + + /*--- Compute roughness in wall units. ---*/ + su2double kPlus = FrictionVel*Roughness_Height*density/laminar_viscosity; + su2double Solution[3]; + Solution[0] = nodes->GetSolution(iPoint, 0); // Intermittency remains unchanged + Solution[1] = nodes->GetSolution(iPoint, 1); // ReThetaT remains unchanged + Solution[2] = 8.0*kPlus; // Set A_r explicitly + + // Apply the solution values at the wall + + nodes->SetSolution_Old(iPoint, Solution); + nodes->SetSolution(iPoint, Solution); + //LinSysRes.SetBlock_Zero(iPoint); + + LinSysRes.AddBlock(iPoint, residual); + + + + + if (implicit) { + Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i); + } + + } else { + /*--- 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 @@ -464,6 +533,9 @@ void CTransLMSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, C su2double Inlet_Vars[MAXNVAR]; Inlet_Vars[0] = Inlet_TurbVars[val_marker][iVertex][0]; Inlet_Vars[1] = Inlet_TurbVars[val_marker][iVertex][1]; + if (options.LMROUGH) { + Inlet_Vars[2] = Inlet_TurbVars[val_marker][iVertex][2]; + } 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)); diff --git a/SU2_CFD/src/variables/CTransLMVariable.cpp b/SU2_CFD/src/variables/CTransLMVariable.cpp index c1c09ecd698..e4521d52df0 100644 --- a/SU2_CFD/src/variables/CTransLMVariable.cpp +++ b/SU2_CFD/src/variables/CTransLMVariable.cpp @@ -28,15 +28,24 @@ #include "../../include/variables/CTransLMVariable.hpp" -CTransLMVariable::CTransLMVariable(su2double Intermittency, su2double ReThetaT, su2double gammaSep, su2double gammaEff, unsigned long npoint, unsigned long ndim, unsigned long nvar, CConfig *config) +CTransLMVariable::CTransLMVariable(su2double Intermittency, su2double ReThetaT, su2double gammaSep, su2double gammaEff, su2double A_r, unsigned long npoint, unsigned long ndim, unsigned long nvar, CConfig *config) : CTurbVariable(npoint, ndim, nvar, config) { + LM_ParsedOptions options = config->GetLMParsedOptions(); + if (options.LMROUGH) { for(unsigned long iPoint=0; iPoint