**3. Numerical discretization procedures**

æ ö - <sup>=</sup> ç ÷ - è ø B ,L C,L d

is the latent heat at *T*B,L and *T*C,L is the critical temperature.

C,L B,L

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,

> = - å( & ) c d

=- + å( & ) c

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

<sup>1</sup> for fuel

( )

*m*

Δ *F i N S F mu*

1

c

*N*

*<sup>k</sup> Y N*

ì ï- = í ï î

molecular weights of the *k*th species and fuel, respectively.

*S V* <sup>d</sup>

Δ

*Q*

d d d,

æ ö

0 for other species

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:

= - combustion, <sup>F</sup>

*k n W S R*

F F *k k*

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

d, d, d d V,sf <sup>1</sup> ( ) Δ 2

*i i*

*u u S Qm h <sup>V</sup>* (33)

*N S m*

1 Δ *<sup>m</sup>*

*T T* (30)

*<sup>V</sup>* (31)

*<sup>V</sup>* (32)

å & (34)

*n W* (35)

*T T*

V V,

*L L*

where *<sup>L</sup>* V,*T*B,L

**2.4. Two-way coupling models**

36 Modeling and Simulation in Engineering Sciences

*T*

A finite difference methodology is used to solve the above governing equations. An explicit Runge-Kutta time-integration methodology is applied, obtaining a third-order time-accurate computation. For example, for a semi-discrete equation of d*Uj* = *L <sup>j</sup>* (*U* )d*t*, the iteration from *n* to iteration *n*+1 is performed as

$$\begin{aligned} \mathbf{U}I\_{j}^{(1)} &= \mathbf{U}I\_{j}^{n} + \Delta t L\_{j} \left(\mathbf{U}^{n}\right) \\ \mathbf{U}I\_{j}^{(2)} &= \frac{3}{4}\mathbf{U}^{n}\_{j} + \frac{1}{4}\mathbf{U}^{(1)}\_{j} + \frac{1}{4}\Delta t L\_{j} \left(\mathbf{U}^{(1)}\_{j}\right) \\ \mathbf{U}I\_{j}^{n+1} &= \frac{1}{3}\mathbf{U}^{n}\_{j} + \frac{2}{3}\mathbf{U}^{(2)}\_{j} + \frac{2}{3}\Delta t L\_{j} \left(\mathbf{U}^{(2)}\_{j}\right) \end{aligned} \tag{36}$$

where *Uj* is the conserved variable at the direction *j*.

The higher-resolution numerical scheme applied for the resolution of supersonic flow is essential to the reliability of simulation. The nonviscous flux *f* ^ *<sup>j</sup>*+1/2, for the interface at *j*+1/2, is evaluated using a fifth-order hybrid Compact-Weighted Essentially Non-Oscillatory scheme [14] for the resolution of turbulent field and shock-capturing calculation in the supersonic flow. The numerical flux of the hybrid scheme is calculated as

$$
\hat{f}\_{\backslash \ast 1/2} = \sigma\_{\backslash \ast 1/2} \hat{f}\_{\backslash \ast 1/2}^{\text{CU}} + \left(\mathbf{1} + \sigma\_{\backslash \ast 1/2}\right) \hat{f}\_{\backslash \ast 1/2}^{\text{WENO}} \tag{37}
$$

which is constructed by hybridizing a fifth-order Compact Upwind (CU) scheme, *f* ^ *j*+1/2 CU ,and a fifth-order WENO one, *f* ^ *j*+1/2 WENO, through a smoothness indicator *rc*. *<sup>σ</sup> <sup>j</sup>*+1/2 is the weight of the numerical flux,

$$\sigma\_{\boldsymbol{\gamma}^{\*1/2}} = \min \left( \mathbf{1}\_{\boldsymbol{\gamma}} \frac{r\_{\boldsymbol{\gamma}^{\*1/2}}}{r\_{\boldsymbol{\varepsilon}}} \right) \tag{38}$$

