**2. Governing equations and models of gas properties**

## **2.1. Governing equations**

There are three numerical simulation methods for compressible reactive flows via solving Navier-Stokes equations, including Reynolds-averaged Navier-Stokes (RANS) equations, large eddy simulation (LES), and direct numerical simulation (DNS). RANS simulates all scales of fluctuations using turbulent models. However, there is still lack of a general model suitable for all flow types, since RANS ignores the details and impacts of turbulent fluctuations. LES divides the flow scales into large and small scales via a low-pass filter. The large-scale flows are solved using the governing equations, and thus LES is able to give the temporal and spatial evolution of large-scale eddies. The small-scale flows are simulated using a sub-grid-scale (SGS) model, but there is also no general SGS model. DNS solves the Navier-Stokes equations for all scales of flows without introducing any models and assumptions, and thus it can accurately predict the full-scale flow characteristics. However, its main limitation in practice

In the early days, RANS or LES was used to simulate the compressible reactive flows consid‐ ering the computing cost. Dauptain et al. [3] simulated the supersonic hydrogen flame using the large eddy simulation, and their numerical results agreed well with the experiment measurements. LES was also used to simulate the supersonic flow and combustion in a model scramjet combustor [4], and the velocity and temperature at different cross-sections were compared with the experimental data. Edwards et al. used the hybrid method of LES/RANS to simulate the supersonic reacting wall-jet flow [5]. There are also other examples using RANS or LES [6–9]. Recently, some researchers have used DNS to simulate the compressible reactive flows. For example, O'Brien et al. [10] simulated the time-developing, H2/air-reacting mixing layers and investigated the dynamics of backscatter of kinetic energy using DNS; Diegelmann et al. [11] did the direct simulations of reactive shock-bubble interaction with detailed chem‐

One of the main difficulties in the simulation of compressible flows is how to distinguish the discontinuities which may widely exist in supersonic flows. Jiang and Shu [12] proposed the weighted essentially non-oscillatory (WENO) scheme to capture shock waves with almost no oscillations. But the WENO scheme is too dissipative in regions with large shear rates. To overcome this drawback, the hybrid scheme was hence developed, which switches the WENO scheme and the linear (high-order) scheme according to a discontinuity detector. Pirozzoli [13] proposed a hybrid compact-WENO scheme to obtain high resolutions for simulation of shockturbulence interaction. Furthermore, Ren et al. [14] proposed a Roe type, characteristic-wise compact-WENO scheme. Considering overcoming the complexity and program paralleliza‐ tion of compact schemes, Hu et al. [15] proposed a simple WENO scheme, which switches between the WENO scheme and its optimal linear scheme. A detailed review of numerical

The spray combustion can occur in the scenario of compressible flows [17–24]. Therefore, the methods to treat the dispersed droplets in the flows will be concerned by different approaches, such as Eulerian and Lagrangian methods. The Eulerian one needs complicated modelling of two-phase correlation terms, which is yet to be done. Even though the high computation cost will be taken by the Lagrangian method to track each particle released to the flow fields, it is the physically clear and mathematically simple method in simulations of two-phase reactive

is the high cost of calculation.

30 Modeling and Simulation in Engineering Sciences

istry under varying initial pressure.

methods used for high-speed flows can be seen in [16].

Since this chapter will extend the simulations to the compressible two-phase turbulent flows considering the dispersed droplets' evaporation and gas-phase combustion, the governing equations will describe the Lagrangian transport of droplets through a continuous, compres‐ sible and chemically reacting carrier gas flow.

The equations for the gas-phase including mass, momentum, and energy exchange between the gas and the dispersed phase read as

$$\frac{\partial \hat{\rho}}{\partial t} + \frac{\partial}{\partial \alpha\_{\\_}} \left( \rho u\_{\\_} \right) = \mathbb{S}\_{m} \tag{1}$$

$$\frac{\partial}{\partial t} \left( \rho u\_i \right) + \frac{\partial}{\partial \mathfrak{x}\_{\neq}} \left( \rho u\_i u\_j + P \delta\_{ij} - \boldsymbol{\tau}\_{\neq} \right) = \mathbf{S}\_F \tag{2}$$

$$
\lambda \frac{\partial}{\partial t} (\rho e\_t) + \frac{\partial}{\partial \mathbf{x}\_j} \left( (\rho e\_t + P)\mu\_j - \lambda \frac{\partial T}{\partial \mathbf{x}\_j} - \mu\_i \boldsymbol{\pi}\_{ij} \right) = \mathbf{S}\_Q \tag{3}
$$

