Direct Numerical Simulations of the Impact of High Turbulence Intensities and Volume Viscosity on Premixed Methane Flames

Abstract

Parametric direct numerical simulations (DNS) of turbulent premixed flames burning methane in the thin reaction zone regime have been performed relying on complex physicochemical models and taking into account volume viscosity (𝜅). The combined effect of increasing turbulence intensities (𝑢^′) and 𝜅 on the resulting flame structure is investigated. The turbulent flame structure is marred with numerous perforations and edge flame structures appearing within the burnt gas mixture at various locations, shapes and sizes. Stepping up 𝑢^′ from 3 to 12 m/s leads to an increase in the scaled integrated heat release rate from 2 to 16. This illustrates the interest of combustion in a highly turbulent medium in order to obtain high volumetric heat release rates in compact burners. Flame thickening is observed to be predominant at high turbulent Reynolds number. Via ensemble averaging, it is shown that both laminar and turbulent flame structures are not modified by 𝜅. These findings are in opposition to previous observations for flames burning hydrogen, where significant modifications induced by 𝜅 were found for both the local and global properties of turbulent flames. Therefore, to save computational resources, we suggest that the volume viscosity transport term be ignored for turbulent combustion DNS at low Mach numbers when burning hydrocarbon fuels.

1. Introduction

Both the availability of electrical energy and all transportation systems are controlled to a large extent by turbulent combustion processes, such as in gas turbines or Internal Combustion engines, burning either fossil or renewable fuels. Optimizing further such well-known systems is only possible with a much better understanding of all relevant processes involved. Detailed quantitative experiments are needed and very useful, but sometimes impossible and usually limited to only a few flame quantities. As a complement, detailed (direct) numerical simulations (DNSs) are increasingly gaining grounds as a reliable tool for detailed investigations towards fundamental understanding of a variety of turbulent combustion phenomena [1]. Progress in numerical techniques as well as computational power now allows quantitative investigations of turbulent reacting flows for increasingly realistic conditions, for which practically relevant fuels are simulated at higher turbulent Reynolds numbers (Re_𝑡), while employing at the same time complex physicochemical models to describe turbulence, molecular transport, and chemistry [2–4].

There are in fact only two possibilities to obtain statistically significant results. Ensemble averaging is the most straightforward way, based on a repetition of the same “numerical experiment,” then averaging the observations [5]. This can sometimes be supported or even replaced by a spatial averaging along statistically homogeneous directions, when available. As an alternative, it is possible to continue the simulation for a long time and employ time averaging between samples separated by a sufficiently long interval. In this last case, turbulence must remain at a constant level in time. Therefore, time-averaging is only possible for spatially evolving turbulence, usually injected through a boundary condition, or by artificially forcing the flow in time-decaying simulations, which is again associated with many theoretical and practical issues. Relying on ergodicity, both averaging methods should deliver the same results. Considering numerical requirements, both methods lead to a considerable increase of the computational cost. The simulation must be either repeated (ensemble averaging) or pursued during a very long physical time (time averaging), leading typically to an order of magnitude increase of the effort compared to a single DNS.

More generally, the need for accurate physical models in practical combustion computations has been demonstrated clearly for various flame configurations (see, e.g., [2, 6–8] and references therein). Unfortunately, the inclusion of the volume (or bulk) viscosity (𝜅) transport term is still controversial in the literature, since it has been traditionally neglected in almost all turbulent combustion studies up to now without any convincing justification [9], even though a few publications have demonstrated that it is wrong to neglect it a priori [10, 11]. Its impact should be indeed small for many applications. On the other hand, the value of 𝜅 for hydrogen, for example, is very high and could lead to noticeable flame modifications. Moreover, studies of the impact of 𝜅 on the flame structure while considering only laminar or mild turbulent conditions might lead to biased conclusions since more realistic, higher values of Re_𝑡 have seldom been assessed.

