Effects of the Equivalence Ratio and Reynolds Number on Turbulence

Article pubs.acs.org/EF

Effects of the Equivalence Ratio and Reynolds Number on Turbulence and Flame Front Interactions by Direct Numerical Simulation Cong Xu, Zhihua Wang,* Wubin Weng, Kaidi Wan, Ronald Whiddon, and Angjian Wu State Key Laboratory of Clean Energy Utilization, Zhejiang University, Hangzhou, Zhejiang 310027, People’s Republic of China ABSTRACT: Direct numerical simulation (DNS) of a two-dimensional premixed turbulent syngas flame (50% H2 and 50% CO) was performed to investigate turbulence−flame interactions. An in-house DNS code, based on the low-Mach-number version of Navier−Stokes equations, was used to predict the influences of the equivalence ratio and Reynolds number on the turbulence−flame interactions using the constant volume enclosure method. The reaction between syngas and air was modeled with the Davis mechanism, which includes 14 species and 38 elemental reactions, and the species had distinct Lewis numbers. Syngas−air flames at combinations of four equivalence ratios and three Reynolds numbers were investigated. The flame structures at turbulence time scales were analyzed to demonstrate the influence of the turbulence on the flame structure, such as localized deformation of the flame front and isolated pockets of flame. The phenomenon was based on real effects detected by planar laser-induced fluorescence imaging of premixed syngas combustion under different turbulence conditions. The mechanism behind flame front disturbance was associated with various flame physical descriptors, e.g., flame thickness, tangential strain rate, curvature, displacement speed, flame length, and OH mass fraction. The flame thickness at the location of a bulge was usually thin, and the flame thickness on the gully was usually thick. A negative tangential strain rate was more likely to appear at lean stoichiometric ratios. Flashback occurred when the flame displacement speed was larger than the laminar flame speed. The flame length, another characteristic variable, had a positive correlation with the Reynolds number. The OH mass fraction, which is an important reaction intermediate, showed a strong negative correlation to flame front curvature. Correlation analysis indicated that the OH mass fraction was a better parameter to characterize the interaction of the turbulence and flame than the flame length.

1. INTRODUCTION With the implementation of increasingly strict pollution limits and unavailability of petroleum fuels, synthetic gas (syngas) is considered as a promising alternative energy source for certain applications, such as industrial process heaters, combined cycle gas turbine systems, etc.1−4 Sourced from pyrolysis and gasification of solid fuels, such as coal, biomass,5 and other carbon-rich materials, syngas provides a significantly cleaner utilization option compared to direct burning of those fuels. Premixed syngas combustion presents several complications when compared to traditionally fueled combustion systems, namely, the likelihood of combustion instabilities as a result of the transport, ignition, and flammability characteristics of hydrogen and carbon monoxide. This provides a strong motivation to develop accurate models of premixed syngas combustion for improving the application in traditional combustion systems. Detailed information regarding syngas combustion dynamics has shown difficulty to be obtained by experimental methods because turbulence produces a highly wrinkled and stretched flame front, which is difficult to be measured in a meaningful way.6 In this situation, computer modeling may be used to predict flame processes. Although large eddy simulation (LES) has been widely used in combustion research,7−9 it does not resolve interactions between the flame front and the full scale of turbulent energy. Therefore, despite direct numerical simulation (DNS) being computationally more expensive than LES, it is a more capable tool when investigating turbulence−flame front interactions. In © XXXX American Chemical Society

contrast to LES, DNS, which solve Navier−Stokes (N−S) equations free of modeling errors, can predict comprehensive information on flow and scalar fields in fundamental flame combustion problems. Thus, DNS is gradually becoming a standard modeling approach with the merits of its inherent high temporal and spatial resolution.10−14 Efforts to study the interactions in the flame front of a premixed turbulent flame have been made previously. One example is the interaction between turbulence and a laminar flame using a homogeneous isotropic turbulent flow (initially proposed by Orszag and Patterson15), which had been investigated using a simplified system, such as a constant Lewis number or one-step mechanism.16,17 Further, Baum et al.18 simulated premixed H2/O2/N2 flames in two dimensions with detailed chemical kinetics coupled with turbulent flow; the model results accurately predicted the complex H2 flame structure and OH field. Because of the accurate description of the reaction, recent studies used DNS with detailed chemistry to study the interaction mechanism in a H2 premixed flame. Day et al.19 suggested that burning in cellular structures was intense and that the burning sections were separated by regions in which the flame was effectively extinguished; the work of Luo et al.20 predicted that the distributions of flame curvature and shape factor of a flame kernel and planar flame indicated similar Received: January 12, 2016 Revised: July 25, 2016

A

DOI: 10.1021/acs.energyfuels.6b00068 Energy Fuels XXXX, XXX, XXX−XXX

Article

Energy & Fuels Table 1. Initial Turbulence Parameters in the Simulation case

ut (m/s)

η (μm)

Lii (mm)

le (mm)

1 2 3 4 5 6

1.8046 3.1562 5.4649 4.8134 2.1859 3.2789

63.0 35.9 20.5 23.1 28.1 26.5

1.023 0.583 0.333 0.375 0.250 0.333

3.847 1.437 1.253 1.411 0.940 1.253

τ (s) 8.1 2.6 8.7 1.1 4.9 8.7

× × × × × ×

10−3 10−3 10−4 10−3 10−4 10−4

Da

Ka

159 159 159 159 89.6 159

0.238 0.238 0.238 0.238 0.127 0.143

Table 2. Initial Chemical Parameters case

Φ

ReLii

X H2

XO2

s0L (m/s)

δl (mm)

1 2 3 4 5 6

0.7 1.0 2.0 3.0 2.0 2.0

100 100 100 100 30 60

0.1136 0.1479 0.2283 0.2788 0.2283 0.2283