$$\frac{\partial}{\partial t} \left( \rho \mathbf{Y}\_k \right) + \frac{\partial}{\partial \mathbf{x}\_j} \left( \rho \mathbf{Y}\_k \left( \mathbf{u}\_j + \mathbf{V}\_j^c \right) - \rho D\_k \frac{\partial \mathbf{Y}\_k}{\partial \mathbf{x}\_j} \right) = \mathbf{S}\_{\text{combustion},k} + \mathbf{S}\_{\mathbf{y}\_k} \tag{4}$$

Here, *ρ* is the gas-phase density, *u<sup>i</sup>* is the velocity, *T* is the temperature, *P* is the pressure, *Yk*is the mass fraction of the *k*th species, and *δij* is the Kronecker delta function. The right-hand side terms of the above equations *Sm*, *S<sup>F</sup>*, and *SQ* describe the two-phase couplings of mass, mo‐ mentum, and energy, respectively. *S*combustion,*<sup>k</sup>* is the source term due to combustion. *τij* is the Newtonian viscous stress tensor

$$\pi\_{ij} = \mu \left( \frac{\partial u\_i}{\partial \mathbf{x}\_j} + \frac{\partial u\_j}{\partial \mathbf{x}\_i} - \frac{2}{3} \frac{\partial u\_k}{\partial \mathbf{x}\_k} \delta\_{ij} \right) \tag{5}$$

A correction velocity *V*c *j*is added to the convection velocity in the species equations to ensure global mass conservation

$$\boldsymbol{V}\_{j}^{c} = \sum\_{k=1}^{N\_{c}} \boldsymbol{D}\_{k} \frac{\boldsymbol{W}\_{k}}{\boldsymbol{W}} \frac{\partial \mathbf{X}\_{k}}{\partial \boldsymbol{\alpha}\_{i}} \tag{6}$$

Here, *NS* is the total number of species. *Wk* and *Xk* are the molecule weight and the mole fraction of the *k*th species, respectively. *e*<sup>t</sup> is the total energy, that is, kinetic energy plus internal (containing chemical) energy, which is defined as follows:

$$e\_t = \sum\_{k=1}^{N\_\xi} Y\_k \left( \int\_{T\_{\text{out}}}^T c\_{p,k} \, \mathbf{d} \, T + h\_{f,k}^0 \right) - \frac{P}{\rho} + \frac{u\_i u\_i}{2} \tag{7}$$

where *cp,k*and *h*0 *f,k* are the specific heat capacity at constant pressure and specific chemical formation enthalpy at the reference temperature, *T*ref, of the *k*th species, respectively. The values of *cp,k*are represented by a fifth-order polynomial, and the coefficients can be found in [25]

$$\frac{c\_{p,k}W\_k}{R} = a\_{1,k} + a\_{2,k}T + a\_{3,k}T^2 + a\_{4,k}T^3 + a\_{5,k}T^4 \tag{8}$$

The equation of state used for ideal multispecies gas is introduced to close the above equations as

$$P = RT \sum\_{k=1}^{N\_\xi} \frac{Y\_k}{W\_k} \tag{9}$$

where *R* is the universal gas constant.

#### **2.2. Models of gas properties**

The kinetic theory for gases is used to compute the transport and thermodynamic properties of the mixture. The Lennard-Jones potentials are used to calculate the inter-molecular forces. The thermal conductivity of each species is determined by the modified Eucken model. The dynamic viscosity of each species and the binary diffusivity are computed with the Chapman-Enskog theory, and their formulas can be found in [26]

$$
\mu\_k = \frac{26.69\sqrt{W\_k T}}{\sigma\_k^2 \Omega\_{v,k}} \tag{10}
$$

$$\frac{\lambda\_k W\_k}{\mu\_k \mathbb{C}\_{v,k}} = 1.32 + \frac{1.77}{\left(\mathbb{C}\_{p,k} / R\right) - 1} \tag{11}$$

$$\mathbf{D}\_{\rm A^{\rm B}} = \frac{0.00266 \mathbf{T}^{3/2}}{\mathbf{P} \mathbf{W}\_{\rm A^{\rm B}}^{1/2} \sigma\_{\rm A^{\rm B}}^2 \Omega\_{\rm D}} \tag{12}$$

The viscosity and thermal conductivity of the mixture are calculated by the expressions of Wake and Wassiljewa [26], respectively,

$$\mu\_{\boldsymbol{w}} = \sum\_{i=1}^{n} \frac{\mathbb{X}\_{i}\mu\_{i}}{\sum\_{j=1}^{n} \mathbb{X}\_{j}\phi\_{ij}} \tag{13}$$

