Combustion Equations

Detailed chemistry formulation

Disclaimer: The equations given here a dimensional!

Selection of conservative variables set and general equations

A set of primitive and conservative variables is defined for two-dimensional flows. Note that the extension to three-dimensional flows is straightforward. The set of conservative variables, namely \(\mathrm{G}_c\), is the non-reacting Navier-Stokes equations conservative variables set expanded by a total of \(N\) additional species mass fraction variables \(Y_k\), required to define the chemical composition of the gas mixture, where the subscript \(k\) denotes the associated species. The set of primitive variables, \(G_p\), is derived from \(G_c\) and the equation of state. In this way, the conservative variables and primitive variables are chosen as $$ \mathrm{G}_c=\left[\rho, \rho \mathbf{U}, \rho E, \rho Y_1, \ldots, \rho Y_N\right] $$

$$ \mathrm{G}_p=\left[\rho, \mathbf{U}, p, Y_1, \ldots, Y_N\right] $$

, where \(\mathbf{U}\) is the velocity vector, \(\rho\) defines the density of the fluid, \(E\) is the total energy and $p$ is the pressure. A total of \((N+4)\) variables are therefore necessary to completely describe the problem and to solve the corresponding equation system for two-dimensional flows.

The introduced mass fraction \(Y_k\) is defined as

\begin{align} Y_k=\frac{m_k}{m}, \end{align}

where \(m_k\) is the mass of species \(k\) in a given volume \(V\), and $m$ is the total mass in the same given volume. The mass fraction is related to the molar fraction \(X_k\), by the following equation

\begin{align} X_k=\frac{\bar{W}}{W_k} Y_k, \end{align}

where \(\bar{W}\) is the mean molar weight of the mixture and \(W_k\) is the species molar weight. The mean molar weight can be computed through

\begin{align} \frac{1}{\bar{W}}=\sum_{k=1}^N \frac{Y_k}{W_k} \quad \text { or } \quad \bar{W}=\sum_{k=1}^N X_k W_k \end{align}

For the energy equation the total non-chemical energy was chosen as the dependent variable. The total energy is related to the temperature field and velocity field through the expression

\begin{align} E=e_s+\frac{1}{2} \mathbf{U U}=\int_{T_0}^T c_v(T) d T-\frac{R T_0}{\bar{W}}+\frac{1}{2} \mathbf{U U}, \end{align}

where \(E\) stands for the total non-chemical energy, $e_s$ is the sensible energy, \(T_0\) is the reference temperature and \(c_v\) is the specific heat of the mixture at constant volume. The term \(R T_0 / \bar{W}\) is the sensible energy level at the reference temperature \(T_0\), which is unequal zero. The sensible energy is related to the other forms of energy through the expressions given in the following table:

Form	Energy	Enthalpy
Sensible	\(e_s=\int_{T_0}^T c_v d T-R T_0 / \bar{W}\)	\(h_s=\int_{T_0}^T c_p d T\)
Sensible \(+\) Chemical	\(e=e_s+\sum_{k=1}^N \Delta h_{f, k}^o Y_k\)	\(h=h_s+\sum_{k=1}^N \Delta h_{f, k}^o Y_k\)
Total Chemical	\(e_t=e+\frac{1}{2} \mathbf{U U}\)	\(h_t=h+\frac{1}{2} \mathbf{U U}\)
Total non Chemical	\(E=e_s+\frac{1}{2} \mathbf{U U}\)	\(H=h_s+\frac{1}{2} \mathbf{U U}\)

A relation can be derived from the first law of thermodynamics to relate both heat capacities \(c_p\) and \(c_v\) for an ideal gas mixture. Similarly, an expression can be derived for the pure species heat capacities \(c_{p, k}\) and \(c_{v, k}\), which are required to compute the pure species sensible enthalpies. Both relations can be written as

$$ c_p-c_v=\frac{R}{\bar{W}} \quad \text { and } \quad c_{p, k}-c_{v, k}=\frac{R}{W_k}, $$

where \(R\) is the universal gas constant and \(c_p\) is the specific heat of the mixture at constant pressure. To completely close the set of equations, the pressure is computed from the temperature and density fields through an equation of state (ideal gas law), expressed as

$$ p=\frac{\rho T R}{\bar{W}} \text {. } $$

Heat capacity computation

The heat capacity is a species thermodynamic property and depends both on the temperature and on the pressure. However, the pressure dependence is usually small compared to the temperature dependence and is usually ignored. No exact, analytical expressions exist to compute the heat capacities as a function of the temperature, however, approximate expressions were developed by the NASA in the 90s. The heat capacity is approximated by a fourth-order polynomial, with coefficients \(a_i, i \in[0,1,2,3,4]\), given as thermodynamic data to the chemistry solver. The polynomial has the form

