Not able to reproduce nasa from Cp values

[dronekill]

Version of pMuTT

• Version 1.2.18

Describe the bug We are trying to calculate the Nasa polynomials from a set of Cp values but are not able to match the reported values.

To Reproduce Steps to reproduce the behavior: Codes and values:

import numpy as np
from pmutt import constants as c
T = np.array([300,400,500,600,800,1000,1500]) # K
Cp = np.array([24.5301147227533,27.4199330783939,29.6099426386233,31.3401051625239,33.9299713193117,35.8200286806883,38.7198852772466]) # cal/mol/K
CpoR = Cp/c.R('cal/mol/K')
T_ref = c.T0('K')
H_ref =-528.858 # kJ/mol
HoRT_ref = H_ref/c.R('kJ/mol/K')/T_ref
S_ref = 351.147 # J/mol/K
SoR_ref = S_ref/c.R('J/mol/K')
T_mid=871.2

from pmutt.empirical.nasa import Nasa

CH3SiCl3_Nasa = Nasa.from_data(name='CH3SiCl3', T=T, CpoR=CpoR,
                                  T_ref=T_ref, HoRT_ref=HoRT_ref,
                                  SoR_ref=SoR_ref,elements=None, T_mid=T_mid)

import json
from pmutt.io.json import pmuttEncoder
with open('CH3OH_Nasa.json', 'w') as f_ptr:
    json.dump(CH3OH_Nasa, f_ptr, cls=pmuttEncoder, indent=True)

Additional context the NASA polynomial are not matching my chemkin string from litreature.

MTS(1)                  H   3Si  1C   1Cl  3G   100.000  5000.000  871.20      1
 8.74632481E+00 1.54050649E-02-7.34893609E-06 1.36409589E-09-9.12369639E-14    2
-6.68166711E+04-1.18594086E+01 1.57446334E+00 5.85801130E-02-9.93279576E-05    3
 8.52490042E-08-2.80368462E-11-6.59558920E+04 1.95201278E+01                   4

[jonlym]

Thanks for your interest in our software.

When I ran the code above, I see the following warnings:

RuntimeWarning: Small set of CpoR data between T_mid and T_high. Fit may not be desirable.
  warn(warn_msg, RuntimeWarning)

RankWarning: Polyfit may be poorly conditioned
  (mse, a_low, a_high) = _get_CpoR_MSE(T=T, CpoR=CpoR, T_mid=T_m)

Essentially, you are trying to fit the 10 polynomial heat capacity coefficients using only 7 data points. To solve this, I would suggest giving more data (the more the better). If you cannot collect any more, I would suggest using a Shomate polynomial instead, which only has 5 polynomial coefficients.

---

Another problem is the literature NASA polynomial has different temperature ranges than your data.

	T_low (K)	T_high (K)
Literature	100	5000
Data	300	1500

Given these, it is unlikely your NASA polynomial will be equivalent.

---

Lastly, I would recommend comparing your thermodynamic quantities (i.e. Cp, H, S) instead of the polynomial. Since it is just a polynomial fit, there can be several combinations of coefficients that describe the data. For example, I would add the following code that plots the thermodynamic quantities on one plot:

# Create NASA polynomial from literature values
# MTS(1) H 3Si 1C 1Cl 3G 100.000 5000.000 871.20 1
# 8.74632481E+00 1.54050649E-02-7.34893609E-06 1.36409589E-09-9.12369639E-14 2
# -6.68166711E+04-1.18594086E+01 1.57446334E+00 5.85801130E-02-9.93279576E-05 3
# 8.52490042E-08-2.80368462E-11-6.59558920E+04 1.95201278E+01 4
MTS_nasa = Nasa(name='MTS(1)', T_low=100., T_mid=871.2, T_high=5000.,
                a_low=np.array([1.57446334E+00, 5.85801130E-02, -9.93279576E-05,
                                8.52490042E-08, -2.80368462E-11,
                                -6.59558920E+04, 1.95201278E+01]),
                a_high=np.array([8.74632481E+00, 1.54050649E-02, -7.34893609E-06,
                                 1.36409589E-09, -9.12369639E-14, 
                                 -6.68166711E+04, -1.18594086E+01]))

# Calculate thermo data
fit_Cp = CH3SiCl3_Nasa.get_Cp(T=T, units='cal/mol/K')
lit_Cp = MTS_nasa.get_Cp(T=T, units='cal/mol/K')

fit_H = CH3SiCl3_Nasa.get_H(T=T, units='kJ/mol')
lit_H = MTS_nasa.get_H(T=T, units='kJ/mol')

fit_S = CH3SiCl3_Nasa.get_S(T=T, units='J/mol/K')
lit_S = MTS_nasa.get_S(T=T, units='J/mol/K')

# Create plots
f, axes = plt.subplots(nrows=3, ncols=1, sharex=True)
axes[0].plot(T, Cp, 'k*')
axes[0].plot(T, fit_Cp, 'r-')
axes[0].plot(T, lit_Cp, 'b--')

axes[1].plot(T_ref, H_ref, 'k*')
axes[1].plot(T, fit_H, 'r-')
axes[1].plot(T, lit_H, 'b--')

axes[2].plot(T_ref, S_ref, 'k*')
axes[2].plot(T, fit_S, 'r-')
axes[2].plot(T, lit_S, 'b--')

axes[0].set_ylabel('Cp (cal/mol/K)')
axes[1].set_ylabel('H (kJ/mol/K)')
axes[2].set_ylabel('S (J/mol/K)')
axes[2].set_xlabel('T (K)')
axes[0].legend(['Data', 'Fit', 'Literature'])

plt.show()

[jonlym]

If you need more help with this, please feel free to reopen the issue. Otherwise, I will assume it's solved.