0.1623 0.1479 0.1141 0.0929 0.1141 0.1141

0.5536 1.0517 1.7923 1.5482 1.7923 1.7923

0.5295 0.4713 0.3815 0.3908 0.3815 0.3815

local geometries. Rocco et al.21 found that flame curvature was associated with local super burning cells and quenching events depending upon its sign. The preceding simulations each focused on homogeneous fuels; simulation-based blended fuels, of which syngas is an example, are somewhat rare. Zhang et al. recently reported the modeling ignition of a CO/H2/air blend in homogeneous charge compression ignition (HCCI)22 and partially premixed combustion23 conditions. Battista et al.,24 using DNS with detailed chemistry, predicted that the addition of CO to H2/air flames reduced both local quenching and local temperature peaks. Understanding the interaction between the turbulence and flame front of a syngas premix flame is of great significance for improving the combustion device design and mitigating instability. In this paper, we presented two-dimensional (2D) DNS modeling with a detailed chemical mechanism of the instantaneous flame structure, variation of combustion characteristics, and Lewis number of syngas flames. Reasonable characterization of syngas premixed combustion was obtained by combining the information garnered in simulation, planar laser-induced fluorescence (PLIF) measurements, and associated statistical analysis. The flame length25 and the OH mass fraction26 were frequently used to characterize the behavior of the flame front in the experiment, and correlation analysis of the simulation statistics was studied. While three-dimensional (3D) DNS can account for the nature of turbulence, because it is computationally expensive, 2D DNS has been used. The statistics of small-scale quantities are different between 2D and 3D DNS; for instance, the smallest scales of motion do not follow the usual Kolmogorov scaling in 2D DNS. Nonetheless, the topology of a propagating flame structure in 3D turbulent flow was found to be primarily 2D, and those surfaces with the highest curvature had 2D characteristics.18 Many features computed with 2D DNS, including scaling of the flame-area-averaged mean tangential strain rate, are in accordance with the 3D findings.16

τf (s) 6.7 3.6 1.8 2.1 1.8 1.8

× × × × × ×

10−4 10−4 10−4 10−4 10−4 10−4

heating, bulk viscosity, body forces, and diffusion by pressure gradients and thermal radiation were not implemented in this model. The governing equations are simplified as below (in Cartesian coordinates)

∂ρ + ∇(ρ u) = 0 ∂t

(1)

∂(ρ u) ∂(ρ uu) + = −∇p + ∇τ ∂t ∂x

(2)

∂(ρYn) + ∇(ρ uYn) = ∇(ρ VnYn) + ρωn ∂t

n = 1, 2, 3, ..., Ns (3)

∂(ρT ) 1 1 + ∇(ρ uT ) = ∇(κ ∇T ) − ∂t cp cp −

1 cp

Ns

∑ ρcp, nYn Vn∇T n=1

Ns

∑ ρhnωn

(4)

n=1

⎛Y ⎞ p0 = ρ ∑ ⎜ n ⎟RT M n=1 ⎝ n ⎠ Ns

(5)

where ρ is the density of the gas mixture (kg/m ), u is the fluid velocity vector (m/s), τ is the fluid viscous vector (kg m−1 s−2), p is the hydrodynamic pressure (Pa) in eqs 2 and 3, ωn is the net formation rate of species n (s−1), Yn is the mass fraction of gas components n, Ns means the number of species, Vn denotes the species diffusion velocity (m/s), and the approximate mixture diffusivities can be calculated using a simple formula reported by Bird.28 In eq 4, κ is the thermal conductivity of the mixture (W m−1 K−1), T is the temperature (K), cp is the average specific heat capacity of mixed species (J kg−1 K−1), cp,n is the specific heat capacity of species n (J kg−1 K−1), and hn is the specific enthalpy of species n (J/kg). p0 is the thermodynamic pressure (Pa), and Mn is the molar weight of species n in eq 5. To improve the stability of the in-house code, a staggered spacetime, conservative discretization of the variable-density transport equations was implemented and a semi-implicit iterative method29 was adopted to integrate the equations. Velocity components were separated with density and other scalars by the staggered space-time mesh.30 The accuracy was enhanced by the fourth central-difference scheme.31 Pressure was solved as a Poisson equation. To provide the most accurate flame combustion simulation, detailed chemical kinetics were implemented, as derived by Davis et al.,32 containing 14 species and 38 elemental reactions and calculated by the DVODE solver33 to provide the source terms of species equations. The combustion environment was defined as a 12 × 12 mm square box, as shown in Figure 1a. No pressure reflection occurred in solving 3

2. NUMERICAL METHOD AND CASE DEFINITION 2.1. Numerical Method. The low-Mach-number series of N−S equations were used as the model cases defined by a low Mach number (Tables 1 and 2). The conservation equations of mass, momentum, species, and energy (as temperature) and the equation of state were derived from Pierce and Moin.27 Sundry key parameters involving acoustic interactions, compressibility, viscous dissipation B


Article

Energy & Fuels

Figure 1. Domain schematic and initial vorticity field generated by the PP spectrum method and initial temperature field: (a) schematic of the domain, (b) initial vorticity field, and (c) initial temperature field.

