air flames at high Karlovitz numbers

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 2 0 2 1 6 e2 0 2 3 2

Available online at www.sciencedirect.com

ScienceDirect journal homepage: www.elsevier.com/locate/he

Direct numerical simulation of lean premixed CH4/air and H2/air flames at high Karlovitz numbers Henning Carlsson*, Rixin Yu, Xue-Song Bai Division of Fluid Mechanics, Lund University, 221 00 Lund, Sweden

article info

abstract

Article history:

Three-dimensional direct numerical simulation with detailed chemical kinetics of lean

Received 5 July 2014

premixed CH4/air and H2/air flames at high Karlovitz numbers (Ka ~ 1800) is carried out. It

Received in revised form

is found that the high intensity turbulence along with differential diffusion result in a

12 September 2014

much more rapid transport of H radicals from the reaction zone to the low temperature

Accepted 22 September 2014

unburned mixtures (~500 K) than that in laminar flamelets. The enhanced concentration of

Available online 24 October 2014

H radicals in the low temperature zone drastically increases the reaction rates of exothermic chain terminating reactions (e.g., H þ O2þM ¼ HO2 þ M in lean H2/air flames),

Keywords:

which results in a significantly enhanced heat release rate at low temperatures. This effect

Turbulent premixed combustion

is observed in both CH4/air and H2/air flames and locally, the heat release rate in the low

Direct numerical simulation

temperature zone can exceed the peak heat release rate of a laminar flamelet. The effects

High Karlovitz number

of chemical kinetics and transport properties on the H2/air flame are investigated, from

Detailed chemistry

which it is concluded that the enhanced heat release rate in the low temperature zone is a

Differential diffusion

convectionediffusion-reaction phenomenon, and to obtain it, detailed chemistry is essential and detailed transport is important. Copyright © 2014, Hydrogen Energy Publications, LLC. Published by Elsevier Ltd. All rights reserved.

Introduction Turbulent premixed flames have for a long time been a focus of research due to their extensive industrial usage, e.g., in gas turbines and internal combustion engines. At low intensity turbulence conditions, the structure of the flame is known to be similar to that of a laminar premixed flame. The flame is characterized by a multiple-zone structure: starting from the unburned mixture, a preheat zone, where there is no significant chemical reaction owing to low temperature, a reaction zone, where fuel is converted to combustion products, and

finally a post-flame zone where temperature is high but the chemical reactions are in equilibrium. For hydrocarbon flames, the reaction zone can be further divided into a fuel consumption layer, which is also known as the inner layer, where fuel is converted to CO and H2, and an oxidation layer downstream the inner layer, where CO and H2 are further oxidized to CO2 and H2O [1,2]. For premixed combustion in a turbulent environment the flame structure is highly depending on the flame/turbulence interaction, which can be classified into different regimes [3] based on non-dimensional parameters such as Karlovitz number (Ka). The Karlovitz number is defined as the ratio

* Corresponding author. Tel.: þ46 46 222 92 76; fax: þ46 46 222 47 17. E-mail addresses: [email protected], [email protected] (H. Carlsson). http://dx.doi.org/10.1016/j.ijhydene.2014.09.173 0360-3199/Copyright © 2014, Hydrogen Energy Publications, LLC. Published by Elsevier Ltd. All rights reserved.

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 2 0 2 1 6 e2 0 2 3 2

between the chemical time scale, tc, and the shortest turbulence time scale, the Kolmogorov time scale, th. At high Karlovitz numbers, the Kolmogorov scale can be smaller than the chemical reaction time scale [3,4] and the small scale turbulence can disrupt the inner thin reaction layer, creating a distributed combustion region of broken and disconnected reaction zones. With the current trend, where many industrial combustion devices are pushed towards leaner operating conditions with subsequent thicker flames and lower laminar flame speed, the Karlovitz number in such applications tends to increase [5]. Very high turbulence intensity can be generated by high-swirl combustor nozzles, which are often used for flame stabilization in gas turbines. Stopper et al. [6] studied an industrial gas turbine burner, where Ka ~ 230 was obtained for a power output of ~1 MW. With increasing demand on power, the conditions of such swirl-operating gas turbines are pushed further into the distributed reaction zone regime. Strakey et al. [7] showed in measurements of an experimental swirl combustor that turbulence intensity of 150 times the laminar flame speed could be obtained (Ka ~ 600). Due to the existence of fine length and time scales in such flames it is difficult to resolve all these fine scales in either experiments or numerical simulations, and as such, the reaction zone structures at high intensity turbulence high Karlovitz number conditions are not well understood. Turbulence/chemistry interaction is typically two-fold: on one hand it is well known that the chemical reaction can affect the flow field through gas expansion from heat release, on the other hand turbulence may also affect chemical reactions both indirectly and directly. Turbulence will typically induce an overall higher effective diffusion rate, which causes the mean structure to develop into a thickened flame brush with an increased mean flame speed; the chemical reactions are then adjusted indirectly to respond to the altered flame profile. Another mechanism of turbulence affecting chemistry may be by more direct disturbance of the chemical pathway through a sequence of elementary reactions. In a turbulent environment, the temperature field and the concentration fields of various major and intermediate species can be significantly different from the manifold of a laminar flame. For high Karlovitz number flames, such scattering may become large enough to cause the reaction to proceed through unconventional pathways. A chemical kinetic mechanism, consisting of detailed elementary reactions, usually contains two groups of important reactions: the chain branching and chain termination reactions [8]. The chemical pathway can be affected by the balance between the chain branching and termination reactions; in a detailed chemistry mechanism, there may exist multiple chain branching/termination reactions. In case turbulence is intense, we may hypothesize (and as we will show later) that some of the species involved in (usually negligible) chain branching/termination reactions can be transported by rapid convection to a level to enhance/suppress those reactions. This effect we refer to as deviation from conventional chemical pathway. A goal of this work is to identify such change of chemical pathway in high Karlovitz number turbulent premixed flame. If the phenomena are verified, it is also important to understand through which elementary reaction such change proceeds and how such change will be

20217

affected by the use of different fuels, chemical kinetic mechanisms, and transport properties. Previous studies of high Karlovitz number combustion systems have been performed with both experimental and computational approaches. Relying on planar laser induced fluorescence (PLIF) technique, several experimental works [9e16] have been performed for Karlovitz number around (and above) 100; most of these studies were aimed at identifying the effect of local quenching and reignition. For numerical studies of fundamental turbulent combustion processes, direct numerical simulation (DNS) [17e21], with the capacity of resolving the full turbulent flow spectrum, has gained increased popularity. Using simplified chemistry (with a model based on first order Arrhenius kinetics and a simplified reactionediffusion model), Poludnenko and Oran [22,23] carried out three-dimensional (3D) DNS of a premixed H2/air flame at Karlovitz number above 100; the work was focused on examination of the flame structures. To describe the details of the multiple reaction zone nature in flames, relatively detailed elementary reactions are needed; it is therefore essential to perform DNS studies with detailed chemistry. Few such studies have however been performed, due to the high demands on computational power and on robust/accurate numerical algorithm. The application of 3D-DNS with detailed chemistry is limited to combustion systems studying autoignition [21], premixed or partially premixed combustion at low Karlovitz number [24e26], and jet or shear layer flames [27,28]. A recent review on this topic is given by Chen [29]. For the study of premixed flames at high Karlovitz number, Aspden et al. [30e32] performed 3D-DNS to examine the effect of differential diffusion in flames of different fuels with detailed chemistry. Savre et al. [33] also performed a DNS study at high Karlovitz number in a two-dimensional (2D) configuration; the study was focused on the phenomenon of quenching/reignition. It should be noted that a detailed analysis of the balance between chain branching and termination reactions is performed in an early 2D-DNS work by Echekki and Chen [34], where the investigated combustion system is of autoignition type. The above-discussed challenge is addressed in this work by performing the state-of-art 3D-DNS of turbulent premixed flames at high Karlovitz number with a detailed chemistry solver. This study focuses on the direct effect of turbulence on chemistry. Two different fuels are investigated in 3D configurations with a similar Karlovitz number: an H2/air flame and a CH4/air flame. While 3D-DNS is time consuming, we also perform additional 2D-DNS to examine the influence of the chemical kinetic mechanisms and transport properties. This article is organized as follow. The numerical method is described briefly in Section Numerical method. In Section Computational cases, the computational setup and the studied cases are stated. Section Results and discussion is divided into three parts, where the flame structure at the high Karlovitz number condition is firstly discussed. In the second part of Section Results and discussion, changes from conventional low Karlovitz number flames are investigated and analyzed in detail and in the last part of Section Results and discussion, amplified heat release rate due to the change in the flame structures is quantified. In Section Conclusions, the main conclusions of this work are given.

