diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 1acfe7a7a797..ad1a80b9573e 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -102,6 +102,7 @@ const su2double EPS = 1.0E-16; /*!< \brief Error scale. */ const su2double TURB_EPS = 1.0E-16; /*!< \brief Turbulent Error scale. */ const su2double ONE2 = 0.5; /*!< \brief One divided by two. */ +const su2double ONE3 = 1.0 / 3.0; /*!< \brief One divided by three. */ const su2double TWO3 = 2.0 / 3.0; /*!< \brief Two divided by three. */ const su2double FOUR3 = 4.0 / 3.0; /*!< \brief Four divided by three. */ diff --git a/SU2_CFD/include/numerics/flow/flow_sources.hpp b/SU2_CFD/include/numerics/flow/flow_sources.hpp index 6834d8a2a380..ee889780c694 100644 --- a/SU2_CFD/include/numerics/flow/flow_sources.hpp +++ b/SU2_CFD/include/numerics/flow/flow_sources.hpp @@ -63,7 +63,16 @@ class CSourceBase_Flow : public CNumerics { * \ingroup SourceDiscr * \author F. Palacios */ -class CSourceAxisymmetric_Flow final : public CSourceBase_Flow { +class CSourceAxisymmetric_Flow : public CSourceBase_Flow { +protected: + bool implicit, viscous; + su2double yinv{0.0}; + + /*! + * \brief Diffusion residual of the axisymmetric source term. + */ + void ResidualDiffusion(); + public: /*! * \brief Constructor of the class. @@ -74,12 +83,31 @@ class CSourceAxisymmetric_Flow final : public CSourceBase_Flow { CSourceAxisymmetric_Flow(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config); /*! - * \brief Residual of the rotational frame source term. + * \brief Residual of the axisymmetric source term. * \param[in] config - Definition of the particular problem. * \return Lightweight const-view of residual and Jacobian. */ ResidualType<> ComputeResidual(const CConfig* config) override; + +}; +/*! + * \class CSourceGeneralAxisymmetric_Flow + * \brief Class for source term for solving axisymmetric problems for a general (non ideal) fluid. + * \ingroup SourceDiscr + * \author F. Dittmann + */ +class CSourceGeneralAxisymmetric_Flow final : public CSourceAxisymmetric_Flow { +public: + + using CSourceAxisymmetric_Flow::CSourceAxisymmetric_Flow; + /*! + * \brief Residual of the general axisymmetric source term. + * \param[in] config - Definition of the particular problem. + * \return Lightweight const-view of residual and Jacobian. + */ + ResidualType<> ComputeResidual(const CConfig* config) override; + }; /*! diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index 6e8530faad5e..16fbd073feb2 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -1865,8 +1865,10 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol else if (config->GetAxisymmetric() == YES) { if (incompressible) numerics[iMGlevel][FLOW_SOL][source_first_term] = new CSourceIncAxisymmetric_Flow(nDim, nVar_Flow, config); - else + else if (ideal_gas) numerics[iMGlevel][FLOW_SOL][source_first_term] = new CSourceAxisymmetric_Flow(nDim, nVar_Flow, config); + else + numerics[iMGlevel][FLOW_SOL][source_first_term] = new CSourceGeneralAxisymmetric_Flow(nDim, nVar_Flow, config); } else if (config->GetGravityForce() == YES) { numerics[iMGlevel][FLOW_SOL][source_first_term] = new CSourceGravity(nDim, nVar_Flow, config); diff --git a/SU2_CFD/src/numerics/flow/flow_sources.cpp b/SU2_CFD/src/numerics/flow/flow_sources.cpp index c871c5a841a4..00055971ff45 100644 --- a/SU2_CFD/src/numerics/flow/flow_sources.cpp +++ b/SU2_CFD/src/numerics/flow/flow_sources.cpp @@ -50,16 +50,17 @@ CSourceAxisymmetric_Flow::CSourceAxisymmetric_Flow(unsigned short val_nDim, unsi Gamma = config->GetGamma(); Gamma_Minus_One = Gamma - 1.0; + + implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); + viscous = config->GetViscous(); } CNumerics::ResidualType<> CSourceAxisymmetric_Flow::ComputeResidual(const CConfig* config) { - su2double yinv, Pressure_i, Enthalpy_i, Velocity_i, sq_vel; + su2double Pressure_i, Enthalpy_i, Velocity_i, sq_vel; unsigned short iDim, iVar, jVar; - - bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); - + if (Coord_i[1] > EPS) { yinv = 1.0/Coord_i[1]; @@ -77,6 +78,8 @@ CNumerics::ResidualType<> CSourceAxisymmetric_Flow::ComputeResidual(const CConfi residual[1] = yinv*Volume*U_i[1]*U_i[2]/U_i[0]; residual[2] = yinv*Volume*(U_i[2]*U_i[2]/U_i[0]); residual[3] = yinv*Volume*Enthalpy_i*U_i[2]; + + /*--- Inviscid component of the source term. ---*/ if (implicit) { jacobian[0][0] = 0.0; @@ -104,6 +107,10 @@ CNumerics::ResidualType<> CSourceAxisymmetric_Flow::ComputeResidual(const CConfi jacobian[iVar][jVar] *= yinv*Volume; } + + /*--- Add the viscous terms if necessary. ---*/ + + if (viscous) ResidualDiffusion(); } @@ -124,6 +131,107 @@ CNumerics::ResidualType<> CSourceAxisymmetric_Flow::ComputeResidual(const CConfi return ResidualType<>(residual, jacobian, nullptr); } +void CSourceAxisymmetric_Flow::ResidualDiffusion(){ + + su2double laminar_viscosity_i = V_i[nDim+5]; + su2double eddy_viscosity_i = V_i[nDim+6]; + su2double thermal_conductivity_i = V_i[nDim+7]; + su2double heat_capacity_cp_i = V_i[nDim+8]; + + su2double total_viscosity_i = laminar_viscosity_i + eddy_viscosity_i; + su2double total_conductivity_i = thermal_conductivity_i + heat_capacity_cp_i*eddy_viscosity_i/Prandtl_Turb; + + su2double u = U_i[1]/U_i[0]; + su2double v = U_i[2]/U_i[0]; + + residual[0] -= 0.0; + residual[1] -= yinv*Volume*total_viscosity_i*(PrimVar_Grad_i[1][1]+ONE3*PrimVar_Grad_i[2][0]); // - 2/3 * v * d(mu)/dx + residual[2] -= yinv*Volume*total_viscosity_i*FOUR3*(PrimVar_Grad_i[2][1]-v*yinv); // - 2/3 * v * d(mu)/dy + residual[3] -= yinv*Volume*(total_viscosity_i*(u*(PrimVar_Grad_i[2][0]+PrimVar_Grad_i[1][1]-TWO3*(PrimVar_Grad_i[2][1]-yinv*v)) + - TWO3*v*(PrimVar_Grad_i[1][1]+PrimVar_Grad_i[1][0])) + + total_conductivity_i*PrimVar_Grad_i[0][1]); // - 2/3 * (uu + uv) * d(mu)/dy +} + + +CNumerics::ResidualType<> CSourceGeneralAxisymmetric_Flow::ComputeResidual(const CConfig* config) { + unsigned short iVar, jVar; + + if (Coord_i[1] > EPS) { + + yinv = 1.0/Coord_i[1]; + + su2double Density_i = U_i[0]; + su2double Velocity1_i = U_i[1]/U_i[0]; + su2double Velocity2_i = U_i[2]/U_i[0]; + su2double Energy_i = U_i[3]/U_i[0]; + + su2double Pressure_i = V_j[3]; + su2double Enthalpy_i = Energy_i + Pressure_i/Density_i; + + /*--- Inviscid component of the source term. ---*/ + + residual[0] = yinv*Volume*U_i[2]; + residual[1] = yinv*Volume*U_i[1]*Velocity2_i; + residual[2] = yinv*Volume*U_i[2]*Velocity2_i; + residual[3] = yinv*Volume*U_i[2]*Enthalpy_i; + + if (implicit) { + + su2double dPdrho_e_i = S_i[0]; + su2double dPde_rho_i = S_i[1]; + + jacobian[0][0] = 0.0; + jacobian[0][1] = 0.0; + jacobian[0][2] = 1.0; + jacobian[0][3] = 0.0; + + jacobian[1][0] = -Velocity1_i*Velocity2_i; + jacobian[1][1] = Velocity2_i; + jacobian[1][2] = Velocity1_i; + jacobian[1][3] = 0.0; + + jacobian[2][0] = -Velocity2_i*Velocity2_i; + jacobian[2][1] = 0.0; + jacobian[2][2] = 2*Velocity2_i; + jacobian[2][3] = 0.0; + + jacobian[3][0] = Velocity2_i*(dPdrho_e_i + dPde_rho_i/Density_i*(Velocity1_i*Velocity1_i + + Velocity2_i*Velocity2_i + - Energy_i) - Enthalpy_i); + jacobian[3][1] = -Velocity1_i*Velocity2_i/Density_i *dPde_rho_i; + jacobian[3][2] = Enthalpy_i - Velocity2_i*Velocity2_i/Density_i *dPde_rho_i; + jacobian[3][3] = Velocity2_i + Velocity2_i/Density_i *dPde_rho_i; + + for (iVar=0; iVar < nVar; iVar++) + for (jVar=0; jVar < nVar; jVar++) + jacobian[iVar][jVar] *= yinv*Volume; + + } + + /*--- Add the viscous terms if necessary. ---*/ + + if (viscous) ResidualDiffusion(); + + } + + else { + + for (iVar=0; iVar < nVar; iVar++) + residual[iVar] = 0.0; + + if (implicit) { + for (iVar=0; iVar < nVar; iVar++) { + for (jVar=0; jVar < nVar; jVar++) + jacobian[iVar][jVar] = 0.0; + } + } + + } + + return ResidualType<>(residual, jacobian, nullptr); +} + + CSourceIncAxisymmetric_Flow::CSourceIncAxisymmetric_Flow(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config) : CSourceBase_Flow(val_nDim, val_nVar, config) { diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 44cfdd1c3b6b..d750249ec66b 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -3046,6 +3046,8 @@ void CEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_contain const bool harmonic_balance = (config->GetTime_Marching() == HARMONIC_BALANCE); const bool windgust = config->GetWind_Gust(); const bool body_force = config->GetBody_Force(); + const bool ideal_gas = (config->GetKind_FluidModel() == STANDARD_AIR) || + (config->GetKind_FluidModel() == IDEAL_GAS); /*--- Pick one numerics object per thread. ---*/ CNumerics* numerics = numerics_container[SOURCE_FIRST_TERM + omp_get_thread_num()*MAX_TERMS]; @@ -3120,6 +3122,17 @@ void CEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_contain /*--- Set y coordinate ---*/ numerics->SetCoord(geometry->nodes->GetCoord(iPoint), geometry->nodes->GetCoord(iPoint)); + /*--- Set primitive variables ---*/ + numerics->SetPrimitive(nodes->GetPrimitive(iPoint), nodes->GetPrimitive(iPoint)); + + /*--- Set gradient of primitive variables ---*/ + numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), nodes->GetGradient_Primitive(iPoint)); + + if (!ideal_gas) { + /*--- Set secondary variables ---*/ + numerics->SetSecondary(nodes->GetSecondary(iPoint), nodes->GetSecondary(iPoint)); + } + /*--- Compute Source term Residual ---*/ auto residual = numerics->ComputeResidual(config);