A sixth-order symmetric compact difference scheme is applied for the viscous diffusion terms,

$$\frac{1}{3}\hat{F}\_{j-1} + \hat{F}\_j + \frac{1}{3}\hat{F}\_{j+1} = \frac{28\left(\hat{f}\_{j+1} - \hat{f}\_{j-1}\right) + \left(\hat{f}\_{j+2} - \hat{f}\_{j-2}\right)}{36\Delta x} \tag{39}$$

where *F* ^ *j* is the difference approximation for ( <sup>∂</sup> *<sup>f</sup>* ^ <sup>∂</sup> *<sup>x</sup>* ) *<sup>j</sup>* at the node *j*. Because of the stiffness of the equation due to the reactions, the time integration of the equations is performed by means of the point-implicit method [27]. For instance, *U* is supposed as the unknown variable at time step *n*+1, and the below equation is solved to obtain *U*n+1=*U*<sup>n</sup> +Δ*U*

$$
\left(I - \left(\frac{\partial \mathcal{S}}{\partial U}\right)^{\
u} \Delta t\right) \Delta U = \left(-\mathcal{C}^{\
u} + V^{\
u} + \mathcal{S}^{\
u}\right) \Delta t \tag{40}
$$

where *C* is the convection term, *V* is the viscous term, *S* is the source term, and *I* is the unity matrix.

In order to obtain the physical quantities of the fluid at a droplet position, the fluid velocities computed on the Eulerian mesh were interpolated at the droplet locations. A fourth-order Lagrangian interpolation was used for this purpose [28]. When the droplet is located at (*x*d, *y*d, *z*d) in a cell whose coordinates are (*x*0, *y*0, *z*0)…(*xl* , *yl* , *zl* ), the fluid velocity at the droplet position is obtained via the Lagrangian interpolation, which is expressed as

$$\mu\left(\mathbf{x}\_{d},\mathbf{y}\_{d},\mathbf{z}\_{d}\right) = \sum\_{l=1-n/2}^{n/2} \sum\_{j=1-n/2}^{n/2} \sum\_{k=1-n/2}^{n/2} \times \left(\prod\_{\substack{l=1-n/2\\l\neq s\neq l}}^{n/2} \frac{\mathbf{x}\_{d}-\mathbf{x}\_{l}}{\mathbf{x}\_{l}-\mathbf{x}\_{l}} \prod\_{\substack{l=1-n/2\\(l\neq s)}}^{n/2} \frac{\mathbf{y}\_{d}-\mathbf{y}\_{l}}{\mathbf{y}\_{l}-\mathbf{y}\_{l}} \prod\_{\substack{l=1-n/2\\(l\neq k)}}^{n/2} \frac{\mathbf{z}\_{d}-\mathbf{z}\_{l}}{\mathbf{z}\_{l}-\mathbf{z}\_{l}}\right) \mathbf{u}\_{i,j,k} \tag{41}$$

where the subscripts *i, j, k*, and *l* represents the cell indexes and *n* = 4.

In the Lagrangian droplet-tracking method, the droplet equation of motion was integrated with the third-order Adams-Bash-forth scheme in direction, *i*, for droplet position,

$$\mathbf{x}\_{d,l}^{n+1} = \mathbf{x}\_{d,l}^{n} + \Delta t \left(\frac{\mathbf{23}}{1\mathbf{2}} \boldsymbol{\mu}\_{d,l}^{n} - \frac{16}{1\mathbf{2}} \boldsymbol{\mu}\_{d,l}^{n-1} + \frac{5}{1\mathbf{2}} \boldsymbol{\mu}\_{d,l}^{n-2}\right) \tag{42}$$

The equations of velocity and temperature of the droplet are solved using the same integration method, as Eq. (42).

#### **4. Numerical examples**

#### **4.1. Taylor-Green vortex (TGV) case**

