**5. Foundation of reaction**

## **5.1. Chemical equilibrium**

A general chemical equilibrium reaction with *v′i,s* and *v′'i,s* representing the stoichiometric coefficient of reaction and products for the chemical species *Mi* 

$$\begin{array}{ccccc}\sum\_{l=1}^{N}\boldsymbol{\upsilon}\prime\_{l,\boldsymbol{s}} & \boldsymbol{M}\_{l} \end{array} \rightleftharpoons \begin{array}{c}\sum\_{l=1}^{N}\boldsymbol{\upsilon}\prime\_{l,\boldsymbol{s}} & \boldsymbol{M}\_{l}\;.\tag{36}$$

Numerical Simulation of Combustion in Porous Media 541

**6.1. The finite volume equations formulation** 

**Figure 2.** Computational molecule in 3D domain (Patankar, 2002)

For more information interested reader can see (Mohammadi, 2010).

*6.2.1. Discretisation of transient convection diffusion equation* 

that govern by equation (40). The fully implicit discretisation equation is:

ap Φp = aw Φw + aE ΦE + as Φs + aN ΦN + aB ΦB + aT ΦT + aº

**6.2. Discretisation and numerical solution of the momentum equation** 

Finally, the momentum equation for the calculation of velocity and pressure by use of continuity equation should be considered. For numerical reasons, it is recommendable to resort to so called staggered grid, i.e. pressure and velocity are calculated on computational grids shifted to each other, the pressure for example in the cells and the velocity on the nodes. The calculations of velocity commonly take place iteratively, for which several algorithms are known (e.g. SIMPLE, PISO, SIMPLER…). In final analysis, all have the fact in common that is first step the momentum equation is solved for the velocities of momentums kept constant. In the second step, pressure corrections are then calculated with the help of a Poisson equation. For pressure with these pressure corrections, new velocities are then calculated again, and that again, until a pre-given break off threshold for the convergence is

Transient three dimensional convection diffusion of a general property *Φ* in a velocity field

<sup>p</sup>Φ<sup>º</sup>

p + Su , (41)

relative to those for all other variables.

reached.

Finite volume equations are derived by the integration of above differential equations over finite control volumes that taken together fully cover the entire domain of interest (Fig. 2). These control volumes are called ''cells'' P, for which the fluid-property value, are regarded as representative of the whole cell. It is surrounded by neighboring nodes which we shall denote by N, S, E, W, B and T. Cells and nodes for velocity components are ''staggered''

State of equilibrium can be interpreted as a situation, in which both the forward as well as reverse reactions progresses with identical speed.

$$K\_p(T) = \prod\_l \left(\frac{p l}{p 0}\right)^{\left(\nu^{\prime\prime} l, s - \nu^{\prime} l, s\right)}.\tag{37}$$

The equilibrium constant *Kp* now contains the information about the equilibrium material composition in term of partial pressure *pi* of the various species *i .* For more information interested reader can see (Mohammadi 2010).

#### **5.2 Reaction kinetics**

A one step chemical reaction of arbitrary complexity can be represented by the following stoichiometric equation:

$$
\Sigma\_{l=1}^{N} \; \upsilon' \; \; M\_l \; \rightarrow \; \Sigma\_{l=1}^{N} \; \upsilon'' \; \; M\_l \; \tag{38}
$$

where *v′s* are the stoichiometric coefficient of reactions and *v′'i* representing the stoichiometric coefficient of products, *Mi* the specification of molecule of *i*th chemical species, and *N* total number of component involved. Usually, they are represented with an Arrhenious formulation form:

$$K = AT^b \exp\left(-\frac{E\_A}{RT}\right). \tag{39}$$

The constant *A* and the exponent *b* as well as the so-called activation energy *EA* are summarized for many chemical reactions in extensive table.

#### **6. The finite volume method**