Figure 2. Validation of the DNS code for simulating syngas premixed combustion: (a) comparison of the laminar flame speed among the experimental data of McLean et al. and simulation from CHEMKIN and the DNS code [(○) experimental data,38 () CHEMKIN results, and (- - -) DNS code results] for 1D code validation and (b) gas component mass fraction and temperature profiles of the flame at the equivalence ratio of 2.0 for 1.0 atm predicted by the DNS code (on the basis of the mechanism by Davis et al.32). X is the length of the simulation domain in one dimension. is the root-mean-square (rms) velocity of the fluctuation. The energetic scale le (le = 2/ke), a second characteristic turbulent time scale τ (τ = ut2/ε),37 and the Kolmogorov scales η of the six cases are also listed (Table 1). The smallest Kolmogorov scale is 20.5 μm, considering that the uniform grid (600 × 600) has a spacing of 20 μm; therefore, all turbulent scales can be fully resolved. Table 2 gives the initial chemical parameters; all variables are previously described other than the equivalence ratio, Φ. The fuel entering the left boundary (inflow) of the domain consists of H2 and CO, at a molar ratio of 1:1, which parallels the experimental work of McLean et al.38 The mole fractions of H2 and O2 for each of the six cases are given in Table 2. The oxidant is 21% O2 + 79% N2. The premixed laminar flame speed calculated by CHEMKIN is used in the Dirichlet boundary condition (Table 2). The flame thickness δl is calculated as δl = (Tb − Tu)/(max(|dT/dx|)) (Tb is the temperature of the burnt side, while Tu is the temperature of the fresh side).18 The flame is fixed at standard operation conditions (temperature of 300 K, pressure of 1 atm, and undiluted fuel). The characteristic chemical time scale τf = δl/s0L is also shown in Table 2. MPICH2 is used to enable parallel computation of the code. The code was run on an 8 CPU IBM blade cluster with Intel quad core 2.0 GHz CPUs and 10 Gb/s Infini-band network. The simulation required 4 weeks to finish. To observe the interaction between turbulence and flame front clearly, the thickness of the flame front should be thinner than the turbulent integral scale but thicker than the Kolmogorov scale. This was the wrinkled or corrugated region of the premixed turbulent combustion diagram as defined by Poinsot et al.39 The Da number used in the six cases was either 159 or 89.6, and the Ka number used was 0.238, 0.127, or 0.143. Because the purpose of this work was to discuss variations in flame−turbulence interactions as a function of Reynolds number and Φ, the flame system was solved for various values of these variables: Φ = 0.7, 1.0, 2.0, and 3.0, and ReLii = 30, 60, and 100.

incompressible equations; therefore, convective boundary conditions were handled with velocity and momentum at the outflow.34 A fifthorder, weighted, essentially non-oscillatory scheme (WENO 5) was used for treating scalars at the boundary35 to improve the accuracy and stability of the code. The top and bottom sections of the domain were regarded as periodic boundaries.15 The inflow boundary adopted Dirichlet boundary conditions with constant velocity and scalars. The model initialization proceeded as follows, which was required to ensure the stability of the simulation. First, a laminar premixed onedimensional (1D) flame was modeled in CHEMKIN 4.1 using the Davis mechanism. The CHEMKIN result was used as the starting condition for the in-house code, which solved the same conditions in 1D. In this way, the agreement between the different discretization schemes used by the in-house code and CHEMKIN could be evaluated. Provided good agreement, the scalars (species and temperature) and velocity fields predicted at the second step were inputted into the 2D DNS model as the initial condition (the initial temperature field was shown in Figure 1c). The initial homogeneous isotropic turbulent flow field was created with the Passot−Pouquet (PP) spectrum method,36 and the velocity field was obtained in the second step (the initial vorticity field is shown in Figure 1b). The turbulence energy spectrum E(k), which is used to obtain the turbulent kinetic energy and the dissipation rate, is defined by

E(k) =

16ut 2 ke

4 ⎛ ⎞ 2 ⎛ k ⎞ −2⎜⎝ kk ⎟⎠ e ⎜ ⎟e π ⎝ ke ⎠

2

(6)

where k is the wavenumber and ke is the most energetic vortex. 2.2. Case Definition. Six cases were defined to model flame− turbulence interactions. The initial turbulence parameters are shown in Table 1, and the chemical parameters are shown in Table 2. The characteristic turbulence intensity ReLii (ReLii = utLii/v) is computed with the autocorrelation integral scale Lii. ReLii = 30, 60, and 100 represent low, medium, and high turbulence intensities, respectively. ut C


Article

Energy & Fuels

Figure 3. Effect of different equivalence ratios and Reynolds numbers on mass fraction fields of O2 and the black solid line of c = 0.5, in which c is a variable defined by eq 8 to characterize the reaction progress, which is also labeled using a black solid line: (a) Φ = 0.7 and ReLii = 100, (b) Φ = 1.0 and ReLii = 100, (c) Φ = 2.0 and ReLii = 100, (d) Φ = 3.0 and ReLii = 100, (e) Φ = 2.0 and ReLii = 30, and (f) Φ = 2.0 and ReLii = 60.

3. RESULTS AND DISCUSSION 3.1. 1D Verification. Laminar flame speeds calculated using CHEMKIN and the in-house DNS code are shown in Figure 2a. The experimental data comes from constant-pressure expanding spherical flame measurements.38 The laminar flame speed from CHEMKIN was calculated from a freely propagating premixed H2/CO/air flame in a 1D domain with an adaptive grid (∼410 grid points) using the PREMIX module. The in-house code was run in a 1D configuration, and the laminar flame speed was calculated from a temperature distribution approach. The integral of the temperature at a spatial coordinate is calculated as

ΔQ =

∫ (Tt+Δt − Tt) dx

(7)

where T is the temperature at a specific position at time t and t + Δt is used to trace the motion of the modeled flame front. A first prediction of the laminar flame speed coming from CHEMKIN is used for the initial inflow velocity; the inflow velocity varies with the value of ΔQ in subsequent iterations. When the simulation has run for more than 20 flow-through times and the position of the flame remains static, the inlet velocity at that time is the laminar flame speed of the current case. D


Article

Energy & Fuels

