InTechOpen uses cookies to offer you the best online experience. By continuing to use our site, you agree to our Privacy Policy.

Computer and Information Science » Numerical Analysis and Scientific Computing » "Modeling and Simulation in Engineering Sciences", book edited by Noreen Sher Akbar and O. Anwar Beg, ISBN 978-953-51-2609-6, Print ISBN 978-953-51-2608-9, Published: August 31, 2016 under CC BY 3.0 license. © The Author(s).

Chapter 2

Numerical Simulation of Compressible Reactive Flows

By Bing Wang, Zhao Xin Ren and Wei Wei
DOI: 10.5772/63906

Article top


Comparison of the total kinetic energy. (a) Re = 30,000; (b) Re = 3000; (c) Re = 300.
Figure 1. Comparison of the total kinetic energy. (a) Re = 30,000; (b) Re = 3000; (c) Re = 300.
Energy spectrum of TGV with Re = 30,000.
Figure 2. Energy spectrum of TGV with Re = 30,000.
Initial distribution of passive scalar f.
Figure 3. Initial distribution of passive scalar f.
Variance curves of passive scalar.
Figure 4. Variance curves of passive scalar.
Passive scalar field in the plane of z = π. (a) t = 0 s; (b) t = 3 s; (c) t = 6 s; (d) t = 10 s.
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.
Schematic of computation model of spatially developing mixing layer.
Figure 6. Schematic of computation model of spatially developing mixing layer.
Instantaneous contour distribution of reactive flows. (a) temperature; (b) H2; (c) H2O.
Figure 7. Instantaneous contour distribution of reactive flows. (a) temperature; (b) H2; (c) H2O.
Flame structure. (a) temperature distribution and (b) water distribution.
Figure 8. Flame structure. (a) temperature distribution and (b) water distribution.
Schematic of a fuel spray in a two-dimensional turbulent shear layer.
Figure 9. Schematic of a fuel spray in a two-dimensional turbulent shear layer.
Flame kernels at (a) time t = tig; (b) time t = tig + 0.25 tL; (c) time t = tig + 0.5 tL; and (d) time t = tig + 0.75 tL. Instantaneous distributions of dimensionless gaseous temperature (left) and reaction rate (right), with red isolines given 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.
Figure 10. Flame kernels at (a) time t = tig; (b) time t = tig + 0.25 tL; (c) time t = tig + 0.5 tL; and (d) time t = tig + 0.75 tL. Instantaneous distributions of dimensionless gaseous temperature (left) and reaction rate (right), with red isolines given 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.

Numerical Simulation of Compressible Reactive Flows

Bing Wang1, Zhao Xin Ren and Wei Wei
Show details


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.

Keywords: compressible reactive flows, spray combustion, numerical simulation, WENO, Lagrangian trajectory model

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 [69]. Recently, some researchers have used DNS to simulate the compressible reactive flows. For example, O’Brien et al. [10] simulated the time-developing, H2/air-reacting mixing layers and investigated the dynamics of backscatter of kinetic energy using DNS; Diegelmann et al. [11] did the direct simulations of reactive shock-bubble interaction with detailed 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 [1724]. 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, ui is the velocity, T is the temperature, P is the pressure, Yk is the mass fraction of the kth species, and δij is the Kronecker delta function. The right-hand side terms of the above equations Sm, SF, and SQ describe the two-phase couplings of mass, momentum, and energy, respectively. Scombustion,k is the source term due to combustion. τij is the Newtonian viscous stress tensor


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


Here, NS is the total number of species. Wk and Xk are the molecule weight and the mole fraction of the kth species, respectively. et is the total energy, that is, kinetic energy plus internal (containing chemical) energy, which is defined as follows:


where cp,k and h0 f,k are the specific heat capacity at constant pressure and specific chemical formation enthalpy at the reference temperature, Tref, of the kth species, respectively. The values of cp,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 φij=(1+(μi/μj)1/2(Wj/Wi)1/4)2(8(1+Wi/Wj))1/2.


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

2.3. Equations of spray droplets

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


where cp is the specific heat of mixture gas, cL is the specific heat of the liquid. The momentum response time, τd, is defined as

where dd 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


f1 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, BM , and is given by


Here, YV is the mass fraction of the vapor on the far-field condition for the droplets and Ysf 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, WV is the molecular weight of the vapor, Patm is the atmospheric pressure, R is the universal gas constant, and TB,L is the liquid boiling temperature at Patm. LV is the latent heat of fuel vapor at liquid temperature, Td, and is given as follows:


where LV,TB,L is the latent heat at TB,L and TC,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 Nc,


where md and m˙d denote the droplet mass and dmd/dt, respectively, ud,i is the droplet velocity, ΔV is the volume of the computational grid for the gas-phase calculation, and hV,sf is the evaporated vapor enthalpy at the droplet surface. Fd and Qd 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, Scombustion,k, is expressed using the combustion reaction rate per unit volume, RF, as follows:


where nk and nF are the molar stoichiometric coefficients of the kth species and the fuel for the one-step global reaction (positive for the production side), respectively. Wk and WF are the molecular weights of the kth 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 dUj=Lj(U)dt, the iteration from n to iteration n+1 is performed as


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


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


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


where F^j is the difference approximation for (f^x)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 Un+1=UnU


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 (xd, yd, zd) in a cell whose coordinates are (x0, y0, z0)…(xl, yl, zl), the fluid velocity at the droplet position 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.

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 [0,2π]×[0,2π]×[0,2π]. The parameters are set as p0=100, ρ0=1, U0=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.

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


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)


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


where δ0 is the inlet layer thickness given as 0.4 mm, and ly is the length of computational domain in the y direction. UA is the hot-air inflow velocity, 2000 m/s, UF is the cold fuel inflow velocity, 1000 m/s. The temperatures of both streams are TA = 1600 K and TF = 390 K. The reactive system pressure is specified as pA = 1 atm.

The stream-wise and transverse domain lengths are lx and ly, respectively. By setting lx = 4ly, 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, fo 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

where WH and WO are the mass fraction of element H and element O, respectively. For the stoichiometric ratio, Zst = 0.111. It can be seen that the peaks of temperature and water are both near Zst. 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 UA and a spray carried by cold air moving at velocity US. 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
Temperature ratio
Pressure ratio
Convective Mach number
US (m/s)TS (K)PS (MPa)Inflow Reynolds number
ud (m/s)Td (K)ρd (kg/m3)Particle Stokes number (St0)
St0=τdτf=(ρddd, 0218μS)/(δ0|UAUS|)

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


where ni is the molar stoichiometric coefficients of each species. The global reaction constant, k (mol/cm3/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 × 1011, E = 1.256 × 105J/mol, a = 0.25, b = 1.5, and n = 0. Therefore, the combustion reaction rate per unit volume (g/cm3/s) is RF = WVk × 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 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 tig defines the moment that the incipient establishment of ignition kernels occurs. We normalize the gaseous temperature by T* = T/Tin and Tin = (TA + TS)/2. The droplet diameters are normalized by d* d = dd/dd0 and dd0 is the initial droplet diameter.

Large eddies rotate and the reaction kernels between two convection vortices are strained at time t = tig + 0.25tL, as shown in Figure 10(b). Here, the time tL 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 tig, the chemical heat release in the reaction kernels at time tig + 0.25tL results in a larger temperature increase and facilitates the self-acceleration of the chemical reaction rate. At time t = tig + 0.5tL, 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 = tig + 0.75tL. Comparison with Figure 10(c) shows that the reaction kernels that have high chemical reaction rates in the high-strain region at time t = tig + 0.5tL disappear. The other reaction kernel entrained in large eddies is distorted due to the turbulence motions.


Figure 10.

Flame kernels at (a) time t = tig; (b) time t = tig + 0.25 tL; (c) time t = tig + 0.5 tL; and (d) time t = tig + 0.75 tL. Instantaneous distributions of dimensionless gaseous temperature (left) and reaction rate (right), with red isolines given 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.

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.


1 - Tomboulides, AG, Lee, JCY and Orszag, SA. Numerical simulation of low Mach number reactive flows. Journal of Scientific Computing. 1997;12(2):139–167.
2 - 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.
3 - 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.
4 - Berglund M, Fureby C. LES of supersonic combustion in a scramjet engine model. Proceedings of the Combustion Institute. 2007;31(2):2497–2504.
5 - 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.
6 - 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.
7 - 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.
8 - 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.
9 - Saghafian A, Terrapon V E, Pitsch H. An efficient flamelet-based combustion model for compressible flows. Combustion and Flame. 2015;162(3):652–667.
10 - 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.
11 - 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.
12 - Jiang G S, Shu C W. Efficient implementation of weighted ENO schemes. Journal of Computational Physics. 1996;126:202–228.
13 - Pirozzoli S. Conservative hybrid compact-WENO schemes for shock-turbulence interaction. Journal of Computational Physics. 2002;178(1):81–117.
14 - 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.
15 - 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.
16 - Pirozzoli S. Numerical methods for high-speed flows. Annual review of fluid mechanics. 2011;43(163–194)
17 - 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.
18 - 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.
19 - Reveillon J, Vervisch L. Analysis of weakly turbulent dilute-spray flames and spray combustion regimes. Journal of Fluid Mechanics. 2005;537:317–347. 4
20 - 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.
21 - 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.
22 - 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.
23 - 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.
24 - 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.
25 - 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.
26 - Reid R C, Prausnitz J M, Poling B E. The properties of gases and liquids. McGraw-Hill, New York; 2001.
27 - 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.
28 - Wang Q, Squires K D, Wu X. Lagrangian statistics in turbulent channel flow. Atmospheric environment. 1995;29(18):2417–2427.
29 - 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.
30 - 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.