Customarily, CFD codes work with the finite volume method. This approach guaranties the numerical preservation of conservative quantities for the incompressible flows. The finite volume (FV) method uses the integral form of the conservation of equations. The solution domain is subdivided into a finite number of control volumes, and the conservation equations are applied to each control volume. As result, an algebraic equation for each *CV* is obtained. The FV method accommodates any type of grid, so it is suitable for complex geometries. However, the computational mesh ideally, be built hexahedratically. The conservation law for transport of a scalar in an unsteady flow has the general form:

$$\frac{\partial}{\partial t}(\rho \,\Phi) + \nabla.\{\rho u \Phi\} = \nabla.\{\Gamma \nabla \Phi\} + \S\_{\Phi} \,\, \, \, \, \tag{40}$$

*(*�uΦ*)* designates convection, (��Φ*)* diffusion flows of and *SΦ* the corresponding local source. For more information interested reader can see (Mohammadi 2010).

## **6.1. The finite volume equations formulation**

540 Numerical Simulation – From Theory to Industry

reverse reactions progresses with identical speed.

interested reader can see (Mohammadi 2010).

**5.2 Reaction kinetics** 

stoichiometric equation:

Arrhenious formulation form:

**6. The finite volume method** 

the general form:

�

∑ �� �,�

��� �� ⇄ ∑ ���

��(�) = ∏ �

∑ �� �

�

summarized for many chemical reactions in extensive table.

State of equilibrium can be interpreted as a situation, in which both the forward as well as

�� ���

The equilibrium constant *Kp* now contains the information about the equilibrium material composition in term of partial pressure *pi* of the various species *i .* For more information

A one step chemical reaction of arbitrary complexity can be represented by the following

where *v′s* are the stoichiometric coefficient of reactions and *v′'i* representing the stoichiometric coefficient of products, *Mi* the specification of molecule of *i*th chemical species, and *N* total number of component involved. Usually, they are represented with an

� = ��� ��� �� ��

The constant *A* and the exponent *b* as well as the so-called activation energy *EA* are

Customarily, CFD codes work with the finite volume method. This approach guaranties the numerical preservation of conservative quantities for the incompressible flows. The finite volume (FV) method uses the integral form of the conservation of equations. The solution domain is subdivided into a finite number of control volumes, and the conservation equations are applied to each control volume. As result, an algebraic equation for each *CV* is obtained. The FV method accommodates any type of grid, so it is suitable for complex geometries. However, the computational mesh ideally, be built hexahedratically. The conservation law for transport of a scalar in an unsteady flow has

*(*�uΦ*)* designates convection, (��Φ*)* diffusion flows of and *SΦ* the corresponding local

source. For more information interested reader can see (Mohammadi 2010).

�

��� �� → ∑ ���

�

�,�

�,��

�����,� � ��

�

�� (� �) � �. (���) = �. (���) � �� , (40)

��� �� . (36)

� . (37)

��� ��, (38)

���. (39)

�

Finite volume equations are derived by the integration of above differential equations over finite control volumes that taken together fully cover the entire domain of interest (Fig. 2). These control volumes are called ''cells'' P, for which the fluid-property value, are regarded as representative of the whole cell. It is surrounded by neighboring nodes which we shall denote by N, S, E, W, B and T. Cells and nodes for velocity components are ''staggered'' relative to those for all other variables.

**Figure 2.** Computational molecule in 3D domain (Patankar, 2002)

For more information interested reader can see (Mohammadi, 2010).

#### **6.2. Discretisation and numerical solution of the momentum equation**

Finally, the momentum equation for the calculation of velocity and pressure by use of continuity equation should be considered. For numerical reasons, it is recommendable to resort to so called staggered grid, i.e. pressure and velocity are calculated on computational grids shifted to each other, the pressure for example in the cells and the velocity on the nodes. The calculations of velocity commonly take place iteratively, for which several algorithms are known (e.g. SIMPLE, PISO, SIMPLER…). In final analysis, all have the fact in common that is first step the momentum equation is solved for the velocities of momentums kept constant. In the second step, pressure corrections are then calculated with the help of a Poisson equation. For pressure with these pressure corrections, new velocities are then calculated again, and that again, until a pre-given break off threshold for the convergence is reached.