20218

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 2 0 2 1 6 e2 0 2 3 2

Numerical method An in-house DNS solver [35] for low Mach number reacting flows is employed taking into account detailed transport properties and stiff chemistry involved in detailed chemical kinetics. The pressure is split into two parts, a thermodynamic pressure, P(t), which in the equation of state determines the mixture density, and a hydrodynamic pressure, p(x,t), which enters into the conservation equations for momentum. x is spatial coordinate and t is time. The turbulent reacting flow is governed by the conservation equations of mass, momentum, species, and energy, as well as the equation of state:

Integration of the chemical reaction rates involved in detailed chemical kinetic mechanisms is performed using an efficient stiff solver, the DVODE solver [40]. In this study, the CH4/air mechanism by Smooke and Giovangigli [41], employing 16 species and 25 reactions (10 reversible and 15 irreversible reactions), and the H2/air mechanism by Li et al. [42], employing 9 species and 21 reversible reactions, are used for detailed chemistry investigations. A single step mechanism by Marinov et al. [43] is also used for comparison. The transport properties, e.g., species diffusion coefficients, thermal conductivity, and viscosity, are mixture averaged based on the individual species obtained from Chemkin thermodynamic database [44].

8   1 vr vr vui > > > þ ui ¼ > > r vt vx vxi > i > >  >     > > vui vui vp v vui vuj 2 vuk > > þ u r þ þ  d ¼  m > j ij > > vxi vxj vt vxj vxj vxi 3 vxk > > >   > > < vYk vrYk Vk;j vYk þ uj þ u_ k r ¼ vt vxj vxj > > > ! >  N     > N > X P > vT vT vT v vT > > _ þ u u rC h  r C Y V þ ¼  l p k k p;k k j k;j > > vt vxj vxj vxj vxj > k¼1 > k¼1 > > > > N Y > P > k > > : P ¼ rℛT k¼1 Wk

where r, m, Cp, and l are the mixture averaged density, dynamic viscosity, heat capacity, and conductivity of the mixture, respectively. ui is velocity component i and dij is the Kronecker delta. Yk, u_ k , hk, Cp,k, and Wk are, respectively, mass fraction, reaction rate, specific enthalpy, heat capacity, and molar mass of species k. T is temperature, N is total number of species, and ℛ is the universal gas constant. Species diffusion, Vk,j, is given by: Yk Vk;j ¼ Dk

N X vYk vYi  Yk Di ; vxj vxj i¼1

(2)

where Dk is mixture averaged diffusion coefficient of species k. The temporal integration of convectionediffusion-reaction (CDR) terms is often performed using an operator splitting technique. The 2nd order symmetrical Strang splitting algorithm [36] is used, where the integration of one full step of the stiff reaction rates is placed in-between two half-step integrations of the diffusion and convection terms. The integration of the diffusion terms is further split into multiple substeps of explicit integrations with the sub-step size satisfying the diffusion stability limit. The discretization of all spatial terms is made using a 6th order central difference scheme except for the convective terms in the species mass and energy equations, where a 5th order weighted essentially nonoscillatory (WENO) scheme [37] is used to avoid unphysical numerical oscillations. The spatial/temporal accuracy of the numerical scheme has been verified with grid/time-step dependency tests [35]. The solver has been applied to studies of various reacting flow systems [21,38,39]. Detailed description of the numerical method is found in Ref. [35].

(1)

Additionally, it is worth mentioning that the code is written in a vector/tensor form, where a switching of dimension is easily done while maintaining an identically performing solver for 1D, 2D, and 3D setups.

Computational cases Two 3D-DNS cases are considered in this study. Initially, a premixed laminar planar flame is placed in an inflow-outflow configuration with a constant mean axial flow speed from an inlet down to an outlet. In both cases the mean streamwise flow velocity is set relatively high as 10 m/s to avoid negative axial velocity at the outlet boundary and to transport turbulent kinetic energy imposed at the inlet downstream to the flame, thus allowing for shorter transition from decaying turbulence to statistically stationary state. Furthermore, the mean flow velocity has to be kept sufficiently low to maintain the flame within the computational domain during the simulation time. Periodic boundary conditions are applied at the two cross-flow boundaries and a convective outflow condition is used at the outlet. Combustion is therefore taking place in a constant pressure system. As shown in Fig. 1, the computational domain is a rectangular prism with a streamwise (x-direction) length of 10 mm and a spanwise (y-direction) width and lateral (z-direction) height of the prism of 5 mm, yielding a prism aspect ratio of 2:1:1. The fluctuating velocity at the inlet plane is taken from a cubic box of periodic turbulence generated by synthesizing Fourier waves [45,46] satisfying a prescribed turbulence energy spectrum [47]. The

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 2 0 2 1 6 e2 0 2 3 2

Fig. 1 e A schematic illustration of the 3D computational domain with the arrow indicating the flow direction, and the two cross-sections showing the instantaneous temperature field for the hydrogen/air flame at t/tl ¼ 1.

same synthesized turbulence is also used to model the initial flow field. An initial planar flame is imposed around x ¼ 4 mm using the profile of species concentrations and temperature computed for a steady one-dimensional (1D) laminar flame with the same chemical kinetics and transport properties as the 3D DNS cases. The 3D domain is discretized on a uniform grid of 1024  512  512 (in x, y, and z directions) cells of cubic shape, yielding a total of approximately 268 million cells and a grid resolution of Dx ¼ 9.8 mm. The synthesized turbulence field is determined by specification of the turbulent intensity, u', and integral length scale, l (or L11), through which the Kolmogorov length scale, h, is obtained. In all cases the integral length scale is set to 625 mm. Table 1 lists all the computational cases studied in this work. In cases 1 and 2, two 3D-DNS cases are listed for two typical fuel/air mixtures, CH4/air (equivalence ratio Ф ¼ 0.6) and H2/air (Ф ¼ 0.4). For all cases, the thermodynamic pressure is 1 atm and the temperature of the unburned mixture is 298 K. The chemical kinetic mechanisms by

20219