$$ W_k \frac{c_{p, k}}{R}=a_{k, 0}+a_{k, 1} T+a_{k, 2} T^2+a_{k, 3} T^3+a_{k, 4} T^4 . $$

A similar expression follows for \(c_{v, k}\) $$ W_k \frac{c_{v, k}}{R}=\left(a_{k, 0}-1\right)+a_{k, 1} T+a_{k, 2} T^2+a_{k, 3} T^3+a_{k, 4} T^4 . $$

Due to limitations in the fitting procedure for the coefficients \(a_i\), two different sets of coefficients exist for the low temperature and high temperature region, \(a_{i, l o w}\) and \(a_{i, \text { high }}\), respectively. In addition, the mixture specific heat capacity is computed as a mass weighted summation of the individual species specific heat capacities as $$ c_p=\sum_{k=1}^N Y_k c_{p, k} . $$

The analytical expression for the total energy after integrating over the low temperature region reads $$ E=\sum_{k=1}^N\left(\frac{Y_k R}{W_k}\left(\sum_{i=0}^4 \frac{a_{k, i, l o w}}{i+1} T^{* i+1}-T^*\right)_{T_0}^T\right)-\frac{R T_0}{\bar{W}}+\frac{1}{2} \mathbf{U U} $$

Transport equations

Species conservation equations

The mass conservation equation for species \(k\) can be written as

\begin{align} \frac{\partial \rho Y_k}{\partial t}+\frac{\partial}{\partial x_i}\left(\rho\left(u_i+V_{k, i}\right) Y_k\right)=\dot{\omega}_k, \end{align}

where \(\dot{\omega}_k\) and \(V_k\) are the reaction rate and the diffusion velocity of species $k$ respectively. The term \(\rho Y_k V_k\) is called species mass flux \(j_k\). The species diffusive velocity \(V_k\) can be computed through different methods.

The following constraints apply for the species conservation equations

\begin{align} \sum_{k=1}^N j_k=0 \quad \text { and } \quad \sum_{k=1}^N Y_k=1 \quad \text { and } \quad \sum_{k=1}^N \dot{\omega}_k=0 . \end{align}

Mass conservation equation

The total mass conservation equation remains unchanged compared to non-reacting flows. This is expected since mass is neither created nor destroyed during the chemical reaction process. The mass conservation equation reads

\begin{align} \frac{\partial \rho}{\partial t}+\frac{\partial\left(\rho u_i\right)}{\partial x_i}=0 . \end{align}

Since the mass conservation equation can be algebraically derived from the species transport equations, it is not an independent equation. Consequently, the set of differential equations is overdetermined and only \((N-1)\) additional species conservation equations need to be solved together with the mass conservation equation.

Alternatively, the \(N\) species conservation equations can be solved. The system's overdetermination poses no additional challenge when dealing with exact expressions for the species diffusion velocity since the mass conservation is ensured. However, this is no longer the case when first-order approximations are used to approximate the species diffusion velocities, and some additional measures are necessary to ensure global mass conservation.

Momentum conservation equation

The momentum conservation equation remains unchanged compared to non-reacting flows:

\begin{align} \frac{\partial}{\partial t} \rho u_j+\frac{\partial}{\partial x_i} \rho u_i u_j=-\frac{\partial p}{\partial x_j}+\frac{\tau_{i j}}{\partial x_i}+\rho \sum_{k=1}^N Y_k f_{k, j} \end{align}

where $f_{k, j}$ are the body forces in direction $i$.

Energy conservation equation

There are multiple forms for the energy conservation equation. The total nonchemical energy \(E\) is chosen as the dependent variable for the formulation found in m-AIA. The energy conservation equation is therefore written as

\begin{align} \frac{\partial \rho E}{\partial t}+\frac{\partial \rho u_i E}{\partial x_i}=\dot{\omega}_T-\frac{\partial q_i}{\partial x_i}+\frac{\partial}{\partial x_j}\left(\sigma_{i j} u_i\right)+\dot{Q}+\rho \sum_{k=1}^N Y_k f_{k, j}\left(u_i+V_{k, i}\right). \end{align}

The heat release term \(\dot{\omega}_T\) is computed as

\begin{align} \dot{\omega}_T=-\sum_{k=1}^N \Delta h_{f, k}^o \dot{\omega}_k, \end{align}