where *φij* = (1 + (*μi* / *μj* )1/2(*Wj* / *Wi* )1/4)2 (8(1 + *Wi* / *Wj* ))1/2 .

mentum, and energy, respectively. *S*combustion,*<sup>k</sup>* is the source term due to combustion. *τij* is the

æ ö ¶ ¶ ¶ = +- ç ÷ ¶¶ ¶ è ø 2 3 *j i k ij ij ji k u u u τ μ δ xx x*

A correction velocity *V*c *j*is added to the convection velocity in the species equations to ensure

*k i W X V D*

Here, *NS* is the total number of species. *Wk* and *Xk* are the molecule weight and the mole fraction

0

where *cp,k*and *h*0 *f,k* are the specific heat capacity at constant pressure and specific chemical formation enthalpy at the reference temperature, *T*ref, of the *k*th species, respectively. The values of *cp,k*are represented by a fifth-order polynomial, and the coefficients can be found in [25]

> =+ + + + , <sup>234</sup> 1, 2, 3, 4, 5,

*kk k k k*

*a aT aT aT aT*

The equation of state used for ideal multispecies gas is introduced to close the above equations

= <sup>=</sup> å1

The kinetic theory for gases is used to compute the transport and thermodynamic properties of the mixture. The Lennard-Jones potentials are used to calculate the inter-molecular forces.

*NS k k k*

*<sup>Y</sup> P RT*

*NS <sup>T</sup> i i*

*k k*

*W x* (6)

*<sup>ρ</sup>* (7)

*<sup>W</sup>* (9)

is the total energy, that is, kinetic energy plus internal

2

*<sup>R</sup>* (8)

= ¶ <sup>=</sup> ¶ <sup>å</sup> <sup>c</sup> 1

( ) <sup>=</sup> <sup>=</sup> å + -+ ò ref

*k pk f k <sup>T</sup> <sup>k</sup> <sup>P</sup> u u e Y c Th*

d

t ,,

*j k*

*Ns*

(5)

Newtonian viscous stress tensor

32 Modeling and Simulation in Engineering Sciences

global mass conservation

of the *k*th species, respectively. *e*<sup>t</sup>

as

(containing chemical) energy, which is defined as follows:

1

*pk k*

*c W*

where *R* is the universal gas constant.

**2.2. Models of gas properties**

$$\lambda = \sum\_{i=1}^{n} \frac{\mathbf{X}\_i \mathbf{\hat{X}}\_i}{\sum\_{j=1}^{n} \mathbf{X}\_j \mathbf{A}\_{ij}} \tag{14}$$

where *Aij* = (1 + (*λ<sup>i</sup>* / *λ<sup>j</sup>* )1/2(*Wj* / *Wi* )1/4)2 (8(1 + *Wi* / *Wj* ))1/2 .

#### **2.3. Equations of spray droplets**

The evaporating fuel droplets are solved individually under the Lagrangian framework. Before describing the Lagrangian model, several usual assumptions should be formulated. The spray of droplets is sparsely dispersed and every droplet is unaware of the existence of the others. The droplet rotations are neglected and an infinite heat conduction coefficient is assumed. Therefore, the inner temperature distribution of each droplet remains uniform. Because of the high-density ratio between the liquid and gas phases, only the drag force acting on the droplets is significant. The Lagrangian equations for the position (in the *i* th direction), *x*d,*<sup>i</sup>* , velocity, *u*d,*<sup>i</sup>* , temperature, *T*d, and mass, *m*d, of a single droplet are given by

$$\frac{d\boldsymbol{u}\_{d,i}}{d\mathbf{t}} = \frac{\boldsymbol{F}\_d}{\boldsymbol{m}\_d} = \left(\frac{\boldsymbol{f}\_1}{\boldsymbol{\tau}\_d}\right) \left(\boldsymbol{u}\_i - \boldsymbol{u}\_{d,i}\right) \tag{15}$$

$$\frac{\mathbf{d}u\_{\mathbf{d},l}}{\mathbf{d}t} = \frac{F\_{\mathbf{d}}}{m\_{\mathbf{d}}} = \left(\frac{f\_{\mathbf{l}}}{\tau\_{\mathbf{d}}}\right) \left(u\_{l} - u\_{\mathbf{d},l}\right) \tag{16}$$

