13. Reacting Flows with Finite-Rate Chemistry#

Reacting-flow calculations couple fluid transport, molecular mixing, thermodynamics, and chemical kinetics. The principal numerical difficulty is the disparity among their characteristic time scales. Convective and diffusive transport normally evolve on flow time scales, whereas individual chemical reactions may be many orders of magnitude faster. A practical finite-rate formulation must therefore preserve the conservation properties of the multicomponent flow equations while integrating a potentially stiff chemical source term [OB01, PV05]. Finite-rate source evaluation and separate kinetics integration have both been used within pressure-based reacting-flow algorithms [RSE93, SCC+94, SCL+95].

The governing equations in Chapter 2 already contain the species transport and multicomponent energy terms. This chapter develops the additional structure needed for finite-rate chemistry:

  • the construction of species production rates from a reaction mechanism,

  • thermochemical coupling between composition and energy,

  • consistent convective and diffusive species fluxes, and

  • source and Strang splitting of transport and chemistry.

The emphasis is on the mathematical relationships that must remain consistent across these pieces. The more general time-integration setting for source and Strang splitting is developed in Higher-Order Temporal Integration and Splitting.

13.1. Species Conservation and Thermochemical Coupling#

For species \(k\), the conservation equation can be written

(13.1.1)#\[\frac{\partial(\rho Y_k)}{\partial t} + \frac{\partial}{\partial x_j} \left(\rho u_jY_k+J_{k,j}\right) = \dot{\omega}_k, \qquad k=1,\ldots,N_s ,\]

where \(\rho\) is the mixture density, \(Y_k\) is the species mass fraction, \(J_{k,j}\) is the diffusive mass flux relative to the mass-averaged mixture velocity, and \(\dot{\omega}_k\) is the volumetric species mass-production rate. The units of \(\dot{\omega}_k\) are mass per unit volume per unit time. It is a species source, not itself a heat-release source.

A gradient-diffusion specialization of the species flux is

(13.1.2)#\[J_{k,j} = -\Gamma_k\frac{\partial Y_k}{\partial x_j}, \qquad \Gamma_k = \frac{\mu}{Sc_k} + \frac{\mu_t}{Sc_{k,t}},\]

where \(Sc_k\) and \(Sc_{k,t}\) are molecular and turbulent Schmidt numbers. More detailed mixture-averaged or multicomponent transport models replace \(\Gamma_k\) with a corresponding diffusion closure; the conservation requirements below are unchanged.

Because the species constitute the mixture,

(13.1.3)#\[\sum_{k=1}^{N_s}Y_k=1, \qquad \sum_{k=1}^{N_s}J_{k,j}=0, \qquad \sum_{k=1}^{N_s}\dot{\omega}_k=0.\]

The first identity defines a complete composition. The second states that diffusion redistributes mass among species without creating a net mixture mass flux relative to the mass-averaged velocity. The third follows from mass conservation by the reaction mechanism. Summing Eq. (13.1.1) over all species then recovers the mixture continuity equation:

(13.1.4)#\[\frac{\partial\rho}{\partial t} + \frac{\partial(\rho u_j)}{\partial x_j} =0.\]

These identities must hold in the discrete equations as well as in the continuous equations. In particular, all species must use the same convective mixture mass flux, their reconstructed face mass fractions must sum to one, and their face diffusion fluxes must sum to zero.

13.1.1. Where Chemical Heat Release Enters#

Chemical reactions change composition. The resulting heat release is the thermochemical consequence of that composition change; it should not be confused with \(\dot{\omega}_k\). In a total-energy or total-enthalpy formulation whose species thermodynamic properties include formation energies, the chemical energy is already part of the transported thermodynamic state. No separate volumetric heat-release term is then added to the total-energy equation. The reaction source changes \(Y_k\), and the equation of state maps the new composition and conserved energy to a new temperature and pressure.

If an alternative formulation transports only sensible energy, the formation-energy contribution must instead appear explicitly. With \(h_{f,k}^{\circ}\) denoting the species standard formation enthalpy, the corresponding sensible-enthalpy source has the form

