Chemical Equilibrium With Applications To Turbulent Reactive Flows

Definition

Chemical Equilibrium

For the purposes of this discussion, lets begin with a definition of chemical equilibrium. Chemical equilibrium is the state in which the composition has no net change over time. The composition $(\phi_\alpha)$ is defined by $\alpha$ species in terms of mole fractions $(X_\alpha)$, mass fractions $(Y_\alpha)$, or even mole numbers $(N_\alpha)$, having the relationships given below:

(1)

\begin{align} \phi_\alpha=X_\alpha=\frac{Y_\alpha M_{mix}}{M_\alpha}=\frac{N_\alpha}{N},~~~~\alpha=1,2,...,N_s \end{align}

where $M_\alpha$ and $M_{mix}$ denote the molecular weight of species $\alpha$ and of the mixture respectively, and $N_s$ denotes the number of species in the system. I will note here, that in our calculations of chemical equilibrium, presented in the next section, we are assuming that the temperature and pressure are specified and of constant value.

Since this discussion is presented with turbulent reactive flows and combustion in mind, we can make the following qualitative assumptions about the systems that we are interested in:

• 1 phase only - gas phase - and this gas phase mixture is an ideal gas
• Approximate temperature range of 300 K < T < 3000 K

These assumptions are very appropriate for most turbulent reactive flows, and will not be discussed further.

Methods to Calculate Chemical Equilibrium

Equilibrium Constant Method

The first method I will present for the calculation of chemical equilibrium is known as the equilibrium constant method. It is the method presented in most relevant texts, and is based on the statistical mechanics approach. I will only present the important results below. A more detailed discussion was given by Mr. Bostandoost Nik and can be found here.

Given a general reaction:

(2)

\begin{align} aA+bB \rightarrow cC+dD \end{align}

we can define the pressure based equilibrium constant as follows:

(3)

\begin{align} K_p=\frac{(P_c/P^o)^c(P_d/P^o)^d}{(P_a/P^o)^a(P_b/P^o)^b} \end{align}

where $P_\alpha=X_\alpha P$ is the partial pressure, and $P^o$ is the standard state pressure. Thus, we can rewrite the above equation in terms of the mole fractions:

(4)

\begin{align} K_p=\frac{(X_c)^c(X_d)^d}{(X_a)^a(X_b)^b}(P/P^o)^{(c+d)-(a+b)} \end{align}

To illustrate how we can calculate chemical equilibrium from this above equation, it is best to look at an example problem.

---

Equilibrium Constant Method Example:

Lets begin with a simple reaction for the dissociation of carbon dioxide into carbon monoxide and oxygen.

(5)

\begin{align} CO_2 \Leftrightarrow CO+\frac{1}{2}O_2 \end{align}

Then we can define the equilibrium constant as follows based on Eq.(4):

(6)

\begin{align} K_p=\frac{(X_{CO})(X_{O_2})^{0.5}}{(X_{CO_2})}(P/P^o)^{0.5} \end{align}

However, at this point we have 3 unknown quantities $(X_\alpha)$ and only one equation. We can close this problem with the introduction of two constraint equations. The first is that the sum of the mole fractions must be equal to one:

(7)

\begin{align} \sum^{N_s}_{\alpha=1} X_\alpha=1 \end{align}

The second constraint equation is that the relative number of carbon atoms to oxygen atoms can never change. That is, before dissociation begins we have 1 carbon atom for every 2 oxygen atoms. This ratio will be constant at any point in the dissociation process. Mathematically we have:

(8)

\begin{align} \frac{1}{2}=\frac{X_{CO}+X_{CO_2}}{X_{CO}+2X_{CO_2}+2X_{O_2}} \end{align}

At this point we have a closed system of 3 equations with 3 unknowns (recall that temperature and pressure are constant). The above set of equations is most easily solved by an iterative method such as the Newton-Raphson method.

---

The equilibrium constant method is theoretically very sound. However, it does have a major drawback. To illustrate this point, lets look at a general equation for the stoichiometric combustion of a hydrocarbon fuel with air:

(9)

\begin{align} C_xH_y+(x+y/4)O_2+zN_2 \rightarrow xCO_2+(y/2)H_2O+zN_2 \end{align}