Figure 4. Local zoom-in images of the temperature gradient magnitude field at different equivalence ratios and Reynolds numbers. The black line is the isoline of c = 0.12, which indicates the maximum temperature gradient magnitude. The red color means a bigger magnitude, and the blue color means a smaller magnitude.

Figure 5. PDF of the flame thickness at the time of t/τ = 1.0: (a) different equivalence ratio effects and (b) different Reynolds number effects.

A total of 10 cases, with Φ = 0.7, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, and 5.0, are considered in this work. The agreement between the experimental data (○), the CHEMKIN results (), and the DNS code results (- - -) demonstrate that the Davis mechanism32 used in concert with our in-house code can render satisfactory predictions for premixed combustion of H2− CO−air mixtures, especially in the equivalence ratio range from 0.7 to 5.0. In addition to laminar flame speed prediction, accurate prediction of the gas composition and temperature are important indicators of model accuracy. The spatial distribution of various reactant/product gas mole fractions and temperature predicted by the in-house code are given in Figure 2b; the predicted distribution is typical of a premixed flame structure. 3.2. 2D Flame Structure. Reaction progress (c) is a useful variable when discussing the behavior of the flame front, e.g., flame front−turbulence interactions; it is determined by the consumption of O2 as expressed by40 c=1−

Using the definition of c and the transport equation of the O2 species, the following equation can be derived: Sd =

S 1 dA = ∇t v + d A dt R

+ ui c = c0

∂c ∂xi

(8)

= Sd |∇c||c = c0 c = c0

(10)

(11)

where v is the velocity vector and R is the effective radius of the flame curvature in eq 11. The first term of the right-hand side (RHS) of the equation is the tangential strain rate caused by hydrodynamic strains acting on the plane parallel to the flame front. The second term on the RHS is the curvature term caused by propagation of the reaction front. This stretch is actually caused by the non-uniform flow field and the propagation of the flame. Initial conditions for the 2D simulation were created from the species profiles, temperature, and velocities solved in the 1D simulation, including the predicted laminar flame speed. The instantaneous snapshot of O2 mass fraction contours for the six cases calculated (Tables 1 and 2) at time t/τ = 1.0 are

According to eq 8, c is bounded by c = 0.0 at the inflow and c = 1.0 at the outflow, providing that complete combustion occurs for Φ > 1.0. The c transport equation for the isoline of c = c0 may be defined as41 ∂c ∂t

c = c0

where Sd is the displacement speed of the isoline of c = c0 and D is the diffusivity. The interaction between the turbulence and the flame front is exemplified by changes in the flame front, which induces changes in the flame area. The change of flame area is defined by eq 1142,43

YO2 max(YO2)

ωc + ∇(ρD∇c) ρ|∇c|

(9) E


Article

Energy & Fuels given in Figure 3. Turbulence distorts an initially planar flame front, which can result in pockets of fresh gas that may separate from the flame front and be transported into the burnt region as a result of the intensified vortical motion (observed in Figure 3d). The average flame front position for Φ = 0.7 and 1.0 (panels a and b of Figure 3) is located far upstream of the initial position shown in Figure 1, while the flame front remains stable at Φ = 2.0 and 3.0 (panels c−f of Figure 3). This flashback in the Φ = 0.7 and 1.0 cases can be attributed to dynamic instability as a result of the local Lewis number being less than 1.0 for lean combustion cases.6 Flame thickness is a variable relating to the flame geometric structure, which is useful in describing the influence of vortical motion and local reaction rate on the flame structure. It is defined as δf = (Tb − Tu)/max( (∂T /∂x)2 + (∂T /∂y)2 )

3.3. Statistics on the Flame Front. The tangential strain rate is used to evaluate the action of vortex stretching and compression on the flame surface. αt* = τf αt = τf

∫c=0.5 αt dl/∫c=0.5 dl

(13)

N

αsd* = τf

∑ (α t i=1

i

− αt)2 /N (14)

where α*t is the normalized flame-area-average tangential strain rate on the isoline of c = 0.5, αt is the tangential strain rate ∇tv in eq 11, and αsd * is the normalized standard deviation of the tangential strain rate. The average tangential strain rate at Φ = 2.0 is slightly larger than that at Φ = 1.0 and 3.0, while the average tangential strain rate at Φ = 0.7 is much smaller than that in the other cases (Figure 7a). In all cases, the average

(12)

where δf is the flame thickness, Tb is the burned gas temperature, Tu is the fresh gas temperature, and ((∂T/∂x)2 + (∂T/∂y)2)1/2 is the temperature gradient magnitude, which may be simplified to |∇T|. The maximum value of |∇T| is found at c = 0.12; this is indicated by a solid black isoline in the calculated |∇T| for the various model cases presented in Figure 4. A bigger |∇T| is usually found on the bulge of a flame front, and a smaller |∇T| is usually on the gully. However, where the flame location is highly variable, a smaller |∇T| is on the bulge. The probability density functions (PDFs) of flame thickness along the isoline of Figure 4 are shown in Figure 5. Among the various equivalence ratios tested (Figure 5a), Φ = 2.0 has the thinnest flame thickness while moving to either richer or leaner equivalence ratio results in broader flame widths. As shown in Figure 5b, the calculated flame thicknesses are more widely distributed as the Reynolds number is increased but the most probable thickness decreases. Stronger vortical motion will make the flame front more complex, but the average comprehensive flame thickness is thinner. Turbulent premixed flame experiments were performed using a fuel composed of 50% H2 and 50% CO in a round jet burner using parameters similar to simulation cases 3 and 5. Singleshot OH images obtained by the PLIF method44 are shown in Figure 6. In the higher turbulence flame, the flame structure becomes more complex, with structures that would lead to pockets of unburned gas separating from the larger flame surface. This phenomenon is consistent with predictions made by simulation.

