Axisymmetric source viscous terms

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. */
 

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;
+
 };
 
 /*!

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);

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) {
 

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);