(13.1.1.1)#\[\dot q_{\mathrm{chem}} = -\sum_{k=1}^{N_s}h_{f,k}^{\circ}\dot{\omega}_k,\]

subject to the sign convention adopted for the energy equation. Adding this term to an energy variable that already includes formation energy would double count the chemical heat release.

Species diffusion also transports energy. The enthalpy-diffusion flux

(13.1.1.2)#\[q_{h,j}^{\mathrm{diff}} = \sum_{k=1}^{N_s}h_kJ_{k,j}\]

must be formed from the same species diffusion fluxes used in Eq. (13.1.1). Correcting the species equation to enforce a zero-sum diffusive mass flux while retaining an uncorrected flux in the energy equation would make the species and energy balances inconsistent. Here \(h_k(T,p)\) is the species enthalpy in the additive ideal-mixture relation \(h=\sum_kY_kh_k(T,p)\), defined in Eq. (16.1.4). The species properties can describe real fluids; the ideal-gas specialization gives \(h_k(T)\), as discussed in Species Enthalpy Diffusion. Combining the enthalpy-diffusion flux with the species balances gives the temperature contribution in Eq. (16.7.3.4). Its face-flux form is given in Eq. (16.7.7.5).

13.2. Finite-Rate Chemical Source Terms#

13.2.1. Reaction Mechanism and Rates of Progress#

Consider a mechanism containing \(N_r\) reversible reactions among \(N_s\) species. Reaction \(r\) is written

(13.2.1.1)#\[\sum_{k=1}^{N_s}\nu_{kr}'\,\mathcal{X}_k \rightleftharpoons \sum_{k=1}^{N_s}\nu_{kr}''\,\mathcal{X}_k,\]

