Numerical Simulation of Compressible Reactive Flows

Abstract

Numerical simulation has been widely employed to investigate the compressible flows since it is difficult to carry out the experimental measurements, especially in the reactive flows. The shock-wave capturing scheme will be necessary for resolving the compressible flows, and moreover the careful treatments of chemical reaction should be considered for proceeding numerical simulation stable and fast. Recently, the present authors have tried to establish a high-resolution numerical solver for computing the compressible reactive flows. This chapter presents the numerical procedures acquired in this solver for computing the fluxes using weighted essentially non-oscillatory (WENO) scheme, dealing with chemical stiffness problems, and tracing droplets and their interaction with the compressible fluids. As examples, the Taylor-Green vortex (TGV) problem considering the passive scalar mixing, the spatially developing reactive mixing layer flows taken gas-phase hydrogen, and liquid n-decane as fuel are simulated, respectively. Many important characteristics of flow, flame, and ignition are analyzed under the supersonic condition.

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₂/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

$\frac{\partial \rho }{\partial t}+\frac{\partial }{\partial x_j}\left(\rho u_j\right)=S_m$E1

$\frac{\partial }{\partial t}\left(\rho u_i\right)+\frac{\partial }{\partial x_j}\left(\rho u_i{u}_j+P{\delta }_{ij}-{\tau }_{ij}\right)=S_F$E2

$\frac{\partial }{\partial t}\left(\rho e_{\text{t}}\right)+\frac{\partial }{\partial x_j}\left(\left(\rho e_{\text{t}}+P\right)u_j-\lambda \frac{\partial T}{\partial x_j}-u_i{\tau }_{ij}\right)=S_Q$E3

$\frac{\partial }{\partial t}\left(\rho Y_k\right)+\frac{\partial }{\partial x_j}\left(\rho Y_k\left(u_j+V_j^{\text{c}}\right)-\rho D_k\frac{\partial Y_k}{\partial x_j}\right)=S_{\text{combustion},k}+S_{Y_k}$E4

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

${\tau }_{ij}=\mu \left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}-\frac{2}{3}\frac{\partial u_k}{\partial x_k}{\delta }_{ij}\right)$E5

A correction velocity Vc jis added to the convection velocity in the species equations to ensure global mass conservation

$V_j^{\text{c}}={\displaystyle {\displaystyle \sum _{k=1}^{N_s}}{D}_k\frac{W_k}{W}}\frac{\partial X_k}{\partial x_i}$E6

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:

$e_{\text{t}}={\displaystyle {\displaystyle \sum _{k=1}^{N_S}}{Y}_k\left({\displaystyle {\int }_{T_{\text{ref}}}^T{c}_{p,k}\text{d}T+h_{f,k}^0}\right)}-\frac{P}{\rho }+\frac{u_i{u}_i}{2}$E7

where c_p,k and h0 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]

$\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$E8

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

$P=\rho RT{\displaystyle {\displaystyle \sum _{k=1}^{N_S}}\frac{Y_k}{W_k}}$E9

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}}$E10

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

$D_{AB}=\frac{0.00266T^{3/2}}{P{W}_{AB}^{1/2}{\sigma }_{AB}^2{\Omega }_D}$E12

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

${\mu }_m={\displaystyle {\displaystyle \sum _{i=1}^n}\frac{X_i{\mu }_i}{{\displaystyle {\displaystyle \sum _{j=1}^n}{X}_j{\phi }_{ij}}}}$E13

where ${\phi }_{ij}=\frac{{\left(1+{\left({\mu }_i/{\mu }_j\right)}^{1/2}{\left(W_j/W_i\right)}^{1/4}\right)}^2}{{\left(8\left(1+W_i/W_j\right)\right)}^{1/2}}$.

$\lambda ={\displaystyle {\displaystyle \sum _{i=1}^n}\frac{X_i{\lambda }_i}{{\displaystyle {\displaystyle \sum _{j=1}^n}{X}_j{A}_{ij}}}}$E14

