Modeling Chemical Reactions in Cantera

Here, we describe how Cantera calculates chemical reaction rates for various reaction types.

Elementary Reactions

The basic reaction type is a homogeneous reaction with a pressure-independent rate coefficient and mass action kinetics. For example:

\begin{equation*} \mathrm{A + B \rightleftharpoons C + D} \end{equation*}

with a forward rate constant \(k_f\) defined as a modified Arrhenius function:

\begin{equation*} k_f = A T^b e^{-E_a / RT} \end{equation*}

where \(A\) is the pre-exponential factor, \(T\) is the temperature, \(b\) is the temperature exponent, \(E_a\) is the activation energy, and \(R\) is the gas constant. The forward reaction rate is then calculated as:

\begin{equation*} R_f = [\mathrm{A}] [\mathrm{B}] k_f \end{equation*}

An elementary reaction can be defined in the YAML format using the elementary reaction type.

In YAML, the reaction type entry can be omitted, as it represents the default. In case the type entry is omitted and a species occurs on both sides, Cantera infers that the reaction type is three-body.

Three-Body Reactions

A three-body reaction is a gas-phase reaction of the form:

\begin{equation*} \mathrm{A + B + M \rightleftharpoons AB + M} \end{equation*}

Here \(\mathrm{M}\) is an unspecified collision partner that carries away excess energy to stabilize the \(\mathrm{AB}\) molecule (forward direction) or supplies energy to break the \(\mathrm{AB}\) bond (reverse direction). In addition to the generic collision partner \(\mathrm{M}\), it is also possible to explicitly specify a colliding species. In this case, the reaction type is automatically inferred by Cantera.

Different species may be more or less effective in acting as the collision partner. A species that is much lighter than \(\mathrm{A}\) and \(\mathrm{B}\) may not be able to transfer much of its kinetic energy, and so would be inefficient as a collision partner. On the other hand, a species with a transition from its ground state that is nearly resonant with one in the \(\mathrm{AB^*}\) activated complex may be much more effective at exchanging energy than would otherwise be expected.

These effects can be accounted for by defining a collision efficiency \(\epsilon\) for each species, defined such that the forward reaction rate is

\begin{equation*} R_f = [\mathrm{A}][\mathrm{B}][\mathrm{M}]k_f(T) \end{equation*}

where

\begin{equation*} [\mathrm{M}] = \sum_{k} \epsilon_k C_k \end{equation*}

where \(C_k\) is the concentration of species \(k\). Since any constant collision efficiency can be absorbed into the rate coefficient \(k_f(T)\), the default collision efficiency is 1.0.

New in Cantera 3.0: The rate coefficient \(k_f(T)\) may be implemented using any rate parameterization supported by Cantera, not just the modified Arrhenius form.

Sometimes, accounting for a particular third body's collision efficiency may require an alternate set of rate parameters entirely. In this case, two reactions are written:

\begin{equation*} \mathrm{A + B + M \rightleftharpoons AB + M \quad (R1)} \end{equation*}
\begin{equation*} \mathrm{A + B + C \rightleftharpoons AB + C \quad (R2)} \end{equation*}

where the third-body efficiency for C in the first reaction is set to zero.

A three-body reaction may be defined in the YAML format using the three-body reaction type.

Changed in Cantera 3.0: The type field of the YAML entry may be omitted. Reactions containing the generic third body M (as in R1) are automatically identified as three-body reactions. Reactions are also identified as three-body reactions if all of the following conditions are met:

  • There is exactly one species appearing as both a reactant and product (as in reaction R2)

  • All reactants and products have integral stoichiometric coefficients

  • The sum of the stoichiometric coefficients for either the reactants or products is 3.

If the third body efficiency of C in R1 is not set to zero, these are considered duplicate reactions and must be marked as such.

Falloff Reactions

A falloff reaction is one that has a rate that is first-order in \([\mathrm{M}]\) at low pressure, like a three-body reaction, but becomes zero-order in \([\mathrm{M}]\) as \([\mathrm{M}]\) increases. Dissociation/association reactions of polyatomic molecules often exhibit this behavior.