Neglecting the volume viscosity term in computational models seems logical at very low Mach number and for boundary layer flows. In practice, the additional computational cost associated with its evaluation in multi-component mixtures also contributes to a larger extent for such exclusions. Even more importantly, an overconfidence in the partly misleading suggestion of Stokes [9] may have contributed to this systematic evasion. Indeed, the absence of volume viscosity effects in dilute monatomic gases such as Ne [12] is an exception and not a rule [9, 13], since 𝜅 arises very noticeably in dense gases and liquids. While monoatomic gases show no evidence of bulk viscosity effects, experimentally measured values of the bulk-to-shear viscosity ratio (𝜅/𝜂) vary widely for most of the diatomic gases used nowadays in both academic and practical experimental and numerical combustion studies [12, 14, 15]. In [15], it is stated that the volume viscosity is at least of the same order of magnitude as the shear viscosity (𝜂) at 300 K, with a ratio 𝜅/𝜂 for hydrogen as high as 52 at a temperature of 1000 K. Our recent findings for turbulent premixed flames burning hydrogen-containing fuels show a nonnegligible influence of the volume viscosity term on the turbulent flame structure, even for average quantities [16]. On the other hand, first DNS computations for flames burning methane show a minor impact of volume viscosity for all considered cases for global flame properties such as the volume-integrated heat release rate and turbulent burning velocity [17]. The resulting picture is hence quite complex. A systematic quantitative analysis of the influence of 𝜅 on turbulent methane flames is therefore necessary to refine previous findings in a definitive manner. This is of high relevance, in particular to the combustion DNS and modeling community.

Since single isolated DNS simulations might sometimes be misleading, systematic DNS computations has been realized in order to carry out ensemble-averaging and obtain a high statistical significance [5] for the observations.

Three main challenges have therefore been finally addressed in the present work: the chemical and transport complexity has been stepped up by considering methane flames and volume viscosity effects, respectively, while increasing simultaneously the turbulence intensity. To the authors knowledge, such a combination of numerical, physical, and flow conditions has never been taken into account simultaneously in the same analysis.

In the next section (Section 2), an outline of the computational procedure is given, where the numerical method is described, followed with the problem configuration and initialization in Section 3. A few comments regarding recently added modules to the employed DNS tool towards fine-grain parallelism are discussed as well. The numerical results are finally presented in Section 4 before the concluding remarks (Section 5).

2. Direct Numerical Simulation

The DNS approach consists in solving as exactly as possible all the physical space and time scales embedded in the representative flow equations, without adding any model for turbulence. A DNS must thus provide an exact solution for both fluid dynamics and flame structure. Even though this method requires prohibitive numerical costs for practical configurations, 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 [18].

The DNS code employed in this work is the massively parallel flame solver, parcomb [19, 20], which solves the full compressible reactive Navier-Stokes system coupled with detailed chemistry and multicomponent transport models (via coupling with the chemkin [21], transport [22] and eglib [23] libraries): 𝜕𝜌+𝜕𝜕𝑡𝜌𝑢_𝑗𝜕𝑥_𝑗𝜕=0,𝜌𝑢_𝑖+𝜕𝜕𝑡𝜌𝑢_𝑗𝑢_𝑖𝜕𝑥_𝑗=−𝜕𝑝𝜕𝑥_𝑖+𝜕𝜏_𝑖𝑗𝜕𝑥_𝑗,𝜕𝜌𝑌_𝑘+𝜕𝜕𝑡𝜌𝑢_𝑗𝑌_𝑘𝜕𝑥_𝑗𝜕=−𝜌𝑌_𝑘𝑉_𝑘𝑗𝜕𝑥_𝑗+𝜔̇_𝑘;𝑘=1,2,…,𝑁_𝑠,𝜕𝜌𝑒_𝑡+𝜕𝜕𝑡𝜌𝑢_𝑗𝑒_𝑡𝜕𝑥_𝑗=−𝜕𝑞_𝑗𝜕𝑥_𝑖+𝜕𝜏_𝑖𝑗𝑢_𝑖𝜕𝑥_𝑗,𝑝𝜌=𝑅𝑊𝑇, (1) where 𝜌 denotes mixture density, 𝑢_𝑗 the components of the hydrodynamic velocity, 𝑝 the pressure, 𝜏_𝑖𝑗 the stress tensor, 𝑁_𝑠 the total number of species, 𝑉_𝑘𝑗 the component of the diffusion velocity of species 𝑘, in the direction 𝑗, 𝜔̇_𝑘 the chemical production rate of species 𝑘 and 𝑞_𝑗 the 𝑗th-component of the heat flux vector. The volume viscosity coefficient 𝜅 is a function of the fluid local properties and appears explicitly within 𝜏_𝑖𝑗 in the momentum and energy equations: 𝜏_𝑖𝑗=−𝑝𝛿_𝑖𝑗+𝜂𝜕𝑢_𝑖𝜕𝑥_𝑗+𝜕𝑢_𝑗𝜕𝑥_𝑖+2𝜅−3𝜂𝜕𝑢_𝑘𝜕𝑥_𝑘𝛿_𝑖𝑗, (2) where 𝛿 is the Kronecker symbol and 𝜂 is the dynamic (or shear) viscosity. The DNS code has been coupled with the eglib transport library [23] for practical evaluation of the volume viscosity transport term. The single input data subroutine EGSKm(), suitable for fine-grain-distributed architectures with large number of processors, has been used. The integer 𝑚 attached to the subroutine name refers to the various methods/models used in evaluating 𝜅 [15, 16]. It varies between 2 and 6. The higher the value of 𝑚, the more expensive the underlying algorithm and the more accurate the expression for 𝜅, hence our choice of 𝑚=6 (EGSK6()) for the computations in the present study. The simplest and cheapest of the models (𝑚=2, EGSK2()) considers only the transport system matrix associated with the internal energy, whereas 𝑚=6 (EGSK6()) corresponds not only to the most accurate but also to the most expensive model since it solves the full matrix system associated with both translational and internal energy and relies on temperature-dependent ratios of collision integrals, which are evaluated as discussed in [15, 24, 25]. At the end, a direct inversion is performed to evaluate the volume viscosity for a given state of the mixture, characterized by the temperature and species mass fractions. We refer the interested reader to the eglib user manual [23] for further details and to [15] for an extensive discussion on the subject.