where $A_{ij}=\frac{{\left(1+{\left({\lambda }_i/{\lambda }_j\right)}^{1/2}{\left(W_j/W_i\right)}^{1/4}\right)}^2}{{\left(8\left(1+W_i/W_j\right)\right)}^{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,i, velocity, u_d,i, temperature, T_d, and mass, m_d, of a single droplet are given by

$\frac{d{x}_{\text{d},i}}{dt}=u_{\text{d},i}$E15

$\frac{\text{d}{u}_{\text{d},i}}{\text{d}t}=\frac{F_{\text{d}}}{m_{\text{d}}}=\left(\frac{f_1}{{\tau }_{\text{d}}}\right)\left(u_i-u_{\text{d},i}\right)$E16

$\frac{\text{d}{T}_d}{\text{d}t}=\frac{Q_{\text{d}}+{{\displaystyle \dot{m}}}_{\text{d}}{L}_{\text{V}}}{m_{\text{d}}{c}_{\text{L}}}=\left(\frac{f_2}{{\tau }_{\text{d}}}\right)\left(\frac{\text{Nu}}{3\mathrm{Pr}}\right)\left(\frac{c_p}{c_{\text{L}}}\right)\left(T-T_{\text{d}}\right)+\left(\frac{{{\displaystyle \dot{m}}}_{\text{d}}}{m_{\text{d}}}\right)\frac{L_{\text{V}}}{c_{\text{L}}}$E17

$\frac{\text{d}{m}_{\text{d}}}{\text{d}t}={{\displaystyle \dot{m}}}_{\text{d}}=-m_{\text{d}}\left(\frac{1}{{\tau }_{\text{d}}}\right)\left(\frac{\text{Sh}}{\text{3Sc}}\right)\text{ln}\left(1+B_{\text{M}}\right)$E18

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

${\tau }_{\text{d}}=\frac{{\rho }_{\text{d}}{d}_{\text{d}}^2}{18\mu }$E19

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

$\mathrm{Pr}=\frac{\mu c_p}{\lambda },\text{Sc}=\frac{\mu }{\rho D}$E20

respectively. The Nusselt and Sherwood numbers are

$\text{Nu}=2+0.552{\mathrm{Re}}_{\text{d}}^{1/2}{\mathrm{Pr}}^{1/3}$E21

$\text{Sh}=2+0.552{\mathrm{Re}}_{\text{d}}^{1/2}{\text{Sc}}^{1/3}$E22

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

${\mathrm{Re}}_{\text{d}}=\frac{\rho \left|u_i-u_{d,i}\right|d_{\text{d}}}{\mu }$E23

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

$f_1=1+0.15{\mathrm{Re}}_{\text{d}}^{0.687}$E24

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

$f_2=\frac{\beta }{e^{\beta }-1}$E25

where the nondimensional evaporation parameter is given by

$\beta =-1.5\mathrm{Pr}{\tau }_{\text{d}}\left(\frac{{{\displaystyle \dot{m}}}_{\text{d}}}{m_{\text{d}}}\right)$E26

The evaporation rate is controlled by the mass transfer number, B_M , and is given by

$B_{\text{M}}=\frac{Y_{\text{sf}}-Y_V}{1-Y_{\text{sf}}}$E27

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

$Y_{\text{sf}}=\frac{{\chi }_{\text{sf}}}{{\chi }_{\text{sf}}+\left(1-{\chi }_{\text{sf}}\right)W/W_{\text{V}}}$E28

${\chi }_{\text{sf}}=\frac{P_{\text{atm}}}{P}\exp \left(\frac{L_{\text{V}}}{R/W_{\text{V}}}\left(\frac{1}{T_{\text{B},\text{L}}}-\frac{1}{T_{\text{d}}}\right)\right)$E29

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:

$L_{\text{V}}=L_{\text{V},T_{\text{B},\text{L}}}\left(\frac{T_{\text{C},\text{L}}-T_{\text{d}}}{T_{\text{C},\text{L}}-T_{\text{B},\text{L}}}\right)$E30

where $L_{\text{V},T_{\text{B},\text{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,

$S_m=-\frac{1}{\Delta V}{\displaystyle {\displaystyle \sum _{N_{\text{c}}}}\left({{\displaystyle \dot{m}}}_{\text{d}}\right)}$E31

$S_F=-\frac{1}{\Delta V}{\displaystyle {\displaystyle \sum _{N_{\text{c}}}}\left(F_{\text{d}}+{{\displaystyle \dot{m}}}_{\text{d}}{u}_{\text{d},}{}_i\right)}$E32

$S_Q=-\frac{1}{\Delta V}{\displaystyle {\displaystyle \sum _{N_{\text{c}}}}\left(Q_{\text{d}}+{{\displaystyle \dot{m}}}_{\text{d}}(\frac{u_{\text{d},i}{u}_{\text{d},i}}{2}+h_{\text{V},\text{sf}})\right)}$E33

$S_{Y_k}=\{\begin{array}{cc}--\frac{1}{\Delta V}{\displaystyle {\displaystyle \sum _N}\left({{\displaystyle \dot{m}}}_{\text{d}}\right)}& \text{for}\text{fuel}\\ 0& \text{for}\text{other}\text{species}\end{array}$E34

where m_d and ${{\displaystyle \dot{m}}}_{\text{d}}$ denote the droplet mass and dm_d/dt, 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:

$S_{\text{combustion},k}=-\frac{n_k}{n_{\text{F}}}\frac{W_k}{W_{\text{F}}}{R}_{\text{F}}$E35

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 $\text{d}{U}_j=L_j(U)\text{d}t$, the iteration from n to iteration n+1 is performed as

$\begin{array}{l}{U}_j^{\left(1\right)}=U_j^n+\Delta t{L}_j\left(U^n\right)\\ U_j^{\left(2\right)}=\frac{3}{4}{U}_j^n+\frac{1}{4}{U}_j^{\left(1\right)}+\frac{1}{4}\Delta t{L}_j\left(U_j^{\left(1\right)}\right)\\ U_j^{n+1}=\frac{1}{3}{U}_j^n+\frac{2}{3}{U}_j^{\left(2\right)}+\frac{2}{3}\Delta t{L}_j\left(U_j^{\left(2\right)}\right)\end{array}$E36

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 ${{\displaystyle \widehat{f}}}_{j+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

${{\displaystyle \widehat{f}}}_{j+1/2}={\sigma }_{j+1/2}{{\displaystyle \widehat{f}}}_{j+1/2}^{\text{CU}}+\left(1+{\sigma }_{j+1/2}\right){{\displaystyle \widehat{f}}}_{j+1/2}^{\text{WENO}}$E37

which is constructed by hybridizing a fifth-order Compact Upwind (CU) scheme, ${{\displaystyle \widehat{f}}}_{j+1/2}^{\text{CU}}$,and a fifth-order WENO one, ${{\displaystyle \widehat{f}}}_{j+1/2}^{\text{WENO}}$, through a smoothness indicator r_c. ${\sigma }_{j+1/2}$ is the weight of the numerical flux,

${\sigma }_{j+1/2}=\min \left(1,\frac{r_{j+1/2}}{r_c}\right)$E38

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

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

where ${{\displaystyle \widehat{F}}}_j$ is the difference approximation for ${\left(\frac{\partial {\displaystyle \widehat{f}}}{\partial x}\right)}_j$ 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ⁿ⁺¹=Uⁿ+ΔU

$\left(I-{\left(\frac{\partial S}{\partial U}\right)}^n\Delta t\right)\Delta U=\left(-C^n+V^n+S^n\right)\Delta t$E40

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₀, y₀, z₀)…(x_l, y_l, z_l), the fluid velocity at the droplet position is obtained via the Lagrangian interpolation, which is expressed as

$u\left(x_{\text{d}},y_{\text{d}},z_{\text{d}}\right)={\displaystyle {\displaystyle \sum _{i=1-n/2}^{n/2}}{\displaystyle {\displaystyle \sum _{j=1-n/2}^{n/2}}{\displaystyle {\displaystyle \sum _{k=1-n/2}^{n/2}}\times }}}\left({\displaystyle {\displaystyle \prod _{\begin{array}{c}l=1-n/2\\ \left(l\ne i\right)\end{array}}^{n/2}}\frac{x_{\text{d}}-x_l}{x_i-x_l}}{\displaystyle {\displaystyle \prod _{\begin{array}{c}l=1-n/2\\ \left(l\ne j\right)\end{array}}^{n/2}}\frac{y_{\text{d}}-y_l}{y_i-y_l}}{\displaystyle {\displaystyle \prod _{\begin{array}{c}l=1-n/2\\ \left(l\ne k\right)\end{array}}^{n/2}}\frac{z_{\text{d}}-z_l}{z_i-z_l}}\right)u_{i,j,k}$E41

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,

$x_{d,i}^{n+1}=x_{d,i}^n+\Delta t\left(\frac{23}{12}{u}_{d,i}^n-\frac{16}{12}{u}_{d,i}^{n-1}+\frac{5}{12}{u}_{d,i}^{n-2}\right)$E42

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{array}{l}u=U_0\sin x\cos y\sin z\\ v=-U_0\cos x\sin y\sin z\\ w=0\\ p=p_0+{\rho }_0{U}_0{}^2(\cos 2z+2)(cos2x+cos2y)/16\end{array}$E43

The computing domain is $[0,2\pi ]\times [0,2\pi ]\times [0,2\pi ]$. The parameters are set as $p_0=100$, ${\rho }_0=1$, $U_0=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 64³ and 32³ 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 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 (\rho f{\displaystyle \overrightarrow{U}})=\nabla \cdot (\rho D\nabla f)$E44

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=\mu /\left(\rho \text{Sc}\right)$

$f(t=0)=\frac{1+erf(C_{0}Q)}{2},C_0=100$E45

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 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 H₂/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

${\displaystyle \overline{U}}=\frac{U_A+U_F}{2}+\frac{U_A-U_F}{2}\tanh (\frac{y-0.5l_y}{2{\delta }_0})$E46

where δ₀ 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 = 4l_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

$v\text{'}=(U_A-U_F)G(y-0.5l_y)A\sin (2\pi f_{\text{o}}t+\xi )$E47

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₂/O₂ reaction mechanism is used, which can be referred in [29].

Figure 7.

Instantaneous contour distribution of reactive flows. (a) temperature; (b) H₂; (c) H₂O.

Instantaneous distributions of temperature and mass fraction of H₂ and H₂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

$Z\equiv \frac{W_{\text{H}}}{W_{\text{H}}+W_{\text{O}}}$E48

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.

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

Figure 9.

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

Velocity ratio; UA/US	Temperature ratio; TA/TS	Pressure ratio; PA/PS	Convective Mach number; ${\text{M}}_{\text{c}}=\frac{U_A-U_S}{a_A+a_S}$
2	4	1	0.4
US (m/s)	TS (K)	PS (MPa)	Inflow Reynolds number; $\mathrm{Re}=\frac{\rho \left|U_A-U_S\right|{\delta }_0}{\mu }$
458.3	363.3	0.1	2731.2
ud (m/s)	Td (K)	ρd (kg/m^3)	Particle Stokes number (St0); ${\text{St}}_0\text{=}\frac{{\tau }_{\text{d}}}{{\tau }_{\text{f}}}=\left(\frac{{\rho }_{\text{d}}{d}_{\text{d, 0}}^2}{18{\mu }_{\text{S}}}\right)/\left(\frac{{\delta }_0}{\left|U_A-U_S\right|}\right)$
458.3	298.2	642.0	10.5

Table 1.

Inflow parameters.

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

${\text{C}}_{10}{\text{H}}_{22}+n_{{\text{O}}_2}{\text{O}}_2\Rightarrow n_{{\text{CO}}_2}{\text{CO}}_2+n_{{\text{H}}_2\text{O}}{\text{H}}_2\text{O}$E49

where n_i is the molar stoichiometric coefficients of each species. The global reaction constant, k (mol/cm³/s), is given by

$k=A\cdot T^n\exp \left(-\frac{E}{RT}\right){\left(\frac{\rho \cdot Y_{\text{F}}}{W_{\text{F}}\times {10}^6}\right)}^a{\left(\frac{\rho \cdot Y_{{\text{O}}_2}}{W_{{\text{O}}_2}\times {10}^6}\right)}^b$E50

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¹¹, E = 1.256 × 10⁵J/mol, a = 0.25, b = 1.5, and n = 0. Therefore, the combustion reaction rate per unit volume (g/cm³/s) is R_F = W_Vk × 10³.

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.25t_L, as shown in Figure 10(b). Here, the time t_L is the large eddy turnover time and is defined as

$t_{\text{L}}=\frac{{\delta }_{\text{L}}}{\left|U_{\text{A}}-U_{\text{S}}\right|}$E51

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.25t_L results in a larger temperature increase and facilitates the self-acceleration of the chemical reaction rate. At time t = t_ig + 0.5t_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.75t_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.5t_L disappear. The other reaction kernel entrained in large eddies is distorted due to the turbulence motions.

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 given by Ω = [11, 15] (kgm⁻³s⁻¹) in (a) and Ω = [15, 24] (kgm⁻³s⁻¹) in (b)–(d). Here, the dots indicate fuel droplets in the figures, with colors corresponding to dimensionless droplet diameters in the right figures.

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.