In reality, for combustion processes, as the temperature increases due to exothermic reaction, the products tend to dissociate. The dissociation of these products leads to a further sequence of reactions, and a realistic representation of a hydrocarbon - air reaction would include many steps and products. The reaction of natural gas ( approx. methane) and air has been studied extensively by researchers. The results of their study is a chemistry reaction mechanism for natural gas known as GRI-Mech (current version 3.0). This mechanism represents the global reaction of natural gas and air by 325 reactions and includes 53 species. It should be clear then, that to calculate chemical equilibrium for this reaction mechanism using the equilibrium constant method, would require the calculation of 325 equilibrium constants, and the solution of a set of 53, mostly nonlinear, equations. Thus, a solution would be computationally expensive.

It is in this regard that we explore other methods for the calculation of chemical equilibrium, namely, the theory of element potentials.

Element Potential Method

For systems at constant temperature and pressure the appropriate fundamental equation describing the system would be that of the Gibbs free energy. We can define the Gibbs free energy as follows:

(10)

\begin{align} G(T,p,\mathbf{N})=U(S,V,\mathbf{N})+pV-TS \end{align}

The differential of the energy $(\mathrm{d}U)$ is readily recognized as the fundamental energy equation. However, for our purposes here we seek an integral form of the fundamental energy equation. This is readily attained by recognizing that the energy is a first-order homogeneous function.

---

First-Order Homogeneous Functions:

A first-order homogeneous function is a function of extensive variables. An extensive property of a system is equivalent to sum of that property of all the subcomponents of a system. The entropy, volume and mole number meet this requirement. Therefore, we can mathematically represent a first-order homogeneous function as follows:

(11)

\begin{align} f(\lambda x_1,\lambda x_2,...,\lambda x_M)=\lambda f(x_1,x_2,...,x_M) \end{align}

The important result of this property is that we can represent the first-order homogeneous function in the following way:

(12)

\begin{align} f=\sum^M_{x=1}x_i\frac{\partial f}{\partial x_i} \end{align}

Applying this result to $U(S,V,\mathbf{N})$ yields:

(13)

\begin{align} U=\left(\frac{\partial U}{\partial S}\right)_{V,\mathbf{N}}S + \left(\frac{\partial U}{\partial V}\right)_{S,\mathbf{N}}V + \sum^{N_s}_{\alpha=1}\left(\frac{\partial U}{\partial N_\alpha}\right)_{S,V,N_{i\neq \alpha}}}N_\alpha \end{align}

Where the partial derivative terms are defined thermodynamic properties, respectively representing the temperature $(T)$, the negative pressure $(-p)$, and the chemical potential for species $\alpha~~(\mu_\alpha)$. Therefore, we can write the integral of the fundamental energy equation as follows:

(14)

\begin{align} U=TS - pV + \sum^{N_s}_{\alpha=1}\mu_\alpha N_\alpha \end{align}

---

Applying the above results to Eq.(10) yields the following:

(15)

\begin{align} G=TS - pV + \sum^{N_s}_{\alpha=1}\mu_\alpha N_\alpha+pV-TS=\sum^{N_s}_{\alpha=1}\mu_\alpha N_\alpha \end{align}

The chemical potential can equivalently be represented by the partial molar Gibbs energy. Thus,

(16)

\begin{align} G=\sum^{N_s}_{\alpha=1}\bar{g}_\alpha N_\alpha \end{align}

where for an ideal gas, we can define the partial molar Gibbs energy as such:

(17)

\begin{align} \bar{g}_\alpha=g_\alpha(T,p)+RT\ln X_\alpha \end{align}

where $g_\alpha(T,p)$ is the partial molar Gibbs energy of pure species $\alpha$ at the specified temperature and pressure. This quantity is treated as known since it can be obtained from the JANAF Thermochemical Tables.

The next step of the element potential method is to define atomic population constraints. These constraints are given by the following equations:

(18)

\begin{align} \sum^{N_s}_{\alpha=1}n_{i\alpha}N_\alpha=P_i,~~~~i=1,2,...,a \end{align}

where $n_{i\alpha}$ is the number of i elements in an $\alpha$ molecule, $P_i$ is the population (moles) of the ith element in the system, and $a$ is the number of elements in the system. It is important to note that $P_i$ can represent absolute or relative populations.

We now have enough information to calculate the equilibrium composition. The solution will be the distribution of $N_\alpha$ that minimizes Eq.(16) according to the constraints represented by Eq.(18). A minimization problem subject to constraints is best solved by the method of Lagrange Multipliers.

