Enable the use of conductivity models different from CONSTANT_PRANDTL for compressible flows. (#2420)

* enable use of conductivity models for compressible flows

* removing unused variables

* setting eddy viscosity in fluid model

* axisymmetric cleaning

* updating residuals

* fix

* fix parallel regression

* fix

* adding test case using CONSTANT_CONDUCTIVITY model

* rename

* update regression test

* additional cleaning

* further cleaning

* making thermal conductivity only laminar

* small fix

* fix heatflux

* fix

* setting reference values

* fix ND heatSolver

* freestream values

* including Cp setHeatFluxVector

* reformulating

* moving up SetHeatFluxVector

* updating residual test case

* fixing missing thermalConductivity

* increasing nPrimVar and removing some variables

* making LambdaVisc independent of prandtl laminar

* addressing review comments

* updating residuals regression cases

* fixing residuals

---------

Co-authored-by: Cristopher-Morales <[email protected]>

Author: Cristopher-Morales

Common/include/CConfig.hpp

@@ -899,6 +899,8 @@ class CConfig {
   ModVel_FreeStreamND,             /*!< \brief Non-dimensional magnitude of the free-stream velocity of the fluid.  */
   Density_FreeStream,              /*!< \brief Free-stream density of the fluid. */
   Viscosity_FreeStream,            /*!< \brief Free-stream viscosity of the fluid.  */
+  ThermalConductivity_FreeStream,  /*!< \brief Free-stream thermal conductivity of the fluid. */
+  SpecificHeatCp_FreeStream,       /*!< \brief Free-stream specific heat capacity at constant pressure of the fluid.  */
   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.  */
@@ -940,6 +942,8 @@ class CConfig {
   Velocity_FreeStreamND[3],   /*!< \brief Farfield velocity values (external flow). */
   Energy_FreeStreamND,        /*!< \brief Farfield energy value (external flow). */
   Viscosity_FreeStreamND,     /*!< \brief Farfield viscosity value (external flow). */
+  ThermalConductivity_FreeStreamND,  /*!< \brief Farfield thermal conductivity value (external flow). */
+  SpecificHeatCp_FreeStreamND,      /*!< \brief Farfield specific heat capacity at constant pressure value (external flow). */
   Tke_FreeStreamND,           /*!< \brief Farfield kinetic energy (external flow). */
   Omega_FreeStreamND,         /*!< \brief Specific dissipation (external flow). */
   Omega_FreeStream;           /*!< \brief Specific dissipation (external flow). */
@@ -1760,6 +1764,18 @@ class CConfig {
    */
   su2double GetViscosity_FreeStream(void) const { return Viscosity_FreeStream; }
 
+  /*!
+   * \brief Get the value of the freestream thermal conductivity.
+   * \return Freestream thermal conductivity.
+   */
+  su2double GetThermalConductivity_FreeStream(void) const { return ThermalConductivity_FreeStream; }
+
+  /*!
+   * \brief Get the value of the freestream heat capacity at constant pressure.
+   * \return Freestream heat capacity at constant pressure.
+   */
+  su2double GetSpecificHeatCp_FreeStream(void) const { return SpecificHeatCp_FreeStream; }
+
   /*!
    * \brief Get the value of the freestream density.
    * \return Freestream density.
@@ -1991,6 +2007,18 @@ class CConfig {
    */
   su2double GetViscosity_FreeStreamND(void) const { return Viscosity_FreeStreamND; }
 
+  /*!
+   * \brief Get the value of the non-dimensionalized freestream thermal conductivity.
+   * \return Non-dimensionalized freestream thermal conductivity.
+   */
+  su2double GetThermalConductivity_FreeStreamND(void) const { return ThermalConductivity_FreeStreamND; }
+
+  /*!
+   * \brief Get the value of the non-dimensionalized freestream heat capacity at constant pressure.
+   * \return Non-dimensionalized freestream heat capacity at constant pressure.
+   */
+  su2double GetSpecificHeatCp_FreeStreamND(void) const { return SpecificHeatCp_FreeStreamND; }
+
   /*!
    * \brief Get the value of the non-dimensionalized freestream viscosity.
    * \return Non-dimensionalized freestream viscosity.
@@ -2610,6 +2638,18 @@ class CConfig {
    */
   void SetViscosity_FreeStream(su2double val_viscosity_freestream) { Viscosity_FreeStream = val_viscosity_freestream; }
 
+  /*!
+   * \brief Set the freestream thermal conductivity.
+   * \param[in] val_thermalconductivity_freestream - Value of the freestream thermal conductivity.
+   */
+  void SetThermalConductivity_FreeStream(su2double val_thermalconductivity_freestream) { ThermalConductivity_FreeStream = val_thermalconductivity_freestream; }
+
+  /*!
+   * \brief Set the freestream specific heat capacity at constant pressure.
+   * \param[in] val_specificheatCp_freestream - Value of the freestream specific heat capacity at constant pressure.
+   */
+  void SetSpecificHeatCp_FreeStream(su2double val_specificheatCp_freestream) { SpecificHeatCp_FreeStream = val_specificheatCp_freestream; }
+
   /*!
    * \brief Set the magnitude of the free-stream velocity.
    * \param[in] val_modvel_freestream - Magnitude of the free-stream velocity.
@@ -2672,6 +2712,18 @@ class CConfig {
    */
   void SetViscosity_FreeStreamND(su2double val_viscosity_freestreamnd) { Viscosity_FreeStreamND = val_viscosity_freestreamnd; }
 
+  /*!
+   * \brief Set the non-dimensional free-stream thermal conductivity.
+   * \param[in] val_thermalconductivity_freestreamnd - Value of the non-dimensional free-stream thermal conductivity.
+   */
+  void SetThermalConductivity_FreeStreamND(su2double val_thermalconductivity_freestreamnd) { ThermalConductivity_FreeStreamND = val_thermalconductivity_freestreamnd; }
+
+  /*!
+   * \brief Set the non-dimensional free-stream specific heat capacity at constant pressure.
+   * \param[in] val_specificheatCp_freestreamnd - Value of the non-dimensional free-stream specific heat capacity at constant pressure.
+   */
+  void SetSpecificHeatCp_FreeStreamND(su2double val_specificheatCp_freestreamnd) { SpecificHeatCp_FreeStreamND = val_specificheatCp_freestreamnd; }
+
   /*!
    * \brief Set the non-dimensional freestream turbulent kinetic energy.
    * \param[in] val_tke_freestreamnd - Value of the non-dimensional freestream turbulent kinetic energy.

Common/src/CConfig.cpp

@@ -4214,10 +4214,6 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
           SU2_MPI::Error("Only SUTHERLAND viscosity model can be used with US Measurement", CURRENT_FUNCTION);
         }
       }
-      if (Kind_ConductivityModel != CONDUCTIVITYMODEL::CONSTANT_PRANDTL) {
-        SU2_MPI::Error("Only CONSTANT_PRANDTL thermal conductivity model can be used with STANDARD_AIR and IDEAL_GAS",
-                       CURRENT_FUNCTION);
-      }
     }
     /*--- Check for Boundary condition option agreement ---*/
   if (Kind_InitOption == REYNOLDS){

Common/src/geometry/CPhysicalGeometryFEM.cpp

@@ -2228,11 +2228,12 @@ void CPhysicalGeometry::DetermineTimeLevelElements(CConfig* config, const vector
           step is based on free stream values at the moment, but this is easy
           to change, if needed. ---*/
     const su2double Gamma = config->GetGamma();
-    const su2double Prandtl = config->GetPrandtl_Lam();
 
     const su2double Density = config->GetDensity_FreeStreamND();
     const su2double* Vel = config->GetVelocity_FreeStreamND();
     const su2double Viscosity = config->GetViscosity_FreeStreamND();
+    const su2double Thermal_Conductivity = config->GetThermalConductivity_FreeStreamND();
+    const su2double Cp = config->GetSpecificHeatCp_FreeStreamND();
 
     su2double VelMag = 0.0;
     for (unsigned short iDim = 0; iDim < nDim; ++iDim) VelMag += Vel[iDim] * Vel[iDim];
@@ -2253,7 +2254,7 @@ void CPhysicalGeometry::DetermineTimeLevelElements(CConfig* config, const vector
     const su2double charVel = sqrt(charVel2);
 
     /* Also the viscous contribution to the time step is constant. Compute it. */
-    const su2double factHeatFlux = Gamma / Prandtl;
+    const su2double factHeatFlux = (Thermal_Conductivity * Gamma) / (Cp * Viscosity);
     const su2double lambdaOverMu = -TWO3;
     const su2double radVisc = max(max(1.0, 2.0 + lambdaOverMu), factHeatFlux) * Viscosity / Density;
 

SU2_CFD/include/numerics/CNumerics.hpp

@@ -53,7 +53,6 @@ class CNumerics {
   su2double Gamma_Minus_One;  /*!< \brief Fluids's Gamma - 1.0  . */
   su2double Minf;             /*!< \brief Free stream Mach number . */
   su2double Gas_Constant;     /*!< \brief Gas constant. */
-  su2double Prandtl_Lam;      /*!< \brief Laminar Prandtl's number. */
   su2double Prandtl_Turb;     /*!< \brief Turbulent Prandtl's number. */
   su2double MassFlux;         /*!< \brief Mass flux across edge. */
   su2double

SU2_CFD/include/numerics/flow/flow_diffusion.hpp

@@ -57,6 +57,8 @@ class CAvgGrad_Base : public CNumerics {
   su2double **Mean_GradPrimVar = nullptr, /*!< \brief Mean value of the gradient. */
   Mean_Laminar_Viscosity,                 /*!< \brief Mean value of the viscosity. */
   Mean_Eddy_Viscosity,                    /*!< \brief Mean value of the eddy viscosity. */
+  Mean_Thermal_Conductivity,              /*!< \brief Mean value of the thermal conductivity. */
+  Mean_Cp,                                /*!< \brief Mean value of the specific heat capacity at constant pressure. */
   Mean_turb_ke,                           /*!< \brief Mean value of the turbulent kinetic energy. */
   Mean_TauWall,                           /*!< \brief Mean wall shear stress (wall functions). */
   TauWall_i, TauWall_j,                   /*!< \brief Wall shear stress at point i and j (wall functions). */
@@ -214,6 +216,16 @@ class CAvgGrad_Base : public CNumerics {
    */
   inline su2double GetStressTensor(unsigned short iDim, unsigned short jDim) const { return tau[iDim][jDim];}
 
+  /*!
+   * \brief Compute the heat flux due to molecular and turbulent diffusivity
+   * \param[in] val_gradprimvar - Gradient of the primitive variables.
+   * \param[in] val_eddy_viscosity - Eddy viscosity.
+   * \param[in] val_thermal_conductivity - Thermal Conductivity.
+   * \param[in] val_heat_capacity_cp - Heat Capacity at constant pressure.
+   */
+  void SetHeatFluxVector(const su2double* const* val_gradprimvar, su2double val_eddy_viscosity,
+                         su2double val_thermal_conductivity, su2double val_heat_capacity_cp);
+
   /*!
    * \brief Get a component of the heat flux vector.
    * \param[in] iDim - The index of the component
@@ -248,16 +260,6 @@ class CAvgGrad_Flow final : public CAvgGrad_Base {
    */
   ResidualType<> ComputeResidual(const CConfig* config) override;
 
-  /*!
-   * \brief Compute the heat flux due to molecular and turbulent diffusivity
-   * \param[in] val_gradprimvar - Gradient of the primitive variables.
-   * \param[in] val_laminar_viscosity - Laminar viscosity.
-   * \param[in] val_eddy_viscosity - Eddy viscosity.
-   */
-  void SetHeatFluxVector(const su2double* const *val_gradprimvar,
-                         su2double val_laminar_viscosity,
-                         su2double val_eddy_viscosity);
-
   /*!
    * \brief Compute the Jacobian of the heat flux vector
    *
@@ -266,13 +268,14 @@ class CAvgGrad_Flow final : public CAvgGrad_Base {
    *
    * \param[in] val_Mean_PrimVar - Mean value of the primitive variables.
    * \param[in] val_gradprimvar - Mean value of the gradient of the primitive variables.
-   * \param[in] val_laminar_viscosity - Value of the laminar viscosity.
+   * \param[in] val_heat_capacity_cp - Value of the heat capacity at constant pressure.
    * \param[in] val_eddy_viscosity - Value of the eddy viscosity.
    * \param[in] val_dist_ij - Distance between the points.
    * \param[in] val_normal - Normal vector, the norm of the vector is the area of the face.
    */
   void SetHeatFluxJacobian(const su2double *val_Mean_PrimVar,
-                           su2double val_laminar_viscosity,
+                           su2double val_heat_capacity_cp,
+                           su2double val_thermal_conductivity,
                            su2double val_eddy_viscosity,
                            su2double val_dist_ij,
                            const su2double *val_normal);
@@ -351,20 +354,6 @@ class CGeneralAvgGrad_Flow final : public CAvgGrad_Base {
   Mean_Thermal_Conductivity, /*!< \brief Mean value of the thermal conductivity. */
   Mean_Cp;                   /*!< \brief Mean value of the Cp. */
 
-  /*!
-   * \brief Compute the heat flux due to molecular and turbulent diffusivity
-   * \param[in] val_gradprimvar - Gradient of the primitive variables.
-   * \param[in] val_laminar_viscosity - Laminar viscosity.
-   * \param[in] val_eddy_viscosity - Eddy viscosity.
-   * \param[in] val_thermal_conductivity - Thermal Conductivity.
-   * \param[in] val_heat_capacity_cp - Heat Capacity at constant pressure.
-   */
-  void SetHeatFluxVector(const su2double* const *val_gradprimvar,
-                         su2double val_laminar_viscosity,
-                         su2double val_eddy_viscosity,
-                         su2double val_thermal_conductivity,
-                         su2double val_heat_capacity_cp);
-
   /*!
    * \brief Compute the Jacobian of the heat flux vector
    *

SU2_CFD/include/numerics_simd/flow/diffusion/viscous_fluxes.hpp

@@ -75,9 +75,7 @@ class CCompressibleViscousFluxBase : public CNumericsSIMD {
 
   const su2double gamma;
   const su2double gasConst;
-  const su2double prandtlLam;
   const su2double prandtlTurb;
-  const su2double cp;
   const bool correct;
   const bool useSA_QCR;
   const bool wallFun;
@@ -97,9 +95,7 @@ class CCompressibleViscousFluxBase : public CNumericsSIMD {
                                const CVariable* turbVars_, Ts&...) :
     gamma(config.GetGamma()),
     gasConst(config.GetGas_ConstantND()),
-    prandtlLam(config.GetPrandtl_Lam()),
     prandtlTurb(config.GetPrandtl_Turb()),
-    cp(gamma * gasConst / (gamma - 1)),
     correct(iMesh == MESH_0),
     useSA_QCR(config.GetSAParsedOptions().qcr2000),
     wallFun(config.GetWall_Functions()),
@@ -256,13 +252,11 @@ class CCompressibleViscousFluxBase : public CNumericsSIMD {
 template<size_t NDIM>
 class CCompressibleViscousFlux : public CCompressibleViscousFluxBase<NDIM, CCompressibleViscousFlux<NDIM> > {
 public:
-  static constexpr size_t nPrimVar = NDIM+7;
+  static constexpr size_t nPrimVar = NDIM+9;
   using Base = CCompressibleViscousFluxBase<NDIM, CCompressibleViscousFlux<NDIM> >;
   using Base::gamma;
   using Base::gasConst;
-  using Base::prandtlLam;
   using Base::prandtlTurb;
-  using Base::cp;
 
   /*!
    * \brief Constructor, initialize constants and booleans.
@@ -275,7 +269,7 @@ class CCompressibleViscousFlux : public CCompressibleViscousFluxBase<NDIM, CComp
    */
   template<class PrimitiveType>
   FORCEINLINE Double thermalConductivity(const PrimitiveType& V) const {
-    return cp * (V.laminarVisc()/prandtlLam + V.eddyVisc()/prandtlTurb);
+    return V.thermalCond() + V.cp() * V.eddyVisc() / prandtlTurb;
   }
 
   /*!

SU2_CFD/include/solvers/CEulerSolver.hpp

@@ -41,7 +41,6 @@ class CEulerSolver : public CFVMFlowSolverBase<CEulerVariable, ENUM_REGIME::COMP
   using BaseClass = CFVMFlowSolverBase<CEulerVariable, ENUM_REGIME::COMPRESSIBLE>;
 
   su2double
-  Prandtl_Lam = 0.0,    /*!< \brief Laminar Prandtl number. */
   Prandtl_Turb = 0.0;   /*!< \brief Turbulent Prandtl number. */
 
   su2double AllBound_CEquivArea_Inv=0.0; /*!< \brief equivalent area coefficient (inviscid contribution) for all the boundaries. */

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -2389,7 +2389,7 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
   unsigned long iVertex, iPoint, iPointNormal;
   unsigned short iMarker, iMarker_Monitoring, iDim, jDim;
   su2double Viscosity = 0.0, Area, Density = 0.0, FrictionVel,
-            UnitNormal[3] = {0.0}, TauElem[3] = {0.0}, Tau[3][3] = {{0.0}}, Cp,
+            UnitNormal[3] = {0.0}, TauElem[3] = {0.0}, Tau[3][3] = {{0.0}},
             thermal_conductivity, MaxNorm = 8.0, Grad_Vel[3][3] = {{0.0}}, Grad_Temp[3] = {0.0},
             Grad_Temp_ve[3] = {0.0}, AxiFactor;
   const su2double *Coord = nullptr, *Coord_Normal = nullptr, *Normal = nullptr;
@@ -2400,10 +2400,8 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
   const su2double RefLength = config->GetRefLength();
   const su2double RefHeatFlux = config->GetHeat_Flux_Ref();
   const su2double RefTemperature = config->GetTemperature_Ref();
-  const su2double Gas_Constant = config->GetGas_ConstantND();
   auto Origin = config->GetRefOriginMoment(0);
 
-  const su2double Prandtl_Lam = config->GetPrandtl_Lam();
   const bool energy = config->GetEnergy_Equation();
   const bool QCR = config->GetSAParsedOptions().qcr2000;
   const bool axisymmetric = config->GetAxisymmetric();
@@ -2544,11 +2542,7 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
       /*--- Compute total and maximum heat flux on the wall ---*/
 
       if (!nemo) {
-        if (FlowRegime == ENUM_REGIME::COMPRESSIBLE) {
-          Cp = (Gamma / Gamma_Minus_One) * Gas_Constant;
-          thermal_conductivity = Cp * Viscosity / Prandtl_Lam;
-        }
-        if (FlowRegime == ENUM_REGIME::INCOMPRESSIBLE) {
+        if ((FlowRegime == ENUM_REGIME::COMPRESSIBLE) || (FlowRegime == ENUM_REGIME::INCOMPRESSIBLE)) {
           thermal_conductivity = nodes->GetThermalConductivity(iPoint);
         }
 

SU2_CFD/include/solvers/CHeatSolver.hpp

@@ -67,7 +67,6 @@ class CHeatSolver final : public CScalarSolver<CHeatVariable> {
     if (!geometry->nodes->GetDomain(iPoint)) return;
 
     const bool implicit = config->GetKind_TimeIntScheme() == EULER_IMPLICIT;
-    const su2double prandtl_lam = config->GetPrandtl_Lam();
     const su2double const_diffusivity = config->GetThermalDiffusivity();
 
     const auto Point_Normal = geometry->vertex[iMarker][iVertex]->GetNormal_Neighbor();
@@ -83,7 +82,8 @@ class CHeatSolver final : public CScalarSolver<CHeatVariable> {
 
     su2double thermal_diffusivity = const_diffusivity;
     if (flow) {
-      thermal_diffusivity = flow_solver->GetNodes()->GetLaminarViscosity(iPoint) / prandtl_lam;
+      thermal_diffusivity =
+          flow_solver->GetNodes()->GetThermalConductivity(iPoint) / flow_solver->GetNodes()->GetSpecificHeatCp(iPoint);
     }
     LinSysRes(iPoint, 0) -= thermal_diffusivity * dTdn * Area;
 
@@ -107,21 +107,19 @@ class CHeatSolver final : public CScalarSolver<CHeatVariable> {
     const CVariable* flow_nodes = flow ? solver_container[FLOW_SOL]->GetNodes() : nullptr;
 
     const su2double const_diffusivity = config->GetThermalDiffusivity();
-    const su2double pr_lam = config->GetPrandtl_Lam();
     const su2double pr_turb = config->GetPrandtl_Turb();
 
     su2double thermal_diffusivity_i{}, thermal_diffusivity_j{};
 
     /*--- Computes the thermal diffusivity to use in the viscous numerics. ---*/
     auto compute_thermal_diffusivity = [&](unsigned long iPoint, unsigned long jPoint) {
       if (flow) {
-        thermal_diffusivity_i = flow_nodes->GetLaminarViscosity(iPoint) / pr_lam +
+        thermal_diffusivity_i = flow_nodes->GetThermalConductivity(iPoint) / flow_nodes->GetSpecificHeatCp(iPoint) +
                                 flow_nodes->GetEddyViscosity(iPoint) / pr_turb;
-        thermal_diffusivity_j = flow_nodes->GetLaminarViscosity(jPoint) / pr_lam +
+        thermal_diffusivity_j = flow_nodes->GetThermalConductivity(jPoint) / flow_nodes->GetSpecificHeatCp(jPoint) +
                                 flow_nodes->GetEddyViscosity(jPoint) / pr_turb;
         numerics->SetDiffusionCoeff(&thermal_diffusivity_i, &thermal_diffusivity_j);
-      }
-      else {
+      } else {
         numerics->SetDiffusionCoeff(&const_diffusivity, &const_diffusivity);
       }
     };

SU2_CFD/src/interfaces/cht/CConjugateHeatInterface.cpp

@@ -75,14 +75,7 @@ void CConjugateHeatInterface::GetDonor_Variable(CSolver *donor_solution, CGeomet
 
   if (compressible_flow) {
 
-    const su2double Gamma = donor_config->GetGamma();
-    const su2double Gas_Constant = donor_config->GetGas_ConstantND();
-    const su2double Cp = (Gamma / (Gamma - 1.0)) * Gas_Constant;
-
-    const su2double Prandtl_Lam = donor_config->GetPrandtl_Lam();
-    const su2double laminar_viscosity = donor_solution->GetNodes()->GetLaminarViscosity(Point_Donor); // TDE check for consistency
-
-    const su2double thermal_conductivityND = Cp*(laminar_viscosity/Prandtl_Lam);
+    const su2double thermal_conductivityND = donor_solution->GetNodes()->GetThermalConductivity(Point_Donor);
     heat_flux_density = thermal_conductivityND*dTdn;
 
     if ((donor_config->GetKind_CHT_Coupling() == CHT_COUPLING::DIRECT_TEMPERATURE_ROBIN_HEATFLUX) ||