Enthalpy at T=Tref(=298.15 K) is not returning zero for species with zero heat of formation

[AlirezaMazaheri]

Here are some data showing some combinations of thermodynamic database with different state model could result in nonzero enthalpies at T=Tref(=298.15 K) for species that have zero heat of formations, such as N2 and O2. The only acceptable outcome is the first row (maybe second row as well) in the table, and the most surprising outcome is the last row.

Combinations	abs(h_s - 0)
RRHO; model: ChemNonEq1T or equil	1e-12
RRHO; model: ChemNonEqTTv	1e-10
NASA-7; model: *	1e-3
NASA-9; model: *	1e-3
NASA-9-New; model: *	N2=1e-3; O2=144

[jbscoggi]

Hey @AlirezaMazaheri,

Thanks for the pointing this out. Since our previous discussion on this earlier (offline), I have been thinking about this a lot and can report some preliminary results. First, I think we can split this problem into 2 separate issues:

1. The error in the NASA polynomial models for all but O_2 in NASA-9-New, and
2. The error in O_2 with NASA-9-New which seems exceptional.

At first glance, I thought that the errors in (1) are due to round-off / IEEE floating-point representation errors. After some estimates I'm now pretty confident that the real culprit is simply the precision of the polynomial coefficients listed in the NASA-7/9 formats. Please see a complete analysis below.

Note that I ruled out other bugs etc by computing the values of h(T) using a the numpy polyval function in python in a standalone file and confirmed that I get exactly what mppequil provides for O2 NASA-7. This provides more confidence in the results below. Specifically, the h_N2(298.15) and h_O2(298.15) for NASA-7/NASA-9 from mppequil give

Database	N_2	O_2
NASA-7	0.000160077	0.000510749
NASA-9	0.000218271	-0.000400241

You can check this for yourself using

mppequil -P 101325 -T 298.15 -s 10 --thermo-db NASA-7 --species-list "N2 O2"

As for (2), I think this deserves some more investigation. I want to reiterate that the NASA-9-New database is not fully validated and should be used with caution and I think you have probably found a good example of why this is the case. I will look more into that issue when I have time, but for now, I suggest to not use that database.

Analysis of possible errors in NASA polynomials

Let's first recall what the NASA polynomials represent. Specifically, I'm going to just look at the NASA-7 format, since the same analysis applies to NASA-9. For NASA-7 we have

c_p = R \sum_{i=1}^5 a_i T^{i-1},

and

h = \int c_p dT + R b.

Note we can therefore write h as the polynomial

h(T) = \sum_i=0^5 c_i T^i,

where c_0 = R b and c_i = R (a_i / i) for i > 0.

Since h(T) is a polynomial, we can apply some standard round-off error analyses:

Error due to round-off of addition/multiplication operations

This is a classical result for evaluation of polynomials using Horner's rule. The round-off error is bounded by

|h(T) - h*(T)| <= 10 \epsilon \sum_{i=0}^5 |c_i| T^i,

where the 10 comes from 2 times the order of the polynomial (5) and \epsilon is the machine precision. Using O_2 from NASA-7 we see that this bound is 10^{-9} J/kg, which does not explain the error we see.

Error due to floating-point representation of T

Another source of error could be the fact that 298.15 is not exactly represented in IEEE double precision floating-point numbers. We can estimate this impact by perturbing T by some error and linearizing the enthalpy around T, such that

h(T+\delta) \approx h(T) + \delta dh/dT(T).

Therefore, we can write a simple bound on this error (assuming exact arithmetic) as

|h(T+\delta) - h(T)| = \delta |dh/dT(T)| = \delta |c_p(T)|.

Again, using O_2 from NASA-7, this error is determined to be roughly 10^{-11}, and so is not important either (but could explain the RRHO results...).

Error due to precision of the coefficients

Finally, I realized that the coefficients in the database have a precision of 9 digits, meaning that the relative error in the coefficients compared to the true coefficients computed by the authors of the database is roughly 10^{-9}. For simplicity, let's compare the error in the exact enthalpy polynomial (with exact coefficients) with the one that is read from the database. We will assume that c*_i are the true coefficients, and c_i are the ones provided, and that

| c_i - c*_i | / |c_i| <= \delta = 10^{-9} for all i.

Then we have

|h(T) - h*(T)| = | \sum_{i=0}^5 (c_i - c*i) T^i | <= \sum{i=0}^5 | c_i - c*i | T^i <= \delta \sum{i=0}^5 |c_i| T^i.

Again, using the values for O2 from NASA-7, this error bound turns out to be about 0.0006 J/kg at 298.15 K, which can completely explain the errors we see.

