Implement NASA polynomial gas model (CNasaGas) and integrate into solvers

Author: Omgiri01

Common/include/CConfig.hpp

@@ -912,6 +912,9 @@ class CConfig {
   array<su2double, N_POLY_COEFFS> CpPolyCoefficientsND{{0.0}};  /*!< \brief Definition of the non-dimensional temperature polynomial coefficients for specific heat Cp. */
   array<su2double, N_POLY_COEFFS> MuPolyCoefficientsND{{0.0}};  /*!< \brief Definition of the non-dimensional temperature polynomial coefficients for viscosity. */
   array<su2double, N_POLY_COEFFS> KtPolyCoefficientsND{{0.0}};  /*!< \brief Definition of the non-dimensional temperature polynomial coefficients for thermal conductivity. */
+  array<su2double, N_NASA_POLY_COEFFS> NASA_CpLowPolyCoefficients{{0.0}};  /*!< \brief Definition of the NASA temperature polynomial coefficients for specific heat Cp (low temperature range). */
+  array<su2double, N_NASA_POLY_COEFFS> NASA_CpHighPolyCoefficients{{0.0}}; /*!< \brief Definition of the NASA temperature polynomial coefficients for specific heat Cp (high temperature range). */
+  su2double NASA_TransitionTemperature; /*!< \brief Transition temperature for the NASA polynomial gas model. */
   su2double TurbIntensityAndViscRatioFreeStream[2]; /*!< \brief Freestream turbulent intensity and viscosity ratio for turbulence and transition models. */
   su2double Energy_FreeStream,     /*!< \brief Free-stream total energy of the fluid.  */
   ModVel_FreeStream,               /*!< \brief Magnitude of the free-stream velocity of the fluid.  */
@@ -1125,6 +1128,9 @@ class CConfig {
   array<su2double, N_POLY_COEFFS> cp_polycoeffs{{0.0}};  /*!< \brief Array for specific heat polynomial coefficients. */
   array<su2double, N_POLY_COEFFS> mu_polycoeffs{{0.0}};  /*!< \brief Array for viscosity polynomial coefficients. */
   array<su2double, N_POLY_COEFFS> kt_polycoeffs{{0.0}};  /*!< \brief Array for thermal conductivity polynomial coefficients. */
+  array<su2double, N_NASA_POLY_COEFFS> nasa_cp_low_polycoeffs{{0.0}};  /*!< \brief Array for NASA high temperature polynomial coefficients. */
+  array<su2double, N_NASA_POLY_COEFFS> nasa_cp_high_polycoeffs{{0.0}}; /*!< \brief Array for NASA high temperature polynomial coefficients. */
+  su2double nasa_transition_temperature; /*!< \brief Transition temperature for the NASA polynomial gas model. */
   bool Body_Force;                      /*!< \brief Flag to know if a body force is included in the formulation. */
 
   struct CMUSCLRampParam {

Common/include/option_structure.hpp

@@ -202,6 +202,7 @@ constexpr passivedouble COLORING_EFF_THRESH = 0.875;  /*!< \brief Below this val
    Cp(T) = b0 + b1*T + b2*T^2 + b3*T^3 + b4*T^4. By default, all coeffs
    are set to zero and will be properly non-dim. in the solver. ---*/
 constexpr int N_POLY_COEFFS = 5; /*!< \brief Number of coefficients in temperature polynomial fits for fluid models. */
+constexpr int N_NASA_POLY_COEFFS = 7; /*!< \brief Number of coefficients in NASA temperature polynomial fits. */
 
 /*!
  * \brief Boolean answers
@@ -551,6 +552,7 @@ enum ENUM_FLUIDMODEL {
   COOLPROP = 10,          /*!< \brief Thermodynamics library. */
   FLUID_FLAMELET = 11,    /*!< \brief lookup table (LUT) method for premixed flamelets. */
   DATADRIVEN_FLUID = 12,  /*!< \brief multi-layer perceptron driven fluid model. */
+  NASA_GAS = 13,          /*!< \brief NASA polynomial ideal gas model. */
 };
 static const MapType<std::string, ENUM_FLUIDMODEL> FluidModel_Map = {
   MakePair("STANDARD_AIR", STANDARD_AIR)
@@ -566,6 +568,7 @@ static const MapType<std::string, ENUM_FLUIDMODEL> FluidModel_Map = {
   MakePair("COOLPROP", COOLPROP)
   MakePair("DATADRIVEN_FLUID", DATADRIVEN_FLUID)
   MakePair("FLUID_FLAMELET", FLUID_FLAMELET)
+  MakePair("NASA_GAS", NASA_GAS)
 };
 
 /*!

Common/src/CConfig.cpp

@@ -1334,6 +1334,12 @@ void CConfig::SetConfig_Options() {
   addDoubleArrayOption("MU_POLYCOEFFS", N_POLY_COEFFS, false, mu_polycoeffs.data());
   /* DESCRIPTION: Definition of the temperature polynomial coefficients for specific heat Cp. */
   addDoubleArrayOption("KT_POLYCOEFFS", N_POLY_COEFFS, false, kt_polycoeffs.data());
+  /* DESCRIPTION: Definition of the NASA temperature polynomial coefficients for specific heat Cp (low temperature range). */
+  addDoubleArrayOption("NASA_CP_LOW_COEFFS", N_NASA_POLY_COEFFS, false, nasa_cp_low_polycoeffs.data());
+  /* DESCRIPTION: Definition of the NASA temperature polynomial coefficients for specific heat Cp (high temperature range). */
+  addDoubleArrayOption("NASA_CP_HIGH_COEFFS", N_NASA_POLY_COEFFS, false, nasa_cp_high_polycoeffs.data());
+  /* DESCRIPTION: Transition temperature for the NASA polynomial gas model. */
+  addDoubleOption("NASA_T_TRANS", nasa_transition_temperature, 1000.0);
 
   /*!\brief REYNOLDS_NUMBER \n DESCRIPTION: Reynolds number (non-dimensional, based on the free-stream values). Needed for viscous solvers. For incompressible solvers the Reynolds length will always be 1.0 \n DEFAULT: 0.0 \ingroup Config */
   addDoubleOption("REYNOLDS_NUMBER", Reynolds, 0.0);

SU2_CFD/include/fluid/CNasaGas.hpp

@@ -0,0 +1,64 @@
+/*!
+ * \file CNasaGas.hpp
+ * \brief Defines the NASA polynomial gas model.
+ * \author Om Giri
+ * \version 8.4.0 "Harrier"
+ */
+
+#pragma once
+
+#include "CIdealGas.hpp"
+#include <array>
+
+/*!
+ * \class CNasaGas
+ * \brief Child class for defining the NASA polynomial gas model.
+ */
+class CNasaGas : public CIdealGas {
+ protected:
+  array<su2double, 7> CpLowPolyCoefficientsND{{0.0}};  /*!< \brief NASA low temp polynomial coefficients. */
+  array<su2double, 7> CpHighPolyCoefficientsND{{0.0}}; /*!< \brief NASA high temp polynomial coefficients. */
+  su2double TransitionTemperatureND{0.0};              /*!< \brief NASA transition temperature. */
+
+ public:
+  /*!
+   * \brief Constructor of the class.
+   */
+  CNasaGas(su2double gamma, su2double R, bool CompEntropy = true);
+
+  /*!
+   * \brief Set specific heat Cp model.
+   */
+  void SetCpModel(const CConfig* config, su2double val_Temperature_Ref) override;
+
+  /*!
+   * \brief Set the Dimensionless State using Density and Internal Energy
+   */
+  void SetTDState_rhoe(su2double rho, su2double e) override;
+
+  /*!
+   * \brief Set the Dimensionless State using Pressure and Temperature
+   */
+  void SetTDState_PT(su2double P, su2double T) override;
+
+  /*!
+   * \brief Set the Dimensionless State using Pressure and Density
+   */
+  void SetTDState_Prho(su2double P, su2double rho) override;
+
+  /*!
+   * \brief Set the Dimensionless State using Density and Temperature
+   */
+  void SetTDState_rhoT(su2double rho, su2double T) override;
+
+  /*!
+   * \brief Set the Dimensionless State using Pressure and Entropy
+   */
+  void SetTDState_Ps(su2double P, su2double s) override;
+
+ protected:
+  /*!
+   * \brief Calculate Cp, Enthalpy and Entropy from Temperature.
+   */
+  void ComputeProperties_T(su2double T);
+};

SU2_CFD/src/fluid/CFluidModel.cpp

@@ -45,6 +45,7 @@
 #include "../../include/fluid/CCoolPropViscosity.hpp"
 #include "../../include/fluid/CConstantLewisDiffusivity.hpp"
 #include "../../include/fluid/CCoolPropConductivity.hpp"
+#include "../../include/fluid/CNasaGas.hpp"
 
 unique_ptr<CViscosityModel> CFluidModel::MakeLaminarViscosityModel(const CConfig* config, unsigned short iSpecies) {
   switch (config->GetKind_ViscosityModel()) {

SU2_CFD/src/fluid/CNasaGas.cpp

@@ -0,0 +1,146 @@
+/*!
+ * \file CNasaGas.cpp
+ * \brief Source of the NASA polynomial gas model.
+ * \author Om Giri
+ * \version 8.4.0 "Harrier"
+ */
+
+#include "../../include/fluid/CNasaGas.hpp"
+
+CNasaGas::CNasaGas(su2double gamma, su2double R, bool CompEntropy) : CIdealGas(gamma, R, CompEntropy) {
+  // Constructor inherited from CIdealGas
+}
+
+void CNasaGas::SetCpModel(const CConfig* config, su2double val_Temperature_Ref) {
+  for (int i = 0; i < 7; ++i) {
+    CpLowPolyCoefficientsND[i] = config->GetNASA_CpLowPolyCoeff(i);
+    CpHighPolyCoefficientsND[i] = config->GetNASA_CpHighPolyCoeff(i);
+  }
+  TransitionTemperatureND = config->GetNASA_TransitionTemperature() / val_Temperature_Ref;
+}
+
+void CNasaGas::ComputeProperties_T(su2double T) {
+  Temperature = T;
+  const array<su2double, 7>* c;
+
+  if (T < TransitionTemperatureND) {
+    c = &CpLowPolyCoefficientsND;
+  } else {
+    c = &CpHighPolyCoefficientsND;
+  }
+
+  /*--- NASA polynomials:
+     Cp/R = a1*T^-2 + a2*T^-1 + a3 + a4*T + a5*T^2 + a6*T^3 + a7*T^4
+     H/RT = -a1*T^-2 + a2*log(T)/T + a3 + a4*T/2 + a5*T^2/3 + a6*T^3/4 + a7*T^4/5 + a8/T
+     S/R  = -a1*T^-2/2 - a2*T^-1 + a3*log(T) + a4*T + a5*T^2/2 + a6*T^3/3 + a7*T^4/4 + a9
+
+     Note: SU2 usually uses a 7-coefficient form:
+     Cp/R = a0 + a1*T + a2*T^2 + a3*T^3 + a4*T^4
+     H/RT = a0 + a1*T/2 + a2*T^2/3 + a3*T^3/4 + a4*T^4/5 + a5/T
+     S/R  = a0*log(T) + a1*T + a2*T^2/2 + a3*T^3/3 + a4*T^4/4 + a6
+
+     Wait, standard NASA 7-coefficient format is:
+     Cp/R = a0 + a1*T + a2*T^2 + a3*T^3 + a4*T^4
+     H/RT = a0 + a1*T/2 + a2*T^2/3 + a3*T^3/4 + a4*T^4/5 + a5/T
+     S/R  = a0*lnT + a1*T + a2*T^2/2 + a3*T^3/3 + a4*T^4/4 + a6
+
+     I will use this standard 7-coefficient format.
+  ---*/
+
+  Cp = ((*c)[0] + (*c)[1]*T + (*c)[2]*T*T + (*c)[3]*T*T*T + (*c)[4]*T*T*T*T) * Gas_Constant;
+  Enthalpy = ((*c)[0] + (*c)[1]*T/2.0 + (*c)[2]*T*T/3.0 + (*c)[3]*T*T*T/4.0 + (*c)[4]*T*T*T*T/5.0 + (*c)[5]/T) * T * Gas_Constant;
+
+  if (ComputeEntropy) {
+    Entropy = ((*c)[0]*log(T) + (*c)[1]*T + (*c)[2]*T*T/2.0 + (*c)[3]*T*T*T/3.0 + (*c)[4]*T*T*T*T/4.0 + (*c)[6]) * Gas_Constant;
+    // Note: This is partial entropy (temperature dependent part). Pressure part added in SetTDState_rhoe if needed.
+    // However, SU2 CIdealGas includes log(1/rho) term.
+  }
+
+  // Cv = Cp - R
+  Cv = Cp - Gas_Constant;
+  // Gamma = Cp / Cv
+  Gamma = Cp / Cv;
+  Gamma_Minus_One = Gamma - 1.0;
+}
+
+void CNasaGas::SetTDState_rhoe(su2double rho, su2double e) {
+  Density = rho;
+  StaticEnergy = e;
+
+  /*--- Solve for T using Newton-Raphson: e(T) = h(T) - R*T ---*/
+  su2double T_iter = Temperature;
+  if (T_iter <= 0.0) T_iter = 1.0; // Initial guess
+
+  for (int i = 0; i < 10; ++i) {
+    ComputeProperties_T(T_iter);
+    su2double e_iter = Enthalpy - Gas_Constant * T_iter;
+    su2double de_dT = Cv; // de/dT = Cv for ideal gas even with variable Cp
+    su2double deltaT = (e - e_iter) / de_dT;
+    T_iter += deltaT;
+    if (abs(deltaT) < 1e-8 * T_iter) break;
+  }
+
+  Temperature = T_iter;
+  ComputeProperties_T(Temperature);
+
+  Pressure = Density * Gas_Constant * Temperature;
+  SoundSpeed2 = Gamma * Pressure / Density;
+
+  /*--- Derivatives ---*/
+  dPde_rho = Density * Gas_Constant / Cv;
+  dPdrho_e = Gas_Constant * Temperature - StaticEnergy * dPde_rho; // From P = rho*R*T and e = e(T)
+
+  dTde_rho = 1.0 / Cv;
+  dTdrho_e = 0.0;
+
+  if (ComputeEntropy) {
+    // Entropy was S(T). Now add pressure/density part: S = S(T) - R*ln(rho*R) + R*ln(R_ref?)
+    // CIdealGas uses: Entropy = (1.0 / Gamma_Minus_One * log(Temperature) + log(1.0 / Density)) * Gas_Constant;
+    // For NASA, S(T) already has the log(T) part.
+    Entropy += log(1.0 / Density) * Gas_Constant;
+  }
+}
+
+void CNasaGas::SetTDState_PT(su2double P, su2double T) {
+  Pressure = P;
+  Temperature = T;
+  ComputeProperties_T(T);
+  Density = P / (Gas_Constant * T);
+  StaticEnergy = Enthalpy - Gas_Constant * T;
+  SetTDState_rhoe(Density, StaticEnergy);
+}
+
+void CNasaGas::SetTDState_Prho(su2double P, su2double rho) {
+  Pressure = P;
+  Density = rho;
+  Temperature = P / (Density * Gas_Constant);
+  ComputeProperties_T(Temperature);
+  StaticEnergy = Enthalpy - Gas_Constant * Temperature;
+  SetTDState_rhoe(Density, StaticEnergy);
+}
+
+void CNasaGas::SetTDState_rhoT(su2double rho, su2double T) {
+  Density = rho;
+  Temperature = T;
+  ComputeProperties_T(T);
+  Pressure = Density * Gas_Constant * T;
+  StaticEnergy = Enthalpy - Gas_Constant * T;
+  SetTDState_rhoe(Density, StaticEnergy);
+}
+
+void CNasaGas::SetTDState_Ps(su2double P, su2double s) {
+  /*--- This would require another nested iteration to find T from P and s.
+     For now, we can use an approximate T or a simple Newton solver. ---*/
+  su2double T_iter = Temperature;
+  if (T_iter <= 0.0) T_iter = 1.0;
+
+  for (int i = 0; i < 10; ++i) {
+    ComputeProperties_T(T_iter);
+    su2double s_iter = Entropy + log(Gas_Constant * T_iter / P) * Gas_Constant;
+    su2double ds_dT = Cp / T_iter;
+    su2double deltaT = (s - s_iter) / ds_dT;
+    T_iter += deltaT;
+    if (abs(deltaT) < 1e-8 * T_iter) break;
+  }
+  SetTDState_PT(P, T_iter);
+}

SU2_CFD/src/meson.build

@@ -13,7 +13,8 @@ su2_cfd_src += files(['fluid/CFluidModel.cpp',
                       'fluid/CNEMOGas.cpp',
                       'fluid/CMutationTCLib.cpp',
                       'fluid/CSU2TCLib.cpp',
-                      'fluid/CDataDrivenFluid.cpp'])
+                      'fluid/CDataDrivenFluid.cpp',
+                      'fluid/CNasaGas.cpp'])
 
 su2_cfd_src += files(['output/COutputFactory.cpp',
                       'output/CAdjElasticityOutput.cpp',

SU2_CFD/src/output/CFlowOutput.cpp

@@ -250,8 +250,19 @@ void CFlowOutput::SetAnalyzeSurface(const CSolver* const*solver, const CGeometry
             Mach              = sqrt(Velocity2)/SoundSpeed;
             Temperature       = Pressure / (Gas_Constant * Density);
             Enthalpy          = flow_nodes->GetEnthalpy(iPoint);
-            TotalTemperature  = Temperature * (1.0 + Mach * Mach * 0.5 * (Gamma - 1.0));
-            TotalPressure     = Pressure * pow( 1.0 + Mach * Mach * 0.5 * (Gamma - 1.0), Gamma / (Gamma - 1.0));
+            if (config->GetKind_FluidModel() == NASA_GAS) {
+              TotalTemperature = Temperature;
+              TotalPressure = Pressure;
+              // For variable Cp, we'd need an iterative or approximate approach for stagnation properties.
+              // For now, use a simplified approach or leave as static if complex.
+              // A common approximation for small Mach: Tt = T + V^2/(2*Cp)
+              TotalTemperature  = Temperature + 0.5 * Velocity2 / flow_nodes->GetSpecificHeatCp(iPoint);
+              // Pt approximation: Pt = P * (Tt/T)^(Cp/R)
+              TotalPressure     = Pressure * pow(TotalTemperature / Temperature, flow_nodes->GetSpecificHeatCp(iPoint) / Gas_Constant);
+            } else {
+              TotalTemperature  = Temperature * (1.0 + Mach * Mach * 0.5 * (Gamma - 1.0));
+              TotalPressure     = Pressure * pow( 1.0 + Mach * Mach * 0.5 * (Gamma - 1.0), Gamma / (Gamma - 1.0));
+            }
           }
 
           /*--- Compute the mass Surface_MassFlow ---*/

SU2_CFD/src/solvers/CEulerSolver.cpp

@@ -881,6 +881,13 @@ void CEulerSolver::SetNondimensionalization(CConfig *config, unsigned short iMes
 
       auxFluidModel = new CDataDrivenFluid(config);
 
+      break;
+
+    case NASA_GAS:
+
+      auxFluidModel = new CNasaGas(Gamma, config->GetGas_Constant());
+      auxFluidModel->SetCpModel(config, 1.0);
+
       break;
     case COOLPROP:
 
@@ -1133,6 +1140,11 @@ void CEulerSolver::SetNondimensionalization(CConfig *config, unsigned short iMes
         FluidModel[thread] = new CDataDrivenFluid(config, false);
         break;
 
+      case NASA_GAS:
+        FluidModel[thread] = new CNasaGas(Gamma, Gas_ConstantND);
+        GetFluidModel()->SetCpModel(config, config->GetTemperature_Ref());
+        break;
+
       case COOLPROP:
         FluidModel[thread] = new CCoolProp(config->GetFluid_Name());
         break;
@@ -1311,6 +1323,9 @@ void CEulerSolver::SetNondimensionalization(CConfig *config, unsigned short iMes
     case PR_GAS:
       ModelTable << "PR_GAS";
       break;
+    case NASA_GAS:
+      ModelTable << "NASA_GAS";
+      break;
     case COOLPROP:
       ModelTable << "CoolProp library";
       break;

SU2_CFD/src/solvers/CFEM_DG_EulerSolver.cpp

@@ -33,6 +33,7 @@
 #include "../../include/fluid/CPengRobinson.hpp"
 #include "../../include/fluid/CCoolProp.hpp"
 #include "../../include/fluid/CDataDrivenFluid.hpp"
+#include "../../include/fluid/CNasaGas.hpp"
 
 enum {
 SIZE_ARR_NORM = 8
@@ -914,6 +915,21 @@ void CFEM_DG_EulerSolver::SetNondimensionalization(CConfig        *config,
       }
 
       break;
+
+    case NASA_GAS:
+      FluidModel = new CNasaGas(Gamma, config->GetGas_Constant(), config->GetCompute_Entropy());
+      FluidModel->SetCpModel(config, 1.0);
+      if (free_stream_temp) {
+        FluidModel->SetTDState_PT(Pressure_FreeStream, Temperature_FreeStream);
+        Density_FreeStream = FluidModel->GetDensity();
+        config->SetDensity_FreeStream(Density_FreeStream);
+      }
+      else {
+        FluidModel->SetTDState_Prho(Pressure_FreeStream, Density_FreeStream );
+        Temperature_FreeStream = FluidModel->GetTemperature();
+        config->SetTemperature_FreeStream(Temperature_FreeStream);
+      }
+      break;
   }
 
   Mach2Vel_FreeStream = FluidModel->GetSoundSpeed();
@@ -1121,6 +1137,12 @@ void CFEM_DG_EulerSolver::SetNondimensionalization(CConfig        *config,
       FluidModel = new CDataDrivenFluid(config);
       FluidModel->SetEnergy_Prho(Pressure_FreeStreamND, Density_FreeStreamND);
       break;
+
+    case NASA_GAS:
+      FluidModel = new CNasaGas(Gamma, Gas_ConstantND, config->GetCompute_Entropy());
+      FluidModel->SetCpModel(config, config->GetTemperature_Ref());
+      FluidModel->SetEnergy_Prho(Pressure_FreeStreamND, Density_FreeStreamND);
+      break;
   }
 
   Energy_FreeStreamND = FluidModel->GetStaticEnergy() + 0.5*ModVel_FreeStreamND*ModVel_FreeStreamND;
@@ -1291,6 +1313,9 @@ void CFEM_DG_EulerSolver::SetNondimensionalization(CConfig        *config,
     case PR_GAS:
       ModelTable << "PR_GAS";
       break;
+    case NASA_GAS:
+      ModelTable << "NASA_GAS";
+      break;
     case COOLPROP:
       ModelTable << "CoolProp library";
       break;