The above system of governing equations is solved in parcomb on a conventional Cartesian grid with high-order accurate and nondissipative numerical schemes. A spatial sixth-order central scheme, progressively reduced to fourth-order near the boundaries, is used. In the course of upgrading the code from the original 2D version to the current 3D state [4, 20], an improved scheme known as the skew-symmetric formulation [26] has been implemented for the convective terms in order to reduce even further numerical dissipation and increase stability. Time integration is performed in an explicit manner with a fourth-order Runge-Kutta scheme. The extended Navier-Stokes Characteristic Boundary Conditions (NSCBCs, [27, 28]) are used, with nonreflecting boundaries and pressure relaxation applied along the outflow faces.

parcomb is parallelized using a three-dimensional domain decomposition and MPI message passing protocol, offering a good peak performance and a near perfect parallel scaling for a full production-scale three-dimensional run for up to 4,096 computing cores on the IBM BlueGene/P. Further details concerning code structure, optimization, application, and recent upgraded/added modules can be found in [4, 29].

In order to access higher values of Re_𝑡 on fine-grain parallel systems, two turbulence generation techniques based on digital filters [30] and random noise diffusion [31] have been hybridized, implemented, and parallelized on massively parallel computers [16]. With this simple, flexible, and accurate approach, the restriction on problem size imposed by the previously implemented generator based on inverse FFT [32] has been alleviated, paving the way for simulations of large domains at considerably higher turbulent Reynolds numbers. For the results presented later, the initial turbulent field has been systematically generated using this hybrid technique.

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.

3. Flame Configuration and Initialization

A stoichiometric spherical premixed methane-air flame is considered in all computations, within a cubic computational domain of sides 𝐿=3.0 up to 4.0 cm and a uniform grid spacing of 26 μm, allowing a sufficient resolution of all spatial scales for all configurations considered. The ignition and subsequent expansion/development of such a premixed flame-kernel under the influence of a turbulent flow field is a very interesting configuration which allows turbulent flames to be studied well away from the influence of external constraints such as walls and artificial boundary conditions. It has furthermore a direct relevance for a number of industrial cases including spark-ignition internal combustion engine and gas turbine reignition as well as safety issues.

Methane oxidation is modeled by a 25-step skeletal scheme comprising 16 species (CH₄, O₂, H₂, H₂O, CH₂O, CO, CO₂, HO₂, OH, H, O, CH₃, HCO, H₂O₂, CH₃O, N₂) [33]. 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 multidimensional direct computations of nonpremixed methane jet flames [34] and, most recently, highly turbulent premixed flames [4, 17]. However, it would be inadequate for methane-rich flames due to the absence of C₂ and higher carbon-chain reactions, the reason why Φ≤1.0 for the present study.

