#### 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 (), 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 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):
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:
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 (*EGSK6()*) for the computations in the present study. The simplest and cheapest of the models (, *EGSK2()*) considers only the transport system matrix associated with the internal energy, whereas (*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 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 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_{4}, O_{2}, H_{2}, H_{2}O, CH_{2}O, CO, CO_{2}, HO_{2}, OH, H, O, CH_{3}, HCO, H_{2}O_{2}, CH_{3}O, N_{2}) [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_{2} and higher carbon-chain reactions, the reason why for the present study.

The initial system is a stationary hot ( K) perfectly spherical laminar flame-kernel of initial radius 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 K. The initial profile of any primitive variable is prescribed according to 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 and at outside the kernel, and and at within. An appropriate nitrogen complement is added everywhere.

To investigate the influence of the turbulent Reynolds number 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 , the eddy turn-over time used to quantify flame/turbulence interaction for time-decaying turbulence [35], and the Karlovitz number are given in Table 1. The characteristic viscosity of the mixture, m^{2}/s, the laminar flame speed m/s, the laminar flame thermal thickness mm, and the integral scale are all constant for all simulations presented here.

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.

The general view of the configuration is illustrated in Figure 2, where the iso-surface of temperature with H_{2}O_{2} slices and CO_{2} 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)**

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 . The longest computation requires about two months. For the results presented below, a nondimensional time , 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 ( and ) 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 ( cm) with the same initial mixture composition as in the turbulent cases. The instantaneous iso-contours of the mass fraction of H_{2}O_{2} at different times (, 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_{2}O, H_{2}O_{2}, and HO_{2} are shown at 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_{2}O_{2} and HO_{2} 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) |

(b) |

**(a)**

**(b)**

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 ( 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 (, , , , , and ). 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 , 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) |

(b) |

(c) |

(d) |

(e) |

(f) |

The iso-contours of the OH radical closely follow those of the temperature and are shown in Figure 7 at the same nondimensional time for *cases 1*–*4* when increasing turbulence intensity (, 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 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) |

(b) |

(c) |

(d) |

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 . As in the laminar comparisons, the differences in both the structure and peak values are negligible.

(a) |

(b) |

###### 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]: 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 after which it increases steadily for about . At 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 ). 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 . 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)**

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_{2}O_{2} 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 (the active flame region) when increasing . Conditional mean profiles for the most sensitive minor radicals (H_{2}O_{2}) 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)**

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