BinarySolutionTabulatedThermo Gibbs/entropy change abruptly at composition limits

[speth]

Problem description

The entropy and Gibbs free energy values for the BinarySolutionTabulatedThermo phase change abruptly to artificial values as the composition approaches the limit of being a pure phase of the "tabulated" species.

Steps to reproduce

# setup
>>> import cantera as ct
>>> import numpy as np
>>> p = ct.Solution('lithium_ion_battery.yaml', 'anode')
>>> X = np.vstack((np.ones(10), np.logspace(-1, -28, 10))).T
>>> s = ct.SolutionArray(p, 10)
>>> s.TPX = 350, 10*101325, X
# A set of compositions approaching pure Li[anode] (tabulated species)
>>> s.X
array([[9.09090909e-01, 9.09090909e-02],
       [9.99900010e-01, 9.99900010e-05],
       [9.99999900e-01, 9.99999900e-08],
       [1.00000000e+00, 1.00000000e-10],
       [1.00000000e+00, 1.00000000e-13],
       [1.00000000e+00, 1.00000000e-16],
       [1.00000000e+00, 1.00000000e-19],
       [1.00000000e+00, 1.00000000e-22],
       [1.00000000e+00, 1.00000000e-25],
       [1.00000000e+00, 1.00000000e-28]])

>>> s.entropy_mole
array([3.75293924e+004, 9.59293272e+004, 1.53364703e+005, 2.10798976e+005,
       2.68230665e+005, 3.24798155e+005, 1.00000000e+300, 1.00000000e+300,
       1.00000000e+300, 1.00000000e+300])

>>> s.gibbs_mole
array([-1.16767237e+007, -3.19710049e+007, -5.20732261e+007,
       -7.21752216e+007, -9.22763127e+007, -1.12074934e+008,
       -3.50000000e+302, -3.50000000e+302, -3.50000000e+302,
       -3.50000000e+302])

# A set of compositions approaching pure secondary species (V[anode])
>>> s.TPX = 350, 10*101325, np.fliplr(X)
>>> s.entropy_mole
array([ 3.70455046e+03,  3.88833900e+00,  3.88868591e-03,  3.88868381e-06,
        3.88267188e-09,  1.96453279e-12, -3.33177996e-14, -3.90612269e-17,
       -4.48046542e-20, -5.05480815e-23])
>>> s.gibbs_mole
array([-2.20297912e+06, -7.76719802e+03, -7.76795943e+00, -7.76795933e-03,
       -7.76585516e-06, -7.09450648e-09,  5.25430986e-12,  7.26450941e-15,
        9.27470897e-18,  1.12849085e-20])

While I think it's perhaps acceptable for the partial molar entropies / chemical potentials to start taking on arbitrary values as their corresponding mole fractions go to zero, the molar entropy and Gibbs free energy should still tend toward a finite value at the limit of a single species.

Another way this is visible is that the standard chemical potentials for the "tabulated" species start to follow an unexpected trend below the last tabulated point (here, at X = 0.00575):

System information

• Cantera version: Cantera 2.6.0 or main at 8fbb4bc

Additional context

This was discovered as part of #1299, which implements Cantera/enhancements#114.

[speth]

Thank you, @dschmider-HSOG for the further investigation of these issues in #1302. I think PR #1442 cleans up most of these issues. The "abrupt" changes at the composition limits were a result of some floating point precision issues. From this starting point, it seems that the main remaining issue is with the calculation of the entropy for the mixture and each species, which has some surprising behavior. The following example, starting from the fixes in #1442, demonstrates this, plotting the mixture and partial molar enthalpies (left) and entropies (right):

import numpy as np
import cantera as ct
import matplotlib.pyplot as plt

p = ct.Solution('lithium_ion_battery.yaml', 'anode')
inpdata = p.input_data
Xt = np.array(inpdata['tabulated-thermo']['mole-fractions'])
St = np.array(inpdata['tabulated-thermo']['entropy'])
Ht = np.array(inpdata['tabulated-thermo']['enthalpy'])

N = 500
X = np.vstack((np.ones(N) - np.linspace(1-1e-5, 1e-5, N),
               np.linspace(1-1e-5, 1e-5, N))).T
s = ct.SolutionArray(p, N)
s.TPX = 350, 10*ct.one_atm, X

f, ax = plt.subplots(1, 2, figsize=(7,4))
R = ct.gas_constant
ax[0].plot(X[:,0], s.partial_molar_enthalpies[:,1] / (R * s.T), label=r'$h_{\rm V[anode]}$')
ax[0].plot(X[:,0], s.partial_molar_enthalpies[:,0] / (R * s.T), label=r'$h_{\rm Li[anode]}$')
ax[0].plot(X[:,0], s.enthalpy_mole / (R * s.T), label=r'$h_{\rm mix}$')
ax[0].plot(Xt, Ht / (ct.gas_constant * s.T[0]), '.', label='tabulated data')
ax[0].set(xlabel=r'$X_{\rm Li[anode]}$', ylabel='$h/RT$')
ax[0].legend()

ax[1].plot(X[:,0], s.partial_molar_entropies[:,1] / R, label=r'$s_{\rm V[anode]}$')
ax[1].plot(X[:,0], s.partial_molar_entropies[:,0] / R, label=r'$s_{\rm Li[anode]}$')
ax[1].plot(X[:,0], s.entropy_mole / R, label=r'$s_{\rm mix}$')
ax[1].plot(Xt, St / R, '.', label='tabulated data')
ax[1].plot(X[:,0], s.standard_entropies_R[:,0], label=r'$s^o_{\rm Li[anode]}$ ')

ax[1].legend()
ax[1].set(xlabel=r'$X_{\rm Li[anode]}$', ylabel='$s/R$')

gives this output:

The following features of this raise suspicions to me:

• The "tabulated entropy" corresponds to neither the partial molar entropy of the tabulated species (Li[anode]) nor the standard-state partial molar entropy ($s^o_{Li}$) of that species. If it isn't supposed to be one of these, what is it?
• $s^o_{Li} \to -\infty$ as $X_{Li} \to 0$, and $s^o_{Li} \to +\infty$ as $X_{Li} \to 1$. I think the standard state entropy should be well-behaved at both limits (true for IdealSolidSolnPhase, at least).
• $s_{Li}$ should go to $+\infty$ as $X_{Li} \to 0$, but seems not to.
• $s_{V}$ approaches to $+\infty$ as $X_{Li} \to 0$, but seems not to.
• The mixture entropy should be tend toward a finite value as $X_{Li} \to 1$.

I think if we can at least answer the question of what $s_{tab}$ is meant to correspond to, it should be possible to find a consistent way of representing these other quantities.