Smooke and Giovangigli [41] and Li et al. [42] together with the corresponding detailed transport properties are used respectively for the two flames. For the above two mixtures of CH4/ air and H2/air, the laminar flame speed (SL) are computed to be 0.12 m/s and 0.22 m/s, respectively, and the thermal thickness of the flame (dL) are 1.12 mm and 0.74 mm, respectively. With the above grid resolution, the thermal thickness of the flames (defined as the distance between the two locations of 10% and 90% of the temperature rise in the corresponding 1D flames) are resolved by ca. 75e115 computational cells and the smallest Kolmogorov length scales are resolved with two or more grid cells. For all cases the Karlovitz number based on the initial turbulence field is set the same, Ka ¼ 1860, which results in a lower turbulence intensity (15 m/s) in the CH4/air flame than the one (32 m/s) in the H2/air flame. The characteristic chemical time scale, tc, is estimated as tc ¼ dL/SL. The characteristic diffusive time scale, tD, is estimated as the diffusive time scale of the fuel: tD ¼ ½Df VðYf =Yf0 Þ$VðYf =Yf0 Þ1 at Yf =Yf0 ¼ 0:5, where Df is the diffusivity of the fuel species (DCH4 z 25$106 m2/s, DH2 z 127$106 m2/s), Yf is the mass fraction of the fuel, and Yf0 is the mass fraction of the fuel at unburned mixture. The characteristic turbulence time scale, tl, is estimated by tl ¼ l=u0. These estimations yield that tc ¼ 9.3 ms and 3.4 ms, tD ¼ 6.8 ms and 2.1 ms, and tl ¼ 42 ms and 20 ms, respectively for cases 1 and 2. The time steps, Dt, employed in the simulations were Dt ¼ 10 ns and 5 ns for cases 1 and 2, respectively. tl is known as the large eddy turnover time, which is resolved here with about 4000 time steps for both cases. The Kolmogorov time scale is 4.98 ms and 1.78 ms, respectively for cases 1 and 2. These are resolved with above 300e500 time steps. The above discussed diffusive and flame time scales are based on the corresponding 1D laminar flames. In flames with Karlovitz number greater than 100, by definition, the time scale of laminar flame is significantly larger than the Kolmogorov time scale. It should be pointed out that while the Kolmogorov scales used here are relevant to those in high Karlovitz number flames found in engineering applications, hence the physics studied are essentially representative for the practical flames. The integral length and time scales allowed in the

Table 1 e Simulation properties. Detailed computation of the diffusion coefficients corresponds to mixture averaging based on the individual species obtained from Chemkin thermodynamic database and unity Le refers to replacing the diffusion coefficients by Dk ¼ l=ðr$Cp Þ. For both cases, the thermodynamic pressure is 1 atm and the temperature of the unburned mixture is 298 K. Case Dim. Comp. Mech. Diff. Coeff. Ф SL [cm/s] dL [mm] u' [m/s] l [mm] tl [ms] Re h [mm] Ka Dx [mm]

1

2

3

4

5

3D CH4/air Smooke [41] Detailed 0.6 12 1.12 15 625 41.7 70 26 1860 9.8

3D H2/air Li [42] Detailed 0.4 22 0.74 32 625 19.5 120 17 1860 9.8

2D H2/air Li [42] Detailed 0.4 22 0.74 32 625 19.5 120 17 1860 9.8

2D H2/air Li [42] Le ¼ 1 0.4 39 0.42 32 625 19.5 120 17 910 9.8

2D H2/air Marinov [43] Detailed 0.4 34 0.40 32 625 19.5 120 17 1140 9.8

20220

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 2 0 2 1 6 e2 0 2 3 2

present 3D DNS domains (and also in pervious 3D DNS studies [30e33]) are smaller than that found in engineering flames. Such DNS result reveals the physics of local flame/turbulence interaction but it may not be used to explain the large scale behavior of an engineering flame. 3D DNS studies are of essential importance to capture the physical energy cascading from large to small scales and to obtain a quantitative capture of the interaction between turbulence and flames. However, for qualitative comparison with the 3D DNS results, three 2D DNS cases are also carried out. For the geometry and grid configurations, the 2D cases 3e5 are similar to that of the 3D H2/air flame by omitting the lateral (z) direction. Furthermore, in case 4 a model species diffusion coefficient is used for the H2/air/products mixture, which assumes that Dk ¼ l/(r$Cp), i.e. without differential diffusion (with unity Lewis number). In case 5 the one-step H2/air chemistry mechanism by Marinov et al. [43] is employed. Cases 4 and 5 are considered to delineate the effect of differential diffusion and the impact of detailed chemical kinetic mechanism on the flames, respectively. This kind of “numerical experiment” is a powerful asset of numerical studies, which allows us to study idealized model situations that may not be realizable in physical experiments. The computational cost for the two 3D-DNS cases is about 2.11 million core-hours, run on 8192 Cray CPU Sandy Bridge clocked at 2.3 GHz (AVX).

Results and discussion 3D flame structures and global statistics For turbulent planar flame simulation in an inflow-outflow configuration, the flame is initially imposed on a synthesized turbulence flow field of uniform spatial turbulence intensity. With no extra turbulence production mechanism the turbulence will decay along the streamwise direction, while at the inlet, turbulence is maintained by using the turbulence inflow condition described in Section Computational cases. At a downstream location, e.g., in front of the flame, the initially imposed turbulence will experience a decay before it reaches its statistically stationary state. Quantification of the actual turbulence intensity of the unburned mixture interacting with the flame is therefore necessary. Fig. 2 shows the temporal evolution of a volume averaged mean turbulence intensity (u0un ), integral length scale (lun), and Kolmogorov length scale

(hun) for the unburned mixture and a corresponding effective Karlovitz number (Kaeff) for cases 1 and 2. For an instantaneous turbulent flow field, u(x,t), u0un is computed as: u0un ðtÞ ¼

1 Xun

X Zun

h i1=2  〈uðx; tÞ  〈uðx; tÞ〉yz 2 〉yz dx;

(3)

0

where Xun is taken to be 3 mm, which represents the flow in the unburned side of the flame. For an arbitrary function j(x), the operator 〈$〉yz , for performing spatial averaging in the y-z plane, is defined as: 〈j〉yz ¼

1 ∬ jðxÞdydz; Ly Lz Ly Lz

(4)

where Ly and Lz are the domain length in the spanwise and lateral directions, respectively. The effective Karlovitz number (Kaeff) is estimated according to: Kaeff ¼

 0 3=2  1=2 uun dL : SL lun

(5)

By assuming statistically homogeneity and isotropy in the y and z directions, an x-varying integral length scale may be estimated as the integration of the longitudinal autocorrelation function from the instantaneous fluctuating velocity field in the y-z plane. This scale is then averaged over 0 mm  x  3 mm to obtain an estimate of lun. It can be seen in Fig. 2 that during 0  t  2tl, lun increases by a factor of ~1.2 for the CH4/air flame and by a factor of ~1.8 for the H2/air flame. hun increases during this time by a factor ~1.7e1.9. It is worth mentioning that the estimated hun is the smallest possible scale; however, since the initial turbulence field is generated from synthesized Fourier waves and is not adapted to physical cascading, it will be shown later that during the very initial stage of the simulation the size of the smallest scale is in fact decreasing before increasing at later stage. The turbulence intensity in the unburned mixture decays in both cases, while the turbulence in the H2/air flame case starts with a larger intensity than in the CH4/air flame, as a consequence of the thinner flame and higher flame speed of the H2/air flame while setting the same Karlovitz number as the CH4/air flame. The effective Karlovitz number for both flames decays during the first two integral time scales from the initial value of 1860 to around 700. In this work most analysis is performed within 2.5 turbulent integral time during which the studied flame is reaching a statistical stationary state while maintaining sufficiently high Ka.

Fig. 2 e (a) Karlovitz number and turbulent intensity of cases 1 and 2 based on unburned quantities with time and (b) integral and Kolmogorov length scales based on unburned quantities.

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 2 0 2 1 6 e2 0 2 3 2

