Finite Rate Chemistry in Aither

As of v0.9.0 Aither has a finite rate reacting flow capability. This post will detail how the chemistry source terms are calculated in Aither. The governing equations for multispecies reacting flow are given below.

In order to calculate the source term due to chemical reactions, some information about each species in the simulation is needed. In Aither the fluid database specifies the enthalpy of formation , molar mass , vibrational temperatures , and the linear component of internal energy . The fluid database also supplies reference properties for pressure , temperature , and entropy . It is possible to add to Aither’s fluid database by specifying these properties for a new species. This information can be found in a variety of sources, including the NIST Chemistry WebBook.

Chemistry Source Terms

Reacting flow simulations differ from nonreacting simulations in that the species mass equations have a source term, . This source term is calculated via the law of mass action. In the equation below represents the stoichiometric coefficient on the product side of the reaction for species in reaction . Similarly, represents the stoichiometric coefficient on the reactant side of the reaction for species in reaction .

The source term depends on the forward and backward reaction rates. These rates are functions of temparature only. In Aither the forward reaction rate is calculated using an Arrhenius curve fit. The backward reaction rate is calculated from thermodynamic properties using the equilibrium rate . The units of the forward rate are and the units of the backward rate are .

Calculation of Equilibrium Rate

The equilibrium rate of the reaction is calculated via the minimization of the Gibbs free energy [1]. For a thermally perfect gas using the vibrational equilibrium model in Aither, the Gibbs free energy for a given species is shown below. In the equations below represents the number of vibrational modes for a given species.

The remaining integral represents the vibrational contribution to the entropy term. It can be solved with integration by parts.

Substituting this expression back into Gibbs free energy equation yields the following.

Where is a constant to account for any difference between the given reference entropy of the species and the calculated entropy at the species’ reference temperature.

From this the equilibrium rate can be calculated as shown below.

Implicit Treatment

For the implicit solvers in Aither the chemistry jacobian and/or spectral radius of the jacobian are needed. The analytical chemistry jacobian is very complicated and somewhat expensive to calculate. For this reason Aither uses a numerical jacobian. The derivative of the chemistry source terms is needed with respect to and . This is done by perturbing the conserved variable of interest by a small amount , while keeping the other conserved variables constant, and recalucuating the chemistry source terms.

For the scalar matrix implicit solvers in Aither (LU-SGS and DPLUR) only the spectral radius of the chemistry jacobian is needed. This is approximated using the procedure in [2] where only the negative terms that would increase the diagonal dominance of the implicit matrix are treated implicitly.

References

[1] Luke, E. A. “A Rule-Based Specification System for Computational Fluid Dynamics”. Ph. D. Thesis. 1999.

[2] Savard, B. et al. “A Computationally-Efficient Semi-Implicit Iterative Method for the Time Integration of Reacting Flows with Stiff Chemistry”. Journal of Computational Physics. 2015.