---

Lagrange Multipliers:

In general, to find the extremum of a function, $f(x_1,x_2,...,x_M)$, subject to the constraints, $g_i(x_1,x_2,...,x_M)=c_i$, where $c_i$ are constants we have:

(19)

\begin{eqnarray} \frac{\partial f}{\partial x_1} - \sum_i \lambda_i \frac{\partial g_i}{\partial x_1}=0 \\ \frac{\partial f}{\partial x_2} - \sum_i \lambda_i \frac{\partial g_i}{\partial x_2}=0 \\ \vdots \\ \frac{\partial f}{\partial x_M} - \sum_i \lambda_i \frac{\partial g_i}{\partial x_M}=0 \\ \end{eqnarray}

An alternative representation of the Lagrange multiplier method is as follows:

(20)

\begin{align} \mathrm{d}(f-\sum_i\lambda_i g_i)=\sum^M_{j=1}\left[\frac{\partial f}{\partial x_j} - \sum_i \lambda_i \frac{\partial g_i}{\partial x_j}\right]\mathrm{d}x_j = 0 \end{align}

where Eq. (19) shows us that each of the bracketed terms above is necessarily equal to 0.

---

Now lets apply the Lagrange multiplier method to the element potential method. For convenience we seek to minimize $G/RT$:

(21)

\begin{align} \frac{G}{RT}=\sum^{N_s}_{\alpha=1}\left[\tilde{g}_\alpha N_\alpha+N_\alpha\ln N_\alpha-N_\alpha\ln N\right] \end{align}

where $\tilde{g}_\alpha=g_\alpha(T,p)/RT$. Differentiating the above equation yields:

(22)

\begin{eqnarray} \mathrm{d}\left(\frac{G}{RT}\right) & = & \sum^{N_s}_{\alpha=1}\left[\tilde{g}_\alpha \mathrm{d}N_\alpha + \mathrm{d}N_\alpha\ln N_\alpha +N_\alpha\frac{\mathrm{d}N_\alpha}{N_\alpha}-\mathrm{d}N_\alpha\ln N\right] \\ ~ & = & \sum^{N_s}_{\alpha=1}\left[\tilde{g}_\alpha \ln X_\alpha\right]\mathrm{d}N_\alpha + \sum^{N_s}_{\alpha=1}N\mathrm{d}X_\alpha \\ ~ & = & \sum^{N_s}_{\alpha=1}\left[\tilde{g}_\alpha \ln X_\alpha\right]\mathrm{d}N_\alpha \end{eqnarray}

where the second summation in the second equality is necessarily 0 because $\sum^{N_s}_{\alpha=1}X_\alpha=1 \Rightarrow \sum^{N_s}_{\alpha=1}\mathrm{d}X_\alpha=0$. Now we must differentiate our constraint equations, Eq.(18):

(23)

\begin{align} \sum^{N_s}_{\alpha=1}n_{i\alpha}\mathrm{d}N_\alpha=0,~~~~i=1,2,...,a \end{align}

Now combining our results in the form of Eq.(20) we have:

(24)

\begin{align} \sum^{N_s}_{\alpha=1}\left[\tilde{g}_\alpha+\ln X_\alpha-\sum^{a}_{i=1}\lambda_i n_{i\alpha}\right]\mathrm{d}N_\alpha=0 \end{align}

where the terms in brackets are necessarily 0 as shown above. Therefore, we can solve for the mole fractions in terms of the Lagrange multipliers:

(25)

\begin{align} X_\alpha=\exp \left[\tilde{g}_\alpha - \sum^{a}_{i=1}\lambda_i n_{i\alpha}\right],~~~~\alpha=1,2,...,N_s \end{align}

and we can solve for the Lagrange multipliers based on the constraint equations, and the fact that we know the sum of the mole fractions must sum to unity:

(26)

\begin{align} \sum^{N_s}_{\alpha=1}n_{i\alpha}NX_\alpha=P_i,~~~~i=1,2,...,a \end{align}

(27)

\begin{align} \sum^{N_s}_{\alpha=1}X_\alpha=1 \end{align}

Equations 25, 26, and 27 represent the solution for the equilibrium composition based on the method of element potentials. Some important results of this solution are presented here:

• From Eq.(25) and Eq.(17) we can gain some physical understanding of what the Lagrange multipliers represent. Namely, it is shown that the Lagrange multipliers represent the contribution of individual elements to the partial molar Gibbs function for a particular species:

(28)

\begin{align} \frac{\bar{g}_\alpha}{RT}=\frac{g_\alpha}{RT}-\frac{g_\alpha}{RT}+\sum^{a}_{i=1}\lambda_i n_{i\alpha}=\sum^{a}_{i=1}\lambda_i n_{i\alpha} \end{align}

• The Lagrange multipliers are thus called "element potentials" (recalling that the partial molar Gibbs function is equivalent to the chemical potential). Thus, the name Element Potential Method.
• For GRI-Mech 3.0, the equilibrium constant method required the knowledge of all 325 reactions, the calculation of each equilibrium constant, and the solution of 53 (mostly nonlinear) equations for the equilibrium solution. The Element Potential method only requires the solution of the element potentials, which amounts to the number of elements of interest. In GRI-Mech 3.0 there are only 5 elements (H,O,C,N,Ar)!
• As alluded to above, the Element Potential method requires no knowledge of reactions or equilibrium constants. However, the species of interest must be specified. Although, it is not necessary to offer any quantitative information about them.

Evidently the Element Potential method is a powerful and efficient method for the solution of chemical equilibrium. This method has been implemented in the interactive program known as STANJAN. Details of this program and its solution algorithm can be found here. STANJAN was such a successful program that it has been implemented into the CHEMKIN software, which is software for calculations relating to complex chemical kinetics.

---

Sample Calculation using STANJAN:

Lets revisit our stoichiometric reaction for hydrocarbon with air from earlier, Eq.(9). We will specify our fuel as methane gas $\mathrm{CH}_4$:

(29)

\begin{align} CH_4+2O_2+7.52N_2\rightarrow ? \end{align}

The question mark appears on the reactant side of the equation because we know from earlier that the reactants will dissociate, and we can't be sure about exactly what we will get. Now we must tell STANJAN what species to calculate equilibrium for. We will use the following: $\mathrm{H_2,~H,~O_2,~OH,~H_2O,~HO_2,~H_2O_2,~CH_3,~CH_4,~CO,~CO_2,~CH_2O,~C_2H_2,~C_2H_4,~C_2H_6,~N_2}$. The specification of these species will be adequate for the purposes of this example. Finally, we must specify our initial conditions:

T	298 K
P	1 atm
$N_{CH_4}$	1
$N_{O_2}$	2
$N_{N_2}$	7.52

We now have enough information to implement STANJAN. The following block of FORTRAN code shows this implementation:

X = 0.
X(I_CH4) = 0.09506
X(I_O2 ) = 0.19011
X(I_N2 ) = 0.71483
T = 298.
P = patm2cgs(1.)
call set_XT(cpr,phi,X,T)

call equil(cpr,phi,P)

A few of the important points of this implementation are presented:

• As I mentioned before, STANJAN has been integrated into CHEMKIN, and CHEMKIN operates by default with the Centimetre Gram Second (cgs) system of units. Therefore, the "patm2cgs" function and "set_XT" subroutine make the appropriate unit conversions to cgs units. Although the cgs system of units is not used or taught for most engineering applications, it is particularly well suited for experiments in chemistry, as the scales of the cgs system match well with the scales observed in experimentation.
• The variable "phi" in the code snippet above is actually an array of size $N_s+1$ elements. The first $N_s$ elements being the mole fractions, and the last element being the temperature.
• The variable "cpr" in the code snippet above is a structure of different data types needed for the communication of this code with the CHEMKIN libraries. For our purposes it is unnecessary to talk about it in any more detail.

The above code snippet will produce the following results:

X_H2   =   3.463047028581901E-003
X_H    =   3.957719341563813E-004
X_O2   =   5.402008329332355E-003
X_OH   =   3.063612303400087E-003
X_H2O  =   0.183569104732570     
X_HO2  =   5.721846447477325E-007
X_H2O2 =   5.027619424992306E-008
X_CH3  =   5.443802177681839E-017
X_CH4  =   2.516705383905423E-017
X_CO   =   8.695652836372861E-003
X_CO2  =   8.568543729085364E-002
X_CH2O =   1.225781883959981E-011
X_C2H2 =   7.532021965068877E-022
X_C2H4 =   5.563498080974568E-027
X_C2H6 =   4.124612875897405E-033
X_N2   =   0.709724743071636     
sum(X) =    1.00000000000000     
T_ad   =    2231.68195644285