Fig. 3 shows the scatter plot and joint probability density function (JPDF) of local equivalence ratio and temperature at various instances of time during the flame development process. The color bar at the bottom row of the figure further displays the JPDF of local equivalence ratio and temperature at t/tl ¼ 1.0. The solid reference lines correspond to the distributions in the corresponding steady 1D laminar flames. The local equivalence ratio, ɸ, is defined as [48]: f¼

YH =ð2$WH Þ þ 2$YC =WC ; YO =WO

(6)

where YH, YC, and YO are mass fractions of hydrogen, carbon, and oxygen elements, respectively, and WH, WC, and WO are molar masses of hydrogen, carbon, and oxygen elements, respectively. Since the initial flame is set as the 1D laminar flame, the local equivalence ratio at t/tl ¼ 0 is a unique function of temperature. Lower local equivalence ratio is seen in the preheat zone and the earlier reaction zone due to differential diffusion effects. For the CH4/air flame the minimal equivalence ratio is at T about 1000 K with deviation of equivalence ratio from the unburned mixture about 0.03; whereas, owing to stronger differential diffusion, the minimal equivalence ratio in the H2/air flame is at temperature about 500 K with a much larger deviation of local

20221

equivalence ratio about 0.15. It is interesting to note that, when the flame evolves in a high Karlovitz turbulent flame, the distribution in the ɸ-T space becomes scattered, primarily towards the direction of lowering the deviation of local equivalence ratio from that in the unburned mixture. This result is consistent with the results of Aspden et al. [30] and it is attributed to that turbulence transports species and energy at the same rate, thereby yielding a unity effective Lewis number. In the study of a fully developed turbulent flame brush, Aspden et al. reported a scattering of cells in the ɸ-T space around mean value of equivalence ratio of that of the unburned mixture (Fig. 4 in Ref. [30]); whereas in the present flame, where a developing turbulent flame is studied, the mean equivalence ratio at a given temperature is close to that in the initial 1D flame. This implies that, even at high Karlovitz number flames, the local structure of a developing turbulent flame is still affected by local differential diffusion, which is owing to transport on the molecular level. Results from a 2D study of CH4/air flame (Fig. 6 in Ref. [33]) are consistent with the present finding. Furthermore, it is well known that in lean H2/air flames thermaldiffusive instability can cause cellular flame wrinkles [4,39,49]. To describe the flame structure, a reaction progress variable, c, based on fuel mass fraction is defined as:

Fig. 3 e Scatter plot (first three and fifth rows from top) and joint probability density function (JPDF) of T and ɸ (fourth row from top) for case 1 (left column) and case 2 (right column) at t/tl ¼ 0.25, 0.5, 0.75, 1.0, and 2.0, respectively, from top row to bottom row. Solid (green) lines correspond to the distributions in 1D flames. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

20222

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 2 0 2 1 6 e2 0 2 3 2

Fig. 4 e HRR* at iso-surfaces conditioned at c ¼ 0.5 for cases 1 (a) and 2 (b).

  . cðx; tÞ ¼ Yf0  Yf ðx; tÞ Yf0  Yfb ;

(7)

