diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index dc96c21ad3b..92746709f4d 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -1121,6 +1121,9 @@ class CConfig { array cp_polycoeffs{{0.0}}; /*!< \brief Array for specific heat polynomial coefficients. */ array mu_polycoeffs{{0.0}}; /*!< \brief Array for viscosity polynomial coefficients. */ array kt_polycoeffs{{0.0}}; /*!< \brief Array for thermal conductivity polynomial coefficients. */ + + bool Cp_NASA_Format{false}; /*!< \brief Flag to use NASA polynomial format for Cp. */ + bool Body_Force; /*!< \brief Flag to know if a body force is included in the formulation. */ struct CMUSCLRampParam { @@ -2592,6 +2595,21 @@ class CConfig { */ void SetEnergy_Ref(su2double val_energy_ref) { Energy_Ref = val_energy_ref; } + /*! + * \brief Get the flag for using NASA polynomial format for Cp. + * \return True if NASA format is used. + */ + bool GetCp_NASA_Format(void) const { return Cp_NASA_Format; } + + /*! + * \brief Set the flag for using NASA polynomial format for Cp. + * \param[in] val - Value of the flag. + */ + void SetCp_NASA_Format(bool val) { Cp_NASA_Format = val; } + + // Setter for TemperatureLimits (for testing) + void SetTemperatureLimits(unsigned short val_index, su2double val) { TemperatureLimits[val_index] = val; } + /*! * \brief Set the reference Omega for nondimensionalization. * \param[in] val_omega_ref - Value of the reference omega. @@ -4255,6 +4273,7 @@ class CConfig { * \param[in] val_index - Index of the array with all polynomial coefficients. */ void SetCp_PolyCoeffND(su2double val_coeff, unsigned short val_index) { CpPolyCoefficientsND[val_index] = val_coeff; } + void SetCp_PolyCoeff(unsigned short val_index, su2double val) { cp_polycoeffs[val_index] = val; } /*! * \brief Set the temperature polynomial coefficient for viscosity. diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 690350fae05..7bfff7fb5c9 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -201,7 +201,7 @@ constexpr passivedouble COLORING_EFF_THRESH = 0.875; /*!< \brief Below this val assume a quartic form (5 coefficients). For example, 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_POLY_COEFFS = 10; /*!< \brief Number of coefficients in temperature polynomial fits for fluid models. */ /*! * \brief Boolean answers diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 8a95237ac3f..0318b7e8803 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -1333,11 +1333,14 @@ void CConfig::SetConfig_Options() { /*--- Options related to temperature polynomial coefficients for fluid models. ---*/ /* DESCRIPTION: Definition of the temperature polynomial coefficients for specific heat Cp. */ - addDoubleArrayOption("CP_POLYCOEFFS", N_POLY_COEFFS, false, cp_polycoeffs.data()); + addDoubleArrayOption("CP_POLYCOEFFS", N_POLY_COEFFS, true, cp_polycoeffs.data()); /* DESCRIPTION: Definition of the temperature polynomial coefficients for specific heat Cp. */ - addDoubleArrayOption("MU_POLYCOEFFS", N_POLY_COEFFS, false, mu_polycoeffs.data()); + addDoubleArrayOption("MU_POLYCOEFFS", N_POLY_COEFFS, true, mu_polycoeffs.data()); /* DESCRIPTION: Definition of the temperature polynomial coefficients for specific heat Cp. */ - addDoubleArrayOption("KT_POLYCOEFFS", N_POLY_COEFFS, false, kt_polycoeffs.data()); + addDoubleArrayOption("KT_POLYCOEFFS", N_POLY_COEFFS, true, kt_polycoeffs.data()); + + /* DESCRIPTION: Flag to use NASA polynomial format for Cp in incompressible ideal gas model. */ + addBoolOption("CP_NASA", Cp_NASA_Format, false); /*!\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); diff --git a/SU2_CFD/include/fluid/CIncIdealGasNASA.hpp b/SU2_CFD/include/fluid/CIncIdealGasNASA.hpp new file mode 100644 index 00000000000..93693f910de --- /dev/null +++ b/SU2_CFD/include/fluid/CIncIdealGasNASA.hpp @@ -0,0 +1,223 @@ +/*! + * \file CIncIdealGasNASA.hpp + * \brief Defines the incompressible Ideal Gas model with NASA polynomials for Cp. + * \author Pratyksh Gupta + * \version 8.4.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#pragma once + +#include +#include + +#include "CFluidModel.hpp" + +/*! + * \class CIncIdealGasNASA + * \brief Child class for defining an incompressible ideal gas model with NASA polynomials. + * \author Pratyksh Gupta + * + * Implements NASA 9-coefficient polynomial format for thermodynamic properties with backward compatibility for NASA-7. + * Ref: McBride, B.J., Zehe, M.J., and Gordon, S., "NASA Glenn Coefficients for Calculating Thermodynamic Properties of Individual Species", NASA/TP-2002-211556, 2002. + * + * NASA-9 format (full): + * 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*ln(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*ln(T) + a4*T + a5*T^2/2 + a6*T^3/3 + a7*T^4/4 + a9 + * + * NASA-7 format (subset): Set a1=a2=0 to recover the traditional 7-coefficient format. + * + * Uses a single temperature range provided via CP_POLYCOEFFS (indices 0-8). + */ +template +class CIncIdealGasNASA final : public CFluidModel { + public: + /*! + * \brief Constructor of the class. + */ + CIncIdealGasNASA(su2double val_gas_constant, su2double val_operating_pressure, su2double val_Temperature_Ref, su2double val_Ref_Temp_Dim = 1.0) { + /*--- In the incompressible ideal gas model, the thermodynamic pressure is decoupled + from the governing equations and held constant. The density is therefore only a + function of temperature variations. The gas is incompressible, so Cp = Cv (gamma = 1). ---*/ + Gas_Constant = val_gas_constant; + Pressure = val_operating_pressure; + Gamma = 1.0; + Std_Ref_Temp_ND = val_Temperature_Ref; + Ref_Temp_Dim = val_Ref_Temp_Dim; + } + + /*! + * \brief Set the NASA polynomial coefficients for variable Cp. + * \param[in] config - configuration container for the problem. + */ + void SetCpModel(const CConfig* config, su2double val_Temperature_Ref) override { + + /*--- Read NASA coefficients from the standard polynomial coefficient array (CP_POLYCOEFFS). + NASA-9 format uses indices 0-8: + Indices 0-1: Inverse temperature terms (a1*T^-2, a2*T^-1) + Indices 2-6: Polynomial terms (a3, a4*T, a5*T^2, a6*T^3, a7*T^4) + Index 7: Enthalpy constant (a8) + Index 8: Entropy constant (a9) + + For NASA-7 compatibility: Set indices 0-1 to zero. ---*/ + for (int i = 0; i < N_COEFFS; ++i) { + if (i < config->GetnPolyCoeffs()) { + coeffs_[i] = config->GetCp_PolyCoeff(i); + } else { + coeffs_[i] = 0.0; + } + } + + Temperature_Min = config->GetTemperatureLimits(0) / Ref_Temp_Dim; + Temperature_Max = config->GetTemperatureLimits(1) / Ref_Temp_Dim; + } + + /*! + * \brief Set the Dimensionless State using Temperature. + * \param[in] t - Temperature value at the point. + */ + void SetTDState_T(su2double t, const su2double *val_scalars = nullptr) override { + using std::pow; + using std::log; + Temperature = t; + Density = Pressure / (Temperature * Gas_Constant); + + const su2double a1 = coeffs_[0]; + const su2double a2 = coeffs_[1]; + const su2double a3 = coeffs_[2]; + const su2double a4 = coeffs_[3]; + const su2double a5 = coeffs_[4]; + const su2double a6 = coeffs_[5]; + const su2double a7 = coeffs_[6]; + const su2double a8 = coeffs_[7]; + const su2double a9 = coeffs_[8]; + + // Convert to dimensional temperature for polynomial evaluation (NASA coeffs expect Kelvin) + const su2double T_dim = t * Ref_Temp_Dim; + + // NASA-9: Cp/R = a1*T^-2 + a2*T^-1 + a3 + a4*T + a5*T^2 + a6*T^3 + a7*T^4 + su2double Cp_over_R = a1 * pow(T_dim, -2.0) + a2 * pow(T_dim, -1.0) + a3 + a4 * T_dim + a5 * pow(T_dim, 2.0) + a6 * pow(T_dim, 3.0) + a7 * pow(T_dim, 4.0); + + Cp = Cp_over_R * Gas_Constant; + + // NASA-9: H/(RT) = -a1*T^-2 + a2*ln(T)/T + a3 + a4*T/2 + a5*T^2/3 + a6*T^3/4 + a7*T^4/5 + a8/T + su2double H_over_RT = -a1 * pow(T_dim, -2.0) + a2 * log(T_dim) / T_dim + a3 + a4 * T_dim / 2.0 + a5 * pow(T_dim, 2.0) / 3.0 + + a6 * pow(T_dim, 3.0) / 4.0 + a7 * pow(T_dim, 4.0) / 5.0 + a8 / T_dim; + + Enthalpy = H_over_RT * Gas_Constant * t; + + // NASA-9: S/R = -a1*T^-2/2 - a2*T^-1 + a3*ln(T) + a4*T + a5*T^2/2 + a6*T^3/3 + a7*T^4/4 + a9 + su2double S_over_R = -a1 * pow(T_dim, -2.0) / 2.0 - a2 * pow(T_dim, -1.0) + a3 * log(T_dim) + a4 * T_dim + + a5 * pow(T_dim, 2.0) / 2.0 + a6 * pow(T_dim, 3.0) / 3.0 + a7 * pow(T_dim, 4.0) / 4.0 + a9; + + Entropy = S_over_R * Gas_Constant; + + Cv = Cp / Gamma; + } + + /*! + * \brief Set the Dimensionless State using enthalpy. + * \param[in] val_enthalpy - Enthalpy value at the point. + */ + void SetTDState_h(su2double val_enthalpy, const su2double* val_scalars = nullptr) override { + using std::pow; + using std::log; + using std::fmin; + using std::fmax; + using std::fabs; + Enthalpy = val_enthalpy; + + const su2double toll = 1e-5; + su2double temp_iter = (Temperature_Min + Temperature_Max) / 2.0; /* Start in middle of allowed range */ + if (temp_iter < 1.0) temp_iter = 300.0; /* Fallback if limits are not set or zero */ + + const su2double a1 = coeffs_[0]; + const su2double a2 = coeffs_[1]; + const su2double a3 = coeffs_[2]; + const su2double a4 = coeffs_[3]; + const su2double a5 = coeffs_[4]; + const su2double a6 = coeffs_[5]; + const su2double a7 = coeffs_[6]; + const su2double a8 = coeffs_[7]; + const su2double a9 = coeffs_[8]; + + su2double Cp_iter = 0.0; + su2double delta_temp_iter = 1e10; + su2double delta_enthalpy_iter; + const int counter_limit = 50; + int counter = 0; + + while ((fabs(delta_temp_iter) > toll) && (counter++ < counter_limit)) { + + const su2double T_dim = temp_iter * Ref_Temp_Dim; + + // NASA-9: Cp/R = a1*T^-2 + a2*T^-1 + a3 + a4*T + a5*T^2 + a6*T^3 + a7*T^4 + su2double Cp_over_R = a1 * pow(T_dim, -2.0) + a2 * pow(T_dim, -1.0) + a3 + a4 * T_dim + a5 * pow(T_dim, 2.0) + + a6 * pow(T_dim, 3.0) + a7 * pow(T_dim, 4.0); + + Cp_iter = Cp_over_R * Gas_Constant; + + // NASA-9: H/(RT) = -a1*T^-2 + a2*ln(T)/T + a3 + a4*T/2 + a5*T^2/3 + a6*T^3/4 + a7*T^4/5 + a8/T + su2double H_over_RT = -a1 * pow(T_dim, -2.0) + a2 * log(T_dim) / T_dim + a3 + a4 * T_dim / 2.0 + + a5 * pow(T_dim, 2.0) / 3.0 + a6 * pow(T_dim, 3.0) / 4.0 + + a7 * pow(T_dim, 4.0) / 5.0 + a8 / T_dim; + + su2double Enthalpy_iter = H_over_RT * Gas_Constant * temp_iter; + + delta_enthalpy_iter = Enthalpy - Enthalpy_iter; + delta_temp_iter = delta_enthalpy_iter / Cp_iter; + temp_iter += delta_temp_iter; + + temp_iter = fmin(fmax(Temperature_Min, temp_iter), Temperature_Max); + } + + Temperature = temp_iter; + Cp = Cp_iter; + + // Evaluate Entropy at final state + const su2double T_dim_final = Temperature * Ref_Temp_Dim; + su2double S_over_R = -a1 * pow(T_dim_final, -2.0) / 2.0 - a2 * pow(T_dim_final, -1.0) + a3 * log(T_dim_final) + + a4 * T_dim_final + a5 * pow(T_dim_final, 2.0) / 2.0 + a6 * pow(T_dim_final, 3.0) / 3.0 + + a7 * pow(T_dim_final, 4.0) / 4.0 + a9; + Entropy = S_over_R * Gas_Constant; + + if (counter == counter_limit) { + if (SU2_MPI::GetRank() == MASTER_NODE) + std::cout << "Warning: Newton-Raphson exceeds max. iterations in temperature computation (NASA Model).\n"; + } + + Density = Pressure / (Temperature * Gas_Constant); + Cv = Cp / Gamma; + } + + private: + su2double Gas_Constant{0.0}; /*!< \brief Specific Gas Constant. */ + su2double Gamma{0.0}; /*!< \brief Ratio of specific heats. */ + su2double Std_Ref_Temp_ND{0.0}; /*!< \brief Nondimensional standard reference temperature for enthalpy. */ + su2double Ref_Temp_Dim{1.0}; /*!< \brief Dimensional reference temperature for evaluating polynomials. */ + + std::array coeffs_; /*!< \brief NASA polynomial coefficients. */ + + su2double Temperature_Min{0.0}; /*!< \brief Minimum temperature allowed. */ + su2double Temperature_Max{1e10}; /*!< \brief Maximum temperature allowed. */ +}; diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index c99f87896c2..03daa5915bd 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -30,6 +30,7 @@ #include "../../include/fluid/CConstantDensity.hpp" #include "../../include/fluid/CIncIdealGas.hpp" #include "../../include/fluid/CIncIdealGasPolynomial.hpp" +#include "../../include/fluid/CIncIdealGasNASA.hpp" #include "../../include/variables/CIncNSVariable.hpp" #include "../../../Common/include/toolboxes/geometry_toolbox.hpp" #include "../../include/fluid/CFluidScalar.hpp" @@ -307,11 +308,21 @@ void CIncEulerSolver::SetNondimensionalization(CConfig *config, unsigned short i config->SetGas_Constant(UNIVERSAL_GAS_CONSTANT/(config->GetMolecular_Weight()/1000.0)); Pressure_Thermodynamic = Density_FreeStream*Temperature_FreeStream*config->GetGas_Constant(); - auxFluidModel = new CIncIdealGasPolynomial(config->GetGas_Constant(), Pressure_Thermodynamic, STD_REF_TEMP); + if (config->GetCp_NASA_Format()) { + /*--- Use the NASA polynomial model (supports NASA-9 format with 9 coefficients, backward compatible with NASA-7) ---*/ + auxFluidModel = new CIncIdealGasNASA<9>(config->GetGas_Constant(), Pressure_Thermodynamic, STD_REF_TEMP); + } else { + /*--- Use the standard polynomial model ---*/ + auxFluidModel = new CIncIdealGasPolynomial(config->GetGas_Constant(), Pressure_Thermodynamic, STD_REF_TEMP); + } + if (viscous) { /*--- Variable Cp model via polynomial. ---*/ - for (iVar = 0; iVar < config->GetnPolyCoeffs(); iVar++) - config->SetCp_PolyCoeffND(config->GetCp_PolyCoeff(iVar), iVar); + if (!config->GetCp_NASA_Format()) { + for (iVar = 0; iVar < config->GetnPolyCoeffs(); iVar++) + config->SetCp_PolyCoeffND(config->GetCp_PolyCoeff(iVar), iVar); + } + /*--- Both models refer to config for coeffs, so setup is handled inside SetCpModel ---*/ auxFluidModel->SetCpModel(config, Temperature_FreeStream); } auxFluidModel->SetTDState_T(Temperature_FreeStream); @@ -499,12 +510,19 @@ void CIncEulerSolver::SetNondimensionalization(CConfig *config, unsigned short i break; case INC_IDEAL_GAS_POLY: - fluidModel = new CIncIdealGasPolynomial(Gas_ConstantND, Pressure_ThermodynamicND, STD_REF_TEMP / config->GetTemperature_Ref()); + if (config->GetCp_NASA_Format()) { + fluidModel = new CIncIdealGasNASA<9>(Gas_ConstantND, Pressure_ThermodynamicND, STD_REF_TEMP / config->GetTemperature_Ref(), config->GetTemperature_Ref()); + } else { + fluidModel = new CIncIdealGasPolynomial(Gas_ConstantND, Pressure_ThermodynamicND, STD_REF_TEMP / config->GetTemperature_Ref()); + } + if (viscous) { /*--- Variable Cp model via polynomial. ---*/ - config->SetCp_PolyCoeffND(config->GetCp_PolyCoeff(0)/Gas_Constant_Ref, 0); - for (iVar = 1; iVar < config->GetnPolyCoeffs(); iVar++) - config->SetCp_PolyCoeffND(config->GetCp_PolyCoeff(iVar)*pow(Temperature_Ref,iVar)/Gas_Constant_Ref, iVar); + if (!config->GetCp_NASA_Format()) { + config->SetCp_PolyCoeffND(config->GetCp_PolyCoeff(0)/Gas_Constant_Ref, 0); + for (iVar = 1; iVar < config->GetnPolyCoeffs(); iVar++) + config->SetCp_PolyCoeffND(config->GetCp_PolyCoeff(iVar)*pow(Temperature_Ref,iVar)/Gas_Constant_Ref, iVar); + } fluidModel->SetCpModel(config, Temperature_FreeStreamND); } fluidModel->SetTDState_T(Temperature_FreeStreamND); diff --git a/UnitTests/SU2_CFD/fluid/CIncIdealGasNASA_tests.cpp b/UnitTests/SU2_CFD/fluid/CIncIdealGasNASA_tests.cpp new file mode 100644 index 00000000000..1c4ea8d7c57 --- /dev/null +++ b/UnitTests/SU2_CFD/fluid/CIncIdealGasNASA_tests.cpp @@ -0,0 +1,153 @@ +/*! + * \file CIncIdealGasNASA_tests.cpp + * \brief Unit tests for the NASA polynomial fluid model class (Single Range). + * \author Pratyksh Gupta + * \version 8.4.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#include "catch.hpp" +#include +#include +#include "../../../SU2_CFD/include/fluid/CIncIdealGasNASA.hpp" + +// Helper to config an already constructed CConfig object +void ConfigureNASA(CConfig& config) { + // Nitrogen (N2) Low Range Coefficients (200K - 1000K) in NASA-7 format + // NASA-9 indexing: a1=0, a2=0 (no inverse T terms), a3-a7 (polynomial), a8-a9 (integration constants) + su2double nasa_coeffs[9] = { + 0.0, // a1 (T^-2 term, not used in NASA-7) + 0.0, // a2 (T^-1 term, not used in NASA-7) + 3.298677, // a3 (constant term) + 1.4082404e-3, // a4 (T term) + -3.963222e-6, // a5 (T^2 term) + 5.641515e-9, // a6 (T^3 term) + -2.444854e-12, // a7 (T^4 term) + -1020.8999, // a8 (H integration constant) + 3.950372 // a9 (S integration constant) + }; + + config.SetCp_NASA_Format(true); + for (int i = 0; i < 9; ++i) { + config.SetCp_PolyCoeff(i, nasa_coeffs[i]); + } + + // Ensure we set the temperature limits that SetCpModel reads + // Note: These limits are usually set via TEMPERATURE_LIMITS option, here manually. + config.SetTemperatureLimits(0, 200.0); + config.SetTemperatureLimits(1, 1000.0); +} + +// Helper to populate the stringstream with minimal required options +void SetupConfigStream(std::stringstream& ss) { + ss << "SOLVER= EULER\n"; + ss << "MATH_PROBLEM= DIRECT\n"; + ss << "KIND_TURB_MODEL= NONE\n"; + ss << "KIND_TRANS_MODEL= NONE\n"; +} + +TEST_CASE("NASA polynomial Cp calculation", "[CIncIdealGasNASA]") { + const su2double R_N2 = 296.8; + const su2double P = 101325.0; + const su2double T_ref = 298.15; + + std::stringstream ss; + SetupConfigStream(ss); + CConfig config(ss, SU2_COMPONENT::SU2_CFD, false); + ConfigureNASA(config); + + CIncIdealGasNASA<9> nasa_model(R_N2, P, T_ref); + nasa_model.SetCpModel(&config, T_ref); + + SECTION("Valid temperature range (300K)") { + nasa_model.SetTDState_T(300.0); + const su2double cp = nasa_model.GetCp(); + // Cp of N2 at 300K is approx 1040 J/kgK. + CHECK(cp == Approx(1040.0).epsilon(0.01)); + } +} + +TEST_CASE("NASA polynomial enthalpy-temperature inversion", "[CIncIdealGasNASA]") { + const su2double R_N2 = 296.8; + const su2double P = 101325.0; + const su2double T_ref = 298.15; + + std::stringstream ss; + SetupConfigStream(ss); + CConfig config(ss, SU2_COMPONENT::SU2_CFD, false); + ConfigureNASA(config); + + CIncIdealGasNASA<9> nasa_model(R_N2, P, T_ref); + nasa_model.SetCpModel(&config, T_ref); + + SECTION("Inversion Consistency") { + su2double T_target = 800.0; + nasa_model.SetTDState_T(T_target); + su2double H_target = nasa_model.GetEnthalpy(); + + // Reset state to something else + nasa_model.SetTDState_T(300.0); + + // Invert from H_target + nasa_model.SetTDState_h(H_target); + su2double T_recovered = nasa_model.GetTemperature(); + + CHECK(T_recovered == Approx(T_target).margin(0.001)); + } +} + +TEST_CASE("NASA polynomial non-dimensional scaling", "[CIncIdealGasNASA]") { + const su2double R_N2_dim = 296.8; + const su2double P_dim = 101325.0; + const su2double T_ref_dim = 288.15; // Reference Temperature + + // Non-dimensional Gas Constant set to 1.0 for simplicity (Cp output will be Cp/R) + const su2double R_nd = 1.0; + const su2double P_nd = 1.0; // Arbitrary for Cp check + const su2double T_ref_nd_for_class = 298.15 / T_ref_dim; // Std_Ref_Temp_ND + + std::stringstream ss; + SetupConfigStream(ss); + CConfig config(ss, SU2_COMPONENT::SU2_CFD, false); + ConfigureNASA(config); + + // Initialize with dimensional T_ref passed as 4th argument + CIncIdealGasNASA<9> nasa_model(R_nd, P_nd, T_ref_nd_for_class, T_ref_dim); + nasa_model.SetCpModel(&config, T_ref_nd_for_class); + + SECTION("Correct scaling for T = 300K") { + const su2double T_target_dim = 300.0; + const su2double T_target_nd = T_target_dim / T_ref_dim; + + nasa_model.SetTDState_T(T_target_nd); + + const su2double cp_nd = nasa_model.GetCp(); + // Expected Cp/R at 300K for N2 from coeffs: + // a1=a2=0, a3=3.298677, a4=1.4082404e-3, a5=-3.963222e-6, a6=5.641515e-9, a7=-2.444854e-12 + // Cp/R = 3.298677 + 1.4082404e-3*300 - 3.963222e-6*300^2 + 5.641515e-9*300^3 - 2.444854e-12*300^4 + // = 3.298677 + 0.422472 - 0.356690 + 0.152321 - 0.019803 + // = ~3.497 + + // Check if Cp output (which is Cp/R * 1.0) matches approx 3.5 + CHECK(cp_nd == Approx(3.497).epsilon(0.01)); + } +} diff --git a/UnitTests/meson.build b/UnitTests/meson.build index 3e64cfb9aa3..6b92257f0b0 100644 --- a/UnitTests/meson.build +++ b/UnitTests/meson.build @@ -15,6 +15,7 @@ su2_cfd_tests = files(['Common/geometry/primal_grid/CPrimalGrid_tests.cpp', 'Common/toolboxes/multilayer_perceptron/CLookUp_ANN_tests.cpp', 'SU2_CFD/numerics/CNumerics_tests.cpp', 'SU2_CFD/fluid/CFluidModel_tests.cpp', + 'SU2_CFD/fluid/CIncIdealGasNASA_tests.cpp', 'SU2_CFD/gradients.cpp', 'SU2_CFD/windowing.cpp']) diff --git a/config_template.cfg b/config_template.cfg index 99eb260e61e..19d4002fbc5 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -426,8 +426,13 @@ MOLECULAR_WEIGHT= 28.96, 16.043 % % Temperature polynomial coefficients (up to quartic) for specific heat Cp. % Format -> Cp(T) : b0 + b1*T + b2*T^2 + b3*T^3 + b4*T^4 +% For NASA 9-coefficient format (see CP_NASA option), use 9 coefficients: +% a1 (T^-2), a2 (T^-1), a3 (T^0), a4 (T), a5 (T^2), a6 (T^3), a7 (T^4), a8 (H const), a9 (S const) CP_POLYCOEFFS= (0.0, 0.0, 0.0, 0.0, 0.0) % +% Flag to use NASA polynomial format for Cp in incompressible ideal gas model. +CP_NASA= NO +% % --- Nonequilibrium fluid options % % Gas model - mixture