Chemical equilibrium for this sample problem is thus represented by the composition given above. Some of the important points of this composition are discussed below:

• It is important to note that the sum of the mole fractions equal one. This was a necessary condition as indicated by Eq.(27)
• Our reaction was stoichiometric. That is, there was enough oxygen supplied for complete combustion of the methane gas. As you can see the methane gas was completely consumed when compared to the other major species in the above snippet.
• "T_ad" is the adiabatic flame temperature, and to this point hasn't been defined yet. It is the topic of the next section.

---

Adiabatic Flame Temperature

For a constant pressure system we can define the adiabatic temperature as follows:

(30)

\begin{equation} H_{reac}(T_i,P)=H_{prod}(T_{ad},P) \end{equation}

That is, the adiabatic flame temperature is the temperature at which the total enthalpy of the products is equivalent to the total enthalpy of the reactants at the initial temperature. We can write the above equation in terms of partial molar quantities:

(31)

\begin{align} \sum_{reac}N_i\bar{h}_i=\sum_{prod}N_i\bar{h}_i \end{align}

with

(32)

\begin{align} \bar{h}_i=\bar{h}^o_{f,i}+\int^T_{298}\bar{c}_{p,i}(T)\mathrm{d}T \end{align}

where $\bar{h}^o_{f,i}$ is known as the enthalpy of formation. The enthalpy of formation, and the functional dependence of the specific heat on the temperature are treated appropriately by the JANAF thermochemical data tables. Often it is assumed that the specific heat is constant, or the specific heat for some average temperature is used in calculation. In these cases Eq.(32) can be easily integrated:

(33)

\begin{align} \bar{h}_i=\bar{h}^o_{f,i}+\bar{c}_{p,i}(T-298) \end{align}

Based on the above definition STANJAN can easily calculate the adiabatic flame temperature in two steps. The first step is to calculate the equilibrium composition from an initial state. Once the equilibrium composition is known, the total enthalpy of the products can easily be calculated in addition to the total enthalpy of the reactants. The adiabatic flame temperature is then the only unknown quantity and can easily be solved for. A sample problem will help to illustrate this.

---

Calculation of Adiabatic Flame Temperature:

For the purpose of comparison, lets use the same reaction we used in our STANJAN example, but for the sake of simplicity lets assume that there is no dissociation of the products. Therefore, we have:

(34)

\begin{align} CH_4+2O_2+7.52N_2 \rightarrow CO_2+2H_2O+7.52N_2 \end{align}

Now, again for simplicity's sake, lets assume the specific heat is constant. We will use an average value for the specific heat, lets say $1200~K=0.5(300~K+2100~K)$. Thus, the appropriate thermochemical data for standard pressure is presented below:

Species	$\bar{h}^o_{f,i}$ (kJ/kmol) @ 298 K	$\bar{c}_{p,i}$ (kJ/kmol K) @ 1200 K
$\mathrm{CH_4}$	-74,831	-
$\mathrm{CO_2}$	-393,546	56.21
$\mathrm{H_2O}$	-241,845	43.87
$\mathrm{N_2}$	0	33.71
$\mathrm{O_2}$	0	-

Now we have all the information necessary to calculate the adiabatic flame temperature. We do so using equations 31 and 33:

(35)

\begin{equation} H_{reac}=(1)(-74,831)+(2)(0)+(7.52)(0)=-74,831 \end{equation}

(36)

\begin{equation} H_{prod}=(1)[-393,546+56.21(T_{ad}-298)]+(2)[-241,845+43.87(T_{ad}-298)]+(7.52)[0+33.71(T_{ad}-298)] \end{equation}

Therefore, we have:

(37)

\begin{align} H_{reac}(T_i,P)=H_{prod}(T_{ad},P)\Rightarrow~T_{ad}=2318~K \end{align}

What is surprising is that for all for the assumptions we made, the adiabatic flame temperature of 2318 K compares rather well to the 2232 K calculated by STANJAN. It turns out that relaxing the assumption for constant specific heat changes the calculated adiabatic flame temperature very little. Therefore, the majority of the discrepancy between these two values is due primarily to the dissociation of the products.

---

This ends our discussion of chemical equilibrium calculation methods. For some applications to turbulent reactive flows see the next tab.

