Impact of Turbulence Intensity and Equivalence Ratio on the Burning Rate of Premixed Methane – Air Flames

Direct Numerical Simulations (DNS) have been conducted to study the response of initially laminar spherical premixed methane–air flame kernels to successively higher turbulence intensities at five different equivalence ratios. The numerical experiments include a 16-species/25-step skeletal mechanism for methane oxidation and a multicomponent molecular transport model. Highly turbulent conditions (with integral Reynolds numbers up to 4 513) have been accessed. The effect of turbulence on the physical properties of the flame, in particular its consumption speed Sc, which is an interesting measure of the turbulent flame speed ST has been investigated. Local quenching events are increasingly observed for highly turbulent conditions, particularly for lean mixtures. The obtained results qualitatively confirm the expected trend regarding correlations between u/SL and the consumption speed: Sc first increases, roughly linearly, with u/SL (low turbulence zone), then levels off (bending zone) before decreasing again (quenching limit) for too intense turbulence. For a fixed value of u/SL, Sc/SL varies with the mixture equivalence ratio, showing that additional parameters should probably enter phenomenological expressions relating these two quantities.


Introduction
Premixed turbulent flames remain a significant research challenge in the combustion community.Their propagation is of significant and practical importance in the design and optimization of efficient low emission burners.One essential property is their turbulent burning velocity S T [1,2].Its influential role in turbulent premixed flame modeling has attracted a wealth of experimental and theoretical investigations in a variety of flow configurations [3].Experimental setups to measure S T differ from one study to another, and many challenges remain-unsteadiness and inhomogeneity of the flow, influence of the ignition source on the flame propagation, difficulties in quantifying heat losses and length scales-leading to a large scattering in the dataset of the latter for a fixed root-mean-square (rms) velocity to laminar flame velocity ratio, u ′ /S L ; an accurate measurement of S T is still difficult (see e.g., [3,4] and references therein).
Considering the correlation of the ratio of turbulent to laminar burning velocities S T /S L and u ′ /S L [5], old experimental data reveal a steady linear increase of S T /S L with u ′ /S L [6] for the higher range of turbulence Reynolds number Re t = u ′ l t /ν, where l t is the integral length scale and ν is the kinematic viscosity.This is contrary to the seminal investigations in, for instance, [2,7] for a lower range of Re t .While most combustion vessel experiments witness a global gas phase extinction at too intense turbulence, this is still considered as an open issue for other configurations and measurement techniques [8][9][10].
Beyond theoretical considerations and complex experimental measurements, numerical simulation, in particular, Direct Numerical Simulation (DNS) is now a widely accepted technique for investigating turbulent combustion problems.DNS provides a unique insight into the physics of flame-turbulence interactions (see e.g., [11] and references therein).The advent of advanced computer technology and the development of efficient computational algorithms have now made it possible to employ DNS with realistic models.Numerical investigation of multi-component multi-dimensional reactive flows are possible in long but affordable computing times, without compromising accuracy.All parameters can be determined in such numerical experiments at the same time, offering a high level of flexibility for determining specific quantities in details under prescribed conditions.This makes DNS a convenient tool for the generation of correlations between S T /S L and u ′ /S L for various flow and burning conditions.
Due to its practical relevance, methane is one of the first fuels for which reduced kinetic models have been developed, since the chemical complexity associated with its direct simulation using full reaction schemes is still often prohibitive.Consequently, a large proportion of numerical studies of hydrocarbon flames is confined either to simple and/or two-dimensional flows with full chemical kinetic schemes (see [12] and references therein) or to complex and/or three-dimensional flows with simple/reduced schemes (see [13] and references therein).However, it is worth noting that significant advances have been made in recent years towards simulating systems burning practically relevant fuels together with detailed and accurate models for both the chemical and transport phenomena [14][15][16].Thanks to the rapid growth of computational capabilities, adequate numerical predictions of the turbulent burning velocity have been provided [17,18] up to laboratory scale configurations [19,20] using high fidelity DNS.Such numerical studies have shown that turbulence enhances the flame speed with a rate faster than the growth in flame surface area, which is consistent with earlier studies conducted with simple models.A recent review is given by Driscoll [3].
In the present work, direct simulations are used to investigate systematically premixed methane-air mixtures under successively higher turbulence intensities using detailed physicochemical models, with the aim of supplying complementary data of turbulent consumption speed S c (which is an interesting measure of the turbulent burning velocity S T ) as a function of turbulence intensity u ′ /S L .It is then possible to determine in particular if and how the turbulent Reynolds number directly affects S c for a given turbulence intensity; and if the still controversial "bending zone" and "quenching limit" on the S c curves versus u ′ /S L can be captured numerically.The setup and range of parameters to be explored mimic combustion vessel experiments (see e.g., [2,7] and references therein).The computational settings are similar to those used in [13] for a single-step reaction mechanism and unity Lewis number assumption.In the present study, both the three-dimensional and turbulent nature of the flame-kernel evolution together with a skeletal reaction scheme and multi-component transport models are considered.
In the next section (Section 2), details of the Direct Numerical Simulations and general computational procedure are given.Two major fine-grain parallelism issues-the initial fluctuating velocity generation and data input/output scenarios-for the employed DNS solver are discussed briefly.The problem configuration and initialization is next outlined in Section 3, directly followed by the numerical results for the various considered cases in Section 4, just before concluding in Section 5.