Taylor-Green vortex problem is an idea test for a numerical scheme whether it has the ability to simulate the transition from laminar to turbulence, since it contains the stages of laminar, transition, and finally turbulence as time develops. The three-dimensional Taylor-Green vortex was first introduced by Taylor and Green, and its initial condition is given as follows:

$$\begin{aligned} u &= \mathcal{U}\_0 \sin x \cos y \sin z \\ v &= -\mathcal{U}\_0 \cos x \sin y \sin z \\ w &= 0 \\ p &= p\_o + \rho\_o \mathcal{U}\_0^{-2} (\cos 2z + 2)(\cos 2x + \cos 2y) / 16 \end{aligned} \tag{43}$$

The computing domain is 0, 2*π* × 0, 2*π* × 0, 2*π* . The parameters are set as *p*<sup>0</sup> =100, *ρ*<sup>0</sup> =1, *U*<sup>0</sup> =1, then the Mach number is about 0.08 and the flow is nearly incompressible. The boundary conditions in the three dimensions are all periodic. Three numerical schemes are used here, WENO [7], WENO-compact upwind [9] (denoted as WENO-CU), and WENO-linear upwind [10] (denoted as WENO-LU), respectively. Two sets of computation grid are taken as 643 and 323 for different flow Reynolds numbers.

Because of the stiffness of the equation due to the reactions, the time integration of the equations is performed by means of the point-implicit method [27]. For instance, *U* is supposed as the unknown variable at time step *n*+1, and the below equation is solved to obtain *U*n+1=*U*<sup>n</sup>

Δ Δ Δ

where *C* is the convection term, *V* is the viscous term, *S* is the source term, and *I* is the unity

In order to obtain the physical quantities of the fluid at a droplet position, the fluid velocities computed on the Eulerian mesh were interpolated at the droplet locations. A fourth-order Lagrangian interpolation was used for this purpose [28]. When the droplet is located at (*x*d, *y*d,

> , *yl* , *zl*

=- =- =- = - = - = -

*nnn nnn*

*xx yy zz ux y z <sup>u</sup>*

with the third-order Adams-Bash-forth scheme in direction, *i*, for droplet position,

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

1 1 2 ,, , , , 23 16 5 <sup>Δ</sup> 12 12 12

*nn n n n*

ååå ÕÕÕ /2 /2 /2 / 2 / 2 / 2

d dd , , 1 /2 1 /2 1 /2 1 /2 1 /2 1 /2

ç ÷ --- <sup>=</sup> ´ ---

*i nj nk n ln ln ln il i l il*

In the Lagrangian droplet-tracking method, the droplet equation of motion was integrated

The equations of velocity and temperature of the droplet are solved using the same integration

Taylor-Green vortex problem is an idea test for a numerical scheme whether it has the ability to simulate the transition from laminar to turbulence, since it contains the stages of laminar, transition, and finally turbulence as time develops. The three-dimensional Taylor-Green vortex was first introduced by Taylor and Green, and its initial condition is given as follows:

è ø

( ) ( ) ( )

*l i l j l k*

æ ö

ddd

*di di di di di x x tu u u* (42)

*l ll*

è ø

¹¹¹

*<sup>U</sup>* (40)

), the fluid velocity at the droplet position

*xx yy zz* (41)

*ijk*

( ) æ ö æ ö ¶ ç ÷ - ç ÷ =- + +

*<sup>S</sup> n nn I tU CV S t*

è ø ¶ è ø