Figure 7. Flame-area-average-tangential-strain-rate-based statistics normalized by the chemical time scale τf versus time normalized by τf: (a) results based on statistics at equivalence ratios of 0.7, 1.0, 2.0, and 3.0 and (b) results based on statistics at Reynolds numbers of 30, 60, and 100. The error bar is the standard deviation of the tangential strain rate.

tangential strain rate remained positive, which is consistent with previous research.18 The strain rate fluctuation (error bars in Figure 7a) trends similarly with the averaged strain rate values, with the largest fluctuations associated with Φ = 2.0 and generally decreasing as the strain rate decreases. The extensive tangential strains on the flame at Φ = 2.0 are greater than the other cases; tangential pulling may explain the decrease in flame thickness found in Figure 4, and the trend is consistent with Figure 5a. The normalized flame-area-average tangential strain rate αt* at normalized time t/τf for different Reynolds numbers is shown in Figure 7b. The error bar is the normalized standard deviation of the tangential strain rate. The average tangential strain rate and the fluctuation become larger as the Reynolds number increases, and the stronger vortical motion can give the flame larger extensive tangential strains and make the flame thickness thinner. The PDFs of the tangential strain rate are used to evaluate variation of the tangential strain rate with the equivalence ratio and Reynolds number. The PDFs of the instantaneous tangential strain rate normalized by the chemical time scale τf at different equivalence ratios are clear in Figure 8a. In the stable flame front cases, Φ = 2.0 and 3.0, the PDF of the tangential strain rate is consistent with previous stud-

Figure 6. PLIF transient OH images of the 50% H2/50% CO premixed flame: (a) low turbulence and (b) high turbulence. F


Article

Energy & Fuels

Figure 8. PDFs of the instantaneous tangential strain rate normalized by the chemical time scale τf at the time of t/τ = 1.0: (a) different equivalence ratio effects and (b) different Reynolds number effects.

ies;16,18,41,45 the PDF of the transient tangential strain rate at Φ = 2.0 is normally distributed over the entire extent of the normalized tangential strain rate, which demonstrates that the flame surface structure is more complex. The PDFs also show that, at Φ = 0.7 and 1.0, the greatest likelihood is to encounter negative tangential strain. This has rarely been reported in previous studies. Previous research on lean fuel premixed flames found positive tangential strain rates with a stable flame and negative tangential strain rates when the flashback occurs. In comparison of the three PDFs in Figure 8b, the probability density of the larger tangential strain rate increases as the Reynolds number increases; therefore, the high turbulent intensity will make separation of the flame front more likely. Curvature is used to evaluate the degree of reaction at the flame surface. Cur* = τf Cur = τf

∫c=0.5 Cur dl/∫c=0.5 dl

Figure 9. Flame-area-average curvature normalized by the chemical time scale τf along time: (a) results based on statistics at equivalence ratios of 0.7, 1.0, 2.0, and 3.0 and (b) results based on statistics at Reynolds numbers of 30, 60, and 100. The error bar is the standard deviation of the curvature.

(15)

N

Cur*sd = τf

∑ (Curi − Cur)2 /N i=1

(16)

The displacement speed defined by eq 10, which characterizes the propagation of the flame front, is used to explain the motion of the flame front. The PDFs of displacement speed Sd on the isoline of c = 0.5 at t = τ0 are shown in Figure 11. There is a high probability in the flames with Φ = 0.7 and 1.0 of high negative speed, where the negative sign indicates that the propagation is in the fresh gas. The universal displacement speed at Φ = 0.7 is much larger than its laminar flame speed, and the universal displacement speed at Φ = 1.0 is slightly larger than its laminar flame speed. In contrast, the universal displacement speeds at Φ = 2.0 and 3.0 are smaller than their laminar flame speeds. It is shown in Figure 12 that the isoline of c = 0.5 coincides with the region of greatest reaction intensity for Φ = 0.7 and 1.0; however, at Φ = 2.0 and 3.0, the isolines are upstream of the most intense reaction, which may explain the small displacement speeds at Φ = 2.0 and 3.0 based on the definition of Sd in eq 10. The PDFs of displacement speed at different Reynolds numbers in Figure 11b show that a Reynolds number of 30 has a narrower range of displacement speed than Reynolds numbers of 60 and 100; however, as the Reynolds number grows, the distribution of the displacement speed is essentially the same. The reason may be that our vortex scales of Re of 60 and 100 are the same, and the turbulence shows little effect on the displacement speeds at Re of 60 and 100. 3.4. Parameter Characterizing the Flame Front. The flame length, which is defined by Lflame = ∫ c = 0.5 dl along the isoline of c = 0.5 and normalized by the initial flame length L0,

where Cur is used to represent Sd/R in eq 11. In a simplified sense, Cur* is the normalized flame-area-average curvature and Cur*sd is the normalized standard deviation of that curvature. The average curvature at Φ = 0.7 is negative, while that at other equivalence ratios is near zero, which is in accordance with the research by Becker et al.46 and Baum et al.18 The fluctuations of curvature at Φ = 0.7, 1.0, and 2.0 have some notable peaks (see the error bars in Figure 9a), each of which corresponds to an instance of flame front separation. Logically, the curvature values are highly influenced by the behavior of the flame front. As for the influence of the Reynolds number shown in Figure 9b, the average curvatures of these three cases under different Reynolds numbers remain zero, while the probability of extreme peaks increases as the Reynolds number increases. Therefore, as the Reynolds number becomes larger, more separated pockets of combustion will appear. The PDFs of instantaneous curvature normalized by the chemical time scale τf at the time of t/τ = 1 for different equivalence ratios are shown in Figure 10a. The PDFs (Figure 10a) are nearly symmetric with the zero curvature line, which is consistent with other studies.16,18,20,41 The broad distribution of curvature also reflects the complex flame front structure. The distribution of the PDF at a high Reynolds number is broader than that at a low Reynolds number (Figure 10b), indicating that the flame surface becomes more complex as the Reynolds number increases. G