Direct Numerical Simulation
A DNS must provide as far as possible an exact solution for both fluid dynamics and flame structure.Even though this method requires prohibitive numerical costs for practical engineering configurations at large scales, it offers an excellent complement to experiments in order to assess the importance of various physical mechanisms, to obtain complementary information on flame structure, and therefore to improve turbulent combustion modeling [21].
In this study, the massively parallel DNS flame solver parcomb [22] is used.It solves the full compressible reactive Navier-Stokes system coupled with detailed physicochemical models [23,24].The solved balance equations are, using classical notations [25] and ignoring all external forces: where ρ denotes mixture density, u j the components of the hydrodynamic velocity, p the pressure, N s the total number of species, V kj the component of the diffusion velocity of species k in the direction j, ωk the chemical production rate of species k, q j the j th -component of the heat flux vector and τ ij the stress tensor: where δ is the Kronecker symbol and η is the dynamic viscosity.
A spatial sixth-order central scheme and an explicit fourth-order Runge-Kutta time integrator are employed.In its recent version [26][27][28], the skew-symmetric formulation [29] has been implemented for the convective terms in order to reduce even further numerical dissipation and increase stability.The extended Navier-Stokes Characteristic Boundary Conditions (NSCBC) [30] are used, with non-reflecting boundaries and pressure relaxation applied along all open faces.The parallelization is in a 3D block structure using the MPI protocol for data exchange.The code offers, for instance on the IBM BlueGene/P 13.5% of the peak performance on each node (139 TeraFlop/s for the full machine) and a near perfect parallel scaling for a real, three-dimensional run using up to 4 096 PowerPC 450 computing cores.Further details concerning code structure, optimization, application and recent upgrades can be found for instance in [27,31].
In order to access higher values of Re t on fine-grain parallel systems, two turbulence generators based on digital filters [32] and random noise diffusion [33] have been hybridized, implemented and parallelized on massively parallel computers.With this simple and flexible approach, the restriction on problem size imposed by the previously implemented generator based on inverse FFT has been removed, paving the way for simulations of larger domains at considerably higher Re t values.For the results presented later, the initial turbulent field has been systematically generated using this hybrid technique (see [27] for further details).
The second major issue that needed attention towards fine-grain parallelism is that of efficient, fully parallel data input/output (I/O).The traditional sequential I/O approaches have been replaced by a fully parallel I/O (via MPI-I/O), where multiple processes of the parallel program access data (for read/write) from a common, shared file.This provides both higher performance (speedup in time needed for writing/reading all files by a factor of at least 3) and single (restart/solution) data files.