Applications to Turbulent Reactive Flows

Time Scales

Before we begin to discuss any relevance that chemical equilibrium might have with turbulent reactive flows, we must first have a discussion about some appropriate time scales. In literature on turbulent reactive flows, the chemistry is often described in one of two ways:

• finite-rate (non-equilibrium)
• infinitely fast (equilibrium)

I will use the figure below to explain what these two regimes are.

In the figure above, on the left we are presented with chemical time scales, and on the right physical time scales of the turbulent flow. When the chemical time scales and physical time scales are comparable (the intermediate and slow zone on the figure above) we say that the chemistry is in non-equilibrium or of a finite-rate. However, when the chemical time scales are much faster than the physical time scales we say that the chemistry is infinitely fast or in equilibrium.

There exist a dimensionless number to describe this effect. It is known as the Damkohler number, and the important results mentioned above can be presented in terms of the Damkohler number:

• infinitely fast (equilibrium) regime

(38)

\begin{align} Da=\frac{\tau_t}{\tau_c}\gg 1 \end{align}

• finite-rate (non-equilibrium) regime

(39)

\begin{align} Da=\frac{\tau_t}{\tau_c}\approx O(1) \end{align}

where $\tau_t$ represents to the turbulent time scales, and $\tau_c$ represents the chemical time scales. From this point on we will only concern ourselves with the infinitely fast regime, that is when the Damkohler number is much greater than one.

Conserved Scalar Approach

When we are in the infinitely fast chemistry regime, we can simplify our calculations by applying the conserved scalar approach. The theory behind this approach is that chemical reaction is completed (carried out until equilibrium) as soon as the reactants are mixed. This would imply that our species composition and temperature anywhere in the domain should be directly related to the "mixedness" of the flow at that location. Therefore, we need to define a variable that adequately describes the mixing, as well as a transport equation for this variable that gives us the evolution of this variable in space and time.

For our purposes, the variable that will adequately describe the "mixedness" of the turbulent flow is the mixture fraction. We will develop the definition of the mixture fraction below. However, we first need to define elemental mass fractions:

(40)

\begin{align} Z_i=\sum^{N_s}_{j=1}\mu_{ij}Y_j,~~~~i=1,2,...,a \end{align}

where $\mu_{ij}$ represents the grams of the ith element in a gram of the jth species. An example calculation is presented below:

---

Calculation of Elemental Mass Fractions:

Suppose we have 1 species $(H_2O)$. Then we clearly have two elements (O,H), and we can calculate the elemental mass fractions as follows:

(41)

\begin{align} Z_O=\frac{16}{18}Y_{H_2O} \end{align}

(42)

\begin{align} Z_H=\frac{2}{18}Y_{H_2O} \end{align}

where the numbers above come from the atomic weights for an oxygen atom, hydrogen atom, and water molecule.

---

Now that we have defined elemental mass fractions, we can define the mixture fraction in terms of the elemental mass fractions:

(43)

\begin{align} \xi=\frac{Z_i-Z^O_i}{Z^F_i-Z^O_i} \end{align}

where $Z^F_i$ represents the elemental mass fraction of the fuel stream, and $Z^O_i$ represents the elemental mass fraction of the oxidizer stream. These values are known a priori, and so long as the fuel supply and oxidizer supply is constant, then so are these values. It follows then that the mixture fraction is effectively a normalization of the elemental mass fraction, and can only take on the values between [0,1]. From the definition above, it is evident that in the fuel stream the mixture fraction is 1, and in the oxidizer stream the mixture fraction is 0. The figure below helps to illustrate the properties of the mixture fraction:

The figure above shows the mixture fraction for a turbulent reacting jet. The center stream is the fuel stream, and the outer regions are a coflow of air (oxidizer). All of the properties discussed above are thus evident in this figure. The intermediate values represent the mixing of fuel and oxidizer.

Now, the next step is to develop a transport equation for the mixture fraction. To do this we will first begin with the species conservation equation:

(44)

\begin{align} \frac{\partial \rho Y_j}{\partial t} + \frac{\partial \rho u_k Y_j}{\partial x_k}=\frac{\partial}{\partial x_k}\left( \rho D \frac{\partial Y_j}{\partial x_k} \right)+\dot{w}_j,~~~~j=1,2,...,N_s \end{align}