[AlirezaMazaheri]

Hi @jbscoggi,

Thank you for the analysis.
If I manually (outside of the code) calculate h(T=298.15) using NASA-9 data set, I get these:


which are better than what we get out of the code (~1e-3). Could you please verify these?

[jbscoggi]

Hi @AlirezaMazaheri, could you confirm the units of your calculations? The results I was reporting were in J/kg. If I multiply your hO2 by Ru/Mw_O2*298.15, I get something close to my values.

[AlirezaMazaheri]

@jbscoggi, Ah units :-) Here are values in J/kg:


Which are close to what you've got. We need more precise constants ;)

[jbscoggi]

@AlirezaMazaheri, OK great. Note that the signs are still different for some reason compared to what I got with mppequil, but I implemented the same in a python code and matched exactly the mppequil output for O2 NASA-7, so I am inclined to agree with it. For reference, here was my code

import numpy as np

def enthalpy_polynomial(data):
    c = np.zeros(6)
    c[0] = data[5]
    c[1:] = [a/(i+1) for i,a in enumerate(data[:5])]
    return np.flip(c)

if __name__ == '__main__':

    Ru = 8.31451070

    # O2 NASA-7
    data = np.array([
        3.78245636E+00,
       -2.99673416E-03,
        9.84730201E-06,
       -9.68129509E-09,
        3.24372837E-12,
       -1.06394356E+03,
        3.65767573E+00,
        0.0
    ])
    Mw = 0.032

    h_poly = enthalpy_polynomial(Ru/Mw*data)
    print(np.polyval(h_poly, 298.15))

Note the value for Ru was taken from the Gordon and McBride paper for the coefficients, it's not what is actually used in M++, but should be close enough for the order of magnitude that we care about.

[AlirezaMazaheri]

@jbscoggi, didn't I mention that my hand-calculation is based on NASA-9?

[jbscoggi]

@AlirezaMazaheri, yes but from my table above, I got a positive enthalpy for N2 and negative for O2 with NASA-9. I think we should get the same sign regardless of differences in Ru and Mw because these are just positive multipliers.

[AlirezaMazaheri]

@jbscoggi NASA-7 and NASA-9 dataset may not necessarily produce the same sign h at Tref. FWIW, here is the dataset along with the calculation. I'm pretty sure these dataset matches the published data from NASA report TP-2002-211556, "NASA Glenn Coefficients for Calculating Thermodynamic Properties of Individual Species," by B.J. McBride, M.J. Zehe, and S. Gordon. September 2002. Am I missing something?

import numpy as np

def h_overRT(T, a):
  return (-a[0]/T/T
          +a[1]*np.log(T)/T
          +a[2]
          +a[3]*T/2
          +a[4]*T**2/3
          +a[5]*T**3/4
          +a[6]*T**4/5
          +a[7]/T #integration constant
          )

if __name__ == '__main__':

  Avogadro  = 6.02214076e+23 # 1/mole
  Boltzmann = 1.380649e-23   #J/K
  Ru = Avogadro * Boltzmann  #J/K-mole

  # NASA-9
  mw_N2 = 0.0280134 #[kg/mol]
  data_N2 = np.array([-3.425563420e4,
                      4.847000970e2,
                      1.119010961,
                      4.293889240e-3,
                      -6.836300520e-7,
                      -2.023372700e-9,
                      1.039040018e-12,
                      -3.391454870e3
                     ])

  mw_O2 = 0.0319988 #[kg/mol]
  data_O2 = np.array([ 2.210371497e4,
                      -3.818461820e2,
                      6.082738360,
                      -8.530914410e-3,
                      1.384646189e-5,
                      -9.625793620e-9,
                      2.519705809e-12,
                      7.108460860e2
                     ])

  Tref = 298.15 #K
  print("h_N2[J/kg]:",  h_overRT(Tref, data_N2) * Ru/mw_N2*Tref,
        "\nh_O2[J/kg]:",h_overRT(Tref, data_O2) * Ru/mw_O2*Tref)

Which produces these results

h_N2[J/kg]: -0.0004571815747510116 
h_O2[J/kg]: 0.00019108540237411262

[jbscoggi]

Hi @AlirezaMazaheri, I think what's different (compared to me) is that you are dividing all the coefficients by T and then multiplying by T again to get the dimensional units. Of course, this will change the relative errors between coefficients and so you could get a different sign then what I got. Anyways, I think we agree that the errors for h(298.15) are totally reasonable in the end. I will close this issue and fill a new one for O2 NASA-9-New by itself, unless you have any objections?

[AlirezaMazaheri]

@jbscoggi Sounds good :-)