*z*d) in a cell whose coordinates are (*x*0, *y*0, *z*0)…(*xl*

38 Modeling and Simulation in Engineering Sciences

( )

method, as Eq. (42).

**4. Numerical examples**

**4.1. Taylor-Green vortex (TGV) case**

, ,

*n*

is obtained via the Lagrangian interpolation, which is expressed as

where the subscripts *i, j, k*, and *l* represents the cell indexes and *n* = 4.

+Δ*U*

matrix.

**Figure 1** shows the results of the total kinetic energy integrated in the entire flow field. For Re = 30,000, the kinetic energy decays very slowly at the beginning and it has almost no change in the first few seconds. However, when the time is at about 4 s, the sub-grid scales are produced and the kinetic energy begins to decay due to the numerical dissipation, which also can be taken as the implicit sub-grid dissipation. The results of Re number equalling 3000 are similar, since in both cases, the numerical viscosity is larger compared to the physical viscosity. For the results of Re = 300, there is only slight difference between the curve of WENO-LU and WENO-CU, since in this case the physical viscosity is dominated.

**Figure 1.** Comparison of the total kinetic energy. (a) Re = 30,000; (b) Re = 3000; (c) Re = 300.

**Figure 2** shows the energy spectrum of the TGV with Re = 30,000. As the time increases, the slope of energy spectrum also changes. If the time is longer than 4, the turbulence develops and the kinetic energy builds up a Kolmogorov inertial range. The slope is then approximately decayed with the −5/3 law. It implies that it can predict the characteristics of turbulence.

**Figure 2.** Energy spectrum of TGV with Re = 30,000.

Since the passive scalar transportation is very important in many fields, especially in the chemical-reacting process, we also solved the conserved passive scalar transport equation

$$\frac{\partial \left(\rho f\right)}{\partial t} + \nabla \cdot \left(\rho f \vec{L}\right) = \nabla \cdot \left(\rho \mathcal{D} \nabla f\right) \tag{44}$$

Here, only the WENO-LU scheme is used since it has less dissipation and costs less computing time. Equation (44) presents the initial condition for *f* and it is a function of the initial *Q* distribution, where *Q* is the value of *q*-criteria. Different Schmidt numbers are taken in the simulations, 1, 0.1, and 0.01, corresponding to different diffusivity coefficients, respected *D* =*μ* / (*ρ*Sc)

$$f(t=0) = \frac{1 + err\{C\_0Q\}}{2}, C\_0 = 100\tag{45}$$

Here, erf is the error function. **Figure 3** shows the initial distribution of *f*. As we can see, in most regions, *f* equals 1 or 0 and there are large gradients of *f* initially.

**Figure 3.** Initial distribution of passive scalar *f*.

**Figure 2** shows the energy spectrum of the TGV with Re = 30,000. As the time increases, the slope of energy spectrum also changes. If the time is longer than 4, the turbulence develops and the kinetic energy builds up a Kolmogorov inertial range. The slope is then approximately decayed with the −5/3 law. It implies that it can predict the characteristics of turbulence.

Since the passive scalar transportation is very important in many fields, especially in the chemical-reacting process, we also solved the conserved passive scalar transport equation

(44)

¶( ) +Ñ× =Ñ× Ñ

<sup>+</sup> = = <sup>=</sup> <sup>0</sup>

1 () ( 0) , 100 <sup>2</sup>

0

*erf C Q f t <sup>C</sup>* (45)

<sup>r</sup> () ( ) *ρf ρfU ρD f <sup>t</sup>*

Here, only the WENO-LU scheme is used since it has less dissipation and costs less computing time. Equation (44) presents the initial condition for *f* and it is a function of the initial *Q* distribution, where *Q* is the value of *q*-criteria. Different Schmidt numbers are taken in the simulations, 1, 0.1, and 0.01, corresponding to different diffusivity coefficients, respected

¶

**Figure 2.** Energy spectrum of TGV with Re = 30,000.

40 Modeling and Simulation in Engineering Sciences

*D* =*μ* / (*ρ*Sc)

**Figure 4** shows the variance of the passive scalar for the Re = 30,000 flow in different Schmidt numbers and different grid levels. Similar to the kinetic energy decaying, the scalar variance decays very slowly in the first 4 s, and it can be concluded that the passive scalar mixing goes on with the laminar diffusion. However, the significant decaying can be observed due to the initial large gradients of passive scalar. For the time longer than four, the turbulent mixing occurs and the decaying of passive scalar variance is larger than that in the beginning stage.

**Figure 4.** Variance curves of passive scalar.

**Figure 5** shows the passive scalar contours for Re = 30,000 and Sc = 1.0 in the plane of *z* = π. During the stage of laminar diffusion, the interface between the high and low *f* region is clear while the interface is wrinkled in the stage of turbulent mixing, which enhances the mixing process.

**Figure 5.** Passive scalar field in the plane of *z* = π. (a) *t* = 0 s; (b) *t* = 3 s; (c) *t* = 6 s; (d) *t* = 10 s.

#### **4.2. Reactive H2/air spatially developing mixing layer**

The schematic of the supersonic mixing layer is shown in **Figure 6**.

**Figure 6.** Schematic of computation model of spatially developing mixing layer.

The inlet conditions are set as follows. The average velocity profile in the stream-wise direction is a hyperbolic tangent, and the average velocity in other directions is set to zero

$$\overline{\mathcal{U}} = \frac{\mathcal{U}\_A + \mathcal{U}\_F}{2} + \frac{\mathcal{U}\_A - \mathcal{U}\_F}{2} \tanh(\frac{y - 0.5l\_y}{2\delta\_0}) \tag{46}$$

where *δ*0 is the inlet layer thickness given as 0.4 mm, and *l*y is the length of computational domain in the *y* direction. *U*<sup>A</sup> is the hot-air inflow velocity, 2000 m/s, *U*<sup>F</sup> is the cold fuel inflow velocity, 1000 m/s. The temperatures of both streams are *T*A = 1600 K and *T*F = 390 K. The reactive system pressure is specified as *p*A = 1 atm.

The stream-wise and transverse domain lengths are *lx* and *ly*, respectively. By setting *lx*= 4*ly*, the transverse domain size is large enough to have minimal influence on the main flow region of the droplet-laden shear layer. The computation grid in the two-dimensional case is taken as 512 × 128 equally spaced computational cells in the stream-wise and the transverse directions, respectively, after a series of grid independence tests.

In order to stimulate the flow instability, the inlet velocity perturbation is given as

while the interface is wrinkled in the stage of turbulent mixing, which enhances the mixing

**Figure 5.** Passive scalar field in the plane of *z* = π. (a) *t* = 0 s; (b) *t* = 3 s; (c) *t* = 6 s; (d) *t* = 10 s.

The schematic of the supersonic mixing layer is shown in **Figure 6**.

**Figure 6.** Schematic of computation model of spatially developing mixing layer.

system pressure is specified as *p*A = 1 atm.

The inlet conditions are set as follows. The average velocity profile in the stream-wise direction

tanh( ) 22 2

where *δ*0 is the inlet layer thickness given as 0.4 mm, and *l*y is the length of computational domain in the *y* direction. *U*<sup>A</sup> is the hot-air inflow velocity, 2000 m/s, *U*<sup>F</sup> is the cold fuel inflow velocity, 1000 m/s. The temperatures of both streams are *T*A = 1600 K and *T*F = 390 K. The reactive

0 0.5

*<sup>δ</sup>* (46)

is a hyperbolic tangent, and the average velocity in other directions is set to zero

+ - - = +

*AF AF <sup>y</sup> UUUU y l <sup>U</sup>*

**4.2. Reactive H2/air spatially developing mixing layer**

process.

42 Modeling and Simulation in Engineering Sciences

$$
\sigma' = (\mathcal{U}\_A - \mathcal{U}\_F) G(y - 0.5 \mathcal{l}\_y) \; A \sin(2\pi f\_\circ t + \xi) \tag{47}
$$

where *G* is the Gauss function, *A* is the fluctuation amplitude, *f*o is the maximum disturbance frequency, and *ξ* is the random phase. The nonreflecting conditions are used in the transverse directions. The outflow is treated by an extrapolation of each primitive variable. A 19-step detailed H2/O2 reaction mechanism is used, which can be referred in [29].

**Figure 7.** Instantaneous contour distribution of reactive flows. (a) temperature; (b) H2; (c) H2O.

Instantaneous distributions of temperature and mass fraction of H2 and H2O at *t* = 2 ms are shown in **Figure 7**. At about *x* = 0.4 m, the mixing layer loses stability, the large-scale vortices begin to form, and efficient mixing occurs inside the vortices. Due to the relatively high temperature of air, fast reaction occurs and a large amount of water concentrates inside the vortices.

**Figure 8** shows the flame structures in the mixing layer flow. The flame structures were obtained by means of plotting all the combustion status in the whole computation domain within the time interval 2<*t*<2.5 ms, that is, recording and plotting the temperature or species fraction at each calculation node in the observing time interval. Here, the mixture fraction is defined as

$$Z \equiv \frac{\mathcal{W}\_{\rm H}}{\mathcal{W}\_{\rm H} + \mathcal{W}\_{\rm O}} \tag{48}$$

where WH and WO are the mass fraction of element H and element O, respectively. For the stoichiometric ratio, *Z*st = 0.111. It can be seen that the peaks of temperature and water are both near *Z*st. The flame structure is similar with that in the flamelet model, since the non-premixed flame dominates in the current case.

**Figure 8.** Flame structure. (a) temperature distribution and (b) water distribution.

#### **4.3. Non-premixed combustion of** *n***-decane droplets in turbulent mixing layers**

The considered flow configuration is that of a spatially developing shear layer, as depicted in **Figure 9**, which is formed between a stream of hot air moving at velocity *U*<sup>A</sup> and a spray carried by cold air moving at velocity *U*S. *x* and *y* refer to the stream-wise and transverse directions, respectively. The inflow parameters of the two streams are listed in **Table 1**.

**Figure 9.** Schematic of a fuel spray in a two-dimensional turbulent shear layer.


**Table 1.** Inflow parameters.

Instantaneous distributions of temperature and mass fraction of H2 and H2O at *t* = 2 ms are shown in **Figure 7**. At about *x* = 0.4 m, the mixing layer loses stability, the large-scale vortices begin to form, and efficient mixing occurs inside the vortices. Due to the relatively high temperature of air, fast reaction occurs and a large amount of water concentrates inside the

**Figure 8** shows the flame structures in the mixing layer flow. The flame structures were obtained by means of plotting all the combustion status in the whole computation domain within the time interval 2<*t*<2.5 ms, that is, recording and plotting the temperature or species fraction at each calculation node in the observing time interval. Here, the mixture fraction is

where WH and WO are the mass fraction of element H and element O, respectively. For the stoichiometric ratio, *Z*st = 0.111. It can be seen that the peaks of temperature and water are both near *Z*st. The flame structure is similar with that in the flamelet model, since the non-premixed

*W W* (48)

<sup>º</sup> <sup>+</sup> H H O

*<sup>W</sup> <sup>Z</sup>*

vortices.

defined as

flame dominates in the current case.

44 Modeling and Simulation in Engineering Sciences

**Figure 8.** Flame structure. (a) temperature distribution and (b) water distribution.

**4.3. Non-premixed combustion of** *n***-decane droplets in turbulent mixing layers**

The considered flow configuration is that of a spatially developing shear layer, as depicted in **Figure 9**, which is formed between a stream of hot air moving at velocity *U*<sup>A</sup> and a spray carried

The initial Mach numbers of the two streams equal 1.2. The cold-air stream is nonuniformly laden with droplets at the cold-air inlet. For the combustion reaction model, a one-step global reaction model of *n*-decane is used [30]. In this model, the chemical reaction is simplified as

$$\rm{C}\_{10}\rm{H}\_{22} + \rm{\eta}\_{\rm{O}\_2}\rm{O}\_2 \Rightarrow \rm{\eta}\_{\rm{CO}\_2}\rm{CO}\_2 + \rm{\eta}\_{\rm{H}\_2\rm{O}}\rm{H}\_2\rm{O} \tag{49}$$

where *ni* is the molar stoichiometric coefficients of each species. The global reaction constant, *k* (mol/cm3 /s), is given by

$$k = A \cdot T^{\
u} \exp\left(-\frac{E}{RT}\right) \left(\frac{\rho \cdot Y\_{\text{p}}}{W\_{\text{F}} \times 10^{6}}\right)^{\circ} \left(\frac{\rho \cdot Y\_{\text{O}\_{2}}}{W\_{\text{O}\_{2}} \times 10^{6}}\right)^{\circ} \tag{50}$$

where *A* is the coefficient of frequency and *E* is the activation energy. In the reaction of *n*decane, these parameters are given as follows: *A* = 3.8 × 1011, *E* = 1.256 × 105 J/mol, *a* = 0.25, *b* = 1.5, and *n* = 0. Therefore, the combustion reaction rate per unit volume (g/cm3 /s) is *R*F = *W*V*k* × 103 .

In all simulation cases, the fuel droplets are initially located in the cold air and cannot vaporize until they traverse into the hot air associated with the rolling-up process of vortices. The hotair stream provides the heat needed for droplet vaporization. The fuel vapor diffuses into the air stream, and it begins to react with the local oxygen as it reaches the high-temperature regions.

Large eddies rotate and entrain the hot air, providing heat for the evaporation of fuel droplets that preferentially accumulate in the high-strain vortex-braid regions adjoining the hot air. The fuel vapor generated by droplet evaporation diffuses toward the hot air and mixes with the surrounding oxygen, forming reactive pockets. Auto-ignition occurs when the fuel mass fraction is high enough, as depicted in **Figure 10(a)**. Under the instantaneous condition, the time *t*ig defines the moment that the incipient establishment of ignition kernels occurs. We normalize the gaseous temperature by *T*\* = *T*/*T*in and *T*in = (*T*A + *T*S)/2. The droplet diameters are normalized by *d*\* d = *d*d/dd0 and dd0 is the initial droplet diameter.

Large eddies rotate and the reaction kernels between two convection vortices are strained at time *t* = *t*ig + 0.25*t*L, as shown in **Figure 10(b)**. Here, the time *t*L is the large eddy turnover time and is defined as

$$t\_{\rm L} = \frac{\delta\_{\rm L}}{\left| \mathcal{U}\_{\rm A} - \mathcal{U}\_{\rm s} \right|} \tag{51}$$

where *δ*<sup>L</sup> is the length of the large eddy. Compared with the reaction rates of the ignition kernels at time *t*ig, the chemical heat release in the reaction kernels at time *t*ig + 0.25*t*L results in a larger temperature increase and facilitates the self-acceleration of the chemical reaction rate. At time *t* = *t*ig + 0.5*t*L, the reaction kernels are separated into two parts due to the stretch effects, as shown by **Figure 10(c)**. One part of the reaction kernels still exists in the high-strain vortexbraid zone and suffers the shear strain continuously. Evaporating fuel droplets segregate in a thin layer along this reaction kernel and feed fuel vapors for the chemical reaction. The other separated part is entrained in the inner of the large cold eddy. It is difficult forthe fuel droplets to penetrate into the high-vorticity core of the large eddy because of the large droplet inertia and fuel droplets surround far from this separated reaction kernel, forming a thick vaporiza‐

tion layer. However, the turbulence motion separates the vaporization layer and the flame kernel. Hence, the fuel vapor diffusion is difficult forfeeding the chemicalreaction in the flame kernel effectively. Here, we call it the weak-fueling mechanism due to the preferential concentration effects of droplets. If we track the separated reaction kernel in the inner of the cold eddy, it is found that the local chemical reaction rate of this flame kernel is decelerated, compared to that ofthe reaction kernel in the high-strain zone accompanied by the rich reactant supplement. **Figure 10(d)** depicts the extinction phenomena at time *t* = *t*ig + 0.75*t*L. Compari‐

where *ni*

regions.

and is defined as

*k* (mol/cm3

/s), is given by

46 Modeling and Simulation in Engineering Sciences

normalize the gaseous temperature by *T*\*

normalized by *d*\* d = *d*d/dd0 and dd0 is the initial droplet diameter.

is the molar stoichiometric coefficients of each species. The global reaction constant,

F O

2

*RT W W* (50)

= *T*/*T*in and *T*in = (*T*A + *T*S)/2. The droplet diameters are

*U U* (51)

J/mol, *a* = 0.25, *b* =

/s) is *R*F = *W*V*k* × 103