Flame Configuration and Initialization
Initially perfectly spherical laminar premixed methane-air flames are considered in all computations, within a cubic computational domain of sides L = 4.0 cm (Figure 1a) and a uniform grid spacing of 35-20 µm for the mild to the most intense turbulent cases, respectively, ensuring resolution of the smallest (Kolmogorov) scales, since their length decreases with increasing Re t .The initialization and subsequent development of a premixed flame-kernel under the influence of a turbulent flow field is an interesting configuration which allows turbulent flames to be studied well away from the influence of external perturbations such as walls and artificial boundary conditions.From a fundamental point of view, it offers the possibility to study complex multi-scale flows involving coupled physicochemical processes.Simultaneously, it has direct practical relevance in a number of industrial cases including spark-ignition internal combustion engine and gas turbine re-ignition, as well as safety issues.
Methane oxidation is modeled by a 25-step skeletal scheme comprising 16 species (CH [34].This reaction mechanism is retained here due to its simplicity and stability, and provides sufficiently accurate results for lean up to stoichiometric conditions.It has been successfully used for large scale multi-dimensional direct computations of non-premixed methane jet flames [35] and most recently, highly turbulent premixed flames [27,36].However, it would be inadequate for methane-rich flames due to the absence of C 2 and higher carbon-chain reactions, the reason why Φ ≤ 1.0 for the present study. The initial mixture composition (Y i ), prescribed burnt (T b ) and unburned (T u ) temperatures, laminar flame speed S L , thermal flame thickness for the various mixture equivalence ratios (Φ) are given in Table 1.The given range for Ka corresponds to the different turbulence intensities, as defined below.The initial system is a hot (T = T b ) perfectly spherical laminar flame-kernel of initial radius r o = 5.0 mm, located at the center of the computational box and surrounded by a fresh premixed atmospheric mixture of methane and air at T u .The initial mass fraction values of Y CH 4 and Y O 2 at T u , and Y CO 2 and Y H 2 O at T b are prescribed outside and within the kernel, respectively.These initial values for any primitive variable ϕ are transformed into smooth profiles according to where ∆ϕ is the difference between the initial values (ϕ o ) in the fresh and burnt gas mixture.The constant s is a measure of the stiffness at the fresh/burnt gas interface and is in the order of a few hundred.In this range, the influence of s is confined to the very early part of the simulation and therefore does not impact the analysis presented below at later times.In all cases, an appropriate nitrogen complement is added everywhere at start.
Table 1.Initial flame and flow parameters.To investigate systematically the influence of Re t on the fuel consumption rate, the calculations for a given Φ were repeated with exactly the same initial composition, but with an initial pseudo-turbulent velocity field at successively higher intensity.The eddy turn-over time τ = l t /u ′ , u ′ and Re t for the various cases are given in Table 2, while l t = 3.2 mm is kept constant and ν ≈ 1.56 × 10 −5 m 2 /s.In Table 2, data for case 1 (lowest turbulence intensity) is omitted due to lack of space.Note also that not all the cases shown here are realized for every Φ given in Table 1, due to the higher sensitivity of the leaner mixtures to increasing turbulence intensity, as will be demonstrated later.For the mixtures with Φ = 0.6 & 0.7, only cases 1-7 are considered, while for Φ = 0.8 and 0.9, only cases 1-8 and 1-9 are performed, respectively.Based on the above turbulence characteristics, it is clear that a wide range of turbulence Reynolds numbers have been explored with more than 90% of the investigated combustion phenomena covering the Thin-Reaction Zone (TRZ) regime, according to the modified regime diagram of Peters [4], as illustrated in Figure 1b.Here, it is believed that small turbulent eddies are capable of penetrating and disrupting the preheat zone but fall short of doing so on the reaction zone because of an order of magnitude disparity in the thickness of these flame layers.For the presented computations, our local Linux-based PC-cluster has been employed in a single-user mode.It is built with Opteron quad-core nodes with 32 GB memory/node and an Infiniband connection, with a peak performance of roughly 5 Tflop/s, reserved exclusively for this project.Additionally, three HPC systems across Europe (BABEL-IDRIS in France, HPCx-and HECToR-EPCC in Scotland) were employed, with up to 4 096 computing cores.Typically, 10 days of computing time are needed for the lowest grid resolution (35 µm spacing) to reach t/τ = 1 for ϕ = 0.6.The longest computation requires about two months.While the fully resolved system was simulated, only a reduced grid was saved for post-processing.For all results presented below, the simulations have been carried out up to a non-dimensional time t ≥ 1.3τ as recommended for DNS relying on time-decaying turbulence [37].

Turbulent Flame Structure
The effect of turbulence on the physical structure of the flame is first investigated visually by examining both the instantaneous and temporal evolution of the temperature and of selected species mass fraction fields, plotted in Figures 2-5.All figures are cut in the middle plane z = 2.0 cm.2a) is progressively being distorted and stretched (Figure 2b) by the very strong turbulent field with time, leading to the creation of islands (in the form of both hot and fresh gas pockets) and edge flame-like structures (Figure 2c,d) [38] at various locations within the computational domain and for various shapes and sizes.For a high turbulence intensity, local flame extinction and flame-flame interactions become important.The temporal evolution of the iso-contours of the minor radical CH 3 O is shown in Figure 3 for the same Φ and Re t values.enThe wrinkling effect induced by turbulence follows the patterns in the temperature field but shows much stronger disruptions at later times (t = 0.9-1.3τ ).Such events are usually not observed when looking at integrated quantities, like temperature, illustrating the importance of a detailed description of chemical processes.Considering CH 3 O as a measure of flame activity, local extinction events can be tracked with a good temporal and spatial resolution.Numerous fresh gas islands and flame pinch-off events are visible (Figure 3c), evidencing local flame extinctions.Looking at the mass fraction values, the peak species mass fraction rises soon after its formation from 1.24 × 10 −5 at t = 0.1 τ to 6.47 × 10 −5 at t = 0.5 τ , after which it drops to 5.76 × 10 −5 at t = 0.9 τ , then down to 3.12 × 10 −5 at t = 1.3 τ .
The iso-contours of most major and minor species exhibit qualitatively similar patterns to those of the temperature and CH 3 O fields, respectively.They are therefore not shown in the interest of space.Snapshots of the iso-surface of the mass fraction of the oxygen atom is shown in Figure 4 for different Re t values at the same time t = 1.3 τ and mixture ratio Φ = 0.9.In the past, DNS easily accessed mild turbulence conditions such as the one shown in Figure 4a-c where the flame is only slightly contorted.Comparing the snapshots in Figure 4, it is observed qualitatively that the amount of wrinkling increases strongly from Figure 4a-c and tends to saturate afterwards (compare in particular Figure 4d with Figure 4e).For higher values of Re t , considerable structural modifications are observed, in particular flame-flame leading to pinch off as evident in the higher Re t snapshots in Figures 4c-e.The resulting turbulent flame structure is then marred with numerous perforations.Further increase in turbulence intensity leads to a further increase of pinch off and mutual annihilation effects, thereby limiting further increase in the flame surface area.Consequently, it drops steadily as evidenced in Figures 4f-i, indicating the growing importance of extinction processes, as also observed experimentally in combustion vessel experiments [2,39].To show the effect of the equivalence ratio on the turbulent flame structure, Figure 5 presents the instantaneous iso-contours of the mass fraction of OH (case 6) for different mixture equivalence ratios Φ = 0.7, 0.8, 0.9 and 1.0 at the same time t = 1.3 τ .The color scale is kept identical for all plots.The lower flame activity when decreasing the equivalence ratio can easily be seen in this figure, explaining again why local flame extinction is systematically observed earlier at lower values of Φ.Since the employed chemical scheme has only been validated for lean to stoichiometric conditions, similar studies for rich conditions cannot be presented.

