search

Cantera / cantera

Chemical kinetics, thermodynamics, and transport tool suite
https://cantera.org
Other
581 stars 342 forks source link

"Piecewise Gibbs" thermo fits generate unrealistic properties #1641

Open speth opened 8 months ago

speth commented 8 months ago

Problem description

The "piecewise-Gibbs" species thermodynamic model, implemented by class Mu0Poly is described as implementing "an interpolation of the Gibbs free energy based on a piecewise constant heat capacity approximation".

If you fit a Mu0Poly object using a set of points for a species, using a different model (e.g. NASA polynomials) to evaluate the standard Gibbs free energy at the interpolation point:

you can see that it does indeed match the specified values of $g^\circ$ at the interpolation points (green circles). However, the interpolated values (orange) oscillate more than you might expect for a fit of a function that's practically a straight line. Things get really wild when you plot the heat capacity, enthalpy, and entropy values:

While the heat capacity is a constant in each interpolation region as specified by the model, the values are bonkers, oscillating between positive and negative values, for values that are about 2 orders of magnitude higher than they should be. This in turn results in the zig-zag behavior of the enthalpy and entropy.

Steps to reproduce

import cantera as ct
import matplotlib.pyplot as plt
import numpy as np
plt.rcParams['figure.constrained_layout.use'] = True

name = 'CH4'
Nfit = 9
Thigh = 800

TT = np.linspace(298.15, Thigh, Nfit)
g0 = []
gas.TP = 298.15, ct.one_atm
h0 = gas[name].standard_enthalpies_RT[0] * ct.gas_constant * T
coeffs = [len(TT), h0]
for T in TT:
    gas.TP = T, ct.one_atm
    g0.append(gas[name].standard_gibbs_RT[0] * ct.gas_constant * T)
    coeffs.extend([T, g0[-1]])

interp_thermo = ct.Mu0Poly(T_low=TT[0], T_high=TT[-1], P_ref=ct.one_atm, coeffs=coeffs)
S = ct.Species(name, gas.species(name).composition)
S.thermo = interp_thermo
gas1 = ct.Solution(thermo='ideal-gas', species=[gas.species(name)])
gas2 = ct.Solution(thermo='ideal-gas', species=[S])

states1 = ct.SolutionArray(gas1, 200)
states1.TP = np.linspace(273, 298.15 + (Thigh - 298.15) * 1.2, 200), ct.one_atm
states2 = ct.SolutionArray(gas2, 200)
states2.TP = np.linspace(273, 298.15 + (Thigh - 298.15) * 1.2, 200), ct.one_atm

fig, ax = plt.subplots(figsize=(5,3))
ax.plot(states1.T, states1.standard_gibbs_RT[:,0] * ct.gas_constant * states1.T, label='NASA polynomial');
ax.plot(states2.T, states2.standard_gibbs_RT[:,0] * ct.gas_constant * states2.T, label='Mu0Poly');
ax.plot(TT, g0, 'o')
ax.set(xlabel='T [K]', ylabel=r'standard gibbs free energy [J/kmol]')
ax.legend()

fig, ax = plt.subplots(1, 3, figsize=(8,3))
ax[0].plot(states1.T, states1.cp_mole);
ax[0].plot(states2.T, states2.cp_mole);
ax[0].set(title='cp [J/kmol/K]')

ax[1].plot(states1.T, states1.enthalpy_mole);
ax[1].plot(states2.T, states2.enthalpy_mole);
ax[1].set(title='enthalpy [J/kmol]')

ax[2].plot(states1.T, states1.entropy_mole);
ax[2].plot(states2.T, states2.entropy_mole);
ax[2].set(title='entropy [J/kmol/K]');

System information