.

2

*b a*

F O 6 6

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

*<sup>n</sup> <sup>E</sup> ρ Y ρ Y k AT*

decane, these parameters are given as follows: *A* = 3.8 × 1011, *E* = 1.256 × 105

1.5, and *n* = 0. Therefore, the combustion reaction rate per unit volume (g/cm3

exp 10 10

where *A* is the coefficient of frequency and *E* is the activation energy. In the reaction of *n*-

In all simulation cases, the fuel droplets are initially located in the cold air and cannot vaporize until they traverse into the hot air associated with the rolling-up process of vortices. The hotair stream provides the heat needed for droplet vaporization. The fuel vapor diffuses into the air stream, and it begins to react with the local oxygen as it reaches the high-temperature

Large eddies rotate and entrain the hot air, providing heat for the evaporation of fuel droplets that preferentially accumulate in the high-strain vortex-braid regions adjoining the hot air. The fuel vapor generated by droplet evaporation diffuses toward the hot air and mixes with the surrounding oxygen, forming reactive pockets. Auto-ignition occurs when the fuel mass fraction is high enough, as depicted in **Figure 10(a)**. Under the instantaneous condition, the time *t*ig defines the moment that the incipient establishment of ignition kernels occurs. We

Large eddies rotate and the reaction kernels between two convection vortices are strained at time *t* = *t*ig + 0.25*t*L, as shown in **Figure 10(b)**. Here, the time *t*L is the large eddy turnover time

> <sup>=</sup> - L

A S *δ*

where *δ*<sup>L</sup> is the length of the large eddy. Compared with the reaction rates of the ignition kernels at time *t*ig, the chemical heat release in the reaction kernels at time *t*ig + 0.25*t*L results in a larger temperature increase and facilitates the self-acceleration of the chemical reaction rate. At time *t* = *t*ig + 0.5*t*L, the reaction kernels are separated into two parts due to the stretch effects, as shown by **Figure 10(c)**. One part of the reaction kernels still exists in the high-strain vortexbraid zone and suffers the shear strain continuously. Evaporating fuel droplets segregate in a thin layer along this reaction kernel and feed fuel vapors for the chemical reaction. The other separated part is entrained in the inner of the large cold eddy. It is difficult forthe fuel droplets to penetrate into the high-vorticity core of the large eddy because of the large droplet inertia and fuel droplets surround far from this separated reaction kernel, forming a thick vaporiza‐

L

*t*

**Figure 10.** Flame kernels at (a) time *t* = *t*ig; (b) time *t* = *t*ig + 0.25 *t*L; (c) time *t* = *t*ig + 0.5 *t*L; and (d) time *t* = *t*ig + 0.75 *t*L. Instantaneous distributions of dimensionless gaseous temperature (left) and reaction rate (right), with red isolines giv‐ en by *Ω* = [11, 15] (kgm−3s−1) in (a) and *Ω* = [15, 24] (kgm−3s−1) in (b)–(d). Here, the dots indicate fuel droplets in the figures, with colors corresponding to dimensionless droplet diameters in the right figures.

son with **Figure 10(c)** shows that the reaction kernels that have high chemical reaction rates in the high-strain region at time *t* = *t*ig + 0.5*t*L disappear. The other reaction kernel entrained in large eddies is distorted due to the turbulence motions.