$$\frac{dT\_d}{\mathbf{d}t} = \frac{Q\_\mathbf{d} + \dot{m}\_\mathbf{d}L\_\mathbf{v}}{m\_\mathbf{d}c\_\mathbf{L}} = \left(\frac{f\_2}{\tau\_\mathbf{d}}\right) \left(\frac{\mathbf{Nu}}{3\mathbf{Pr}}\right) \left(\frac{c\_p}{c\_\mathbf{L}}\right) (T - T\_\mathbf{d}) + \left(\frac{\dot{m}\_\mathbf{d}}{m\_\mathbf{d}}\right) \frac{L\_\mathbf{v}}{c\_\mathbf{L}} \tag{17}$$

$$\frac{\mathbf{d}m\_{\rm d}}{\mathbf{d}t} = \dot{m}\_{\rm d} = -m\_{\rm d} \left(\frac{1}{\tau\_{\rm d}}\right) \left(\frac{\mathbf{Sh}}{3\mathbf{Sc}}\right) \text{In} \left(1 + B\_{\rm M}\right) \tag{18}$$

where *cp* is the specific heat of mixture gas, *c*<sup>L</sup> is the specific heat of the liquid. The momentum response time, *τ*d, is defined as

$$
\tau\_d = \frac{\rho\_d d\_d^2}{18\mu} \tag{19}
$$

where *d*d is the droplet diameter. The gas-phase Prandtl and Schmidt numbers in gaseous phase are given by

$$\mathbf{Pr} = \frac{\mu \mathbf{c}\_p}{\lambda}, \mathbf{Sc} = \frac{\mu}{\rho D} \tag{20}$$

respectively. The Nusselt and Sherwood numbers are

$$\mathbf{Nu} = \mathbf{2} + 0.552 \,\mathrm{Re}\_d^{1/2} \,\mathrm{Pr}^{1/3} \tag{21}$$

$$\text{Sh} = 2 + 0.552 \,\text{Re}\_d^{1/2} \,\text{Sc}^{1/3} \tag{22}$$

respectively. Here, the droplet Reynolds number based on the slip velocity is defined as

Numerical Simulation of Compressible Reactive Flows http://dx.doi.org/10.5772/63906 35

$$\text{Re}\_d = \frac{\rho |u\_i - u\_{d,i}| d\_d}{\mu} \tag{23}$$

*f*1 is an empirical correction to the Stokes drag law and is given as

Because of the high-density ratio between the liquid and gas phases, only the drag force acting

, temperature, *T*d, and mass, *m*d, of a single droplet are given by

*i di*

d,i

d

(17)

*λ ρD* (20)

= + 1/2 1/3 Nu 2 0.552Re Pr <sup>d</sup> (21)

<sup>d</sup> Sh 2 0.552Re Sc (22)

( ) *<sup>d</sup> <sup>d</sup>*

( ) *<sup>i</sup> i*

d L d L d L

*d d du <sup>F</sup> <sup>f</sup> u u*

*<sup>u</sup> <sup>F</sup> <sup>f</sup> u u*

æ ö = = - ç ÷ ç ÷ è ø t

d d

( ) *<sup>p</sup> d V dT Q m L <sup>f</sup> <sup>c</sup> m L T T t mc c m c* d d 2 d V

Nu

<sup>+</sup> æ öæ öæ ö æ ö <sup>=</sup> <sup>=</sup> ç ÷ç ÷ç ÷ - + ç ÷ ç ÷ ç ÷ ç ÷ è øè øè ø è ø & &

> ( ) *<sup>m</sup> m m B*

æ öæ ö <sup>=</sup> <sup>=</sup> - + ç ÷ç ÷ ç ÷ è øè ø

t

=

Pr ,Sc = = *<sup>p</sup> μc <sup>μ</sup>*

= + 1/2 1/3

respectively. Here, the droplet Reynolds number based on the slip velocity is defined as

d d M d <sup>d</sup> 1 Sh In 1 d 3Sc

where *cp* is the specific heat of mixture gas, *c*<sup>L</sup> is the specific heat of the liquid. The momentum

2 d d <sup>d</sup> 18 *ρ d <sup>τ</sup> μ*

where *d*d is the droplet diameter. The gas-phase Prandtl and Schmidt numbers in gaseous

*d m*

*t m* d, d 1

t

d d

d 3Pr

&

*t* d

respectively. The Nusselt and Sherwood numbers are

response time, *τ*d, is defined as

phase are given by

th direction),

(15)

(16)

(18)

(19)