where \(\Delta h_{f, k}^o\) is the standard formation enthalpy of species \(k\) and \(\dot{\omega}_k\) is the species reaction rate. The standard formation enthalpy is defined as the standard-state enthalpy \(h^o\) minus the standard-state enthalpy of the elements forming the molecular compound. The superscript indicates standard-state conditions, which are taken to be a pressure of 1 bar. The standard temperature can be chosen freely. However, a single reference temperature has to be defined to ensure the consistency of the results. The reference temperature is \(T_0=298 \mathrm{~K}\). From its definition, it follows that for a pure species molecule at the standard pressure and temperature \(\Delta h_f^o=0\). The energy flux \(q_i\) is defined as

\begin{align} q_i=-\lambda \frac{\partial T}{\partial x_i}+\rho \sum_{k=1}^N h_{s, k} Y_k V_{k, i}-\sum_{k=1}^N \frac{R T}{W_k X_k} D_k^T d_{k, i} \end{align}

where \(\lambda\) is the heat transfer coefficient, \(h_{s, k}\) is the species sensible enthalpy, \(D_k^T\) is the thermal diffusivity of species \(k\) and \(d_k\) is the species diffusion driving force. The species' sensible enthalpy is computed as

\begin{align} h_{s, k}=\int_{T_0}^T c_{p, k}(T) d T . \end{align}

Special attention is given to the additional transport terms that appear in the conservation equations. These transport terms are all diffusive in nature, as they involve the second derivative of a field variable. This, in turn, implies that the convective flux remains unchanged, and only the diffusive flux is affected by the additional terms introduced by the species transport.

Species mass flux

Before proceeding with the discussion of the different models for the computation of the transport coefficients, the species mass diffusion velocity \(V_k\) is reintroduced, defined as

\begin{align} \mathbf{V}_{\mathbf{k}}=\mathbf{V}-\tilde{\mathbf{V}}_{\mathbf{k}} \end{align}

where \(\mathbf{V}\) is the mass average velocity and \(\tilde{\mathbf{V}}_{\mathbf{k}}\) is the average velocity of species \(k\) relative to the fixed frame of reference. The mass average velocity is defined as

\begin{align} \mathbf{V}=\frac{1}{\rho} \sum_{k=1}^N \rho_k \tilde{\mathbf{V}}_{\mathbf{k}}=\sum_{k=1}^N Y_k \tilde{\mathbf{V}}_{\mathbf{k}} . \end{align}

Once the mass diffusion velocity \(V_k\) is computed, the species mass flux \(j_k\) can be computed as

\begin{align} j_k=\rho_k\left(\tilde{\mathbf{V}}_{\mathbf{k}}-\mathbf{V}\right)=\rho_k V_k=\rho Y_k V_k . \end{align}

Similar expressions can be derived for molar quantities. The molar flux of species \(k\) related to the molar averaged velocity, defined as \(V^*\), is given by

\begin{align} J_k^*=\left[X_k\right]\left(\tilde{\mathbf{V}}_{\mathbf{k}}-\mathbf{V}^*\right)=\left[X_k\right] \mathbf{V}_{\mathbf{k}}^* \end{align}

where \(\left[X_k\right]\) is the molar concentration of species \(k\) and \(\mathbf{V}_{\mathbf{k}}^*\) its molar diffusion velocity. The computation of the diffusion velocity \(V_k\) is a central aspect of each transport model, and depends on the computed diffusion coefficients.

Mixture-averaged transport model

The mixture-averaged model employs averaging rules on the pure species transport and binary diffusion coefficients to compute transport coefficients of species \(\mathrm{k}\) into the mixture. The model offers a practical compromise between accuracy and computational expense.

For the computation of the mixture-averaged viscosity, the semi-empirical expression developed by Wilke and modified by Bird et al. is used

\begin{align} \mu=\sum_{k=1}^N \frac{X_k \mu_k}{\sum_{k=1}^N X_j \Phi_{k j}}, \end{align}

where \(\Phi_{k j}\) is defined by

\begin{align} \Phi_{k j}=\frac{1}{\sqrt{8}}\left(1+\frac{W_k}{W_j}\right)^{-1 / 2}\left[1+\left(\frac{\mu_k}{\mu_j}\right)^{1 / 2}\left(\frac{W_j}{W_k}\right)^{1 / 4}\right]^2 . \end{align}

For the mixture-averaged thermal conductivity, a similar expression is derived

\begin{align} \lambda=\frac{1}{2}\left(\sum_{k=1}^N X_k \lambda_k+\frac{1}{\sum_{k=1}^N X_k / \lambda_k}\right) . \end{align}

