Skip to content

Commit eca925a

Browse files
guptapratykshhguptapratykshh
committed
Implement full NASA-9 polynomial format with NASA-7 backward compatibility
1 parent 3701aa2 commit eca925a

3 files changed

Lines changed: 58 additions & 24 deletions

File tree

‎SU2_CFD/include/fluid/CIncIdealGasNASA.hpp‎

Lines changed: 39 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -37,15 +37,19 @@
3737
* \brief Child class for defining an incompressible ideal gas model with NASA polynomials.
3838
* \author Pratyksh Gupta
3939
*
40-
* Implements NASA 7-coefficient polynomial format (NASA SP-273) for thermodynamic properties:
40+
* Implements NASA 9-coefficient polynomial format for thermodynamic properties with backward compatibility for NASA-7.
4141
* 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.
42-
* Cp/R = a1 + a2*T + a3*T^2 + a4*T^3 + a5*T^4
43-
* H/(R*T) = a1 + a2*T/2 + a3*T^2/3 + a4*T^3/4 + a5*T^4/5 + a6/T
44-
* S/R = a1*ln(T) + a2*T + a3*T^2/2 + a4*T^3/3 + a5*T^4/4 + a7
4542
*
46-
* Uses a single temperature range provided via CP_POLYCOEFFS (indices 0-6).
43+
* NASA-9 format (full):
44+
* Cp/R = a1*T^-2 + a2*T^-1 + a3 + a4*T + a5*T^2 + a6*T^3 + a7*T^4
45+
* 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
46+
* 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
47+
*
48+
* NASA-7 format (subset): Set a1=a2=0 to recover the traditional 7-coefficient format.
49+
*
50+
* Uses a single temperature range provided via CP_POLYCOEFFS (indices 0-8).
4751
*/
48-
template <int N_COEFFS = 7>
52+
template <int N_COEFFS = 9>
4953
class CIncIdealGasNASA final : public CFluidModel {
5054
public:
5155
/*!
@@ -68,9 +72,13 @@ class CIncIdealGasNASA final : public CFluidModel {
6872
void SetCpModel(const CConfig* config, su2double val_Temperature_Ref) override {
6973

7074
/*--- Read NASA coefficients from the standard polynomial coefficient array (CP_POLYCOEFFS).
71-
Indices 0-4: Cp coefficients (a1-a5)
72-
Index 5: Enthalpy constant (a6)
73-
Index 6: Entropy constant (a7) ---*/
75+
NASA-9 format uses indices 0-8:
76+
Indices 0-1: Inverse temperature terms (a1*T^-2, a2*T^-1)
77+
Indices 2-6: Polynomial terms (a3, a4*T, a5*T^2, a6*T^3, a7*T^4)
78+
Index 7: Enthalpy constant (a8)
79+
Index 8: Entropy constant (a9)
80+
81+
For NASA-7 compatibility: Set indices 0-1 to zero. ---*/
7482
for (int i = 0; i < N_COEFFS; ++i) {
7583
if (i < config->GetnPolyCoeffs()) {
7684
coeffs_[i] = config->GetCp_PolyCoeff(i);
@@ -97,12 +105,20 @@ class CIncIdealGasNASA final : public CFluidModel {
97105
const su2double a4 = coeffs_[3];
98106
const su2double a5 = coeffs_[4];
99107
const su2double a6 = coeffs_[5];
108+
const su2double a7 = coeffs_[6];
109+
const su2double a8 = coeffs_[7];
100110

101-
su2double Cp_over_R = a1 + a2*t + a3*t*t + a4*t*t*t + a5*t*t*t*t;
111+
const su2double t_inv = 1.0 / t;
112+
const su2double t_inv2 = t_inv * t_inv;
113+
114+
// NASA-9: Cp/R = a1*T^-2 + a2*T^-1 + a3 + a4*T + a5*T^2 + a6*T^3 + a7*T^4
115+
su2double Cp_over_R = a1*t_inv2 + a2*t_inv + a3 + a4*t + a5*t*t + a6*t*t*t + a7*t*t*t*t;
102116

103117
Cp = Cp_over_R * Gas_Constant;
104118

105-
su2double H_over_RT = a1 + a2*t/2.0 + a3*t*t/3.0 + a4*t*t*t/4.0 + a5*t*t*t*t/5.0 + a6/t;
119+
// 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
120+
su2double H_over_RT = -a1*t_inv2 + a2*std::log(t)*t_inv + a3 + a4*t/2.0 + a5*t*t/3.0 +
121+
a6*t*t*t/4.0 + a7*t*t*t*t/5.0 + a8*t_inv;
106122

107123
Enthalpy = H_over_RT * Gas_Constant * t;
108124
Cv = Cp / Gamma;
@@ -125,6 +141,8 @@ class CIncIdealGasNASA final : public CFluidModel {
125141
const su2double a4 = coeffs_[3];
126142
const su2double a5 = coeffs_[4];
127143
const su2double a6 = coeffs_[5];
144+
const su2double a7 = coeffs_[6];
145+
const su2double a8 = coeffs_[7];
128146

129147
su2double Cp_iter = 0.0;
130148
su2double delta_temp_iter = 1e10;
@@ -134,14 +152,19 @@ class CIncIdealGasNASA final : public CFluidModel {
134152

135153
while ((abs(delta_temp_iter) > toll) && (counter++ < counter_limit)) {
136154

137-
su2double Cp_over_R = a1 + a2*temp_iter + a3*temp_iter*temp_iter +
138-
a4*temp_iter*temp_iter*temp_iter + a5*temp_iter*temp_iter*temp_iter*temp_iter;
155+
const su2double t_inv = 1.0 / temp_iter;
156+
const su2double t_inv2 = t_inv * t_inv;
157+
158+
// NASA-9: Cp/R = a1*T^-2 + a2*T^-1 + a3 + a4*T + a5*T^2 + a6*T^3 + a7*T^4
159+
su2double Cp_over_R = a1*t_inv2 + a2*t_inv + a3 + a4*temp_iter + a5*temp_iter*temp_iter +
160+
a6*temp_iter*temp_iter*temp_iter + a7*temp_iter*temp_iter*temp_iter*temp_iter;
139161

140162
Cp_iter = Cp_over_R * Gas_Constant;
141163

142-
su2double H_over_RT = a1 + a2*temp_iter/2.0 + a3*temp_iter*temp_iter/3.0 +
143-
a4*temp_iter*temp_iter*temp_iter/4.0 +
144-
a5*temp_iter*temp_iter*temp_iter*temp_iter/5.0 + a6/temp_iter;
164+
// 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
165+
su2double H_over_RT = -a1*t_inv2 + a2*std::log(temp_iter)*t_inv + a3 + a4*temp_iter/2.0 +
166+
a5*temp_iter*temp_iter/3.0 + a6*temp_iter*temp_iter*temp_iter/4.0 +
167+
a7*temp_iter*temp_iter*temp_iter*temp_iter/5.0 + a8*t_inv;
145168

146169
su2double Enthalpy_iter = H_over_RT * Gas_Constant * temp_iter;
147170

‎SU2_CFD/src/solvers/CIncEulerSolver.cpp‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -309,8 +309,8 @@ void CIncEulerSolver::SetNondimensionalization(CConfig *config, unsigned short i
309309
config->SetGas_Constant(UNIVERSAL_GAS_CONSTANT/(config->GetMolecular_Weight()/1000.0));
310310
Pressure_Thermodynamic = Density_FreeStream*Temperature_FreeStream*config->GetGas_Constant();
311311
if (config->GetCp_NASA_Format()) {
312-
/*--- Use the NASA polynomial model (requires 7 coefficients) ---*/
313-
auxFluidModel = new CIncIdealGasNASA<7>(config->GetGas_Constant(), Pressure_Thermodynamic, STD_REF_TEMP);
312+
/*--- Use the NASA polynomial model (supports NASA-9 format with 9 coefficients, backward compatible with NASA-7) ---*/
313+
auxFluidModel = new CIncIdealGasNASA<9>(config->GetGas_Constant(), Pressure_Thermodynamic, STD_REF_TEMP);
314314
} else {
315315
/*--- Use the standard polynomial model ---*/
316316
auxFluidModel = new CIncIdealGasPolynomial<N_POLY_COEFFS>(config->GetGas_Constant(), Pressure_Thermodynamic, STD_REF_TEMP);
@@ -511,7 +511,7 @@ void CIncEulerSolver::SetNondimensionalization(CConfig *config, unsigned short i
511511

512512
case INC_IDEAL_GAS_POLY:
513513
if (config->GetCp_NASA_Format()) {
514-
fluidModel = new CIncIdealGasNASA<7>(Gas_ConstantND, Pressure_ThermodynamicND, STD_REF_TEMP / config->GetTemperature_Ref());
514+
fluidModel = new CIncIdealGasNASA<9>(Gas_ConstantND, Pressure_ThermodynamicND, STD_REF_TEMP / config->GetTemperature_Ref());
515515
} else {
516516
fluidModel = new CIncIdealGasPolynomial<N_POLY_COEFFS>(Gas_ConstantND, Pressure_ThermodynamicND, STD_REF_TEMP / config->GetTemperature_Ref());
517517
}

‎UnitTests/SU2_CFD/fluid/CIncIdealGasNASA_tests.cpp‎

Lines changed: 16 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -32,11 +32,22 @@
3232

3333
// Helper to config an already constructed CConfig object
3434
void ConfigureNASA(CConfig& config) {
35-
// Nitrogen (N2) Low Range Coefficients (200K - 1000K)
36-
su2double nasa_coeffs[7] = {3.298677, 1.4082404e-3, -3.963222e-6, 5.641515e-9, -2.444854e-12, -1020.8999, 3.950372};
35+
// Nitrogen (N2) Low Range Coefficients (200K - 1000K) in NASA-7 format
36+
// NASA-9 indexing: a1=0, a2=0 (no inverse T terms), a3-a7 (polynomial), a8-a9 (integration constants)
37+
su2double nasa_coeffs[9] = {
38+
0.0, // a1 (T^-2 term, not used in NASA-7)
39+
0.0, // a2 (T^-1 term, not used in NASA-7)
40+
3.298677, // a3 (constant term)
41+
1.4082404e-3, // a4 (T term)
42+
-3.963222e-6, // a5 (T^2 term)
43+
5.641515e-9, // a6 (T^3 term)
44+
-2.444854e-12, // a7 (T^4 term)
45+
-1020.8999, // a8 (H integration constant)
46+
3.950372 // a9 (S integration constant)
47+
};
3748

3849
config.SetCp_NASA_Format(true);
39-
for (int i = 0; i < 7; ++i) {
50+
for (int i = 0; i < 9; ++i) {
4051
config.SetCp_PolyCoeff(i, nasa_coeffs[i]);
4152
}
4253

@@ -64,7 +75,7 @@ TEST_CASE("NASA polynomial Cp calculation", "[CIncIdealGasNASA]") {
6475
CConfig config(ss, SU2_COMPONENT::SU2_CFD, false);
6576
ConfigureNASA(config);
6677

67-
CIncIdealGasNASA<7> nasa_model(R_N2, P, T_ref);
78+
CIncIdealGasNASA<9> nasa_model(R_N2, P, T_ref);
6879
nasa_model.SetCpModel(&config, T_ref);
6980

7081
SECTION("Valid temperature range (300K)") {
@@ -85,7 +96,7 @@ TEST_CASE("NASA polynomial enthalpy-temperature inversion", "[CIncIdealGasNASA]"
8596
CConfig config(ss, SU2_COMPONENT::SU2_CFD, false);
8697
ConfigureNASA(config);
8798

88-
CIncIdealGasNASA<7> nasa_model(R_N2, P, T_ref);
99+
CIncIdealGasNASA<9> nasa_model(R_N2, P, T_ref);
89100
nasa_model.SetCpModel(&config, T_ref);
90101

91102
SECTION("Inversion Consistency") {

0 commit comments

Comments
 (0)