## 1. Introduction

Numerical simulation has become one of important approaches in academia and industry to research-complicated reactive flows in recent years, though it is still a fast developing subject. For combustion in low-speed flows, usually the low-Mach-number assumption is used, which assumes that the change of temperature is unrelated to the hydrodynamic pressure [1,2]. While for transonic or supersonic reactive flows, such as in the scramjet combustor, the compressible effects could not be ignored and the set of fully compressible Navier-Stokes equations should be solved.

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 is the high cost of calculation.

In the early days, RANS or LES was used to simulate the compressible reactive flows considering 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, H_{2}/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 chemistry under varying initial pressure.

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 shock-turbulence interaction. Furthermore, Ren et al. [14] proposed a Roe type, characteristic-wise compact-WENO scheme. Considering overcoming the complexity and program parallelization 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 methods used for high-speed flows can be seen in [16].

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 flows. This chapter will also consider the spray combustion in mixing layer flows as a computation example.

This chapter tries to contribute the numerical simulation in the direction of compressible reactive flows. It is organized as follows. In section 2, the governing equations are presented; the methods to calculate the multispecies physical parameters and to deal with the chemical stiffness are also presented. In section 3, the discretization schemes and time integration methods are discussed, including the discretization of inviscid/viscous terms. In section 4, three examples of numerical simulations of compressible flows are presented. In section 5, the conclusions are drawn finally.

## 2. Governing equations and models of gas properties

### 2.1. Governing equations

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, compressible 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

Here, *ρ* is the gas-phase density, *u*_{i} is the velocity, *T* is the temperature, *P* is the pressure, *Y*_{k }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 *S*_{m}, *S*_{F,} and *S*_{Q} describe the two-phase couplings of mass, momentum, and energy, respectively. *S*_{combustion,k} is the source term due to combustion. *τ*_{ij} is the Newtonian viscous stress tensor

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

Here, *N*_{S} is the total number of species. *W*_{k} and *X*_{k} are the molecule weight and the mole fraction of the *k*^{th} species, respectively. *e*_{t} is the total energy, that is, kinetic energy plus internal (containing chemical) energy, which is defined as follows:

where *c*_{p,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 *c*_{p,k }are represented by a fifth-order polynomial, and the coefficients can be found in [25]

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

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]

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

where

where

### 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,i}, velocity, *u*_{d,i}, temperature, *T*_{d}, and mass, *m*_{d}, of a single droplet are given by

where *c*_{p} is the specific heat of mixture gas, *c*_{L} is the specific heat of the liquid. The momentum response time, *τ*_{d}, is defined as

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

respectively. The Nusselt and Sherwood numbers are

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

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

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

where the nondimensional evaporation parameter is given by

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

Here, *Y*_{V} 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

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*_{V} is the latent heat of fuel vapor at liquid temperature, *T*_{d}, and is given as follows:

where *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},

where *m*_{d} and *m*_{d}/d*t*, respectively, *u*_{d,i} 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*_{d} 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,k}, is expressed using the combustion reaction rate per unit volume, *R*_{F}, as follows:

where *n*_{k} 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. *W*_{k} and *W*_{F} are the molecular weights of the *k*^{th} species and fuel, respectively.

## 3. Numerical discretization procedures

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 *n* to iteration *n*+1 is performed as

where *U*_{j} 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 *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

which is constructed by hybridizing a fifth-order Compact Upwind (CU) scheme, *r*_{c}.

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

where *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*^{n}+Δ*U*

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})…(*x*_{l}, *y*_{l}, *z*_{l}), the fluid velocity at the droplet position is obtained via the Lagrangian interpolation, which is expressed as

(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,

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:

The computing domain is ^{3} and 32^{3} for different flow Reynolds numbers.

**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 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

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

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 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 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.

### 4.2. Reactive H_{2}/air spatially developing mixing layer

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

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

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*_{A} is the hot-air inflow velocity, 2000 m/s, *U*_{F} 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 *l*_{x} and *l*_{y}, respectively. By setting *l*_{x }= 4*l*_{y}, 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

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 H_{2}/O_{2} reaction mechanism is used, which can be referred in [29].

Instantaneous distributions of temperature and mass fraction of H_{2} and H_{2}O 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

where W_{H} and W_{O} 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.

### 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*_{A} 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**.

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

where *n*_{i} is the molar stoichiometric coefficients of each species. The global reaction constant, *k* (mol/cm^{3}/s), is given by

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 × 10^{11}, *E* = 1.256 × 10^{5}J/mol, *a* = 0.25, *b* = 1.5, and *n* = 0. Therefore, the combustion reaction rate per unit volume (g/cm^{3}/s) is *R*_{F }= *W*_{V}*k* × 10^{3}.

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 hot-air 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}/d_{d0} and d_{d0} 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

where *δ*_{L} 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 vortex-braid 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 for the 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 vaporization layer. However, the turbulence motion separates the vaporization layer and the flame kernel. Hence, the fuel vapor diffusion is difficult for feeding the chemical reaction 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 of the 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}. Comparison 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.

## 5. Conclusions

This chapter introduced the state of art of the numerical simulation on compressible flows. The high-order hybrid WENO scheme, as one of the important issues for simulating the supersonic flows, is required to compute the numerical fluxes, in which the WENO scheme is used near the discontinuities and linear scheme is used in the smooth regions of the flows. In the solver developed by the present authors, the point-implicit method is utilized to treat the chemical stiffness problem.

A TGV problem, defined to investigate the transports of passive scalar, is simulated by the present numerical procedures. The decaying of the kinetic energy and variance of passive scalar due to initial gradients are analyzed based on the numerical results. The spatially developing mixing layer flows are simulated, taking the hydrogen and liquid *n*-decane droplets as fuel and air as oxidizer in the scenario of supersonic inflows. The stable combustion of hydrogen/air and the flame structure are analyzed by the numerical results. The ignition considering the eddy turnover effects is carefully analyzed, since the large-scale structures and droplet dispersion can be captured by the present numerical methods.

This chapter tries to contribute the field of numerically analyzing compressible reactive flows. The new high-order schemes and numerical strategy will be developed in the future, so that the stabilities and robustness will be enhanced by the solver based on the Navier-Stokes equations. Thus, deeply understanding the flow physics of compressible reactive flows will be benefited from results obtained by numerical simulations.