The simplest expression for the rate coefficient for a falloff reaction is the Lindemann form 1:

\begin{equation*} k_f(T, [{\mathrm{M}}]) = \frac{k_0[{ \mathrm{M}}]}{1 + \frac{k_0{ [\mathrm{M}]}}{k_\infty}} \end{equation*}

In the low-pressure limit, this approaches \(k_0{[\mathrm{M}]}\), and in the high-pressure limit it approaches \(k_\infty\).

Defining the non-dimensional reduced pressure:

\begin{equation*} P_r = \frac{k_0 [\mathrm{M}]}{k_\infty} \end{equation*}

The rate constant may be written as

\begin{equation*} k_f(T, P_r) = k_\infty \left(\frac{P_r}{1 + P_r}\right) \end{equation*}

More accurate models for unimolecular processes lead to other, more complex, forms for the dependence on reduced pressure. These can be accounted for by multiplying the Lindemann expression by a function \(F(T, P_r)\):

\begin{equation*} k_f(T, P_r) = k_\infty \left(\frac{P_r}{1 + P_r}\right) F(T, P_r) \end{equation*}

This expression is used to compute the rate coefficient for falloff reactions. The function \(F(T, P_r)\) is the falloff function.

A falloff reaction may be defined in the YAML format using the falloff reaction type.

The Troe Falloff Function

A widely-used falloff function is the one proposed by Gilbert et al. 2:

\begin{equation*} \log_{10} F(T, P_r) = \frac{\log_{10} F_{cent}(T)}{1 + f_1^2} \end{equation*}
\begin{equation*} F_{cent}(T) = (1-A) \exp(-T/T_3) + A \exp (-T/T_1) + \exp(-T_2/T) \end{equation*}
\begin{equation*} f_1 = (\log_{10} P_r + C) / (N - 0.14 (\log_{10} P_r + C)) \end{equation*}
\begin{equation*} C = -0.4 - 0.67\; \log_{10} F_{cent} \end{equation*}
\begin{equation*} N = 0.75 - 1.27\; \log_{10} F_{cent} \end{equation*}

A Troe falloff function may be specified in the YAML format using the Troe field in the reaction entry. The first three parameters, \((A, T_3, T_1)\), are required. The fourth parameter, \(T_2\), is optional; if omitted, the last term of the falloff function is not used.

Tsang's Approximation to \(F_{cent}\)

Wing Tsang presented approximations for the value of \(F_{cent}\) for Troe falloff in databases of reactions, for example, Tsang and Herron 3. Tsang's approximations are linear in temperature:

\begin{equation*} F_{cent} = A + BT \end{equation*}

where \(A\) and \(B\) are constants. The remaining equations for \(C\), \(N\), \(f_1\), and \(F\) from Troe falloff are not affected:

\begin{equation*} \log_{10} F(T, P_r) = \frac{\log_{10} F_{cent}(T)}{1 + f_1^2} \end{equation*}
\begin{equation*} f_1 = (\log_{10} P_r + C) / (N - 0.14 (\log_{10} P_r + C)) \end{equation*}
\begin{equation*} C = -0.4 - 0.67\; \log_{10} F_{cent} \end{equation*}
\begin{equation*} N = 0.75 - 1.27\; \log_{10} F_{cent} \end{equation*}

A Tsang falloff function may be specified in the YAML format using the Tsang field in the reaction entry. (New in Cantera 2.6)

The SRI Falloff Function

This falloff function is based on the one originally due to Stewart et al. 4, which required three parameters \(a\), \(b\), and \(c\). Kee et al. 5 generalized this function slightly by adding two more parameters \(d\) and \(e\). (The original form corresponds to \(d = 1\) and \(e = 0\).) Cantera supports the extended 5-parameter form, given by:

\begin{equation*} F(T, P_r) = d \bigl[a \exp(-b/T) + \exp(-T/c)\bigr]^{1/(1+\log_{10}^2 P_r )} T^e \end{equation*}