where $\dot{w}_j$ represents the rate of production/destruction of species j. It is easy to write this equation in terms of the elemental mass fractions. First multiply by $\mu_{ij}$, then sum over all j. Doing so yields the following:

(45)

\begin{align} \frac{\partial \rho Z_i}{\partial t} + \frac{\partial \rho u_k Z_i}{\partial x_k}=\frac{\partial}{\partial x_k}\left( \rho D \frac{\partial Z_i}{\partial x_k} \right),~~~~i=1,2,...,a \end{align}

where

(46)

\begin{align} \sum^{N_s}_{j=1}\mu_{ij}\dot{w}_j=0,~~~~i=1,2,...,a \end{align}

Equation 46 is easily understood. Since $\dot{w}_j$ represents the rate of production/destruction of species j, then $\sum^{N_s}_{j=1}\mu_{ij}\dot{w}_j$ represents the rate of production/destruction of element i. We know that the number of elements cannot change, and thus the quantity is necessarily zero for each element.

The final step is to substitute for the elemental mass fractions with the mixture fraction according to Eq.(43). This step is less obvious than the previous steps so I will present it in more detail:

(47)

\begin{align} \frac{\partial \rho (\xi(Z^F_i-Z^O_i)+Z^O_i)}{\partial t} + \frac{\partial \rho u_k (\xi(Z^F_i-Z^O_i)+Z^O_i)}{\partial x_k}=\frac{\partial}{\partial x_k}\left( \rho D \frac{\partial (\xi(Z^F_i-Z^O_i)+Z^O_i)}{\partial x_k} \right) \end{align}

Next we can pull all constants out of the differential terms:

(48)

\begin{align} (Z^F_i-Z^O_i)\frac{\partial \rho \xi}{\partial t} + Z^O_i\frac{\partial \rho}{\partial t} + (Z^F_i-Z^O_i)\frac{\partial \rho u_k \xi}{\partial x_k} + Z^O_i\frac{\partial \rho u_k}{\partial x_k} = (Z^F_i-Z^O_i)\frac{\partial}{\partial x_k}\left( \rho D \frac{\partial \xi}{\partial x_k} \right) \end{align}

Now we can simplify:

(49)

\begin{align} (Z^F_i-Z^O_i)\left[\frac{\partial \rho \xi}{\partial t} + \frac{\partial \rho u_k \xi}{\partial x_k}\right] + Z^O_i\left[\frac{\partial \rho}{\partial t} + \frac{\partial \rho u_k}{\partial x_k}\right] = (Z^F_i-Z^O_i)\frac{\partial}{\partial x_k}\left( \rho D \frac{\partial \xi}{\partial x_k} \right) \end{align}

where the second bracketed term is the continuity equation, and thus is necessarily zero. Therefore, the remaining coefficients cancel out, and we are left with the conserved scalar transport equation:

(50)

\begin{align} \frac{\partial \rho \xi}{\partial t} + \frac{\partial \rho u_k \xi}{\partial x_k} = \frac{\partial}{\partial x_k}\left( \rho D \frac{\partial \xi}{\partial x_k} \right) \end{align}

Now that we have defined the mixture fraction, and the transport equation for the mixture fraction, we can effectively describe the species compositions and temperature in terms of the mixture fraction. That is,

(51)

\begin{align} Y_j=Y^e_j(\xi),~~~~j=1,2,...,N_s \end{align}

(52)

\begin{align} T=T^e(\xi) \end{align}

where the superscript e indicates the equilibrium composition. From Eq.(43) we can define the elemental mass fractions in terms of the mixture fraction. It should be apparent that the elemental mass fractions are akin to the mole populations (Eq.(18)) that make up the constraint equations for the element potential method. Therefore, knowledge of the mixture fraction everywhere in the domain is enough information for calculation of the species compositions and temperature throughout the domain.

I should note here that the functional dependence of the species compositions and temperature on the mixture fraction is not always as simple as the relationships given in Eq.(51) and Eq.(52). In practice the relationship could be much more complicated. The compositions and temperature can be described by more than one variable (i.e. mixture fraction, progress variable, strain rate). A multitude of research exists on this topic. If the reader is curious, a few of the more prominent and successful methods are offered below:

• Intrinsic Low Dimensional Manifold (ILDM) Approach
  • Maas and Pope
  • Maas and Pope