Turbulent Burning Velocity
Now, the time evolution of the turbulent burning velocity as a function of the turbulence intensity is investigated in moderate to intense turbulent flows for various Φ.Following [3], it is emphasized that different definitions of S T may be useful in different contexts.In the present simulations, the consumption speed S c is defined following [40] and computed as the volume integrated overall rate of fuel combustion per unit flame area (assuming complete fuel consumption on the burned gas side), in a similar manner as employed, for instance, in [41,42]: where W f is the fuel molecular weight and ρ f and Y f are the density and mass fraction of the fuel in the fresh gas mixture, respectively.The associated scaling area has been obtained in a post-processing step using AnaFlame [43].It is chosen as the area of a suitable iso-contour of the reaction progress variable c, (defined following [44] in terms of reduced temperature) representing the wrinkled flame front.For instance, a choice of c = 0.6, corresponding to the position of maximum heat release in the laminar flame could be adopted.Thus the value of S c is not unique.We note also that ensemble averaging would be needed to provide a mean isosurface scaling area, which could account further for flame wrinkling effects.The computed fuel consumption rate is then directly related to the overall burning rate (by scaling with this area), and constitutes an interesting measure of the burning velocity.In all results shown, S c is scaled by its initial value S c,0 (t = 0) (before starting the interaction with the prescribed field of turbulence) at the corresponding equivalence ratio.This global flame parameter is of primary importance to burner designers as far as the characterization of premixed turbulent combustion is concerned.Its prediction remains a challenge for turbulent premixed combustion modelers [45].Most existing models are still phenomenological in nature and often fail beyond the flammability limits, whenever the pre-mixture is either too lean or too rich to sustain a flame in the mean.
The temporal evolution of the consumption speed is shown in Figure 6 for the various cases.Initially, it remains nearly constant up to t ≈ 0.3 ms and t ≈ 0.6 ms for Φ = 1.0 and Φ = 0.6, respectively, after which it increases steadily.Under the same flow conditions, the time required by the leaner mixtures to sustain and propagate a flame front is delayed relative to the richer ones, as expected.After a successful take-off, the profiles for the various cases take different courses depending on the turbulence intensity u ′ , mostly with a steady increase in slope.For a given value of u ′ /S L (which is the parameter mostly used for modeling purposes), the gradient of the profiles increases with increasing Φ as shown in Figure 6a, where all the profiles for the five mixture compositions are plotted together for a nearly constant value of u ′ /S L .Also, it is obvious that the leaner the mixture, the lower the fuel consumption rate and vice versa.Note the progressive exponential nature of the curves associated with the richer mixtures (Φ ≥ 0.9) under highly turbulent conditions.However, as the turbulence intensity further increases, a maximum in S c is observed.This critical (or saturation) rms velocity will be denoted u ′ s .As found in the DNS, the richer the mixture, the higher the scaled saturation velocity: u ′ s /S L is about 29.0 for Φ = 0.6 and up to 39.0 at stoichiometry.Increasing u ′ beyond u ′ s , the obtained temporal profiles of S c first cluster together (e.g., Figures 6b & c for Φ = 1.0 & 0.9) before declining steadily (e.g., Figures 6d & e for Φ = 0.8 & 0.7).This scenario is a well known phenomena observed in several combustion vessel experiments and termed "bending effect" [44,46].Whether such a decline persist until the flame reaches global quenching (the so called "quenching limit", also observed in many experiments) for too intense turbulence, will be the subject of further studies.Figure 6f includes a case for which, due to a high initial turbulence level (u ′ s /S L = 88), the pre-mixture is unable to propagate and slowly fades out.The present DNS results thus confirm the theoretical and experimental assertions that turbulence cannot increase burning rates indefinitely [7,9].Furthermore, leaner methane-air mixtures (e.g., Φ = 0.6 at u ′ /S L ≥ 74.0) exhibit a greater decline in S c with increasing u ′ /S L than do richer ones.
Since high volumetric heat release rates in compact burners at low stoichiometries are often required, combustion in a highly turbulent medium is often necessary [47].Intense turbulence stirring usually results in flames that burn and/or propagate faster and release higher amounts of energy in the form of heat.However, the obtained results show that optimum conditions exist, depending on mixture stoichiometry, above which the resulting flames will be weakened by a further increase in turbulent stirring.