In keeping with the nomenclature of Kee et al. 5, we will refer to this as the SRI falloff function.

An SRI falloff function may be specified in the YAML format using the SRI field in the entry.

Chemically-Activated Reactions

For these reactions, the rate falls off as the pressure increases, due to collisional stabilization of a reaction intermediate. Example:

\begin{equation*} \mathrm{Si + SiH_4 (+M) \leftrightarrow Si_2H_2 + H_2 (+M)} \end{equation*}

which competes with:

\begin{equation*} \mathrm{Si + SiH_4 (+M) \leftrightarrow Si_2H_4 (+M)} \end{equation*}

Like falloff reactions, chemically-activated reactions are described by blending between a low-pressure and a high-pressure rate expression. The difference is that the forward rate constant is written as being proportional to the low-pressure rate constant:

\begin{equation*} k_f(T, P_r) = k_0 \left(\frac{1}{1 + P_r}\right) F(T, P_r) \end{equation*}

and the optional blending function \(F\) may described by any of the parameterizations allowed for falloff reactions.

Chemically-activated reactions can be defined in the YAML format using the chemically-activated reaction type.

Pressure-Dependent Arrhenius Rate Expressions (P-Log)

This parameterization represents pressure-dependent reaction rates by logarithmically interpolating between Arrhenius rate expressions at various pressures. Given two rate expressions at two specific pressures:

\begin{equation*} P_1: k_1(T) = A_1 T^{b_1} e^{-E_1 / RT} \end{equation*}
\begin{equation*} P_2: k_2(T) = A_2 T^{b_2} e^{-E_2 / RT} \end{equation*}

The rate at an intermediate pressure \(P_1 < P < P_2\) is computed as

\begin{equation*} \log k(T,P) = \log k_1(T) + \bigl(\log k_2(T) - \log k_1(T)\bigr) \frac{\log P - \log P_1}{\log P_2 - \log P_1} \end{equation*}

Multiple rate expressions may be given at the same pressure, in which case the rate used in the interpolation formula is the sum of all the rates given at that pressure. For pressures outside the given range, the rate expression at the nearest pressure is used.

Negative A-factors can be used for any of the rate expressions at a given pressure. However, the sum of all of the rates at a given pressure must be positive, due to the logarithmic interpolation of the rate for intermediate pressures. When a P-log type reaction is initialized, Cantera does a validation check for a range of temperatures that the sum of the reaction rates at each pressure is positive. Unfortunately, if these checks fail, the only options are to remove the reaction or contact the author of the reaction/mechanism in question, because the reaction is mathematically unsound.

P-log reactions can be defined in the YAML format using the pressure-dependent-Arrhenius reaction type.

Chebyshev Reaction Rate Expressions

Chebyshev rate expressions represent a phenomenological rate coefficient \(k(T,P)\) in terms of a bivariate Chebyshev polynomial. The rate constant can be written as:

\begin{equation*} \log k(T,P) = \sum_{t=1}^{N_T} \sum_{p=1}^{N_P} \alpha_{tp} \phi_t(\tilde{T}) \phi_p(\tilde{P}) \end{equation*}

where \(\alpha_{tp}\) are the constants defining the rate, \(\phi_n(x)\) is the Chebyshev polynomial of the first kind of degree \(n\) evaluated at \(x\), and

\begin{equation*} \tilde{T} \equiv \frac{2T^{-1} - T_\mathrm{min}^{-1} - T_\mathrm{max}^{-1}} {T_\mathrm{max}^{-1} - T_\mathrm{min}^{-1}} \end{equation*}
\begin{equation*} \tilde{P} \equiv \frac{2 \log P - \log P_\mathrm{min} - \log P_\mathrm{max}} {\log P_\mathrm{max} - \log P_\mathrm{min}} \end{equation*}

are reduced temperatures and reduced pressures which map the ranges \((T_\mathrm{min}, T_\mathrm{max})\) and \((P_\mathrm{min}, P_\mathrm{max})\) to \((-1, 1)\).