Article

Energy & Fuels

Figure 10. PDFs of the instantaneous curvature normalized by the chemical time scale τf at the time of t/τ = 1.0: (a) different equivalence ratio effects and (b) different Reynolds number effects.

Figure 11. PDFs of the displacement speed at the time of t/τ = 1.0: (a) different equivalence ratio effects and (b) different Reynolds number effects.

Figure 12. Local images of the O2 reaction rate field at different equivalence ratios. The black line is c = 0.5.

flame length. However, the Reynolds number is found to be positively correlated with the flame length (Figure 13b). The increase in flame length is due to turbulence transiting the flame front. The turbulence formed at high Reynolds number is much stronger; thus, the interaction is stronger.

along the normalized time at different equivalence ratios and different Re is displayed in panels a and b of Figure 13, respectively. From Figure 13a, no monotonic trend can be observed from different cases, which means there is no significant correlation between the equivalence ratio and the H


Article

Energy & Fuels

Figure 13. Flame length normalized by the initial flame length L0 along time: (a) different equivalence ratio effects and (b) different Reynolds number effects.

Figure 14. Scatter plots of statistics normalized by the chemical time scale τf versus the mass fraction of OH normalized by the initial mass fraction YOH,0: (a) statistics of the tangential strain rate and (b) statistics of curvature.

Table 3. Correlation Coefficient between the Area Change and the Parameters (Lf and OH) case Lf OH

p r p r

1

2

3

4

5

6

0.000 −0.855 0.000 −0.718

0.060 0.268 0.000 −0.900

0.004 0.400 0.000 −0.620

0.000 −0.502 0.000 −0.676

0.000 −0.647 0.000 −0.501

0.000 −0.758 0.000 −0.851

As an important combustion intermediate product,47,48 OH can be a significant intermediate product in syngas combustion and can be used as a flame front marker. Figure 14 presents scatter plots of the tangential strain rate normalized by the chemical time scale τf versus the mass fraction of OH normalized by the initial mass fraction YOH,0 (Figure 14a) and scatter plots of normalized tangential strain rate versus the normalized mass fraction of OH (Figure 14b). The distribution of the OH mass fraction is generally in the positive region of the tangential strain rate, but there is no apparent correlation

between the OH mass fraction and the tangential strain rate. On the other hand, the OH mass fraction shows a strong negative correlation with the curvature. The OH mass fraction is higher at the bulges of the flame front and smaller at the gullies of the flame front. The OH mass fraction grows with the decrease of the radius of the bulges. This trend in fuel-lean conditions is consistent with the research by Rocco et al. on lean-H2 premixed flames21 and the research by Battista et al. on lean hydrogen−carbon monoxide turbulent premixed flames.24 The temperature gradient magnitude on the bulge of a flame I


Article

Energy & Fuels

flame length. Some considerations regarding changes of the flame front and the occurrence of flashback and separated flame pockets were presented, which could improve the understanding of stable syngas combustion. The DNS presented should prove useful for better understanding the flame− turbulence interaction and development of strategies for stable turbulent premixed combustion.

surface is bigger, as shown in Figure 4, which is due to the higher chemical activity at the bulge, which, in turn, causes the increase of the OH mass fraction. This same trend also exists in the fuel-rich conditions. The temperature gradient on the bulge is bigger on the fuel-rich cases, and it causes an increase of the OH mass fraction. The tangential strain rate and curvature, which are the direct parameters that characterize the change of the flame front, could not be obtained directly from experiments. However, the flame length and OH mass fraction can be obtained. Thus, a partial correlation analysis was made between the change of the flame front area (dA/dt/A), sum of the tangential strain rate, and curvature, defined in eq 11), and the above parameters (flame length Lf and OH mass fraction) using the statistics software SPSS. The result is shown in Table 3, with p representing the significance in the two-tailed significance test and r representing the Pearson correlation coefficient49 r=

cov(X , Y ) σXσY

■

AUTHOR INFORMATION

Corresponding Author

*Telephone: +86-181-67-102328. E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by the National Natural Science Foundation of China (51422605) and the National Basic Research Program of China (2012CB214906).

■

(17)

where X means dA/dt/A and Y means Lf or OH mass fraction in eq 17, cov(X,Y) is the covariance of X and Y, and σX and σY are the standard deviations of X and Y, respectively. Significant linear correlation exists when p < 0.05, and on that basis, the correlation is stronger when |r| is approaching 1. It can be concluded that the correlation between the flame front area change and the OH mass fraction is more significant by comparison of the p values of Lf and the OH mass fraction in the six cases because the p value of Lf is bigger than 0.05 at case 2. r describes the degree of linear correlation between the two variables. Generally speaking, the Pearson correlation coefficient of the OH mass fraction is bigger than the flame length, indicating the OH mass fraction being a better variable to characterize the change of the flame front than the flame length. The Pearson correlation coefficient of the OH mass fraction is negative, indicating that the larger the value of the OH mass fraction, the smaller the change of the flame front.

REFERENCES