Concluding Remarks
Direct Numerical Simulations have been realized in a parametric study in order to investigate the response of turbulent consumption speed described by the fuel consumption rate for premixed methane-air combustion.Five different mixture equivalence ratios Φ have been considered, each for increasing values of the integral Reynolds number.Highly turbulent conditions (Re t up to 4 513) have been accessed for each Φ and the effect of turbulence on the physical structure of the flame has been investigated.Local flame quenching has been observed within the computational domain.Quenching effects appear earlier with the leaner mixtures, as expected.Profiles of the fuel consumption rate rise steadily with increasing turbulence intensity till a saturation rms velocity is reached, beyond which the time evolutions collapse together before decreasing steadily.These DNS results confirm qualitatively theoretical and experimental findings from the literature regarding the correlation between scaled turbulence intensity u ′ /S L and turbulent flame speed quantified by S c [44]: S c first increases, roughly linearly, with u ′ /S L (low turbulence zone) and then levels off (bending zone).For a fixed value of u ′ /S L , which already incorporates compositional variations, values of S c /S L still appear to depend on the mixture equivalence ratio.As a consequence, simplified phenomenological expressions relating S c /S L simply to u ′ /S L should probably also take into account the influence of further parameters [3].
Assuming a generic relation in the form S c /S L = [1 + C(u ′ /S L ) n ] 1/n [48], the model constants C and n are found to be ≈ 0.25 and ≈ 1.16, respectively, from fits of the DNS data of S c /S L for single runs, in comparison with experimental estimations of 0.45 and 1.0 [47].The DNS data are presently extended in order to propose quantitative correlations that might be used in turbulent combustion models with a sufficient level of generality.Ensemble averaging should also be conducted in order to obtain ultimately a statistically significant analysis of the results [49].This ensure that a mean isosurface area of c is employed in conjunction with Equation 7 in order to account for flame wrinkling effects.

Figure 1 .
Figure 1.(a) Exemplary view of the configuration showing the heavily wrinkled iso-surface of the mass fraction of CO 2 colored by the HO 2 flame radical and (b) modified combustion diagram of Peters [4] with positions of current DNS.

Figure 4 .
Figure 4. Instantaneous iso-surface of the mass fraction of O for different integral Reynolds number Re t at the same mixture equivalence ratio Φ = 0.9 and same time t = 1.3 τ .

Figure 5 .
Figure 5. Instantaneous iso-contours of the mass fraction of OH (case 6, Re t = 2 051) for different mixture equivalence ratios Φ at the same time t = 1.3 τ .

Figure 6 .
Figure 6.Temporal evolution of the volume-integrated fuel (CH 4 ) consumption rate S c (normalized by their laminar flame value at the corresponding equivalence ratio and at t = 0) for different turbulence intensities u ′ /S L and mixture equivalence ratios Φ .