Feature: NASA polynomial gas model implementation (CNasaGas)

Common/include/CConfig.hpp

@@ -326,6 +326,7 @@ class CConfig {
   su2double *Engine_Area;                    /*!< \brief Specified engine area for nacelle boundaries. */
   su2double *Outlet_Pressure;                /*!< \brief Specified back pressures (static) for outlet boundaries. */
   su2double *Isothermal_Temperature;         /*!< \brief Specified isothermal wall temperatures (static). */
+  mutable std::map<unsigned short, std::map<unsigned long, su2double>> Marker_Custom_Temperature; /*!< \brief Spatially varying isothermal wall temperatures. */
   su2double *HeatTransfer_Coeff;             /*!< \brief Specified heat transfer coefficients. */
   su2double *HeatTransfer_WallTemp;          /*!< \brief Specified temperatures at infinity alongside heat transfer coefficients. */
   su2double *Heat_Flux;                      /*!< \brief Specified wall heat fluxes. */
@@ -911,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.  */
@@ -1124,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 {
@@ -7426,6 +7433,22 @@ class CConfig {
    */
   su2double GetIsothermal_Temperature(const string& val_index) const;
 
+  /*!
+   * \brief Get the custom wall temperature (static) at an isothermal boundary node.
+   * \param[in] iMarker - Marker index.
+   * \param[in] iVertex - Vertex index on the marker.
+   * \return The custom wall temperature (returns -1.0 if not set).
+   */
+  su2double GetCustom_Temperature(unsigned short iMarker, unsigned long iVertex) const;
+
+  /*!
+   * \brief Set the custom wall temperature (static) at an isothermal boundary node.
+   * \param[in] iMarker - Marker index.
+   * \param[in] iVertex - Vertex index on the marker.
+   * \param[in] val - The custom wall temperature.
+   */
+  void SetCustom_Temperature(unsigned short iMarker, unsigned long iVertex, su2double val);
+
   /*!
    * \brief Get the wall heat flux on a constant heat flux boundary.
    * \param[in] val_index - Index corresponding to the constant heat flux boundary.

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);
@@ -9608,6 +9614,21 @@ su2double CConfig::GetIsothermal_Temperature(const string& val_marker) const {
   return Isothermal_Temperature[0];
 }
 
+su2double CConfig::GetCustom_Temperature(unsigned short iMarker, unsigned long iVertex) const {
+  auto it_marker = Marker_Custom_Temperature.find(iMarker);
+  if (it_marker != Marker_Custom_Temperature.end()) {
+    auto it_vertex = it_marker->second.find(iVertex);
+    if (it_vertex != it_marker->second.end()) {
+      return it_vertex->second;
+    }
+  }
+  return -1.0;
+}
+
+void CConfig::SetCustom_Temperature(unsigned short iMarker, unsigned long iVertex, su2double val) {
+  Marker_Custom_Temperature[iMarker][iVertex] = val;
+}
+
 su2double CConfig::GetWall_HeatFlux(const string& val_marker) const {
 
   for (unsigned short iMarker_HeatFlux = 0; iMarker_HeatFlux < nMarker_HeatFlux; iMarker_HeatFlux++)

SU2_CFD/include/drivers/CDriver.hpp

@@ -555,6 +555,30 @@ class CDriver : public CDriverBase {
    */
   void SetMarkerTranslationRate(unsigned short iMarker, passivedouble vel_x, passivedouble vel_y, passivedouble vel_z);
 
+  /*!
+   * \brief Set a custom temperature at a boundary vertex.
+   * \param[in] iMarker - Index of the boundary marker.
+   * \param[in] iVertex - Index of the boundary vertex.
+   * \param[in] Temperature - Temperature value.
+   */
+  void SetMarkerCustomTemperature(unsigned short iMarker, unsigned long iVertex, su2double Temperature);
+
+  /*!
+   * \brief Get the local speed of sound at a boundary vertex.
+   * \param[in] iMarker - Index of the boundary marker.
+   * \param[in] iVertex - Index of the boundary vertex.
+   * \return Local speed of sound.
+   */
+  passivedouble GetMarkerLocalSpeedOfSound(unsigned short iMarker, unsigned long iVertex);
+
+  /*!
+   * \brief Get the Mach number at a boundary vertex.
+   * \param[in] iMarker - Index of the boundary marker.
+   * \param[in] iVertex - Index of the boundary vertex.
+   * \return Mach number.
+   */
+  passivedouble GetMarkerMachNumber(unsigned short iMarker, unsigned long iVertex);
+
   /*!
    * \brief Get the Freestream Density for nondimensionalization
    * \return Freestream Density

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',