on the droplets is significant. The Lagrangian equations for the position (in the *i*

,i 1 , t æ ö = = - ç ÷ ç ÷ <sup>t</sup> è ø

*x*d,*<sup>i</sup>*

, velocity, *u*d,*<sup>i</sup>*

34 Modeling and Simulation in Engineering Sciences

$$f\_1 = 1 + 0.15 \operatorname{Re}\_d^{0.687} \tag{24}$$

The corrections of heat transfer for an evaporating droplet are given as

$$\,\_{1}f\_{2} = \frac{\beta}{e^{\beta} - 1} \tag{25}$$

where the nondimensional evaporation parameter is given by

$$\beta = -1.5 \Pr \, \pi\_d \left(\frac{\dot{m}\_d}{m\_d}\right) \tag{26}$$

The evaporation rate is controlled by the mass transfer number, *B*M , and is given by

$$B\_{\rm M} = \frac{Y\_{\rm sf} - Y\_{\rm V}}{1 - Y\_{\rm sf}} \tag{27}$$

Here, *Y*<sup>V</sup> is the mass fraction of the vapor on the far-field condition for the droplets and *Y*sf is the vapor surface mass fraction calculated directly from the surface molar fraction (*χ*sf), which is obtained using the equilibrium assumption

$$Y\_{sl} = \frac{X\_{sl}}{\chi\_{sl} + \left(1 - \chi\_{sl}\right)W / W\_{\rm V}}\tag{28}$$

$$\chi\_{\rm sd} = \frac{P\_{\rm atm}}{P} \exp\left(\frac{L\_{\rm v}}{R/W\_{\rm v}} \left(\frac{1}{T\_{\rm B,L}} - \frac{1}{T\_{\rm d}}\right)\right) \tag{29}$$

where *W* is the mean molecular weight of mixture gas, *W*V is the molecular weight of the vapor, *P*atm is the atmospheric pressure, *R* is the universal gas constant, and *T*B,L is the liquid boiling temperature at *P*atm. *L*<sup>V</sup> is the latent heat of fuel vapor at liquid temperature, *T*d, and is given as follows:

$$L\_{\rm V} = L\_{\rm V, T\_{\rm B,L}} \left( \frac{T\_{\rm C,L} - T\_{\rm d}}{T\_{\rm C,L} - T\_{\rm B,L}} \right) \tag{30}$$

where *<sup>L</sup>* V,*T*B,L is the latent heat at *T*B,L and *T*C,L is the critical temperature.

#### **2.4. Two-way coupling models**

The source terms *S*, due to interactions between continuous phase and dispersed phase, are expressed by summating total droplets existing in one computation cell, denoting as *N*c,

$$\mathcal{S}\_{\rm m} = -\frac{1}{\Delta V} \sum\_{N\_{\rm c}} \left( \dot{m}\_{\rm d} \right) \tag{31}$$

$$S\_{\rm f} = -\frac{1}{\Delta V} \sum\_{N\_c} \left( F\_{\rm d} + \dot{m}\_{\rm d} \mu\_{\rm d,i} \right) \tag{32}$$

$$S\_Q = -\frac{1}{\Delta V} \sum\_{N\_\epsilon} \left( Q\_\mathbf{d} + \dot{m}\_\mathbf{d} \left( \frac{\mu\_{\mathbf{d},l} u\_{\mathbf{d},l}}{2} + h\_{\mathbf{V},\mathbf{st}} \right) \right) \tag{33}$$

$$\mathcal{S}\_{\mathbf{y}\_{\mathbf{i}}} = \begin{cases} -\frac{1}{\Delta V} \sum\_{\mathbf{N}} \left( \dot{m}\_{\mathbf{d}} \right) & \text{for fuel} \\ 0 & \text{for other species} \end{cases} \tag{34}$$

where *m*d and *m*˙ <sup>d</sup> denote the droplet mass and d*m*d/d*t*, respectively, *u*d,*<sup>i</sup>* is the droplet velocity, Δ*V* is the volume of the computational grid for the gas-phase calculation, and *h*V,sf is the evaporated vapor enthalpy at the droplet surface. *F*d and *Q*<sup>d</sup> are the drag force and the convective heat transfer, respectively, described later. For the combustion reaction model, a one-step global reaction is adopted for the reaction of hydrocarbon-air mixtures. The source term, *S*combustion,*<sup>k</sup>*, is expressed using the combustion reaction rate per unit volume, *R*F, as follows:

$$\mathcal{S}\_{\text{conclusion},k} = -\frac{\boldsymbol{\mu}\_k}{\boldsymbol{\mu}\_\text{F}} \frac{\boldsymbol{W}\_k}{\boldsymbol{W}\_\text{F}} \boldsymbol{R}\_\text{p} \tag{35}$$

where *nk* and *n*F are the molar stoichiometric coefficients of the *k*th species and the fuel for the one-step global reaction (positive for the production side), respectively. *Wk* and *W*F are the molecular weights of the *k*th species and fuel, respectively.