(1) Casleton, K. H.; Breault, R. W.; Richards, G. A. System issues and tradeoffs associated with syngas production and combustion. Combust. Sci. Technol. 2008, 180 (6), 1013−1052. (2) Basshuysen, R. V.; Schaefer, F. Internal Combustion Engine Handbook, Basics, Components, Systems, and Perpectives; SAE International: Warrendale, PA, 2004; DOI: 10.4271/R-345. (3) Reitz, R. D. Directions in internal combustion engine research. Combust. Flame 2013, 160 (1), 1−8. (4) Gas Turbine Performance; Walsh, P. P., Fletcher, P., Eds.; Blackwell Science: Oxford, U.K., 2004; DOI: 10.1002/ 9780470774533. (5) McKendry, P. Energy production from biomass (part 1): Overview of biomass. Bioresour. Technol. 2002, 83 (1), 37−46. (6) Poinsot, T.; Veynante, D. Theoretical and Numerical Combustion; RT Edwards, Inc.: Philadelphia, PA, 2005. (7) Veynante, D.; Vervisch, L. Turbulent combustion modeling. Prog. Energy Combust. Sci. 2002, 28 (3), 193−266. (8) Dodoulas, I. A.; Navarro-Martinez, S. Large Eddy Simulation of Premixed Turbulent Flames Using the Probability Density Function Approach. Flow, Turbul. Combust. 2013, 90 (3), 645−678. (9) Schmitt, T.; Boileau, M.; Veynante, D. Flame Wrinkling Factor Dynamic Modeling for Large Eddy Simulations of Turbulent Premixed Combustion. Flow, Turbul. Combust. 2015, 94 (1), 199−217. (10) Wang, Z.; Fan, J.; Zhou, J.; Cen, K. Direct numerical simulation of hydrogen turbulent lifted jet flame in a vitiated coflow. Chin. Sci. Bull. 2007, 52 (15), 2147−2156. (11) Li, X.; Fu, D.; Ma, Y. DNS of compressible turbulent boundary layer around a sharp cone. Sci. China, Ser. G: Phys., Mech. Astron. 2008, 51 (6), 699−714. (12) Luo, K.; Fan, J.; Cen, K. Transitional phenomenon of particle dispersion in gas-solid two-phase flows. Chin. Sci. Bull. 2007, 52 (3), 408−417. (13) Huang, Z.; Zhou, H.; Luo, J. The investigation of coherent structures in the wall region of a supersonic turbulent boundary layer based on DNS database. Sci. China, Ser. G: Phys., Mech. Astron. 2007, 50 (3), 348−356. (14) Huang, Z.; Zhou, H. Inflow conditions for spatial direct numerical simulation of turbulent boundary layers. Sci. China, Ser. G: Phys., Mech. Astron. 2008, 51 (8), 1106−1115. (15) Orszag, S.; Patterson, G. S., Jr. Numerical simulation of turbulence. In Statistical Models and Turbulence; Rosenblatt, M., Van Atta, C., Eds.; Springer: Berlin, Germany, 1972; Lecture Notes in Physics, Vol. 12, pp 127−147, DOI: 10.1007/3-540-05716-1_8. (16) Haworth, D.; Poinsot, T. Numerical simulations of Lewis number effects in turbulent premixed flames. J. Fluid Mech. 1992, 244, 405−436.

4. CONCLUSION Six cases of premixed turbulent syngas combustion, comprising four different equivalence ratios and three different Reynolds numbers, were simulated by DNS using a detailed chemical mechanism and various Lewis numbers. Two obvious phenomena were observed in the instantaneous snapshot of the O2 field. The flame transitioned upstream (flashback) when the equivalence ratio was less than 1.0. Alternatively, at higher equivalence ratios, combustion pockets separated from the flame and were transported into the burnt region. It was found that a negative tangential strain rate was more likely when Φ < 1.0, which was rarely noted in previous research and could indicate the presence of flashback. A premixed flame with a high Reynolds number would make an increase in the tangential strain rate and curvature, leading to a more complex flame structure. The fluctuation of the tangential strain rate under a high Reynolds number was also stronger than that under a low Reynolds number, which was also substantiated by PLIF images. The mass fraction of OH exhibited a strong negative correlation with flame curvature. It was found that the trend existed not only on the fuel-lean side, which had been found in previous research, but also on the fuel-rich side. The correlation analysis suggested that the OH mass fraction was found to be a better marker for flame front change than the J


Article