The initial system is a stationary hot (𝑇_𝑏≈2200 K) perfectly spherical laminar flame-kernel of initial radius 𝑟_𝑜=5.0 mm, located at the center of the computational box and surrounded by a fresh atmospheric premixed mixture of methane and air at an unburnt temperature 𝑇_𝑢=300 K. The initial profile of any primitive variable Φ is prescribed according to Φ=Φ_𝑜+ΔΦ21−tanh𝑠⋅𝑟−𝑟_𝑜𝑟_𝑜, (3) where ΔΦ is the difference between the initial values (Φ_𝑜) in the fresh and burnt gas mixture. The constant 𝑠 is a measure of the stiffness of the fresh/burnt gas interface. The initial mass fraction values are 𝑌_CH4=0.055 and 𝑌_O2=0.220 at 𝑇_𝑢 outside the kernel, and 𝑌_CO2=0.120 and 𝑌_H2O=0.137 at 𝑇_𝑏 within. An appropriate nitrogen complement is added everywhere.

To investigate the influence of the turbulent Reynolds number Re_𝑡 based on the integral scale on the flame structure, the same calculations were repeated with exactly the same initial pseudoturbulent structures and initial composition, but with successively higher turbulent velocity fluctuations. The rms velocity fluctuation 𝑢′, the turbulent Reynolds number Re_𝑡, the eddy turn-over time 𝜏 used to quantify flame/turbulence interaction for time-decaying turbulence [35], and the Karlovitz number Ka≈[(𝑢^′/𝑠_𝐿)³(𝑙_𝑡/𝛿_𝑡ℎ)⁻¹]^1/2 are given in Table 1. The characteristic viscosity of the mixture, 𝜈=1.56×10⁻⁵ m²/s, the laminar flame speed 𝑠_𝐿=0.504 m/s, the laminar flame thermal thickness 𝛿_𝑡ℎ=(𝑇_𝑏−𝑇_𝑢)/max|∇𝑇|=0.36 mm, and the integral scale 𝑙_𝑡/𝛿_𝑡ℎ=8.89 are all constant for all simulations presented here.

Table 1 

Initial turbulence parameters.

Based on these turbulence characteristics, the flames considered here all fall within the thin reaction zone (TRZ) regime according to the modified regime diagram of Peters [36, 37], as illustrated in Figure 1. It is the regime in which 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 characteristic layers.

Figure 1 

Modified combustion diagram of Peters [18, 19] with positions of current DNS.

The general view of the configuration is illustrated in Figure 2, where the iso-surface of temperature with H₂O₂ slices and CO₂ slices with stream lines (thin white lines) is shown, revealing the heavily wrinkled flame front after interacting with the three-dimensional time-decaying homogeneous isotropic synthetic turbulent velocity field for 0.8 millisecond.

(a)

(b)

Figure 2 

General view of the configuration for the direct simulation, showing (a) the iso-surface of temperature with H2O2 slices and (b) CO2 slices with stream lines (thin white lines), revealing the heavily wrinkled flame front after interacting with the three-dimensional time-decaying homogeneous isotropic synthetic turbulent velocity field for 0.8 millisecond.

A total number of computing cores ranging between 512 and 4 096 were employed to solve the problems on the IBM POWER5 cluster (HPCx) system at EPCC (which delivers a total peak performance of 6 Tflop/s sustain). Typically, 10 days of computing time are needed to reach 𝑡/𝜏=1. The longest computation requires about two months. For the results presented below, a nondimensional time 𝑡=1.2𝜏, needed to obtain appropriate equilibrium conditions between turbulence and chemistry [35], is systematically retained for the analysis.

4. Results and Discussion

The obtained DNS datasets are processed using an in-house postprocessing library called AnaFlame [38, 39]. For cases 3 and 4, each experiment is performed twice with exactly the same initial and boundary conditions including turbulent features, except for the fact that the volume viscosity terms are deactivated in one but activated in the other. Furthermore, ensemble averaging is used in order to improve the statistical significance of the results by repeating each pair (𝜅≠0 and 𝜅=0) of case 3 with different turbulence fields initialized with different random seeds but having the same initial intensity (𝑢′) and then averaging the observations. Since a random number generator is used to determine the velocity phases [40], each DNS pair is associated with the same global properties of turbulence (spectrum, correlations, fluctuations, Reynolds number, etc.) but corresponds to different initial spatial features, and thus to a different evolution (i.e., realization) of the flame. These further realizations will provide statistical significance to the results obtained here [5]. In what follows, all profiles shown for case 3 are therefore the average of the six realizations.

