Databases

Combustion Toolbox generates its own databases using an up-to-date version of NASA’s 9-coefficient polynomial fits [McBride, 2002] that incorporates the Third Millennium database [Burcat and Ruscic, 2005], including the available values from Active Thermochemical Tables. This fitting allows to evaluate the dimensionless thermodynamic functions of the species for the specific heat at constant pressure, enthalpy, and entropy as a function of temperature, namely

\begin{eqnarray} c_p^\circ/R &=& a_1T^{-2} + a_2T^{-1} + a_3 + a_4T + a_5T^2 + a_6T^3 + a_7T^4,\\ h^\circ/RT &=& -a_1T^{-2} + a_2T^{-1} \ln{T} + a_3 + a_4T/2 + a_5T^2/3 + a_6T^3/4 \\ &\phantom{{}={}}& + a_7T^4/5 + a_8/T,\\ s^\circ/R &=& -a_1T^{-2}/2 - a_2T^{-1} + a_3\ln{T} + a_4T + a_5T^2/2 + a_6T^3/3 \\ &\phantom{{}={}}& + a_7T^4/4 + a_9, \end{eqnarray}

where \(a_i\) from \(i=1, \dots, 7\) are the temperature coefficients and \(i =8, 9\) are the integration constants, respectively. Depending of the species the polynomials fit up to 20000 K [McBride, 2002]. These values are available in the source code and can be also obtained from NASA’s thermo build website.

To compute the dimensionless Gibbs energy, \(g_i^\circ (T) / RT\), from NASA’s polynomials we use the next expression

\begin{equation} g_i^\circ/RT = h^\circ/RT - s^\circ/R, \end{equation}

or equivalently

\begin{eqnarray} g_i^\circ/RT &=& -a_1T^{-2}/2 + a_2T^{-1} (1 + \ln{T}) + a_3(1 - \ln{T}) - a_4T/2 - a_5T^2/6 - a_6T^3/12 \\ &\phantom{{}={}}& - a_7T^4/20 + a_8/T - a_9. \end{eqnarray}

This data is collected from the thermo_CT.inp file and next formatted into a more accessible structure (\textit{DB_master}) with the built-in function generate_DB_master.m. Then, for faster data access, we generate a new database DB using griddedInterpolant objects that contain piecewise cubic Hermite interpolating polynomials (PCHIP) [Fritsch and Carlson, 1980]. This allows the evaluation of the thermodynamic functions at a given temperature with ease. The use of piecewise cubic Hermite interpolating polynomials increments the performance of Combustion Toolbox in approximate 200% as shown in Figure 1 obtaining the same results as seen in Figure 2.

Note

For temperatures outside the bounds, we avoid the higher order terms of the polynomials by linear extrapolation, similar to Stock et al. [2018], extending the range of validity of the thermodynamic data available. It should be emphasized that this extension is limited to a narrow temperature range and may not apply to temperatures significantly outside of this range.

To evaluate the thermodynamic functions, e.g., the Gibbs energy [kJ/mol] function, or the thermal enthalpy [kJ/mol] of \(\text{CO}_2\) at \(T = 2000 \text{ K}\) is as simple as using these callbacks

g0_CO2  = species_g0('CO2', 2000, DB)
DhT_CO2 = species_DhT('CO2', 2000, DB) 

Figure 1: Performance test, execution times for over $10^5$ calculations of the specific heat at constant pressure, enthalpy, Gibbs energy, and entropy, denoted as $c_p$, $h_0$, $g_0$, and $s_0$, respectively, using the NASA’s 9 coefficient polynomials (dark blue) and the piecewise cubic Hermite interpolating polynomials (teal). The test has been carried out with an Intel(R) Core(TM) i7-8700 CPU @ 3.20GHz. Note: lower is better.

Figure 2: Comparison of entropy [kJ/(mol-K)] as a function of temperature [K] obtained using the piecewise cubic Hermite interpolating polynomials (lines) and using the NASA’s 9 coefficient polynomials (symbols) for a set of species.

Another important parameter comes from the conservation of mass, which is the stoichiometric matrix \(A_0\), by generalizing this constraint condition we have

\begin{equation} \sum\limits_{j = 1}^{\text{NS}} a_{ij} n_j - b_i^\circ = 0, \end{equation}

or in matricial form

\begin{equation} \underbrace{\left(\begin{array}{c c c c} a_{11} & a_{21} & \cdots & a_{\text{NS}1} \\ a_{12} & a_{22} & \cdots & a_{\text{NS}2} \\ \vdots & \vdots & & \vdots \\ a_{1\text{NE}} & a_{2\text{NE}} & \cdots & a_{\text{NS}\text{NE}} \\ \end{array}\right)}_{\mathbf{A^T}} \underbrace{\left(\begin{matrix} \ce{n_1}\\ \ce{n_2}\\ \ce{\vdots}\\ \ce{n_{\text{NS}}} \end{matrix}\right)}_{\mathbf{N}} - \underbrace{\left(\begin{matrix} \ce{b_1}\\ \ce{b_2}\\ \ce{\vdots}\\ \ce{b_{\text{NE}}} \end{matrix}\right)}_{\mathbf{b^\circ}} = 0, \end{equation}

where \(a_{ij}\) are the stoichiometric coefficients of the species, which represent the number of atoms of element \(i\) per mole of species \(j\). The number of moles and the total number of atoms of the \(j\) species and \(i\) element reads \(n_j\) and \(b_i\), respectively. This is computed during the initialization of the variable self. Using one of the predefined list of species, e.g., soot formation, this can be initializate as

self = App('soot formation')

A simple test can be performed by considering the global exothermic reaction of hydrogen bromide with \(n_j\) moles

\begin{equation} \text{H}_2 + \text{Br}_2 \rightleftharpoons 2\text{HBr} \end{equation}

which have only three species involve. The system obtained is

\begin{equation} \left(\begin{array}{c c c} 0 & 2 & 1\\ 2 & 0 & 1\\ \end{array}\right) \left(\begin{matrix} n_{\ce{H_2}}\\ n_{\ce{Br_2}}\\ n_{\ce{HBr}}\\ \end{matrix}\right) - \left(\begin{matrix} b_{\ce{Br}}^\circ\\ b_{\ce{H}}^\circ\\ \end{matrix}\right) = 0. \end{equation}

A quick check using Combustion Toolbox:

self = App({'H2', 'Br2', 'HBr'});
print_stoichiometric_matrix(self, 'transpose')

Transpose stoichiometric matrix:

         H2    Br2    HBr
         __    ___    ___

   BR    0      2      1 
   H     2      0      1