Energy & Fuels (17) Poinsot, T.; Haworth, D.; Bruneaux, G. Direct simulation and modeling of flame-wall interaction for premixed turbulent combustion. Combust. Flame 1993, 95 (1), 118−132. (18) Baum, M.; Poinsot, T.; Haworth, D.; Darabiha, N. Direct numerical simulation of H2/O2/N2 flames with complex chemistry in two-dimensional turbulent flows. J. Fluid Mech. 1994, 281, 1−32. (19) Day, M.; Bell, J.; Bremer, P.-T.; Pascucci, V.; Beckner, V.; Lijewski, M. Turbulence effects on cellular burning structures in lean premixed hydrogen flames. Combust. Flame 2009, 156 (5), 1035− 1045. (20) Luo, K.; Wang, H.; Bushe, W. K.; Fan, J. Direct numerical simulation and reaction rate modelling of premixed turbulent flames. Int. J. Hydrogen Energy 2014, 39 (23), 12158−12165. (21) Rocco, G.; Battista, F.; Picano, F.; Troiani, G.; Casciola, C. Curvature effects in turbulent premixed flames of H2/Air: A DNS study with reduced chemistry. Flow, Turbul. Combust. 2015, 94 (2), 359−379. (22) Zhang, F.; Yu, J.; Yu, R.; Jangi, M.; Bai, X.-S. DNS of CO/H2/air ignition in a constant volume enclosure relevant to direct injection HCCI engine. Proceedings of the 7th Mediterranean Combustion Symposium, MCS7; Chia Laguna, Cagliari, Sardinia, Italy, Sept 11− 15, 2011. (23) Zhang, F.; Yu, R.; Bai, X. S. Detailed numerical simulation of syngas combustion under partially premixed combustion engine conditions. Int. J. Hydrogen Energy 2012, 37 (22), 17285−17293. (24) Battista, F.; Picano, F.; Troiani, G.; Casciola, C. Direct Numerical Simulation of Hydrogen−Carbon Monoxide Turbulent Premixed Flame. In Direct and Large-Eddy Simulation IX; Springer: Cham, Switzerland, 2015; ERCOFTAC Series, Vol. 20, pp 541−546, DOI: 10.1007/978-3-319-14448-1_69. (25) Yamamoto, K.; Isii, S.; Ohnishi, M. Local flame structure and turbulent burning velocity by joint PLIF imaging. Proc. Combust. Inst. 2011, 33 (1), 1285−1292. (26) Yang, L.; Wang, Z.; Zhu, Y.; Li, Z.; Zhou, J.; Huang, Z.; Cen, K. Premixed jet flame characteristics of syngas using OH planar laser induced fluorescence. Chin. Sci. Bull. 2011, 56 (26), 2862−2868. (27) Pierce, C. D.; Moin, P. Progress-variable approach for large-eddy simulation of non-premixed turbulent combustion. J. Fluid Mech. 1999, 504 (504), 73−97. (28) Bird, R. B. Transport phenomena. Appl. Mech. Rev. 2002, 55 (1), R1−R4. (29) Najm, H. N.; Wyckoff, P. S.; Knio, O. M. A semi-implicit numerical scheme for reacting flow: I. Stiff chemistry. J. Comput. Phys. 1998, 143 (2), 381−402. (30) Harlow, F. H.; Welch, J. E. Numerical calculation of timedependent viscous incompressible flow of fluid with free surface. Phys. Fluids 1965, 8 (12), 2182. (31) Morinishi, Y.; Lund, T.; Vasilyev, O.; Moin, P. Fully conservative higher order finite difference schemes for incompressible flow. J. Comput. Phys. 1998, 143 (1), 90−124. (32) Davis, S. G.; Joshi, A. V.; Wang, H.; Egolfopoulos, F. An optimized kinetic model of H 2/CO combustion. Proc. Combust. Inst. 2005, 30 (1), 1283−1292. (33) Brown, P. N.; Byrne, G. D.; Hindmarsh, A. C. VODE: A variable-coefficient ODE solver. SIAM J. Sci. Stat. Comput. 1989, 10 (5), 1038−1051. (34) Orlanski, I. A simple boundary condition for unbounded hyperbolic flows. J. Comput. Phys. 1976, 21 (3), 251−269. (35) Jiang, G.-S.; Shu, C.-W. Efficient Implementation of Weighted ENO Schemes. J. Comput. Phys. 1996, 126, 202−228. (36) Passot, T.; Pouquet, A. Numerical simulation of compressible homogeneous flows in the turbulent regime. J. Fluid Mech. 1987, 181, 441−466. (37) Pope, S. B. Turbulent Flows; Cambridge University Press: Cambridge, U.K., 2000; DOI: 10.1017/CBO9780511840531. (38) McLean, I. C.; Smith, D. B.; Taylor, S. C. The use of carbon monoxide/hydrogen burning velocities to examine the rate of the CO + OH reaction. Symp. (Int.) Combust., [Proc.] 1994, 25, 749−757.

(39) Poinsot, T.; Veynante, D.; Candel, S. Quenching processes and premixed turbulent combustion diagrams. J. Fluid Mech. 1991, 228, 561−606. (40) Vreman, A.; Van Oijen, J.; De Goey, L.; Bastiaans, R. Direct numerical simulation of hydrogen addition in turbulent premixed Bunsen flames using flamelet-generated manifold reduction. Int. J. Hydrogen Energy 2009, 34 (6), 2778−2788. (41) Echekki, T.; Chen, J. H. Unsteady strain rate and curvature effects in turbulent premixed methane-air flames. Combust. Flame 1996, 106 (1), 184−202. (42) Poinsot, T.; Echekki, T.; Mungal, M. A study of the laminar flame tip and implications for premixed turbulent combustion. Combust. Sci. Technol. 1992, 81 (1−3), 45−73. (43) Candel, S. M.; Poinsot, T. J. Flame stretch and the balance equation for the flame area. Combust. Sci. Technol. 1990, 70 (1−3), 1− 15. (44) Weng, W.; Nilsson, E.; Ehn, A.; Zhu, J.; Zhou, Y.; Wang, Z.; Li, Z.; Aldén, M.; Cen, K. Investigation of formaldehyde enhancement by ozone addition in CH4/air premixed flames. Combust. Flame 2015, 162 (4), 1284−1293. (45) Chakraborty, N.; Klein, M.; Swaminathan, N. Effects of Lewis number on the reactive scalar gradient alignment with local strain rate in turbulent premixed flames. Proc. Combust. Inst. 2009, 32 (1), 1409− 1417. (46) Becker, H.; Monkhouse, P.; Wolfrum, J.; Cant, R.; Bray, K.; Maly, R.; Pfister, W.; Stahl, G.; Warnatz, J. Investigation of extinction in unsteady flames in turbulent combustion by 2D-LIF of OH radials and flamelet analysis. Symp. (Int.) Combust., [Proc.] 1991, 23, 817− 823. (47) Mueller, M.; Kim, T.; Yetter, R.; Dryer, F. Flow reactor studies and kinetic modeling of the H2/O2 reaction. Int. J. Chem. Kinet. 1999, 31 (2), 113−125. (48) Li, J.; Zhao, Z.; Kazakov, A.; Chaos, M.; Dryer, F. L.; Scire, J. J. A comprehensive kinetic mechanism for CO, CH2O, and CH3OH combustion. Int. J. Chem. Kinet. 2007, 39 (3), 109−136. (49) Pearson, K. Note on Regression and Inheritance in the Case of Two Parents. Proc. R. Soc. London 1895, 58, 240−242.

K


Effects of the Equivalence Ratio and Reynolds Number on Turbulence

Recommend Documents