Instantaneous solutions are analyzed in terms of conditional statistics of quantities relevant for modeling such as temperature, heat release, and selected mass fractions illustrating the turbulence-impaired flame structure. The impact of volume viscosity is assessed by comparing the above profiles for each of the twin computations of cases 3 and 4.

4.1. Laminar Flame Structure

First, a laminar case (referred hereafter as case 0) was computed in a smaller computational box (𝐿=3.0 cm) with the same initial mixture composition as in the turbulent cases. The instantaneous iso-contours of the mass fraction of H₂O₂ at different times (𝑡=0.1, 0.25, and 0.5 millisecond) with and without volume viscosity are shown in Figure 3. The concentric circles correspond, respectively, to the three time instances (with increasing radii). As a complement, the flame structure is shown in Figure 4 against the reaction progress, defined in terms of reduced temperature: 𝐶=(𝑇−𝑇_𝑢)/(𝑇_𝑏−𝑇_𝑢) [1]. Corresponding temporal evolution of the fuel consumption (𝑆_𝐶) and integrated heat release (𝐻_𝑟) rates (defined later) are shown in Figure 5. The instantaneous profiles of the temperature, heat release rate, and the mass fractions of OH, O, H₂O, H₂O₂, and HO₂ are shown at 𝑡=1.0 millisecond with (circled points) and without (solid line) taking into account the volume viscosity transport term. The same plotting style applies for the temporal profiles in Figure 5. It is absolutely impossible to differentiate the two numerical results in Figures 3–5 in both physical and progress variable space, even for flame radicals like H₂O₂ and HO₂ as well as global flame quantities like 𝑆_𝐶 and 𝐻_𝑟. All fields are qualitatively and quantitatively identical with relative differences well below 1% at all points. All other analyzed quantities (figures not shown) are identical as well. Hence, volume viscosity has no effect on the laminar premixed methane-air flame structure. This is due to the fact that the dilatation term is approximately zero in this computation, with a peak Mach number below 0.001.

(a) 𝜅≠0

(b) 𝜅=0

Figure 3 

Instantaneous contours of the mass fraction of H2O2 (along the 𝑧=1.5 cm plane) at different times (𝑡=0.1, 0.25, and 0.5 millisecond) with (a) and without (b) volume viscosity. The concentric circles correspond, respectively, to three time instances with increasing radii.

(a)

(b)

Figure 4 

Instantaneous profiles of laminar premixed methane-air flame with (circle points) and without (solid line) volume viscosity in the reaction progress variable space, showing the (a) heat release rate and H2O mass fraction and (b) HO2, H2O2, O, and OH radicals at 𝑡=1.0 millisecond.

Figure 5 

Temporal evolution of the laminar fuel consumption and integrated heat release rates for Case 0 with (circle points) and without (solid line) volume viscosity.

The above observations are not surprising and fully confirm the findings in [9], stating that when numerical simulations are performed with identical boundary conditions and at low Mach numbers, the quantitative differences between simulations are extremely small whether or not the volume viscosity is included. Nevertheless, it is interesting to check now the long-time influence of the chaotic fluctuations induced by flow turbulence on this finding.

4.2. Turbulent Flame Structure

4.2.1. Physical Flame Structure

Highly turbulent conditions (Re_𝑡 up to 2460) have been accessed during this project. The effect of turbulence on the physical structure of the flame is first investigated through examination of primitive variables. The temporal evolution of the iso-contours of temperature is shown in Figure 6 for the most turbulent realization (case 4) without accounting for volume viscosity effects, where a cut through the instantaneous iso-contour of the variable illustrates the physical flame structure at different nondimensional times (𝑡=0.005𝜏, 0.25𝜏, 0.5𝜏, 0.75𝜏, 1.0𝜏, and 1.2𝜏). The initially perfectly spherical laminar flame kernel is progressively and heavily distorted and stretched with time by the very strong turbulent field to the extent that local extinction becomes important. From 𝑡>0.5𝜏, the turbulent flame structure is marred with numerous perforations (in the form of burnt gas pockets in the fresh mixture and fresh gas pockets within the burnt gas) and edge flame-like structures [41] appearing within the burnt gas mixture at various locations, shapes, and sizes. Mutual flame annihilation, extinction, and reignition processes lead to a very complex flame topology.

(a) 𝑡=0.005𝜏

(b) 𝑡=0.25𝜏

(c) 𝑡=0.5𝜏

(d) 𝑡=0.75𝜏

