From d90e8e14ddda3582f782043fba9a00e3d3f4edcf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Corentin=20Thom=C3=A9e?= Date: Fri, 12 Jan 2024 12:21:17 +0100 Subject: [PATCH 1/3] Implemented SA axisymmetric source terms --- Common/src/CConfig.cpp | 5 -- .../numerics/turbulent/turb_sources.hpp | 49 +++++++++++++++++++ SU2_CFD/src/solvers/CTurbSASolver.cpp | 7 +++ 3 files changed, 56 insertions(+), 5 deletions(-) diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index f660ec5633d..f8b2e15dc78 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -3448,11 +3448,6 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i saParsedOptions = ParseSAOptions(SA_Options, nSA_Options, rank); } - /*--- Check if turbulence model can be used for AXISYMMETRIC case---*/ - if (Axisymmetric && Kind_Turb_Model != TURB_MODEL::NONE && Kind_Turb_Model != TURB_MODEL::SST){ - SU2_MPI::Error("Axisymmetry is currently only supported for KIND_TURB_MODEL chosen as SST", CURRENT_FUNCTION); - } - /*--- Postprocess LM_OPTIONS into structure. ---*/ if (Kind_Trans_Model == TURB_TRANS_MODEL::LM) { lmParsedOptions = ParseLMOptions(LM_Options, nLM_Options, rank, Kind_Turb_Model); diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp index af920cd34c1..2382f84388d 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp @@ -77,9 +77,41 @@ class CSourceBase_TurbSA : public CNumerics { const FlowIndices idx; /*!< \brief Object to manage the access to the flow primitives. */ const SA_ParsedOptions options; /*!< \brief Struct with SA options. */ + const bool axisymmetric = false; bool transition_LM; + /*! + * \brief Add contribution from diffusion due to axisymmetric formulation to 2D residual + */ + inline void ResidualAxisymmetricDiffusion(su2double sigma) { + if (Coord_i[1] < EPS) return; + + const su2double yinv = 1.0 / Coord_i[1]; + const su2double& nue = ScalarVar_i[0]; + + const auto& density = V_i[idx.Density()]; + const auto& laminar_viscosity = V_i[idx.LaminarViscosity()]; + + const su2double nu = laminar_viscosity/density; + + su2double nu_e; + + if (options.version == SA_OPTIONS::NEG && nue < 0.0) { + const su2double cn1 = 16.0; + const su2double Xi = nue / nu; + const su2double fn = (cn1 + Xi*Xi*Xi) / (cn1 - Xi*Xi*Xi); + nu_e = nu + fn * nue; + } else { + nu_e = nu + nue; + } + + /* Diffusion source term */ + const su2double dv_axi = (1.0/sigma)*nu_e*ScalarVar_Grad_i[0][1]; + + Residual += yinv * dv_axi * Volume; + } + public: /*! * \brief Constructor of the class. @@ -89,6 +121,7 @@ class CSourceBase_TurbSA : public CNumerics { CSourceBase_TurbSA(unsigned short nDim, const CConfig* config) : CNumerics(nDim, 1, config), idx(nDim, config->GetnSpecies()), + axisymmetric(config->GetAxisymmetric()), options(config->GetSAParsedOptions()), transition_LM(config->GetKind_Trans_Model() == TURB_TRANS_MODEL::LM) { /*--- Setup the Jacobian pointer, we need to return su2double** but we know @@ -217,6 +250,9 @@ class CSourceBase_TurbSA : public CNumerics { SourceTerms::get(ScalarVar_i[0], var, Production, Destruction, CrossProduction, Jacobian_i[0]); Residual = (Production - Destruction + CrossProduction) * Volume; + + if (axisymmetric) ResidualAxisymmetricDiffusion(var.sigma); + Jacobian_i[0] *= Volume; } @@ -520,6 +556,19 @@ class CCompressibilityCorrection final : public ParentClass { const su2double d_CompCorrection = 2.0 * c5 * ScalarVar_i[0] / pow(sound_speed, 2) * aux_cc * Volume; const su2double CompCorrection = 0.5 * ScalarVar_i[0] * d_CompCorrection; + /* Axisymmetric contribution */ + if (config->GetAxisymmetric() && this->Coord_i[1] > EPS) { + const su2double yinv = 1.0 / this->Coord_i[1]; + const su2double nue = ScalarVar_i[0]; + const su2double v = V_i[idx.Velocity() + 1]; + + const su2double d_axiCorrection = 2.0 * c5 * nue * pow(v * yinv / sound_speed, 2) * Volume; + const su2double axiCorrection = 0.5 * nue * d_axiCorrection; + + this->Residual -= axiCorrection; + this->Jacobian_i[0] -= d_axiCorrection; + } + this->Residual -= CompCorrection; this->Jacobian_i[0] -= d_CompCorrection; diff --git a/SU2_CFD/src/solvers/CTurbSASolver.cpp b/SU2_CFD/src/solvers/CTurbSASolver.cpp index 2d301dc7722..143a09323d6 100644 --- a/SU2_CFD/src/solvers/CTurbSASolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSASolver.cpp @@ -308,6 +308,8 @@ void CTurbSASolver::Viscous_Residual(unsigned long iEdge, CGeometry* geometry, C void CTurbSASolver::Source_Residual(CGeometry *geometry, CSolver **solver_container, CNumerics **numerics_container, CConfig *config, unsigned short iMesh) { + bool axisymmetric = config->GetAxisymmetric(); + const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); const bool harmonic_balance = (config->GetTime_Marching() == TIME_MARCHING::HARMONIC_BALANCE); const bool transition_BC = config->GetSAParsedOptions().bc; @@ -383,6 +385,11 @@ void CTurbSASolver::Source_Residual(CGeometry *geometry, CSolver **solver_contai numerics->SetIntermittency(solver_container[TRANS_SOL]->GetNodes()->GetSolution(iPoint, 0)); } + if (axisymmetric) { + /*--- Set y coordinate ---*/ + numerics->SetCoord(geometry->nodes->GetCoord(iPoint), geometry->nodes->GetCoord(iPoint)); + } + /*--- Compute the source term ---*/ auto residual = numerics->ComputeResidual(config); From 32144a20e7274015690b82bf180d65d736459843 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Corentin=20Thom=C3=A9e?= Date: Fri, 12 Jan 2024 13:52:37 +0100 Subject: [PATCH 2/3] Fix reorder warning --- SU2_CFD/include/numerics/turbulent/turb_sources.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp index 2382f84388d..23e51975504 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp @@ -121,8 +121,8 @@ class CSourceBase_TurbSA : public CNumerics { CSourceBase_TurbSA(unsigned short nDim, const CConfig* config) : CNumerics(nDim, 1, config), idx(nDim, config->GetnSpecies()), - axisymmetric(config->GetAxisymmetric()), options(config->GetSAParsedOptions()), + axisymmetric(config->GetAxisymmetric()), transition_LM(config->GetKind_Trans_Model() == TURB_TRANS_MODEL::LM) { /*--- 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. ---*/ @@ -556,7 +556,7 @@ class CCompressibilityCorrection final : public ParentClass { const su2double d_CompCorrection = 2.0 * c5 * ScalarVar_i[0] / pow(sound_speed, 2) * aux_cc * Volume; const su2double CompCorrection = 0.5 * ScalarVar_i[0] * d_CompCorrection; - /* Axisymmetric contribution */ + /*--- Axisymmetric contribution ---*/ if (config->GetAxisymmetric() && this->Coord_i[1] > EPS) { const su2double yinv = 1.0 / this->Coord_i[1]; const su2double nue = ScalarVar_i[0]; From d9f6fa2eb879df1a82f94de522e7c33ce21e50bc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Corentin=20Thom=C3=A9e?= Date: Fri, 12 Jan 2024 14:33:30 +0100 Subject: [PATCH 3/3] Remove useless call to getAxisymmetric() --- SU2_CFD/include/numerics/turbulent/turb_sources.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp index 23e51975504..956979af846 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp @@ -557,7 +557,7 @@ class CCompressibilityCorrection final : public ParentClass { const su2double CompCorrection = 0.5 * ScalarVar_i[0] * d_CompCorrection; /*--- Axisymmetric contribution ---*/ - if (config->GetAxisymmetric() && this->Coord_i[1] > EPS) { + if (this->axisymmetric && this->Coord_i[1] > EPS) { const su2double yinv = 1.0 / this->Coord_i[1]; const su2double nue = ScalarVar_i[0]; const su2double v = V_i[idx.Velocity() + 1];