where Yfb is the fuel mass fraction in the burned mixture at equilibrium. According to the above definition, c ¼ 0.1 corresponds to the initial part of the preheat zone; c ¼ 0.5 corresponds to the central part of the preheat zone, and c ¼ 0.96 corresponds to a position close to peak heat release rate position in a 1D laminar flame solution for both case 1 and case 2. 3D iso-surfaces conditioned on these values are therefore used for the following discussion of the flame structure. Fig. 4 shows normalized heat release rate per unit mass, HRR*, at iso-surfaces corresponding to c ¼ 0.5 for cases 1 and 2 extracted at times t/tl ¼ 0.5, 1.0, 1.5, 2.0, and 2.5. The heat release rates per unit mass, HRR, are normalized by their corresponding peak HRR value (HRRpL) in the initially imposed 1D laminar flame profiles. For the CH4/air flame, HRRpL ¼ 2.15 GJ kg1 s1 and for the H2/air flame, 2.20 GJ kg1 s1; for both flames the peak heat release rates are found at about c ¼ 0.96. The flames are developing during 0  t/ tl ( 1, where the flame wrinkling is increasing. At t/tl > 1.5 the flame wrinkling approaches a nearly statistically stationary level along with the initially decaying turbulence approaching a nearly statistically stationary state (cf. Fig. 2). It appears that

heat release rate on the c ¼ 0.5 iso-surface is not homogeneous, with most of the surface showing a lower heat release rate compared with the peak heat release rate in the corresponding laminar flame (HRR* < 1), whereas certain parts of the surfaces is having an enhanced heat release as compared with the laminar flame (i.e. HRR* > 1), in particular in the H2/ air flame. To further examine the reaction layer structure, the HRR* fields together with iso-vorticity lines at t/tl ¼ 0, 0.2, 0.4, 0.6, 1.0, and 2.0 for cases 1 and 2 are shown in Figs. 5a and b, respectively, on a 2D x-y plane (in the region 2 mm
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 2 0 2 1 6 e2 0 2 3 2

20223

Fig. 5 e Snapshots of HRR* along with vorticity lines at t/tl ¼ 0, 0.2, 0.4, 0.6, 1.0, 2.0 for case 1 (a) and case 2 (b). The vorticity lines are given at 105, 2∙105 (solid lines) and at ¡105, ¡2∙105 (dashed lines) for case 1 and at 4∙105, 8∙105 (solid lines) and at ¡4∙105, ¡8∙105 (dashed lines) for case 2.

energy cascading of the synthesized initial turbulence field, which was based on the prescribed energy spectrum [47]. This cascading tends to change the shape of energy spectrum to a more physical one with decreasing Kolmogorov scales during the initial stage of the simulation. Owing to the change of turbulent energy cascade, in particular the increase of the length of the Kolmogorov scale with time (cf. Fig. 2b), the HRR zone at t/tl ¼ 2 is smoother than that at t/tl ¼ 1. For t/tl < 1 the mean HRR zone is seen within x ¼ 0.4e0.5 mm, indicating a fuel consumption speed of the flame on the order of the mean flow speed (10 m/s); later when t/tl > 1 the HRR zone is pushed downstream by about 1 mm, indicating a lower fuel consumption speed than that of the mean flow velocity. Due to the great challenge in fully resolving the fine scales in high Karlovitz number flames in both experimental studies and numerical simulations, there is a lack of understanding of the structures of high Karlovitz number flames. Turbulence in high Karlovitz number flames have been assumed to distort the structures of the flames in such a way that there will be very low spatial variation of composition, temperature and reaction rates in the reaction zones. Hence, high Karlovitz number flames may be considered as a well-stirred reactor, cf. the discussions of Aspden et al. [32] and the references therein. Subsequently, one may also expect that differential diffusion and thermo-diffusive effects on high Karlovitz number flames are lower than those on low Karlovitz number flames. To improve the understanding of high Karlovitz number flames we look into the details of the reaction zones of the H2 flames, i.e. a zoomed region indicated by the square

in the lower right corner in Fig. 5b at t/tl ¼ 1.0, which is displayed in Fig. 6 where the HRR* is shown along with the local equivalence ratio, the normalized species mass reaction, Yi* , of H2, H, OH, HO2, H2O2, and reaction rates of species H2 and H. The species mass fractions are normalized as: . p Yi* ¼ Yi YiL ; p

(8)

where YiL is the peak mass fraction of species i in the initial 1D laminar flame profile. The dashed line in Fig. 6 represents a relative low temperature of 1000 K and the solid lines represent HRR* ¼ 1. The H and OH distributions show that there are fine structures in the reaction zone; a peninsula of H and OH radicals (marked as A in the H and H2 field) is protruding to the unburned side of the flame, whereas a peninsula of the fuel/ combustion intermediates (e.g. HO2 and H2O2) is protruding to the burned side of the flame (marked as B in the H and H2 field). In the center of peninsula A the reaction rate of H2 is large whereas at the edge of this structure the consumption rate of H is large where the mass fraction HO2 is also high. In the peninsula of B the reaction rate of H2 is low but the consumption rate of H radicals is high, which corresponds again to the high mass fraction of HO2. The spatial variations of the species and reaction rates are large, which is far from that described in the idealized well-stirred reactor model. Fig. 6 shows that a significantly enhanced HRR (corresponding to HRR* > 1) is found locally and that the T ¼ 1000 K line is positioned on the burned side of these positions,

20224

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 2 0 2 1 6 e2 0 2 3 2

* * * * * Fig. 6 e Snapshots of HRR*, YH , YH , YOH , u_ H2 , u_ H , YHO , YH , and ɸ at t/tl ¼ 1.0 in the H2/air flame (case 2), in the zoomed 2 2 2 O2 region of the snapshot as indicated in Fig. 5. The dashed line corresponds to T ¼ 1000 K, and the solid lines correspond to HRR* ¼ 1.

indicating a high heat release rate at rather low temperature regions. The distributions of H and OH radicals show an enhanced concentration of H locally, whereas the OH concentration does not show any significant increase compared

with the 1D peak reference value. Furthermore, the regions where the concentration of H radicals exceeds its peak 1D value do not show a clear correlation with the high HRR regions. In 1D H2/air flame profile the intermediate species HO2

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 2 0 2 1 6 e2 0 2 3 2

and H2O2 can serve as “markers” of preheat zone. The distributions of these two species are shown in Fig. 6. It is interesting to note that the concentration of HO2 is locally enhanced, whereas the concentration of the subsequent species H2O2 is generally lowered. In conventional flamelet combustion regime, lean H2/air flames are well known for being strongly susceptible to differential diffusion effects through the classical thermal-diffusive instability, which results in “focusing” of the hydrogen concentration in convex regions [49]. Such effect is equivalent to an enhanced local equivalence ratio, which may result in an overall enhancement of chemical reaction and consequently heat release rate. Fig. 6 indicates that indeed locally high equivalence ratio is found in the convex region of the flame (protruding to the unburned mixture), and local low equivalence ratio is seen in the concave region of the flame. However, contrary to the flames in the laminar flamelet regime, the regions of very large heat release (HRR* > 2) do not show a correlation with high values of the local equivalence ratio. In fact, the highest HRR* shown in Fig. 6 is at the lowest local equivalence ratio regions where the temperature is below 1000 K.

20225

To quantify the relative effect of different transport and reaction terms within a high Karlovitz number flame, convection, C ¼ uj $vYk =vxj , diffusion, D ¼ r1 $vðrYk Vk;j Þ=vxj , and reaction, u_ k =r terms in the transport equations for the mass fraction of the fuel species (k ¼ CH4 and H2, respectively) are presented in terms of JPDFs verses temperature in Fig. 7 for cases 1 and 2 at t/ tl ¼ 1.0. It is shown in Fig. 7 that at high Karlovitz number conditions, the convection and diffusion terms are about one order of magnitude lager than the reaction terms, whereas the convection term is in general larger than the diffusion term. The different terms are corresponding to the inverse of the effective respective time scales for the convective and diffusive transport and chemical reactions of reactive species. It is therefore confirmed that the scalar transport time scales are (locally) both much smaller than the reaction time scales. This is likely responsible the effect responsible for the existence of the long and narrow protruding structures of unburned fuel/air/combustion intermediate in the flames, c.f. Fig. 6, peninsula B. The relatively fast convection by small structures of unburned mixture into the burned side of the flame causes the formation of long and narrow protruding structures into the flame. This, in

Fig. 7 e Transport and reaction terms of the fuel species verses temperature for cases 1 and 2. C: convection term, D: diffusion terms, and R: reaction term.

20226

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 2 0 2 1 6 e2 0 2 3 2

turn, results in large gradients of composition and therefore molecular diffusion of the species. The relatively low reaction rates allow for the survival of such structures deep into the burned side of the flame. By carefully inspecting the HRR zone in Fig. 5 at t/tl ¼ 1 one can find long and narrow “jet-like” highly reactive mixtures with high heat release rate protruding to the unburned mixture, e.g., at x ¼ 3.4 mm and y ¼ 3.1 mm.

Identification of change of chemical reaction path at high Karlovitz number From examination of the flame structure in Fig. 6 we notice an interesting occurrence of high release rate in a relatively low

temperature region. To further examine this phenomenon, in Fig. 8 we show JPDFs of the instantaneous HRR*, YH, u_ H =r, YHO2 , and u_ HO2 =r with temperature for cases 1 and 2 at t/tl ¼ 1, respectively. The solid lines represent the 1D reference results. For the CH4/air and the H2/air flames, the 1D laminar flame profiles show a negligible HRR* in the low temperature region. e.g., T < 800 K for the CH4/air flame and T < 600 K the H2/air flames, which may be referred to as the preheat zones in the corresponding laminar flames. In the high Karlovitz number turbulent flames, however, significantly high HRR* can be found in these low temperature regions. For both flames, HRR* in the low temperature range can even exceed unity

Fig. 8 e Joint PDFs of HRR*, YH, u_ H =r, YHO2 , and u_ HO2 =r for case 1 (column 1) and case 2 (column 2). Solid lines correspond to the distributions in a 1D flame and dashed lines correspond to the mean quantities conditioned on the given temperature.

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 2 0 2 1 6 e2 0 2 3 2

20227

(especially in the H2/air flame, in which HRR* can reach a value as high as 3 at T ¼ 500 K). This is an interesting phenomenon for high Karlovitz number combustion that deserves further investigation. Different species concentrations are examined to locate the source of heat release rate in the low temperature zone. In Fig. 8 it can be seen that the mass fractions of both H and HO2 in the low temperature zone are significantly increased. An increased concentration of H-radicals in the low temperature zone has also previously been reported for a 2D CH4/air flame [33]. Fig. 8 further shows that the reaction rates of these species are drastically enhanced in the low temperature zone, with a large production rate of HO2 and a large consumption rate of H. In typical laminar flame these rates are negligible in the low temperature zone, cf. Fig. 8. The consumption of H radicals along with the formation of HO2 at low temperature is associated with the exothermal chain-terminating (radical recombination) reaction: H þ O2 þ M/HO2 þ M;

(9)

where M is a third body species. To examine the relative importance of the above elementary reaction within the total heat release, we define a normalized heat release rate for an elementary reaction, HRR*ER , as:  p (10) HRR*ER ¼ ðRER $Q=rÞ HRRL ; where subscript ER denotes elementary reaction, RER is the reaction rate of the given elementary reaction [mole$m3$s1] obtained through an Arrhenius expression, and Q is the heat released in the reaction [J$mole1]. The normalized heat release rate from reaction (9) in various computational cells is computed and plotted against its corresponding HRR* in Fig. 9 for cases 1 and 2 at t/tl ¼ 1.0 for T < 700 K. The solid line represents a reference condition that the total heat release rate in the flame at the given cell is equal to that from reaction (9), i.e. HRR*R(9) ¼ HRR*. For the H2/air flame it is evident that the major part of heat release rate in the low temperature zone is originated from the terminating reaction (9). The heat released in this reaction at low temperature (e.g., T < 700 K) is negligible in a 1D laminar flame (Fig. 8), or low Karlovitz number turbulent flames in the laminar flamelet regime. This is owing to the competition between the diffusion time for the H radicals from the high temperature reaction zone (cf. Fig. 8 in the region with T > 800 K where H radicals are formed through various reactions) to the low temperature zone and the local chemical reaction time during which H radicals are consumed through reaction (9). In laminar flames or laminar flamelets, the transport of H radicals within the flame is mainly due to molecular diffusion, which is slower than the chemical consumption of H radicals through reaction (9); as such, H radicals can exist only in the relatively high temperature zone near where they are formed (cf. Fig. 8). In high Karlovitz number turbulent flames, the effective diffusion time is significantly shortened due to turbulence convection. When the turbulent convective transport time of H radicals from the high temperature zone to the low temperature zone is comparable to the local chemical reaction time of reaction (9), it is clear that H radicals will be found in the low temperature region farther away from where it is generated.

Fig. 9 e Normalized heat release rate caused by elementary reaction (9) against HRR* for cases 1 and 2 at T < 700 K. Solid (red) lines represent one-to-one correspondence. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

For the CH4/air flame, the correlation between HRR*R(9) and HRR* is not as strong as that for the H2/air flame, cf. Fig. 9. This is likely due to the lower concentration of H radical and that the CH4/air chemistry contains several other important low temperature reactions, where the species HCO, CH2O, and CH3O play an important role. An exothermic reaction that shows a significant reaction rate at low temperature in the CH4/air flame is: H þ CH2 O/HCO þ H2 :

(11)

The normalized heat release rate caused by reaction (11) is plotted against HRR* in Fig. 10 for case 1 at t/tl ¼ 1.0 for T < 700 K. It can be seen that the low temperature reaction (11) is releasing a substantial amount of heat in the low temperature zone. The reaction (11) for the CH4/air flame is influenced by the high Karlovitz number conditions in the same way as reaction (9), where the transport of H radicals from the reaction zone to the low temperature zone has a significant effect on the reaction rate. It can be seen that the heat release rate in (11) can be slightly higher than the total net heat release rate due to the existence of endothermic reactions in the system. To further investigate the transport of H radicals and the heat release process in the low temperature region, we examine the combustion process in 2D configurations and the effect of chemistry and transport properties. Three additional 2D DNS cases are considered (cases 3e5, Table 1). Case 3 is

20228

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 2 0 2 1 6 e2 0 2 3 2

Fig. 10 e Normalized heat release rate caused by elementary reaction (11) against HRR* for case 1 at T < 700 K. Solid (red) line represents one-to-one correspondence. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

identical to case 2, except that the computational domain is 2D in case 3, which is to verify the above discussed low temperature zone behavior in 2D. Case 4 is a case with simplified transport properties (assuming a unity Lewis number for all species) to scrutinize the effect of differential diffusion and turbulence transport to the heat release rate in the low temperature zone. Case 5 is a case of single-step chemistry to investigate if the enhanced heat release rate in the low temperature zone can be captured with simplified chemistry. Fig. 11 shows the scatter-plots of HRR* in various computational cells as a function of local temperature for cases 3e5, respectively, at t/tl ¼ 1.0. Fig. 11a shows that HRR* obtained in a 2D domain is qualitatively similar to the corresponding 3D results shown for case 2 in Fig. 8, verifying that the enhanced heat release rate in the low temperature zone can also be obtained in 2D simulations. Cases 4 and 5 are therefore valid test cases. Fig. 11b shows that with unity Lewis number, and hence without the effect of rapid molecular diffusion of H radicals, the heat release rate in the low temperature zone is also enhanced; however, it is not enhanced to the same extent as that in Fig. 11a. From this result it is clear that the enhanced heat release rate in the low temperature zone is a combined, closely coupled process with both turbulence transport and molecular diffusion, where differential diffusion plays an important role. By a close examination of HRR* in Fig. 11a and b one can notice that molecular diffusion can affect also the high temperature zone with T > 1100 K. With unity Lewis number, hence neglected rapid diffusion of H and H2, the HRR* in the high temperature is essentially decreased to be lower than the laminar 1D HRR* profile. In the scatter plot the density of the scattering indicates the probability of the various mixture points achieving a given HRR* at a given temperature. From Fig. 11a and b it appears that, in the low temperature zone, a significant portion of the mixture points lie in the same domain in the HRR*-T space, which indicates that turbulence is responsible for the major transport of H radicals from the reaction zone to the low temperature zone. Fig. 11c shows that, when the detailed chemistry is simplified to a single step mechanism, no enhanced heat release rate in the preheat zone can be found, owing to the

Fig. 11 e Scatter plots of HRR* with temperature at t/tl ¼ 1 for cases 3e5 in aec, respectively. Solid (red) lines correspond to the distributions in a 1D flame and dashed (green) lines correspond to the mean quantities conditioned on the given temperature. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

absence of H radicals in the system. This result indicates that in high Karlovitz number flames, omit of detailed chemistry can lead to a significant loss of information of the flame. Furthermore, for different chemical kinetic mechanisms published in the literature, the constants in the Arrhenius rate formulations (the activation energy, pre-exponential factor, etc) for the same elementary reaction can be significantly different. The published chemical mechanisms were typically validated by matching certain global combustion properties, such as laminar flame speed or ignition delay time at certain operation ranges of temperature and pressure. A study is carried out on the elementary reaction rate of reaction (9) in different chemical kinetic mechanisms. Fig. 12 shows HRR*R(9) computed through five different detailed chemical kinetic mechanisms vs. temperature, based on a 1D steady flame solution obtained through the Li et al. mechanism. In addition to the Li et al. mechanism [42], results from the mechanism by Conaire et al. [50], SmookeeGiovangigli [41], GRI 3.0 [51], and the San Diego mechanism (UCSD) [52] are shown. It can be seen that HRR*R(9) is similar for the Li et al. and Conaire et al. mechanism, as they are both dedicated H2/air mechanisms. The UCSD mechanism show similar HRR*R(9) as the Li et al. and Conaire et al. mechanisms, whereas the SmookeeGiovangigli and GRI 3.0 mechanisms show higher heat release rate from

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 2 0 2 1 6 e2 0 2 3 2

20229

Fig. 12 e Normalized heat release rate from elementary reaction (9) in a 1D flame profile as a function of T, calculated employing different chemical kinetic mechanisms.

reaction (9) than the Li et al., Conaire et al., and UCSD mechanisms do. It is worth pointing out that although the three constants (temperature exponent n, the pre-exponential factor A, and the activation energy Ea) in the Li et al. and Conaire et al. H2/O2 mechanisms are significantly different, both the Li et al. and Conaire et al. mechanisms give similar predictions of HRR*R(9). This similarity suggests that the phenomenon of high heat release at low temperature may be observable if the Conaire et al. mechanism is adopted. For the CH4/O2 mechanisms the rates of reaction (9) from different mechanisms differ by a factor of 2. This may affect the quantitative prediction of HRR*R(9). Further validation of elementary reaction rate constants at low temperature might be needed to quantitatively predict the enhanced heat release phenomenon at low temperatures in high Karlovitz number flames. From the various 2D and 3D results one may conclude that for high Karlovitz number flames, the effect of turbulence motion is to broaden the reaction zone of the radical recombination reactions (e.g., reactions 9 and 11), resulting in deviation from the typical chemical pathway in a laminar or low Karlovitz number flame.

Characterization of heat release and length scales of the reaction zone structures at low temperature To quantify the phenomenon of enhanced heat release rate in the low temperature zone, heat release rate at iso-surfaces of c ¼ 0.1 for case 1 and c ¼ 0.5 for case 2 are extracted. These two different values of c are chosen since they correspond to a similar low temperature value of T z 440 K. It should be noted that conditioning on the same value of c for both flames would therefore not provide a consistent comparison. An area ratio, A*, is defined as: *

Z

A ða; c* Þ ¼



H HRRðxÞ  Ujc*

p a$HRRL dS

, Z dS;

Fig. 13 e A* with time for cases 1 (a) and 2 (b) at different a.

flame, the area of enhanced heat release rate is increasing. For the CH4/air flame, as indicated already in Fig. 8, the heat release rate at the given low temperature (440 K) is rather small. For the H2/air flame, it can be seen that A* almost reaches a value of unity for a ¼ 102 and above 70% for a ¼ 0.1. This indicates that on the majority of the iso-surfaces the local heat release rate can reach as high as 10% of the peak heat release rate in a corresponding 1D laminar flame. To quantify a length scale of the enhanced reaction rates at low temperature, Fourier analysis is performed in the spanwiseelateral planes. These planes have periodic boundary conditions and 2D Fourier transform is therefore applicable. A two-point correlation function, Rj, in the 2D yz plane for variable j is defined as: Rj ðrÞ ¼ 〈jðx þ rÞjðxÞ〉yz ;

(13)

where r represents a 2D vector in the y-z plane. The Fourier b j , are obcoefficients of the two-point correlation function, R tained from the Fourier transform: b j ðkÞ ¼ F ½Rj ðrÞ; R

(14)

where k is a 2D wavenumber vector and F is the 2D Fourier transform operator. By assuming isotropy in the y-z plane, a 1D energy spectrum, Ej(k), can be defined as:

(12)

Ujc*

where H is the Heaviside step function, a is a chosen constant, dS is a surface element, and U is the entire iso-surface of c ¼ c*. A* is plotted against time for cases 1 and 2 in Fig. 13a and b, respectively. During the initial development of the turbulent

Ej ðkÞ ¼

1 b j ðkÞdðjkj  kÞdk; ∬R 2p$k

(15)

where d(k) is the Dirac delta function. Fig. 14a shows an x-z plane along with a y-z plane conditioned on cyz ¼ 0.1 at t/ tl ¼ 1.0, at which reaction rate of H radicals is plotted. Fourier

20230

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 2 0 2 1 6 e2 0 2 3 2

Fig. 14 e (a) u_ H =r at an x-z plane and a y-z plane conditioned on cyz ¼ 0.1 at t/tl ¼ 1.0. The solid line is iso-line of temperature of 800 K and the dashed line for temperature of 1000 K. (b) Energy spectrum of u_ H =r and u_ H2 =r at cyz ¼ 0.1 and cyz ¼ 0.95. kh is the wavenumber of the Kolmogorov length scale.

transform analysis is performed for such periodic cross-flow planes. In Fig. 14a, a set of larger and smaller scales are observable, where the chain branching reactions (with positive reaction rate at high T) form H radicals and the chain terminating reactions (with negative reaction rate at low temperature) consume H radicals. Fig. 14b shows energy spectrums of the reaction rates of H and H2 conditioned on cyz ¼ 0.1 and cyz ¼ 0.95. The difference between the scales of the fuel consumption (u_ H2 =r) and scales of u_ H =r is evident at cyz ¼ 0.1 as two distinct scales are obtained for u_ H =r, whereas u_ H2 =r shows a monotonically decreasing energy with wavenumber. The small scale peak for u_ H =r is corresponding to a length scale of ~60 mm, corresponding to ~6Dx, which is just marginally larger than the Kolmogorov length scale. At cyz ¼ 0.95 the reaction rates energy spectrum of u_ H =r, and u_ H2 =r show similar behavior, where no small scale of amplified energy is obtained. The existence of small scales in the spectrum of the u_ H =r is consistent with the flame structure shown in Fig. 6, e.g., the structure of peninsula B. Careful inspection of the thickness of the reaction rate of H in this structure (more clearly seen as the region of Y*HO2 > 1) shows that the thickness of the high H reaction rate layer is about 58 mm, which is very close to the length scale of 60 mm identified in the spectrum of u_ H =r.

Conclusions Three-dimensional direct numerical simulations are performed to study turbulent premixed CH4/air and H2/air flames at high Karlovitz number and detailed flame structures are presented. The finest Kolmogorov flow scales and the thinnest layers of intermediate species are fully resolved and their interaction is captured. The main conclusions of this paper can be summarized as:  High heat release rate is identified at low temperature in both H2/air flame and CH4/air flame, which is induced by the high Karlovitz number employed. This effect is associated with a chain terminating reaction

(H þ O2 þ M / HO2 þ M), in which the elementary reaction rate is magnified by rapid H radical convective transport from the high temperature zone to the low temperature zone in the presented cases.  It is shown that the observed phenomenon will not appear if single step chemistry is employed and that differential diffusion has a significant effect on the presented phenomenon. The phenomenon is also quantified and proven significant in a major part of the flame. By the use of Fourier transform analysis significant fluctuation energy of H radical reaction rate is shown at a scale which is on the same order as the Kolmogorov length scale.  It is found that due to the relatively large convective transport at high Karlovitz number conditions long and narrow unburned mixture structures can protrude deep into the burned side of the flame and similarly long and narrow structures of reactive mixtures with high heat release rate can penetrate deep into the unburned mixtures. The reaction zones of the flames are thickened by turbulence, but remains connected and not distributed.  The composition gradient is high in the reaction zone of high Karlovitz number flames but the chemical reaction rate is an order of magnitude lower than the convective transport. This is the reason that the fine structures of unburned mixtures can penetrate deep into the flame without being consumed or the reacting mixtures to the unburned mixtures without being quenched. Corresponding to the existence of high composition gradient the diffusive transport is on the same order of that of the convective transport. Consequently, the global effect of differential diffusion is shown to be significant at these high Karlovitz number conditions.

Acknowledgments This work was sponsored by the Swedish Research Council (VR) and the National Centre for Combustion Science and Technology (CeCOST). We acknowledge PRACE for awarding us access to resource CURIE TN based in France at TGCC.

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 2 0 2 1 6 e2 0 2 3 2

Computations were also performed using the computer facilities provided by the Centre for Scientific and Technical Computing at Lund University (LUNARC), High Performance Computing Center North (HPC2N) and Center for Parallel Computers (PDC).

references

[1] Peters N, Williams FA. Premixed combustion in a vortex. Proc Combust Inst 1989;22:495e503. [2] Seshadri K, Peters N. The inner structure of methane-air flames. Combust Flame 1990;81:96e118. [3] Peters N. Turbulent combustion. Cambridge University press; 2000. p. 79. [4] Williams FA. Combustion theory. 2nd ed. California: The Benjamin/Cummings Publishing Company, Inc; 1985. [5] Swaminathan N, Bray KNC. Turbulent premixed flames. Cambridge University press; 2011. p. 15. € hr M, Aigner M, [6] Stopper U, Meier W, Sadanandan R, Sto Bulat G. Experimental study of industrial gas turbine flames including quantification of pressure influence on flow field, fuel/air premixing and flame shape. Combust Flame 2013;160:2103e18. [7] Strakey P, Sidwell T, Ontko J. Investigation of the effects of hydrogen addition on lean extinction in a swirl stabilized combustor. Proc Combust Inst 2007;31:3173e80. [8] Turns SR. An introduction to combustion. McGraw-Hill; 2006. [9] Mansour MS, Peters N, Chen YC. Investigation of scalar mixing in the thin reaction zones regime using a simultaneous CH-LIF/Rayleigh laser technique. Proc Combust Inst 1998;27:767e73. n M. [10] Li ZS, Kiefer J, Zetterberg J, Linvin M, Bai XS, Alde Development of improved PLIF CH detection using an Alexandrite laser for single-shot investigation of turbulent and lean flames. Proc Combust Inst 2007;31:727e35. n M. Investigation of [11] Kiefer J, Li ZS, Zetterberg J, Bai XS, Alde local flame structures and statistics in partially premixed turbulent jet flames using simultaneous single-shot CH and OH planar laser-induced fluorescence imaging. Combust Flame 2008;154:802e18. n M. Turbulence and [12] Li ZS, Li B, Sun ZW, Bai XS, Alde combustion interaction: high resolution local flame front structure visualization using simultaneous single-shot PLIF imaging of CH, OH, and CH2O in a piloted premixed jet flame. Combust Flame 2010;157:1087e96. € holm J, Rosell J, Li B, Richter M, Li ZS, Bai XS, et al. [13] Sjo Simultaneous visualization of OH, CH, CH2O and toluene PLIF in a methane jet flame with varying degrees of turbulence. Proc Combust Inst 2013;34:1475e82. [14] Dunn MJ, Masri AR, Bilger RW. A new piloted premixed jet burner to study strong finite-rate chemistry effects. Combust Flame 2007;151:46e60. [15] Dunn MJ, Masri AR, Bilger RW, Barlow RS, Wang G-H. The compositional structure of highly turbulent piloted premixed flames issuing into a hot coflow. Proc Combust Inst 2009;32:1779e86. [16] Dunn MJ, Masri AR, Bilger RW, Barlow RS. Finite rate chemistry effects in highly sheared turbulent premixed flames. Flow Turb Combust 2010;85:621e48. [17] Sankaran R, Hawkes ER, Chen JH, Lu T, Law CH. Structure of a spatially developing turbulent lean methaneeair Bunsen flame. Proc Combust Inst 2007;31:1291e8. [18] Chakraborty N, Cant RS. Influence of Lewis number on strain rate effects in turbulent premixed flame propagation. Int J Heat Mass Transf 2006;46:2158e72.

20231

[19] Han I, Huh KY. Effects of the Karlovitz number on the evolution of the flame surface density in turbulent premixed flames. Proc Combust Inst 2009;32:1419e25. [20] Hawkes ER, Chen JH. Evaluation of models for flame stretch due to curvature in the thin reaction zones regime. Proc Combust Inst 2005;30:647e55. [21] Yu R, Bai XS. Direct numerical simulation of lean hydrogen/ air auto-ignition in a constant volume enclosure. Combust Flame 2013;160:1706e16. [22] Poludnenko AY, Oran ES. The interaction of high-speed turbulence with flames: global properties and internal flame structure. Combust Flame 2010;157:995e1011. [23] Poludnenko AY, Oran ES. The interaction of high-speed turbulence with flames: turbulent flame speed. Combust Flame 2011;158:301e26. [24] Thevenin D. Three-dimensional direct simulations and structure of expanding turbulent methane flames. Proc Combust Inst 2005;30:629e37. [25] Thevenin D, Gicquel O, De Charentenay J, Hilbert R, Veynante D. Two- versus three-dimensional direct simulations of turbulent methane flame kernels using realistic chemistry. Proc Combust Inst 2002;29:2031e9. [26] Day M, Tachibana S, Bell J, Lijewski M, Beckner V, Cheng RK. Combust Flame 2012;159:275e90. [27] Yoo CS, Richardsson ES, Sankaran R, Chen JH. A DNS study on the stabilization mechanism of a turbulent lifted ethylene jet flame in highly-heated coflow. Proc Combust Inst 2011;33:1619e27. [28] Hawkes ER, Sankaran R, Cutherland JC, Chen JH. Scalar mixing in direct numerical simulations of temporally evolving plane jet flames with skeletal CO/H2 kinetics. Proc Combust Inst 2007;31:1633e40. [29] Chen JH. Petascale direct numerical simulation of turbulent combustiondfundamental insights towards predictive models. Proc Combust Inst 2011;33:99e123. [30] Aspden AJ, Day MS, Bell JB. Lewis number effects in distributed flames. Proc Combust Inst 2011;33:1473e80. [31] Aspden AJ, Day MS, Bell JB. Characterization of low Lewis number flames. Proc Combust Inst 2011;33:1463e71. [32] Aspden AJ, Day MS, Bell JB. Turbulence-flame interactions in lean premixed hydrogen: transition to the distributed burning regime. J Fluid Mech 2011;680:287e320. [33] Savre J, Carlsson H, Bai XS. Turbulent methane/air premixed flame structure at high Karlovitz numbers. Flow Turb. Combust 2013;90:325e41. [34] Echekki T, Chen JH. Direct numerical simulation of autoignition in non-homogeneous hydrogen-air mixtures. Combust Flame 2003;134:169e91. [35] Yu R, Yu JF, Bai XS. An improved high-order scheme for DNS of low Mach number turbulent reacting flows based on stiff chemistry solver. J Comp Phys 2012;231:5504e21. [36] Strang G. On the construction and comparison of difference schemes. SIAM J Num Anal 1968;5:506e17. [37] Jiang GS, Shu CW. Efficient implementation of weighted ENO schemes. J Comp Phys 1996;126:202e28. [38] Zhang F, Yu R, Bai XS. Detailed numerical simulation of syngas combustion under partially premixed combustion engine conditions. Int J Hydrogen Energy 2012;37:17285e93. [39] Yu JF, Yu R, Fan XQ, Christensen M, Konnov AA, Bai XS. Onset of cellular flame instability in adiabatic CH4/O2/CO2 and CH4/air laminar premixed flames stabilized on a flatflame burner. Combust Flame 2013;160:1276e86. [40] Brown PN, Byrne GD, Hindmarsh AC. VODE: a variablecoefficient ODE solver. SIAM J Sci Stat Comp 1989;10:1038e51. [41] Smooke MD, Giovangigli V. Formulation of the premixed and nonpremixed test problems. In: Reduced kinetic mechanisms and asymptotic approximation for methane-air flames. Berlin: Springer Verlag; 1991.

20232

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 2 0 2 1 6 e2 0 2 3 2

[42] Li J, Zhao Z, Kazakov A, Dryer FL. An updated comprehensive kinetic model of hydrogen combustion. Int J Chem Kinet 2004;36:566e75. [43] Marinov NM, Westbrook CK, Pitz WJ. In: Chan el SH, editor. Transport phenomena in combustion. Washington, DC: Taylor and Francis; 1996. p. 118e29. [44] Kee KJ, Rupley FM, Miller JA, Coltrin ME, Grcar JF, Meeks E, et al. CHEMKIN collection. Release 3.6. San Diego, CA: Reaction Design, Inc.; 2000. [45] Kraichnan RH. Diffusion by a random velocity field. Phys Fluids 1970;13:22e31. [46] Yu R, Bai XS. A fully divergence-free method for generation of inhomogeneous and anisotropic turbulence with large spatial variation. J Comp Phys 2014;256:234e53. [47] Hinze JO. Turbulence. 2nd ed. New York: McGraw-Hill; 1975.

[48] Pope CJ, Shandross RA, Howard JB. Variation of equivalence ratio and element ratios with distance from burner in premixed one-dimensional flames. Combust Flame 1999;116:605e14. [49] Linan A, Williams FA. Fundamental aspects of combustion. Oxford university press; 1993.  Curran HJ, Simmie JH, Pitz WJ, Westbrook CK. A [50] Conaire MO, comprehensive modeling study of hydrogen oxidation. Int J Chem Kin 2004;36:603e22. [51] Smith GP, Golden DM, Frenklach M, Moriarty NW, Eiteneer B, Goldenberg M, et al. . [52] San Diego Mechanism web page Chemical-kinetic mechanisms for combustion applications. Mechanical and Aerospace Engineering (Combustion Research), University of California at San Diego; 2012-09-07., http://combustion.ucsd.edu.