(e) 𝑡=1.0𝜏

(f) 𝑡=1.2𝜏

Figure 6 

Time evolution of the iso-contours of temperature (along the 𝑧=2.0 cm plane) for case 4 (without volume viscosity).

The iso-contours of the OH radical closely follow those of the temperature and are shown in Figure 7 at the same nondimensional time 𝑡=1.2𝜏 for cases 1–4 when increasing turbulence intensity (Re_𝑡=615, 1 230, 1 845, and 2 460) for computations where the volume viscosity effects are not accounted for. The OH radical is a widely employed flame marker, often used to define the location of the turbulent flame front [42]. Previously, DNS mostly accessed mild turbulence conditions such as the one shown in Figure 7(a), where the flame is only slightly contorted. When increasing turbulent stirring, flame-flame interactions begin to appear, eventually initiating the formation of fresh gas pockets within the burnt gas mixture as exemplified in Figure 7(b). Moving on to the most turbulent cases shown in Figures 7(c) and 7(d), the amount of wrinkling increases strongly, with considerable structural modifications observed in the form of discontinuous flame fronts, which were hard to figure out from the temperature iso-contours plots on Figure 6. These changes are a first indication that Re_𝑡 should ultimately exceed noticeably 1 000, in order to reach realistic conditions, depending of course on the application and on the corresponding regime for turbulent combustion. The obtained flame topology results from local flame extinctions induced by intense turbulence straining and the combined effect of fresh gas islands creation and flame-flame interactions, which eventually lead to flame pinch-off. At high turbulence, pinch-off and mutual annihilation of flame surface are found to be a dominant mechanism limiting the flame surface area generated by wrinkling due to turbulence. Such flame-flame interactions are very important for combustion device designers, being, for instance, partly responsible for combustion-induced noise [43, 44]. The obtained nonlinear relation between turbulent flame speed and turbulence intensity will be quantified in details in the near future by postprocessing systematically corresponding DNS. The contours of other major and minor species exhibit similar patterns to those of the temperature and OH fields and are therefore omitted in the interest of space.

(a) Re𝑡=615

(b) Re𝑡=1230

(c) Re𝑡=1845

(d) Re𝑡=2460

Figure 7 

Instantaneous iso-contours of OH mass fraction (along the 𝑧=2.0 cm plane) when increasing turbulent Reynolds number Re𝑡 at 𝑡=1.2𝜏. Identical colors are used for all plots.

For a first assessment of the impact of the volume viscosity on the turbulent flame structure, Figure 8 depicts the instantaneous flame front as defined by the iso-contours of the mass fraction of OH for case 3, with and without volume viscosity effects at the same time instance 𝑡=1.2𝜏. As in the laminar comparisons, the differences in both the structure and peak values are negligible.

(a) 𝜅≠0

(b) 𝜅=0

Figure 8 

Instantaneous flame front defined by the iso-contours of the mass fraction of OH (along the 𝑧=2.0 cm plane) for case 3 at 𝑡=1.2𝜏 with (a) and without (b) volume viscosity effects.

4.2.2. Temporal Flame Evolution

Two particularly important global flame quantities—the burning rate 𝑆_𝑐 and the volume-integrated heat release rate 𝐻_𝑟 (introduced earlier for the laminar case)—have been tracked throughout each of the experiments and will be systematically compared to check the impact of volume viscosity on the evolution of the premixed methane flame kernel. 𝑆_𝑐 is an interesting measure of the turbulent burning speed [45]—a parameter of central importance to burner designers relying on premixed turbulent combustion. Here, the burning rate is defined following [46] as the volumetric rate per unit flame area at which the fuel is consumed 𝜔̇_𝑐, and when suitably scaled, serves as a measure of the burning velocity of the flame [47, 48]: 𝑆_𝑐𝑊≡−^𝑓𝜌_𝑓𝑌_𝑓𝜔̇_𝑐𝑑𝑉, (4) where 𝑊_𝑓 is the fuel molecular weight, and 𝜌_𝑓 and 𝑌_𝑓 are the density and mass fraction of the fuel in the fresh gas mixture, respectively. Both quantities are scaled by their respective initial laminar values.