For the formulation of the mixture-averaged diffusion coefficients, an equivalent diffusion coefficient for each species \(k\) into the mixture of gases is computed. This mixture averaged diffusion coefficient is denoted by the subscript \(m\). When computing the mixture-averaged diffusion coefficients, different expressions are available. An overview of the relations, as well as the exact expressions to compute the mixture averaged diffusion coefficients, is presented in the following table:

Flux	Diffusion Coefficient	Diffusion Velocity
\(\mathbf{J}_k^*=-c D_{k m}^* \nabla X_k\)	\(D_{k m}^*=\frac{1-X_k}{\sum_{j \neq k}^N X_j / \mathcal{D}_{k j}}\)	\(\mathbf{V}_{\mathbf{k}}^*=-\frac{c}{X_k} D_{k m}^* \nabla X_k\)
\(\mathbf{j}_k=-\rho D_{k m} \nabla Y_k\)	\(\frac{1}{D_{k m}}=\sum_{j \neq k}^N \frac{X_j}{\mathcal{D}_{k j}}+\frac{X_k}{1-Y_k} \sum_{j \neq k}^N \frac{Y_j}{\mathcal{D}_{k j}}\)	\(\mathbf{V}_k=-\frac{D_{k m}}{Y_k} \nabla Y_k\)
\(\mathbf{j}_k=-\rho \frac{Y_k}{X_k} D_{k m}^{\prime} \nabla X_k\)	\(D_{k m}^{\prime}=\frac{1-Y_k}{\sum_{j \neq k}^N X_j / \mathcal{D}_{k j}}\)	\(\mathbf{V}_k=-\frac{1}{X_k} D_{k m}^{\prime} \nabla X_k\)

The mixture-averaged expressions are first-order approximations, and unlike in the Stefan-Maxwell equations and in the multicomponent model, no constraint exists to ensure that the net species diffusion flux is equal to zero, e.g.

\begin{align} \sum_{k=1}^N \rho Y_k V_k=0 . \end{align}

Multicomponent transport model

The Chapman-Enskog theory provides the basis for the multicomponent model. This model is obtained after a rigorous approach to the derivation of the transport coefficients. The thermal diffusion coefficient, the thermal conductivity and the multicomponent diffusion coefficients are evaluated from the solution of a system of equations, which are defined by the matrix \(L\). The matrix \(L\) can be further decomposed into nine block submatrices. The equation system takes the form

\begin{align} \left(\begin{array}{ccc} L^{(00,00)} & L^{(00,10)} & 0 \\ L^{(10,00)} & L^{(10,10)} & L^{(10,01)} \\ 0 & L^{(01,10)} & L^{(01,01)} \end{array}\right)\left(\begin{array}{c} a_{00}^1 \\ a_{10}^1 \\ a_{01}^1 \end{array}\right)=\left(\begin{array}{c} 0 \\ \mathbf{X} \\ \mathbf{X} \end{array}\right) . \end{align}

where the vector \(\mathbf{X}\) is composed of the single mole fraction components \(X_k\). The multicomponent diffusion coefficients are calculated from the inverse of the submatrix \(L^{(00,00)}\) as

\begin{align} D_{j k}=X_j \frac{16 T}{25 p m_k} \bar{W}\left(P_{j k}-P j j\right), \end{align}

with

\begin{align} (P)=\left(L^{(00,00)}\right)^{-1} \end{align}

The thermal conductivity is found from the solution of the system of equations as

\begin{align} \begin{aligned} \lambda_{0, \text { trans }} & =-4 \sum_{k=1}^N X_k a_{k 10}^1 \\ \lambda_{0, \text { int }} & =-4 \sum_{k=1}^N X_k a_{k 01}^1 \\ \lambda_0 & =\lambda_{0, \text { trans }}+\lambda_{0, \text { int }} \end{aligned} \end{align}

and the thermal diffusion coefficient is computed as

\begin{align} D_k^T=\frac{8 m_k X_k}{5 R} a_{k 00}^1 . \end{align}

Similarly to the mixture-averaged model, the diffusion velocity can be computed through the equation

\begin{align} V_k=\underbrace{\frac{1}{X_k \bar{W}} \sum_{j \neq k}^N W_j D_{k j} d_j}_{\text {Ficksian diffusion }}-\underbrace{\frac{D_k^T}{\rho Y_k} \frac{1}{T} \nabla T}_{\text {Soret diffusion }}, \end{align}

where \(q_j\) is the diffusion driving force, which can be computed as

\begin{align} d_j=\underbrace{\nabla X_j}_{\text {pecies gradient diffusion }}+\overbrace{\left(X_j-Y_j\right) \frac{1}{p} \nabla p}^{\text {Pressure gradient diffusion }} . \end{align}

