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.
- compressible reactive flows
- spray combustion
- numerical simulation
- Lagrangian trajectory model
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.  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 , 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 . 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.  simulated the time-developing, H2/air-reacting mixing layers and investigated the dynamics of backscatter of kinetic energy using DNS; Diegelmann et al.  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  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  proposed a hybrid compact-WENO scheme to obtain high resolutions for simulation of shock-turbulence interaction. Furthermore, Ren et al.  proposed a Roe type, characteristic-wise compact-WENO scheme. Considering overcoming the complexity and program parallelization of compact schemes, Hu et al.  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 .
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
A correction velocity
The equation of state used for ideal multispecies gas is introduced to close the above equations as
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 
The viscosity and thermal conductivity of the mixture are calculated by the expressions of Wake and Wassiljewa , respectively,
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
respectively. The Nusselt and Sherwood numbers are
respectively. Here, the droplet Reynolds number based on the slip velocity is defined 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,
where is the latent heat at
2.4. Two-way coupling models
The source terms
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 , the iteration from
The higher-resolution numerical scheme applied for the resolution of supersonic flow is essential to the reliability of simulation. The nonviscous flux , for the interface at
which is constructed by hybridizing a fifth-order Compact Upwind (CU) scheme, ,and a fifth-order WENO one, , through a smoothness indicator
A sixth-order symmetric compact difference scheme is applied for the viscous diffusion terms,
where is the difference approximation for at the node
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 . For instance,
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 . When the droplet is located at (
where the subscripts
In the Lagrangian droplet-tracking method, the droplet equation of motion was integrated with the third-order Adams-Bash-forth scheme in direction,
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 . The parameters are set as , , , 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 , WENO-compact upwind  (denoted as WENO-CU), and WENO-linear upwind  (denoted as WENO-LU), respectively. Two sets of computation grid are taken as 643 and 323 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
Here, erf is the error function. Figure 3 shows the initial distribution of
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
4.2. Reactive H2/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
The stream-wise and transverse domain lengths are
In order to stimulate the flow instability, the inlet velocity perturbation is given as
Instantaneous distributions of temperature and mass fraction of H2 and H2O at
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<
where WH and WO are the mass fraction of element H and element O, respectively. For the stoichiometric ratio,
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
||Convective Mach number
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
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
Large eddies rotate and the reaction kernels between two convection vortices are strained at time
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
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.
Tomboulides, AG, Lee, JCY and Orszag, SA. Numerical simulation of low Mach number reactive flows. Journal of Scientific Computing. 1997; 12(2):139–167.
Desjardins O, Blanquart G, Balarac G, et al. High order conservative finite difference scheme for variable density low Mach number turbulent flows. Journal of Computational Physics. 2008; 227(15):7125–7159.
Dauptain A, Cuenot B, Poinsot T. Large eddy simulation of a supersonic hydrogen-air diffusion flame. Complex Effects in Large Eddy Simulation. 2005; 71:98.
Berglund M, Fureby C. LES of supersonic combustion in a scramjet engine model. Proceedings of the Combustion Institute. 2007; 31(2):2497–2504.
Edwards J R, Boles J A, Baurle R A. Large-eddy/Reynolds-averaged Navier--Stokes simulation of a supersonic reacting wall jet. Combustion and Flame. 2012; 159(3):1127–1138.
Selle L, Lartigue G, Poinsot T, et al. Compressible large eddy simulation of turbulent combustion in complex geometry on unstructured meshes. Combustion and Flame. 2004; 137(4):489–505.
Lacaze G, Cuenot B, Poinsot T, et al. Large eddy simulation of laser ignition and compressible reacting flow in a rocket-like configuration . Combustion and Flame. 2009; 156(6):1166–1180.
Boudier G, Gicquel L Y M and Poinsot T, et al. Comparison of LES, RANS and experiments in an aeronautical gas turbine combustion chamber. Proceedings of the Combustion Institute. 2007; 31(2):3075–3082.
Saghafian A, Terrapon V E, Pitsch H. An efficient flamelet-based combustion model for compressible flows. Combustion and Flame. 2015; 162(3):652–667.
O’Brien J, Urzay J, Ihme M, et al. Subgrid-scale backscatter in reacting and inert supersonic hydrogen–air turbulent mixing layers. Journal of Fluid Mechanics. 2014; 743:554–584.
Diegelmann F, Tritschler V, Hickel S, et al. On the pressure dependence of ignition and mixing in two-dimensional reactive shock-bubble interaction. . Combustion and Flame. 2016; 163:414–426.
Jiang G S, Shu C W. Efficient implementation of weighted ENO schemes. Journal of Computational Physics. 1996; 126:202–228.
Pirozzoli S. Conservative hybrid compact-WENO schemes for shock-turbulence interaction. Journal of Computational Physics. 2002; 178(1):81–117.
Ren Y X, Zhang H. A characteristic-wise hybrid compact-WENO scheme for solving hyperbolic conservation laws. Journal of Computational Physics. 2003; 192(2):365--386.
Hu X Y, Wang B, Adams N A. An efficient low-dissipation hybrid weighted essentially non-oscillatory scheme. Journal of Computational Physics. 2015; 301:415–424.
Pirozzoli S. Numerical methods for high-speed flows. Annual review of fluid mechanics. 2011; 43(163–194)
Miller R S, Bellan J. Direct numerical simulation of a confined three-dimensional gas mixing layer with one evaporating hydrocarbon-droplet-laden stream. Journal of Fluid Mechanics. 1999; 384:293–338.
Miller R S. Effects of nonreacting solid particle and liquid droplet loading on an exothermic reacting mixing layer. Physics of Fluids. 2001; 13(11):3303–3320.
Reveillon J, Vervisch L. Analysis of weakly turbulent dilute-spray flames and spray combustion regimes. Journal of Fluid Mechanics. 2005; 537:317–347. 4
Wang Y, Rutland C J. Direct numerical simulation of ignition in turbulent n-heptane liquid-fuel spray jets. Combustion and Flame. 2007; 149(4):353–365.
Schroll P, Wandel A P, Cant R S, et al. Direct numerical simulations of autoignition in turbulent two-phase flows. Proceedings of the Combustion Institute. 2009; 32(2):2275–2282.
Neophytou A, Mastorakos E, Cant R S. The internal structure of igniting turbulent sprays as revealed by complex chemistry DNS. Combustion and Flame. 2012; 159(2):641–664.
Borghesi G, Mastorakos E, Cant R S. Complex chemistry DNS of n-heptane spray autoignition at high pressure and intermediate temperature conditions. Combustion and Flame. 2013; 160(7):1254–1275.
Xia J, Zhao H, Megaritis A, et al. Inert-droplet and combustion effects on turbulence in a diluted diffusion flame. Combustion and Flame. 2013; 160(2):366–383.
Burcat A, Ruscic B. Third millenium ideal gas and condensed phase thermochemical database for combustion with updates from active thermochemical tables. Argonne National Laboratory Argonne, IL; 2005.
Reid R C, Prausnitz J M, Poling B E. The properties of gases and liquids. McGraw-Hill, New York; 2001.
Eberhardt S, Imlay S. Diagonal implicit scheme for computing flows with finite rate chemistry. Journal of Thermophysics and Heat Transfer. 1992; 6(2):208–216.
Wang Q, Squires K D, Wu X. Lagrangian statistics in turbulent channel flow. Atmospheric environment. 1995; 29(18):2417–2427.
Nishioka M, Law C K. A numerical study of ignition in the supersonic hydrogen/air laminar mixing layer. Combustion and flame. 1997; 108(1):199–219.
Westbrook C K, Dryer F L. Simplified reaction mechanisms for the oxidation of hydrocarbon fuels in flames. Combustion science and technology. 1981; 27(1–2):31–43.