The temporal evolution of these two global quantities is shown in Figure 9 for cases 1–4. The computations with volume viscosity for case 3 and 4 are also shown on the same figure with circled symbols. Initially, 𝐻_𝑟 is approximately constant up to about 𝑡≈0.2𝜏 after which it increases steadily for about 0.1𝜏. At 𝑡≥0.3𝜏 the profiles for the various cases take noticeably different courses depending on the level of turbulence stirring. The turbulent profiles progressively turn exponential with time (at a pace increasing rapidly with increasing Re_𝑡). Increasing the initial turbulent velocity fluctuation from 3 to 12 m/s leads to an increase in the total heat release normalized by its laminar value from 2 to 16 at 𝑡/𝜏=1.4. The observation is similar for the burning rate, as expected for this combustion regime. A clear saturation effect is not yet observed on these DNS results, highlighting the need for simulations at even higher Reynolds numbers and possibly at other mixture equivalence ratios.

(a)

(b)

Figure 9 

Temporal evolution of the scaled (a) fuel consumption and (b) integrated heat release rates at different turbulence intensities for stoichiometric premixed methane-air flames with (circle points: cases 3 and 4) and without (solid line: cases 1–4) volume viscosity effects.

Considering now the effect of volume viscosity, the temporal profiles show no noticeable impact of this term. For case 4 (single DNS realization), slight instantaneous differences appear, probably due to insufficient statistics. This is confirmed by the fact that, for case 3 (for which the simulations have been averaged over six realizations), the profiles are identical in the presence or absence of 𝜅.

4.2.3. Conditional Analysis

Conditional analysis is of central importance for many turbulent combustion models. Hence, conditional mean values have been computed for different variables characterizing flame behavior as a function of turbulence intensity (Figure 10). The conditional mean of the progress variable gradient |∇𝐶| is of particular interest, since its inverse gives a measure of the flame thickness in analogy to the thermal flame thickness [3]. The conditional profiles of mean |∇𝐶| and H₂O₂ mass fraction show a noticeable dependency on turbulence intensity, with progressively lowered peaks. The decrease of the conditional mean |∇𝐶| reveals that flame thickening is predominant, especially around 0.1≤𝐶≤0.7 (the active flame region) when increasing Re_𝑡. Conditional mean profiles for the most sensitive minor radicals (H₂O₂) as well as major species (not shown) are indifferent in the presence of volume viscosity, irrespective of ensemble averaging. Volume viscosity effects are noticeable on the mean |∇𝐶| profile for case 4 (single DNS realization). Again, this appears to be a statistical artifact, since all conditional mean profiles for case 3 (averages over several realizations) show no influence of 𝜅, whatsoever.

(a)

(b)

Figure 10 

Conditional profiles of mean (a) progress variable gradient, |∇𝐶| and (b) mass fraction of H2O2, conditioned on the progress variable 𝐶 at 𝑡=1.2𝜏.

5. Conclusions

Progress in numerical techniques as well as computational power now allows quantitative investigations of turbulent reacting flows for increasingly realistic conditions including detailed physicochemical models. Direct Numerical Simulations of stoichiometric premixed methane flames have been considered in a parametric study including a case with six pairs of computations for Re_𝑡=615 up to 2 460 (factor 4) by accessing high-end parallel computers at European level. After solving issues associated with I/O and initial turbulence generation, a detailed analysis has been performed using the dedicated Matlab-based library AnaFlame [39].

The flame speed is found to rapidly increase with increasing turbulence stirring, associated with a noticeable flame thickening. Simultaneously, peak conditional profiles of heat release and major and minor species mass fractions are systematically lowered with increasing turbulence intensity. In all considered cases, premixed methane flames show a negligible impact of volume viscosity. No differences at all are found in laminar computations, confirming theoretical findings for low Mach number conditions. Previous observations for turbulent hydrogen flames, revealing a clear impact of volume viscosity, do not apply for the present methane flames at all Reynolds numbers considered in the study. To save computing resources, the inclusion of volume viscosity effects in multicomponent multidimensional turbulent flame computations burning methane and probably higher hydrocarbons as well is therefore discouraged, as long as low Mach numbers are considered.

This study demonstrates also the importance of repeating DNS realizations in order to obtain statistically significant data. Single realizations might lead to spurious discrepancies, rapidly smoothed out when averaging over several results, as observed here when comparing the findings for case 3 (several DNS realizations) and case 4 (single DNS realization).

Acknowledgments

The computational resources have been provided by the DEISA Extreme Computing Initiative (DECI), as part of the 7th Framework program financed by the European Union. The authors are grateful for the support of the DEISA help desk.