A Chebyshev rate expression is specified in terms of the coefficient matrix \(\alpha\) and the temperature and pressure ranges.

Note that the Chebyshev polynomials are not defined outside the interval \((-1,1)\), and therefore extrapolation of rates outside the range of temperatures and pressure for which they are defined is strongly discouraged.

Chebyshev reactions can be defined in the YAML format using the Chebyshev reaction type.

Blowers-Masel Reactions

In some circumstances like thermodynamic sensitivity analysis, or modeling heterogeneous reactions from one catalyst surface to another, the enthalpy change of a reaction (\(\Delta H\)) can be modified. Due to the change in \(\Delta H\), the activation energy of the reaction must be adjusted accordingly to provide accurate simulation results. To adjust the activation energy due to changes in the reaction enthalpy, the Blowers-Masel rate expression is available. This approximation was proposed by Blowers and Masel 6 to automatically scale activation energy as the reaction enthalpy is changed. The activation energy estimation can be written as:

\begin{equation*} E_a = \begin{cases} 0 & \text{if } \Delta H \leq -4 E_a^0 \\ \Delta H & \text{if } \Delta H \geq 4 E_a^0 \\ \frac{\left( w + \frac{\Delta H }{2} \right) (V_P - 2 w + \Delta H) ^2} {V_P^2 - 4 w^2 + \Delta H^2} & \text{Otherwise} \end{cases} \end{equation*}

where

\begin{equation*} V_P = 2 w \frac{w + E_a^0}{w - E_a^0}, \end{equation*}

\(w\) is the average of the bond dissociation energy of the bond breaking and that being formed, \(E_a^0\) is the intrinsic activation energy, and \(\Delta H\) is the enthalpy change of the reaction. Note that the expression is insensitive to \(w\) as long as \(w \ge 2 E_a^0\), so we can use an arbitrarily high value of \(w = 1000\text{ kJ/mol}\).

After \(E_a\) is evaluated, the reaction rate can be calculated using the modified Arrhenius expression

\begin{equation*} k_f = A T^b e^{-E_a / RT}. \end{equation*}

Blowers Masel reactions can be defined in the YAML format using the Blowers-Masel reaction type. (New in Cantera 2.6)

Surface Reactions

Heterogeneous reactions on surfaces are represented by an extended Arrhenius- like rate expression, which combines the modified Arrhenius rate expression with further corrections dependent on the fractional surface coverages \(\theta_{k}\) of one or more surface species. The forward rate constant for a reaction of this type is:

\begin{equation*} k_f = A T^b \exp \left( - \frac{E_a}{RT} \right) \prod_k 10^{a_k \theta_k} \theta_k^{m_k} \exp \left( \frac{- E_k \theta_k}{RT} \right) \end{equation*}

where \(A\), \(b\), and \(E_a\) are the modified Arrhenius parameters and \(a_k\), \(m_k\), and \(E_k\) are the coverage dependencies from species \(k\).

In the YAML format, surface reactions are identified by the presence of surface species and support several additional options.

In YAML, the surface reaction type defaults to interface-Arrhenius, where the rate expression uses the Arrhenius parameterization (see YAML documentation). As an alternative, Cantera also supports the interface-Blowers-Masel surface reaction type, which uses the Blowers-Masel parameterization (see YAML documentation; New in Cantera 2.6).

Sticking Reactions

Sticking reactions represent a special case of surface reactions, where collisions between gas-phase molecules and surfaces result in the gas-phase molecule sticking to the surface. This process can be described as a reaction which is parameterized by a sticking coefficient:

\begin{equation*} \gamma = a T^b e^{-c/RT} \end{equation*}

where \(a\), \(b\), and \(c\) are constants specific to the reaction. The values of these constants must be specified so that the sticking coefficient \(\gamma\) is between 0 and 1 for all temperatures.

The sticking coefficient is related to the forward rate constant by the formula:

\begin{equation*} k_f = \frac{\gamma}{\Gamma_\mathrm{tot}^m} \sqrt{\frac{RT}{2 \pi W}} \end{equation*}

where \(\Gamma_\mathrm{tot}\) is the total molar site density, \(m\) is the sum of all the surface reactant stoichiometric coefficients, and \(W\) is the molecular weight of the gas phase species.

Sticking reactions can be defined in the YAML format by specifying the rate constant in the reaction's sticking-coefficient field.

The sticking reaction type defaults to sticking-Arrhenius, where the rate expression uses the Arrhenius parameterization (see YAML documentation). As an alternative, Cantera also supports the sticking-Blowers-Masel surface reaction type, which uses the Blowers-Masel parameterization (see YAML documentation; New in Cantera 2.6).

Two-Temperature-Plasma Reactions

The two-temperature-plasma reaction is commonly used for non-equilibrium plasmas. The reaction rate of a two-temperature-plasma reaction depends on both gas and electron temperature 8, and can be expressed as:

\begin{equation*} k_f = A T_e^b \exp \left( - \frac{E_{a,g}}{RT} \right) \exp \left(\frac{E_{a,e}(T_e - T)}{R T T_e}\right), \end{equation*}

where \(A\) is the pre-exponential factor, \(T\) is the temperature, \(T_e\) is the electron temperature, \(b\) is the electron temperature exponent, \(E_{a,g}\) is the activation energy for gas, \(E_{a,e}\) is the activation energy for electron and \(R\) is the gas constant. (New in Cantera 2.6)

Additional Options

Reaction Orders

Explicit reaction orders different from the stoichiometric coefficients are sometimes used for non-elementary reactions. For example, consider the global reaction:

\begin{equation*} \mathrm{C_8H_{18} + 12.5 O_2 \rightarrow 8 CO_2 + 9 H_2O} \end{equation*}

the forward rate constant might be given as 7:

\begin{equation*} k_f = 4.6 \times 10^{11} [\mathrm{C_8H_{18}}]^{0.25} [\mathrm{O_2}]^{1.5} \exp\left(\frac{30.0\,\mathrm{kcal/mol}}{RT}\right) \end{equation*}

Special care is required in this case since the units of the pre-exponential factor depend on the sum of the reaction orders, which may not be an integer.

Note that you can change reaction orders only for irreversible reactions.

Normally, reaction orders are required to be positive. However, in some cases negative reaction orders are found to be better fits for experimental data. In these cases, the default behavior may be overridden in the input file.

References

1

  1. Lindemann. Trans. Faraday Soc., 17:598, 1922.

2

R. G. Gilbert, K. Luther, and J. Troe. Ber. Bunsenges. Phys. Chem., 87:169, 1983.

3

  1. Tsang and J. Herron. Journal of Physical and Chemical Reference Data, 20:4, 1991.

4

P. H. Stewart, C. W. Larson, and D. Golden. Combustion and Flame, 75:25, 1989.

5(1,2)

R. J. Kee, F. M. Rupley, and J. A. Miller. Chemkin-II: A Fortran chemical kinetics package for the analysis of gas-phase chemical kinetics. Technical Report SAND89-8009, Sandia National Laboratories, 1989.

6

Blowers, P., & Masel, R. (2000). Engineering approximations for activation energies in hydrogen transfer reactions. AIChE Journal, 46(10), 2041-2052. https://doi.org/10.1002/aic.690461015

7

C. K. Westbrook and F. L. Dryer. Simplified reaction mechanisms for the oxidation of hydrocarbon fuels in flames. Combustion Science and Technology 27, pp. 31--43. 1981.

8

I. A. Kossyi, A. Y. Kostinsky, A. A. Matveyev. and V. P. Kinetic scheme of the non-equilibrium discharge in nitrogen-oxygen mixtures. mechanisms for the oxidation of hydrocarbon fuels in flames. Plasma Sources Science and Technology 1, no. 3, pp. 207. 1992. DOI: https://doi.org/10.1088/0963-0252/1/3/011