#### *6.2.1. Discretisation of transient convection diffusion equation*

Transient three dimensional convection diffusion of a general property *Φ* in a velocity field that govern by equation (40). The fully implicit discretisation equation is:

$$\mathbf{a}\_{\mathbf{P}} \bullet \mathbf{b}\_{\mathbf{P}} = \mathbf{a}\_{\mathbf{W}} \bullet \mathbf{b}\_{\mathbf{W}} + \mathbf{a}\_{\mathbf{E}} \bullet \mathbf{b}\_{\mathbf{E}} + \mathbf{a}\_{\mathbf{S}} \bullet \mathbf{b}\_{\mathbf{S}} + \mathbf{a}\_{\mathbf{N}} \bullet \mathbf{b}\_{\mathbf{N}} + \mathbf{a}\_{\mathbf{I}} \bullet \mathbf{b}\_{\mathbf{I}} + \mathbf{a}\_{\mathbf{I}} \bullet \mathbf{b}\_{\mathbf{I}} + \mathbf{a}\_{\mathbf{P}} \bullet \mathbf{b}\_{\mathbf{P}} + \mathbf{S}\_{\mathbf{U}} \tag{41}$$

where:

$$\mathbf{a}\_{P} = \mathbf{a}\_{\text{av}} + \mathbf{a}\varepsilon + \mathbf{a}\_{\text{v}} + \mathbf{a}\mathbf{v} + \mathbf{a}\mathbf{v} + \mathbf{a}\mathbf{v} + \mathbf{a}\_{P}^{\circ} + \Delta\text{F-S}\_{P} \tag{42}$$

Numerical Simulation of Combustion in Porous Media 543

*\** , where t lies on the z-

+ w′, (44)

� �� (45)

dimensional ones. Similar equations can be written for *vn\** and *wt*

+ p′, u = u\*

�� � ��

**Figure 4.** Control volume for driving the pressure correction equation (Patankar, 1980)

5. Calculate *u, v* and w from *u\*, v\** and *w\** using velocity correction equations.

6. Treat the corrected pressure *p* as the new guess *p\** and iterate the preceding procedure

The simple algorithm has the following main steps:

2. Solve the momentum equation to obtain *u\*,v\*,* and *w\**. 3. Solve the pressure correction equation to obtain *p'*. 4. Add *p'* to *p\** to obtain the corrected pressure *p*.

1. Guess the pressure field *p\*.* 

to convergence.

correct pressure and *p'* the pressure correction, we may write:

p = p\*

correction formula becomes:

equation in then developed as:

�� � ��

direction grid line between grid points P and T. if p is the correct pressure and *p* is the

+ u′, v = v\*

where the prime indicate corrections needed to reach the correct values that satisfy the continuity equation. Omitting the correction terms due to the neighbors, an iterative solution may be developed to solve for the pressure and the velocity field. Then, the velocity

> <sup>∗</sup> <sup>+</sup>�� �� ��� � � �� � �.

<sup>∗</sup> <sup>+</sup>�� �� ��� � � ��

And similarly for *wt*. From the time dependent continuity equation, the pressure correction

where *b* is a mass source which must be eliminated through pressure correction so that continuity is satisfied. Here, T and B are neighboring grid points on the z direction grid line.

+ v′, w = w\*

ap p'p = aE p'e + aw p'w + aN p'N + as p's + aT p'T + aB p'B +b, (46)

$$\text{with } a\_p^\circ = \frac{\rho\_{p\text{ AV}}^\circ}{4} \text{ and } \overline{\mathcal{S}}\Delta V = \mathcal{S}\_u + \mathcal{S}\_p\Phi\_p$$

For more information interested reader can see (Mohammadi, 2010).