• Flamelet Generate Manifold (FGM)
  • van Oijen

However, the general idea of all such infinitely fast models is to generate look-up tables for Eq.(51) and Eq.(52) no matter what the complexity of the functional relationship is, since all the information is known a priori (i.e. the mixture fraction is bound by [0,1]). In this regard, the figure presented below offers an example of the information that such a look-up table would contain.

From the figure above, it is obvious that knowledge of the mixture fraction will give values for the mass fractions of various species.

Comparison of Calculations - Equilibrium Vs. Non-equilibrium

In order to judge the success of a numerical method for the simulation of a turbulent reacting flow, we must validate our numerical methods by comparison with experiments. A popular set of experiments for validation purposes were conducted at Sandia National Laboratories. They developed 6 experiments, labeled A-F, for reacting flows, ranging from laminar to fully turbulent. I will present some results below for Sandia Flame D, a turbulent reacting flow for which the infinitely fast chemistry assumption is valid. All the details about Sandia Flame D can be found here.

For the results presented below, the calculations on the left were done by me. They represent the non-equilibrium approach. That is, I used no assumption of fast chemistry. For this approach I used an augmented reduced mechanism based on GRI-Mech 3.0. It contains 16 species and 12 reaction steps. Therefore, in theory, I solved Eq.(44) for all 16 species, with the chemical source terms coming from the 12 reaction steps in the reduced mechanism. The results presented on the right are the equilibrium calculations. They make use of the infinitely fast chemistry assumption. They were done in 2005 by a former student of our lab. I present both results here to show how accurate the infinitely fast chemistry assumption can be when applied to an appropriate flow.

The figure above shows the radial distribution of the time averaged temperature at various locations downstream. The points on the above figure represent the experimental values as measured by Sandia National Laboratories.

The figure above shows the radial distribution of the time averaged mass fraction of $CH_4$ at various locations downstream. The points on the above figure represent the experimental values as measured by Sandia National Laboratories.

The figure above shows the radial distribution of the time averaged mass fraction of $CO$ at various locations downstream. The points on the above figure represent the experimental values as measured by Sandia National Laboratories.

As you can see, the general agreement of the above results is very good. However, there is one major distinction to make, and that is in regard to the computation time. The results on the left, for non-equilibrium kinetics, were generated on 16 processors on a supercomputer at the Pittsburgh Supercomputing Center (PSC). It took about 750 hours of computation time to generate these results. The results on the right, for equilibrium kinetics, were generated on 6 processors of a Dell Xenon PC. It took about 48 hours to generate those results! As you can see, the infinitely fast chemistry assumption is extremely efficient, and should be employed when appropriate. Therefore, what I have outlined for you here, is a valuable tool for people interested in the research of turbulent reactive flows.

References

Bibliography

1. Carey, Van P. Statistical Thermodynamics and Microscale Thermophysics. Cambridge University Press, Cambridge, UK, 1999.

2. Dill, Ken A.; Bromberg, Sarina. Molecular Driving Forces: Statistical Thermodynamics in Chemistry and Biology. Garland Science, New York, NY, 2003.

3. Poinsot, T.; Veynante, D. Theoretical and Numerical Combustion. R. T. Edwards, Philadelphia, PA, 2001.

4. Turns, Stephen R. An Introduction to Combustion: Concepts and Applications. McGraw-Hill, New York, NY, 2000.

5. Abbott, M. M.; Van Ness, H. G. Thermodynamics with Chemical Applications. McGraw-Hill, New York, NY, second edition, 1989.

6. Reynolds, W. C. The Element Potential Method for Chemical Equilibrium Analysis: Implementation in the Interactive Program STANJAN. Technical Report, Stanford University, 1986.

7. Pope, S. B. CEQ: A Fortran Library to Compute Equilibrium Compositions Using Gibbs Function Continuation. http://eccentric.mae.cornell.edu/~tcg/CEQ, 2003.

8. Pope, S. B. The Computation of Constrained and Unconstrained Equilibrium Compositions of Ideal Gas Mixtures using Gibbs Function Continuation. Cornell University Report, 2003.

9. Pope, S. B. Gibbs Function Continuation for the Stable Computation of Chemical Equilibrium. Combustion and Flame, 139:222-226, 2004.

10. Centimetre Gram Second System of Units. http://en.wikipedia.org/wiki/Centimetre_gram_second_system_of_units, 2008.