where \(\mathcal{X}_k\) denotes species \(k\), and \(\nu_{kr}'\) and \(\nu_{kr}''\) are the reactant and product stoichiometric coefficients. The net stoichiometric coefficient is

(13.2.1.2)#\[\nu_{kr}=\nu_{kr}''-\nu_{kr}'.\]

The molar concentration of species \(k\) is

(13.2.1.3)#\[C_k=\frac{\rho Y_k}{W_k},\]

where \(W_k\) is its molecular weight. For an elementary reaction governed by mass action, the net molar rate of progress is

(13.2.1.4)#\[q_r = k_{f,r}\prod_{k=1}^{N_s}C_k^{\nu_{kr}'} - k_{b,r}\prod_{k=1}^{N_s}C_k^{\nu_{kr}''}.\]

Third-body efficiencies, falloff functions, and other pressure-dependent modifiers may multiply the appropriate forward and reverse contributions. They alter the effective rate of progress but not the stoichiometric conversion from reaction progress to species production.

The mass-production rate of species \(k\) is

(13.2.1.5)#\[\dot{\omega}_k = W_k\sum_{r=1}^{N_r}\nu_{kr}q_r.\]

A valid reaction mechanism conserves every chemical element. If \(a_{\ell k}\) is the number of atoms of element \(\ell\) in species \(k\), each reaction satisfies

(13.2.1.6)#\[\sum_{k=1}^{N_s}a_{\ell k}\nu_{kr}=0.\]

Together with the molecular weights, elemental conservation gives

(13.2.1.7)#\[\sum_{k=1}^{N_s}W_k\nu_{kr}=0 \quad\Longrightarrow\quad \sum_{k=1}^{N_s}\dot{\omega}_k=0.\]

Thus chemistry can redistribute mass among species without appearing as a source in the mixture continuity equation. Mechanism construction and the evaluation of elementary, third-body, and falloff reactions are described in standard combustion texts [KCG03, PV05].

13.2.2. Forward and Reverse Rate Constants#

A common forward rate constant is the modified Arrhenius expression

(13.2.2.1)#\[k_{f,r} = A_rT^{\beta_r} \exp\left(-\frac{E_{a,r}}{R_uT}\right),\]

where \(A_r\) is the pre-exponential factor, \(\beta_r\) is the temperature exponent, \(E_{a,r}\) is the activation energy, and \(R_u\) is the universal gas constant. The units of \(A_r\) depend on the overall reaction order.

For a thermodynamically reversible reaction, the reverse rate should be consistent with the equilibrium constant:

(13.2.2.2)#\[k_{b,r}=\frac{k_{f,r}}{K_{c,r}}.\]

Let

(13.2.2.3)#\[\Delta G_r^\circ = \sum_{k=1}^{N_s}\nu_{kr}g_k^\circ, \qquad \Delta\nu_r = \sum_{k=1}^{N_s}\nu_{kr},\]

where \(g_k^\circ\) is the standard molar Gibbs free energy. The dimensionless pressure-based equilibrium constant is

(13.2.2.4)#\[K_{p,r} = \exp\left(-\frac{\Delta G_r^\circ}{R_uT}\right).\]

When the concentration products in Eq. (13.2.1.4) use dimensional molar concentrations, the corresponding concentration equilibrium constant is

(13.2.2.5)#\[K_{c,r} = K_{p,r} \left(\frac{p^\circ}{R_uT}\right)^{\Delta\nu_r},\]

where \(p^\circ\) is the standard-state pressure. If nondimensional activities are used instead, the standard-state factors are absorbed into their definition. A kinetics implementation must use one convention consistently; mixing a dimensionless equilibrium constant with dimensional concentration products gives an incorrect reverse rate.

13.2.3. Cell-Local Chemistry Problem#

During a reaction-only substep, face transport is absent and each cell solves

(13.2.3.1)#\[\frac{d(\rho Y_k)}{dt} = \dot{\omega}_k(\rho,\boldsymbol{Y},T,p), \qquad k=1,\ldots,N_s.\]

The mixture density is fixed during this cell-local reaction solve because the reaction mechanism conserves mass. Temperature and pressure must still be supplied to the rate expressions. Two useful thermodynamic closures are:

  1. Frozen thermodynamic state. The incoming temperature and pressure are held fixed throughout the chemistry integration.

  2. EOS-refreshed thermodynamic state. The equation of state is reevaluated as the composition changes, using the fixed density and the selected conserved energy or enthalpy for the chemistry substep.

The second approach preserves thermodynamic coupling within the ODE solve and is suitable when an ideal-gas specialization is not appropriate. In either case, the chemistry integrator must return a composition compatible with the fixed mixture mass.

Detailed mechanisms produce stiff ODE systems because their reaction time scales can differ greatly. Backward-differentiation-formula methods with adaptive step size and order are therefore commonly used. Stream uses the CVODE stiff integrator for this cell-local problem [HBG+05]. CVODE controls the internal chemistry steps independently of the outer flow time step; this improves the accuracy and robustness of the reaction solve but does not by itself remove the splitting error between chemistry and transport.

13.3. Consistent Species Fluxes#

The continuous constraints in Eq. (13.1.3) do not automatically survive independent finite-volume discretizations of each species. Small inconsistencies can accumulate because species are limited and diffused separately. The convective reconstruction and diffusive face fluxes must therefore be constructed so that the discrete sum of the species equations remains the discrete continuity equation.

13.3.1. Convective Face Reconstruction#

Let \(P\) denote a cell and \(f\) one of its faces. A second-order reconstruction has the form

(13.3.1.1)#\[Y_{k,f}^{\,*} = Y_{k,P}+\delta Y_{k,f}, \qquad \delta Y_{k,f} = \phi_{k,f}\nabla Y_{k,P}\cdot \left(\boldsymbol{x}_f-\boldsymbol{x}_P\right),\]

where \(\phi_{k,f}\) is a limiter. Even if \(\sum_kY_{k,P}=1\), independently limited increments need not satisfy \(\sum_k\delta Y_{k,f}=0\). Consequently, \(\sum_kY_{k,f}^{\,*}\) may differ from one, and summing the species convective fluxes no longer recovers the mixture mass flux.

A conservative correction can be constructed without reversing any reconstruction increment. Define the positive and negative increment magnitudes

(13.3.1.2)#\[P_f = \sum_{k=1}^{N_s}\max(\delta Y_{k,f},0), \qquad N_f = \sum_{k=1}^{N_s}\max(-\delta Y_{k,f},0).\]

The larger signed group is reduced to match the smaller one. Equivalently,

(13.3.1.3)#\[\widehat{\delta Y}_{k,f} = \alpha_f^+\max(\delta Y_{k,f},0) - \alpha_f^-\max(-\delta Y_{k,f},0),\]

with

(13.3.1.4)#\[\begin{split}(\alpha_f^+,\alpha_f^-) = \begin{cases} \left(N_f/P_f,\,1\right), & P_f>N_f,\\ \left(1,\,P_f/N_f\right), & N_f>P_f,\\ \left(1,\,1\right), & P_f=N_f. \end{cases}\end{split}\]

If one signed sum vanishes while the other does not, the nonzero group is scaled to zero. The corrected face value is

(13.3.1.5)#\[Y_{k,f} = Y_{k,P}+\widehat{\delta Y}_{k,f}.\]

It then follows that

(13.3.1.6)#\[\sum_{k=1}^{N_s}\widehat{\delta Y}_{k,f}=0, \qquad \sum_{k=1}^{N_s}Y_{k,f}=1.\]

The rescaling factors lie between zero and one, so the operation only reduces the magnitude of an already limited increment. It does not create a new extremum beyond that candidate reconstruction. This construction is closely related to mass-fraction-consistent methods developed for compressible multicomponent flow [Lar91].

With one shared face mass flux \(\dot m_f\), the species convective fluxes are \(\dot m_fY_{k,f}\) and therefore satisfy

(13.3.1.7)#\[\sum_{k=1}^{N_s}\dot m_fY_{k,f} = \dot m_f.\]

This is the discrete convective identity required when the species equations are summed.

13.3.2. Diffusive Face Fluxes#

Let \(F_{k,f}^{\,*}\) be the provisional diffusive mass flow rate of species \(k\) through face \(f\), positive along the stored face normal. The cell-outward contribution is \(\sigma_{Pf}F_{k,f}^{\,*}\), with the incidence sign defined in Loci Face Orientation and Cell Incidence. Species-dependent transport coefficients and independent gradient approximations can make \(\sum_kF_{k,f}^{\,*}\ne0\), even though the continuous mixture-relative diffusion flux must sum to zero.

The same signed-group rescaling can be applied to the face flow rates. Define

(13.3.2.1)#\[P_f^J = \sum_{k=1}^{N_s}\max(F_{k,f}^{\,*},0), \qquad N_f^J = \sum_{k=1}^{N_s}\max(-F_{k,f}^{\,*},0).\]

Reducing the larger group to the magnitude of the smaller gives corrected fluxes

(13.3.2.2)#\[F_{k,f} = \alpha_{J,f}^+\max(F_{k,f}^{\,*},0) - \alpha_{J,f}^-\max(-F_{k,f}^{\,*},0),\]

where \(\alpha_{J,f}^+\) and \(\alpha_{J,f}^-\) are defined by the same cases as Eq. (13.3.1.4), with \(P_f\) and \(N_f\) replaced by \(P_f^J\) and \(N_f^J\). Hence

(13.3.2.3)#\[\sum_{k=1}^{N_s}F_{k,f}=0.\]

The corrected flux \(F_{k,f}\) must be used consistently in the species residual and in the species-enthalpy diffusion term. This preserves mixture mass while retaining the energetic effect of counter-diffusing species.

13.3.3. Conservation, Boundedness, and Renormalization#

Flux consistency and positivity are related but distinct requirements:

  • reconstruction consistency enforces \(\sum_kY_{k,f}=1\) at a face;

  • diffusion consistency enforces \(\sum_kF_{k,f}=0\) at a face;

  • a spatial limiter can keep candidate face values within desired bounds;

  • a positivity-preserving time update is needed to guarantee \(Y_{k,P}^{n+1}\ge0\) in a cell; and

  • post-update clipping or renormalization can repair an inadmissible composition, but is a separate operation with its own conservation and thermodynamic consequences.

Thus the zero-sum correction should not be described as a general proof of positivity. It supplies the algebraic identities needed for species-to-mixture consistency. Small negative mass fractions remain nonphysical even when they arise from roundoff or solver tolerance. If clipping or renormalization is used, the resulting composition should be returned to the simplex

(13.3.3.1)#\[Y_k\ge0, \qquad \sum_{k=1}^{N_s}Y_k=1,\]

and the thermodynamic state should be made consistent with that corrected composition.

13.4. Chemistry–Transport Operator Splitting#

The semidiscrete reacting-flow equations can be written schematically as

(13.4.1)#\[\frac{dQ}{dt} = F_T(Q)+F_R(Q),\]

where \(F_T\) contains convection, diffusion, pressure–velocity coupling, and the remaining flow physics, while \(F_R\) contains the finite-rate chemical source. Let \(\mathcal{T}_{\Delta t}\) and \(\mathcal{R}_{\Delta t}\) denote numerical propagators for the transport and reaction subproblems. Operator splitting allows the stiff, cell-local reaction system to use a specialized ODE integrator without placing every chemical time scale directly inside the coupled flow solve [VS98, Yan71].

13.4.1. Sequential Source Splitting#

The simplest split advances one operator and then the other. A reaction-first composition is

(13.4.1.1)#\[Q^{(a)} = \mathcal{R}_{\Delta t}(Q^n), \qquad Q^{n+1} = \mathcal{T}_{\Delta t}(Q^{(a)}).\]

The transport-first ordering is

(13.4.1.2)#\[Q^{(a)} = \mathcal{T}_{\Delta t}(Q^n), \qquad Q^{n+1} = \mathcal{R}_{\Delta t}(Q^{(a)}).\]

These two orderings are not equivalent because transport and reaction generally do not commute. Both are first-order accurate splitting approximations.

A source-coupled implementation of the reaction-first form can integrate the cell chemistry over \(\Delta t\) and convert the result to a time-averaged species source:

(13.4.1.3)#\[\dot{\omega}_{k,\mathrm{eff}} = \frac{(\rho Y_k)^{(a)}-(\rho Y_k)^n}{\Delta t}.\]

The flow transport equations are then solved with \(\dot{\omega}_{k,\mathrm{eff}}\) held fixed during that outer update. The dual transport-first construction can similarly form an effective transport source from the transport-only change and include it during the reaction integration. These effective sources reproduce their respective substep changes, but they do not make the one-sided composition second order.

Sequential splitting is attractive because it is simple and because a stiff chemistry integrator can subcycle within each flow step. Its splitting error can nevertheless become important where transport and chemistry interact strongly, such as flames whose structure is controlled by simultaneous diffusion and reaction [BV00].

13.4.2. Symmetric Strang Splitting#

Strang splitting restores formal second-order accuracy by composing the operators symmetrically [Str68]. A transport-centered step is

(13.4.2.1)#\[Q^{(a)} = \mathcal{T}_{\Delta t/2}(Q^n), \qquad Q^{(b)} = \mathcal{R}_{\Delta t}(Q^{(a)}), \qquad Q^{n+1} = \mathcal{T}_{\Delta t/2}(Q^{(b)}).\]

A reaction-centered step is equally valid:

(13.4.2.2)#\[Q^{(a)} = \mathcal{R}_{\Delta t/2}(Q^n), \qquad Q^{(b)} = \mathcal{T}_{\Delta t}(Q^{(a)}), \qquad Q^{n+1} = \mathcal{R}_{\Delta t/2}(Q^{(b)}).\]

The choice determines which operator occupies the full middle substep; it does not change the symmetry argument. For a constant time step, adjacent half steps of the same operator at successive step boundaries may be combined into one full substep. With a variable time step, their durations must instead be combined using the two neighboring step sizes.

Second-order behavior requires more than writing a symmetric sequence. Each subproblem must be integrated to sufficient accuracy, the exchange state between subproblems must be thermodynamically consistent, and boundary conditions must be applied at the appropriate substep times. Stiff reaction systems can also exhibit order reduction when these assumptions are not met. The numerical convergence of a split reacting-flow calculation should therefore be checked rather than inferred only from the formal composition [BV00, VHB98].

13.4.3. Coupling Interpretation in Stream#

The two reacting-flow paths in Stream follow the preceding distinction:

  • the source-splitting path converts a cell-local chemistry update into an average species source for the transport equations; and

  • the Strang-splitting path constructs a reacted intermediate composition and uses that state in the species transient update so that chemistry and transport can be composed symmetrically.

In both paths, the chemistry problem advances species densities at fixed mixture density, while the chosen thermodynamic closure supplies the temperature and pressure needed by the rates. The transport problem retains the pressure–velocity coupling and face fluxes. This separation is numerical, not physical: the complete time step must still represent one coupled reacting-flow state.

The detailed relationship between the reacted intermediate state and the one-step transport integrator is given in Higher-Order Temporal Integration and Splitting. The essential consistency requirements are:

  1. the chemical mechanism must conserve mass and elements;

  2. the reaction solve and EOS must describe the same thermodynamic state;

  3. species convection and diffusion must sum to the mixture mass equation;

  4. the species diffusion used in energy must match that used in species transport; and

  5. the splitting sequence and substep time intervals must match the intended temporal order.

13.5. References#

[BV00] (1,2)

J.G. Blom and J.G. Verwer. A comparison of integration methods for atmospheric transport-chemistry problems. Journal of Computational and Applied Mathematics, 126:381–396, 2000.

[HBG+05]

Alan C. Hindmarsh, Peter N. Brown, Keith E. Grant, Steven L. Lee, Radu Serban, Dan E. Shumaker, and Carol S. Woodward. SUNDIALS: suite of nonlinear and differential/algebraic equation solvers. ACM Transactions on Mathematical Software, 31(3):363–396, 2005. doi:10.1145/1089014.1089020.

[KCG03]

Robert J. Kee, Michael E. Coltrin, and Peter Glarborg. Chemically Reacting Flow: Theory and Practice. Wiley-Interscience, 2003.

[Lar91]

B. Larrouturou. How to preserve the mass fractions positivity when computing compressible multi-component flows. Journal of Computational Physics, 95(1):59–84, 1991. doi:10.1016/0021-9991(91)90253-H.

[OB01]

E.S. Oran and J.P. Boris. Numerical Simulation of Reactive Flow. Cambridge University Press, 2nd edition, 2001.

[PV05] (1,2)

T. Poinsot and D. Veynante. Theoretical and Numerical Combustion. R.T. Edwards, 2nd edition, 2005.

[RSE93]

C.M. Rhie, S.T. Stowers, and H.B. Ebrahimi. Numerical analysis of reacting flows using finite rate chemistry models. Journal of Propulsion and Power, 9:119–126, 1993.

[SCC+94]

H.M. Shang, Y.S. Chen, Z.J. Chen, C.P. Chen, and T.S. Wang. Numerical studies of finite rate reacting flows with a pressure-based method. In 30th AIAA/ASME/SAE/ASEE Joint Propulsion Conference. June 1994.

[SCL+95]

H.M. Shang, Y.S. Chen, P. Liaw, C.P. Chen, and T.S. Wang. Investigation of chemical kinetics integration algorithms for reacting flow. In 33rd Aerospace Sciences Meeting and Exhibit. January 1995.

[Str68]

G. Strang. On the construction and comparison of difference schemes. SIAM Journal on Numerical Analysis, 5:506–517, 1968.

[VHB98]

J.G. Verwer, W.H. Hundsdorfer, and J.G. Blom. Numerical time integration for air pollution models. Technical Report, CWI, 1998.

[VS98]

J.G. Verwer and B. Sportisse. A note on operator splitting analysis in the stiff linear case. Technical Report, CWI, 1998.

[Yan71]

N.N. Yanenko. The Method of Fractional Steps. Springer-Verlag, 1971.