The necessary components of the \(L\) matrix are

\begin{align} \begin{aligned} & L_{j k}^{00,00}=\frac{16 T}{25 p} \sum_{l=1}^K \frac{X_l}{m_j \mathcal{D}_{j l}}\left\{m_k X_k\left(1-\delta_{j l}\right)-m_j X_j\left(\delta_{j k}-\delta_{k l}\right)\right\}, \\ \hspace 2mm & L_{j k}^{00,10}=\frac{8 T}{5 p} \sum_{l=1}^K X_k X_l\left(\delta_{j k}-\delta_{j l}\right) \frac{m_l\left(1.2 C_{k l}^*-1\right)}{\left(m_k+m_l\right) \mathcal{D}_{k l}}, \\ \hspace 2mm & L_{j k}^{10,00}=L_{k j}^{00,10}, \\ \hspace 2mm & L_{j k}^{01,00}=L_{k j}^{00,01}=0, \\ \hspace 2mm & L_{j k}^{10,10}=\frac{16 T}{25 p} \sum_{l=1}^K \frac{m_j}{m_k} \frac{X_j X_l}{\left(m_j+m_l\right)^2 \mathcal{D}_{j l}}\left\{\left(\delta_{k l}-\delta_{j k}\right)\left[\frac{15}{2} m_k^2+\frac{25}{4} m_l^2-3 m_l^2 B_{j l}^*\right]\right. \\ & \left.-4 m_k m_l A_{j l}^*\left(\delta_{k l}+\delta_{j k}\right)\left[1+\frac{5}{3 \pi}\left(\frac{C_{j, \text { rot }}}{R \xi_{j l}}+\frac{C_{l, \text { rot }}}{R \xi_{l j}}\right)\right]\right\} \text {, } \\ \hspace 2mm & L_{j j}^{10,10}=-\frac{16 m_j X_j^2}{R \mu_j}\left(1+\frac{10 C_{j, \mathrm{rot}}}{R \xi_{j j}}\right)-\frac{16 T}{25 p} \sum_{l \neq j}^K \frac{X_j X_l}{\left(m_j+m_l\right)^2 \mathcal{D}_{j l}} \\ \hspace 2mm & \times\left\{\frac{15}{2} m_j^2+\frac{25}{4} m_l^2-3 m_l^2 B_{j l}^*+4 m_j m_l A_{j l}^*\left[1+\frac{5}{3 \pi}\left(\frac{C_{j, \mathrm{rot}}}{R \xi_{j l}}+\frac{C_{l, \mathrm{rot}}}{R \xi_{l j}}\right)\right]\right\}, \\ \hspace 2mm & L_{j k}^{10,01}=\frac{32 T}{5 \pi p C_{k, \text { int }}} \sum_{l=1}^K \frac{m_k A_{k l}^*}{\left(m_k+m_l\right) \mathcal{D}_{k l}}\left(\delta_{j l}+\delta_{j k}\right) X_k X_l \frac{C_{k, \mathrm{rot}}}{R \xi_{k l}}, \\ \hspace 2mm & L_{j j}^{10,01}=\frac{16}{3 \pi} \frac{m_j X_j^2}{\mu_j C_{j, \text { int }}} \frac{C_{j, \text { rot }}}{R \xi_{j j}}+\frac{32 T R}{5 \pi p C_{j, \text { int }}} \sum_{l \neq j}^K \frac{m_j A_{j l}^*}{\left(m_j+m_l\right) \mathcal{D}_{j l}} X_j X_l \frac{C_{j, \text { rot }}}{R \xi_{j l}}, \\ & L_{j k}^{01,10}=L_{k j}^{10,01}, \\ \hspace 2mm & L_{j j}^{01,01}=-\frac{8 R^2}{\pi C_{j, \text { int }}^2} \frac{m_j X_j^2}{R \mu_j} \frac{C_{j, \text { rot }}}{R \xi_{j j}}-\frac{4 R T}{C_{j, \text { int }} p}\left\{\sum_{l=1}^K \frac{X_j X_l}{\left[\mathcal{D}_{j \mathrm{int}, l}\right]}+\sum_{l \neq j}^K \frac{12 X_j X_l}{5 \pi C_{j, \text { int }}} \frac{m_j}{m_l} \frac{A_{j l}^*}{\mathcal{D}_{j l}} \frac{C_{j, \mathrm{rot}}}{\xi_{j j}}\right\}, \\ \hspace 2mm & L_{j k}^{01,01}=0(j \neq k) . \\ \end{aligned} \end{align}

For a closer look at the computation and terms, please refer to [Kee]