Premixed Turbulent Flames

1 Department of Applied Mechanics, Chalmers University of Technology, 41296 Gothenburg, Sweden 2 School of Mechanical Engineering, The University of Leeds, Woodhouse Lane, Leeds LS2 9JT, UK 3 Center for Interdisciplinary Studies on Resource Recovery and Refinery in Asia (CISRA), EcoTopia Science Institute, F3-4 (670), Nagoya University, Nagoya 464-8603, Japan 4 ONERA/DEFA, Centre de Palaiseau, Chemin de la Hunière, 91761 Palaiseau, France

Turbulent combustion of premixed gases is widely used for energy conversion in transportation, for example, piston engines in cars and light aircrafts, aero-engine afterburners, and stationary power generation, for example, gas turbines. Current environmental and efficiency challenges call urgently for new solutions that satisfy stringent requirements for ultra low emissions as well as highly efficient combustion technology. To make any significant progress in this area, both fundamental and applied studies of numerous unresolved issues relevant to premixed turbulent flames are strongly necessary. This special issue contains seven papers aimed at contributing to fundamental understanding of premixed turbulent combustion as well as to more efficient use of this process in future engines.
Since the early days of research into premixed turbulent flames, turbulent burning velocity or flame speed was in the focus of numerous experimental, theoretical, and numerical studies. Despite long-term extensive research, certain basic issues have not yet been resolved and some of them are addressed in the first three contributions. In the paper entitled "Determination of the burning velocity domain of a statistically stationary turbulent premixed flame in presence of counter-gradient transport," the classical academic problem of the burning velocity in a fully developed statistically stationary, planar, and one-dimensional turbulent flame is theoretically analyzed in order to investigate the influence of countergradient scalar transport (i.e., transport of combustion products towards the burned side of the flame, caused by pressure gradient induced by the heat release) on the burning velocity. In particular, it is demonstrated that, even if turbulent scalar transport shows the countergradient behavior at the cold boundary, "there still exists a possibility of observing a steady regime of propagation." The paper entitled "Burning rate in impinging jet flames" deals also with the problem of countergradient turbulent scalar transport by invoking a recently developed simple model of this phenomenon in order to validate a method for evaluating turbulent burning velocity proposed earlier based on a theoretical consideration. The influence of hydrodynamic effects caused by heat release on turbulent burning velocity is also discussed in the paper entitled "Influence of turbulent scalar mixing physics on premixed flame propagation," which highlights another effect of that kind, that is, dilatation, rather than countergradient transport.
Direct numerical simulation (DNS) is a rapidly developing powerful numerical tool for fundamental research into turbulent flows and turbulent flames, in particular. The present special issue contains three contributions relevant to DNS of premixed turbulent combustion. In the paper entitled "Effects of turbulent Reynolds number on the displacement speed statistics in the thin reaction zones regime of turbulent premixed combustion," a DNS database obtained from premixed turbulent flames characterized by different Reynolds numbers is analyzed in order to investigate statistical characteristics of the instantaneous flame front (its local displacement speed, surface density, curvature, and strain rate) as well as relations between these statistical characteristics and the influence of the Reynolds number on them. In the paper entitled "Statistics of conditional fluid velocity in the corrugated flamelets regime of turbulent premixed combustion: a direct numerical simulation study," a single DNS database is considered in order to study the behavior of the first and second moments of turbulent velocity field, conditioned either to unburned mixture, or to burned mixture, or to the instantaneous flame front. In particular, it is shown that flow velocities conditioned within unburned or unburned side of the front are approximately equal in the largest part of turbulent flame brush with the exception of its leading edge. In the paper entitled "Direct numerical simulations of the impact of high turbulence intensities and volume viscosity on premixed methane flames," DNS data obtained from highintensity turbulent premixed flames are reported, and it is argued that, contrary to hydrogen-air flames, volume viscosity weakly affects methane-air flames and may be ignored in DNS research into turbulent combustion of hydrocarbon-air mixtures.
While the above six papers address basic issues of premixed turbulent combustion either theoretically or numerically, the paper entitled "Effects of turbulence on the combustion properties of partially premixed flames of canola methyl ester and diesel blends" is an example of an applied experimental study aimed at developing clean combustion technology by utilizing an alternative biofuel canola methyl ester (CME). The paper reports measured data on the influence of turbulence on emissions from partially premixed flames of blends of CME and a Diesel fuel.

Introduction
Mean fluid velocities conditional in reactants, products, and instantaneous flame surface often play a key role in the theoretical analysis and modelling of turbulent premixed combustion [1][2][3][4][5][6][7]. Theoretical analysis of Bray et al. [1] demonstrated the relation between the turbulent scalar fluxes ρu i c of reaction progress variable c and the difference in the mean conditional velocities in products and reactants ((u i ) P − (u i ) R ), where the overbar signifies a Reynolds averaging operation, q = q − q, q = q − q, and q = ρq/ρ are the Reynolds fluctuation, Favre fluctuation, and Favre mean of a general quantity q, u i is the ith component of velocity, ρ is the gas density, and the subscripts R and P are used to refer to conditional mean quantities in reactants and products, respectively. Domingo and Bray [4] showed the relation between the slip velocity (u i ) P − (u i ) R and the pressure transport terms in the Reynolds stress ρu i u j transport equation. In addition to the conditional mean velocities (i.e. (u i ) R and (u i ) P ), the surface averaged quantities play a pivotal role in the closure of turbulent premixed flames using the Flame Surface Density (FSD) [3,5,8,9] or conditional mean equations [6,7,[10][11][12][13].
In the case of an asymptotically high Damköhler number Da, which characterises a ratio of the large-scale turbulent time scale to the chemical time scale, the conditional mean velocities may be extracted from the Favre mean velocities, which are readily available in conventional Reynolds Averaged Navier Stokes (RANS) simulations. Under Da → ∞, the probability density function (pdf) of reaction progress variable c could be approximated by a bimodal distribution and analytical expressions for conditional-velocity-related quantities could easily be derived based on a presumed bimodal pdf approach [1,2]. In many practical flames, the underlying combustion process can be characterised by a high but finite Da within the corrugated flamelets regime 2 Journal of Combustion on the regime diagram [14,15], where the flame thickness remains smaller than the Kolmogorov length scale and the pdf of reaction progress variable c is expected to be bimodal.
However, in spite of the key importance of conditional velocities in turbulent premixed combustion modelling, there is yet to be a detailed Direct Numerical Simulation (DNS) based assessment of the validity of the expressions of the mean conditional velocities, which can be derived based on a presumed bi-modal pdf of c, in the corrugated flamelets regime combustion.
In this respect the main objectives of this study are as follows.
(1) To analyse the statistical behaviours of the conditional mean velocity and their relation to the Favre mean velocity in the context of RANS for the corrugated flamelets regime combustion.
(2) To analyse the statistical behaviours of the conditional surface averaged velocities and their relation to the conditional mean velocities for the corrugated flamelets regime combustion.
The rest of the paper will be organised as follows. The necessary mathematical background will be provided in the next section of the paper. This will be followed by a brief discussion on numerical implementation. Following this, results will be presented and subsequently discussed. Finally the main findings will be summarised and conclusions will be drawn.

Mathematical Background
DNS of turbulent combustion should account for both threedimensionality of turbulence and detailed chemical mechanism. However, the limitation of computer storage capacity until recently restricted DNS either to two dimensions with detailed chemistry or to three dimensions with simplified chemistry. As turbulent velocity field is inherently threedimensional in nature and vortex stretching mechanism is absent in two dimensions, the second approach takes precedence in the present study, which is based on a single-step irreversible Arrhenius-type chemistry.
In premixed combustion the species field is often represented in terms of a reaction progress variable c which rises monotonically from zero in pure reactants to unity in fully burned products. The reaction progress variable c can be defined in terms of the mass fraction of a suitable reactant Y R in the following manner: where subscripts 0 and ∝ indicate values in pure reactants and in fully burned products, respectively. According to BML analysis [1,2], the pdf of reaction progress variable c and the joint pdf of velocity vector u and c are given by P c; x = α c x δ(c) + β c x δ(1 − c) P u, c; x = α c x P R u; 0; x δ(c) + β c x P P u; 1; x δ(1 − c) where α c , β c , and γ c are the coefficients for the progress variable pdf, P R ( u; x) and P P ( u; x) refer to pdfs of u in unburned gases and fully burned products, respectively, and the functions f 1 and f 2 originate due to burning gases at the interior of the flame. According to BML analysis [1,2], it is assumed that P R ( u; x), P P ( u; x), and f 2 ( u, c; x) are normalised in such a manner that they individually integrate to unity. The last term on the right-hand side of (2a) and (2b) refers to the contribution of burning fluid. Based on this pdf, one obtains 1 0 P(c)dc = α c + β c + O γ c = 1, where ρ b is the mean burned gas density, ρ 0 is the unburned gas density, and ρ c is the mean value of density conditional on c. For low Mach number adiabatic unity Lewis number flames ρ c and ρ b are given by ρ c = ρ 0 (1 + τc) , where τ = (T ad − T 0 )/T 0 is the heat release parameter with T ad and T 0 being the adiabatic flame temperature and unburned gas temperature respectively. Under high Damköhler number (i.e., Da 1) the contribution of O(γ c ) becomes negligible in (3), which yields the following expressions of α c and β c : Based on (2b) and (5) one obtains the following relation for the Favre mean velocity components u i [1,2]: The conditional mean values in unburned reactants and fully burned products (i.e., (q) R and (q) P ) for a general quantity q dependent on u are evaluated as where ε is a small number (i.e., 0 < ε 1). It is possible to obtain the following relation for turbulent scalar flux ρu i c using (2b), (5), and (6a) [1,2]: Using (6a)-(6b) and (7) one obtains the following relations for (u i ) R and (u i ) P : Using (2b) it is possible to express ρu i u j and ρu i u j c in the following manner: The expressions for Reynolds stresses conditional in reactants and products (i.e., (u i u j ) R and (u i u j ) P ) can be obtained using (7) and (9) in the following manner: The surface averaged value of fluid velocity component based on the generalised Flame Surface Density (i.e. Σ gen = |∇c|) [8,9] can be expressed in the following manner: As proposed in [10] based on the DNS data by Im et al. [6], the simplest possible model for the surface averaged velocity components conditional on unburned gas side of flamelets (i.e. (u i ) Rs ) can be given as follows: In the reference frame attached to the flame, the conditional mean velocity (u 1 ) R approaches the turbulent flame speed S T (i.e., (u 1 ) R → u 1 = S T ) at the leading edge for a hypothetical fully developed, constant-density, statistically planar, one-dimensional flame that propagates from right to left. Under this condition, the surface averaged velocity component (u 1 ) Rs approaches S L (i.e., (u 1 ) Rs → S L < S T ) at the leading edge because the instantaneous propagation velocity (u 1 − S L ) of a flamelet should be equal to zero in the selected reference frame under steady state (otherwise the flamelet would cross the leading edge, and that is impossible by definition). Therefore, (11b) yields wrong result at c → 0 in this simple case. For these reasons, the following linear interpolation has also been proposed in [10]: Here, σ = (1 + τ) and M i = ∂ c /∂x i /|∇ c | in which c is either equal to c or equal to c. This model effectively ensures that (u i ) Rs → (u i ) P /σ when either c → 0 or c → 0. Similarly (u i ) Rs → (u i ) R when either c → 1 or c → 1 according to (11c). In deriving (11c), it is assumed [10] that, at the leading and trailing edges, (i) the flamelet is parallel to the flame brush (i.e., M i = N i = ∂c/∂x i /|∇c|), (ii) the velocity component tangential to the flamelet is assumed to be unaffected (i.e. (u i ) R − N i N j (u j ) R = (u i ) P − N i N j (u j ) P ), and (iii) the velocity component normal to the flamelet is taken to be equal to laminar burning velocity (i.e., (u i ) R N i → S L and (u i ) P N i → σS L ). The simplest possible model for the conditional surfaceweighted mean of the correlation between velocity components (i.e., (u i u j ) Rs ) can be proposed as [10] (u i u j ) Rs = (u i u j ) R , whereas the following linear interpolation [10] ensures the correct behaviour of (u i u j ) Rs → (u i u j ) P at c → 0 for a constant density "flame":

Numerical Implementation
For the present study, a decaying turbulence DNS database of statistically planar freely propagating turbulent premixed flame representing the corrugated flamelets regime combustion has been considered. The initial values of normalised root-mean-square value of turbulent velocity fluctuation u /S L , integral length scale normalised by the thermal flame thickness l/δ th , heat release parameter τ = (T ad − T 0 )/T 0 , global Lewis number Le, Damköhler number Da = l · S L /u δ th , Karlovitz number Ka = (u /S L ) 3/2 (l/δ th ) −1/2 , and 4

Journal of Combustion
Re t = ρ 0 u l/μ 0 are given in Table 1, where S L is the unstrained planar laminar burning velocity, and δ th = (T ad − T 0 )/ Max |∇ T| L is the thermal laminar flame thickness, ρ 0 and μ 0 are the unburned gas density and viscosity respectively. The subscript L is used to denote the quantities in the unstrained planar steady laminar flame. It is evident from Table 1 that the Karlovitz number Ka remains smaller than unity in the case considered here which suggests that combustion in this case belongs to the corrugated flamelets regime [14,15]. Standard values are chosen for Zel'dovich number β = T ac (T ad − T 0 )/T 2 ad , Prandtl number Pr, and the ratio of specific heats (i.e., β = 6.0, Pr = 0.7, and γ = C P /C V = 1.4) for the case considered in the present study, where T ac is the activation temperature.
The DNS database considered here is taken from [16] where the simulation is carried out using low Mach number assumption. A domain of size 118.64α T0 /S L × 131.65α T0 /S L × 131.65α T0 /S L is discretised by a Cartesian grid of 261 × 128 × 128 with uniform grid spacing in each direction where α T0 is the thermal diffusivity in the unburned gas. In this case, inlet and outlet boundaries are specified in the mean direction of flame propagation. Realistic turbulent velocity fluctuations are specified at the inlet boundary using Taylor's hypothesis by scanning a plane through a precomputed box of frozen turbulence. The inlet mean flow velocity is taken to be equal to S L . The outlet boundary is taken to be partially nonreflecting in nature and is specified according to Navier Stokes Characteristic Boundary Condition (NSCBC) formulation [17]. The transverse directions are assumed to be periodic.
A sixth-order finite difference scheme is used for evaluating spatial derivatives in the direction of flame propagation. The spatial derivatives in the transverse direction are evaluated using a pseudospectral method. The density change is accounted for by the relation between density and reaction progress variable according to the BML formulation for the unity Lewis number flames [1,2]. In this case, the initial velocity field is specified by an initially homogeneous isotropic turbulence using a pseudo-spectral method following the Batchelor-Townsend spectrum [18]. The time advancement for viscous and diffusive terms in this case is carried out using an implicit solver based on Crank-Nicholson scheme whereas the convective terms are time advanced with the help of a third-order low storage Runge-Kutta method [19]. The flame is initialised using a steady unstrained planar laminar flame solution. The grid resolution is determined by the resolution of the flame structure. About 10 grid points are placed within the thermal flame thickness δ th .
In this case, flame-turbulence interaction takes place under decaying turbulence. The simulation in this case was run for about 4 initial eddy turn over times (t sim ∼ 4τ f = 4l/u ) which is about 27-fold greater than the chemical time scale t c = δ th /S L . It is recognised that the simulation time remains small but previous studies with similar or smaller simulation time provided useful information in the past [8,9,[20][21][22][23][24][25][26]. The values of u /S L in the unburned reactants ahead of the flame at the time when the statistics were extracted decreased by about 78%. The value of l/δ th has increased from its initial value by a factor of about 1.62, but there are still enough turbulent eddies on each side of the computational domain. For the final values of u /S L and l/δ th the underlying combustion situation belongs to the corrugated flamelets regime according to the combustion regime diagram [14,15]. For the purpose of postprocessing DNS data, the mean quantities are assumed to be a function of the coordinate in the direction of mean flame propagation (x 1 direction in the present case). The Reynolds/Favre averaged quantities are evaluated based on averaging the relevant quantities in transverse directions (x 2 − x 3 planes). The statistical convergence of the Reynolds/Favre averaged quantities is checked by comparing the corresponding values obtained using half of the sample size in the transverse directions with those obtained based on full available sample size. Both the qualitative and quantitative agreements between the values obtained based on full and half sample sizes are found to be satisfactory. In the next section the results obtained based on full sample size available in transverse directions will be presented for the sake of brevity.
In the present study the conditional mean quantities are evaluated using (6b), where ε is taken to be 0.1 following the previous analysis by Domingo and Bray [4]. This essentially suggests that the mean quantities conditional in reactants are evaluated using the data corresponding to 0 < c < 0.1. Similarly the mean quantities conditional in products are evaluated using the data corresponding to 0.9 < c < 1. A smaller value of ε yields qualitatively similar results to those obtained using ε = 0.1 but the sample size for evaluating the conditional mean values decreases with decreasing ε. By contrast, increasing ε value (e.g., 0.15 ≥ ε ≥ 0.1) gives rise to similar qualitative trends to obtained using ε = 0.1 but increasing the value of ε by a large margin is likely to lead to q R and q P values which are not representative of conditional means in reactants and products, respectively. The quantitative agreement between the conditional velocity statistics obtained for 0.05 ≤ ε ≤ 0.15 is indeed found to be excellent.

Results and Discussion
The contours of the reaction progress variable for the flame considered here are plotted in Figure 1(a) which shows that c contours remain mostly parallel to each other because the flame thickness remains smaller than the Kolmogorov length scale in this case and thus turbulent eddies cannot penetrate into the flame structure while background turbulence only wrinkles the flame. This is consistent with the behaviour expected in the corrugated flamelets regime [15]. However, a careful examination of Figure 1 of progress variables for small values of c (i.e., c < 0.2) are not perfectly parallel to the contours representing the reaction zone (i.e. 0.7 ≤ c ≤ 0.9) for the present thermochemistry. This suggests that the preheat zone is relatively more affected by turbulence than the reaction zone, indicating that some eddies of the order of the Kolmogorov length scale may perturb the preheat zone which is likely to give rise to minor departure from pure bimodal distribution of P(c; x) and P( u, c; x).
The nature of pdf of c can be characterised in terms of the variance of reaction progress variable c 2 . For a bimodal distribution of c, the variance of reaction progress variable c 2 is given by The variations of c 2 with c are shown in Figure 1(b) which indicates that c 2 remains close to c(1 − c) for this case. This implies that the assumption of bimodal distribution of c, as done in BML analysis [1,2], holds reasonably well for the case considered here. The pdf of c at the location corresponding to c = 0.5 is shown in Figure 1(c), (a)  which clearly demonstrates that the bimodal distribution of c satisfactorily approximates the pdf of c in the flame considered here. However, Figure 1(b) clearly indicates that 7 The increase in u 1 /S L due to thermal expansion and the effects of flame propagation are responsible for the increase in (u 1 ) R /S L with increasing c. However, the comparison between (u 1 ) R /S L and u 1 /S L in Figure 2(a) reveals that there is a noticeable difference between (u 1 ) R /S L and u 1 /S L . The variations of (u 2 ) R /S L and (u 3 ) R /S L are not shown because these components are identically equal to zero for statistically planar flames. It is evident from Figure 2(a) that the relation given by (8a) satisfactorily captures the variation of u R with c for the flame considered here. As (8a) is derived based on the assumption of the bimodal pdf of c, this expression satisfactorily captures (u i ) R /S L for this flame where P(c) is roughly bimodal in nature.
The variation of (u 1 ) P /S L with c is shown in Figure 2(b). In this case (u 1 ) P /S L also increases from unburned gas side to the burned gas side of the flame brush and there is a noticeable difference between (u 1 ) P /S L and u 1 /S L in this flame. The effects of increase in u 1 /S L due to thermal expansion and flame propagation are also responsible for the increase of (u 1 ) P /S L from unburned gas side to the burned gas side of the flame brush. The variations of (u 2 ) P /S L and (u 3 ) P /S L are not shown because these components are identically equal to zero for statistically planar flames. In this case the pdf of c remains roughly bimodal and thus (8b) satisfactorily predicts (u 1 ) P /S L obtained in the DNS (see Figure 2(b)).
Comparing Figures 2(a) and 2(b) it can be seen that (u 1 ) P /S L remains greater than (u 1 ) R /S L for this case. This behaviour can be observed clearly from the variation of [(u 1 ) P − (u 1 ) R ]/S L with c across the flame brush as shown in Figure 2(c). According to (7) a positive (negative) value of [(u 1 ) P − (u 1 ) R ]/S L leads to a countergradient (gradient) transport of turbulent scalar flux ρu 1 c , which can be substantiated from the variations of ρu 1 c × ∂ c/∂x 1 shown in Figure 2(d), where a positive (negative) value implies countergradient (gradient) type transport. In the present case, the sign of the normalised slip velocity [(u 1 ) P − (u 1 ) R ]/S L indicates the direction of turbulent scalar flux. A similar observation was made earlier by Veynante et al. [24] and Chakraborty and Cant [27]. The modelling of ρu i c has been discussed in [24,[27][28][29] and thus will not be further addressed in this paper for the sake of conciseness. It can be seen from Figure 2(c) that the value of [(u 1 ) P − (u 1 ) R ]/S L remains close to, but smaller than, τ(∼ 0.8τ) for a major portion of the flame brush (i.e., 0.2 ≤ c ≤ 0.8). This behaviour is found to be consistent with the experimental and DNS based findings reviewed in [30,Section 3.3]. Veynante et al. [24] proposed a simple estimate for the normalised slip velocity as [(u 1 ) P − (u 1 ) R ]/S L = [−2α 2 k/3 + τS L ]/S L , where α is an appropriate efficiency function which assumes a magnitude of the order of unity. The second term in the square bracket on the right-hand side is associated with the velocity jump induced by the pressure gradient along the normal to the laminar flame. For the present database the velocity jump τS L remains much greater than 2 k/3 (i.e., τS L 2 k/3), which ultimately gives rise to [(u 1 ) P − (u 1 ) R ]/τS L ∼ O(1) for the major part of the flame brush (see Figure 2(c)). Because τS L 2 k/3, the fact that the ratio of [(u 1 ) P − (u 1 ) R ]/τS L ∼ O(1) is smaller than unity in the present case is associated with (N 1 ) < 1 due to random orientations of the instantaneous flame normal vectors with respect to the x 1 -axis. As ρu 1 c remains positive for the major part of the flame brush (see Figure 2(c)), (u 1 ) R /S L ((u 1 ) P /S L ) assumes smaller (greater) values than u 1 /S L in the current case (see Figures 2(a) and 2(b)).
The variations of (u 1 u 1 ) R /S 2 L with c are shown in Figure 3(a). It is evident from Figure 3(a) that (u 1 u 1 ) R /S 2 L decays from unburned gas side to burned gas side in the present case. The variations of ρu 1 u 1 /ρS 2 L and the prediction of (10a) are also shown in Figure 3(a). The difference between ρu 1 u 1 /ρS 2 L and (u 1 u 1 ) R /S 2 L is significant for this case, which suggests that the contribution The variations of (u 2 u 2 ) R /S 2 L and ρu 2 u 2 /ρS 2 L with c are shown in Figure 3(b). As (u 2 u 2 ) R /S 2 L and ρu 2 u 2 /ρS 2 L are statistically similar to (u 3 u 3 ) R /S 2 L and ρu 3 u 3 /ρS 2 L , respectively, the variations of (u 3 u 3 ) R /S 2 L and ρu 3 u 3 /ρS 2 L are not shown here for the sake of conciseness.
The evolution of ρu 1 u 1 and ρu 2 u 2 within a statistically planar one-dimensional turbulent flame brush (that propagates from right to left) is affected by several physical mechanisms [30][31][32], which can be evident from the following transport equations: The transport equation for ρu 3 u 3 is not shown here because it takes the same form as that of (14b) because x 2 and x 3 directions are statistically similar. The statistical behaviours of the various terms of (14a) and (14b) for this database have been discussed elsewhere [31,32] in detail and will be not be repeated here for the sake of brevity. The main behavioural trends will be summarised here. The mean dilatation term

Equation 10a
(b) Figure 3: (a) Variations of (u 1 u 1 ) R /S 2 L and ρu 1 u 1 /ρS 2 L with c across the flame brush along with the predictions of (10a). (b) Variations of (u 2 u 2 ) R /S 2 L and ρu 2 u 2 /ρS 2 L with c across the flame brush along with the predictions of (10a).
T 1 is a sink in (14a) and acts to reduce ρu 1 u 1 but this term is identically equal to zero for (14b) in the context of statistically planar flames. The mean pressure gradient T 2 acts as a source term in (14a) for the present database due to countergradient transport (u 1 > 0), which serves to increase ρu 1 u 1 . However, the term T 2 vanishes for (14b) in the case of statistically planar flames. The pressure dilatation term T 3 acts as a source term in both (14a) and (14b). The term T 4 acts as a sink due to viscous dissipation which serves to decrease ρu 1 u 1 and ρu 2 u 2 . The terms T 5 and T 6 denote the dispersion due to pressure fluctuation and velocity fluctuation, respectively, and these terms remain small in comparison to the magnitude of the contributions of T 2 , T 3 , and T 4 (T 3 and T 4 ) in the context of ρu 1 u 1 (ρu 2 u 2 ) transport. The generation mechanisms associated with T 2 and T 3 in the context of ρu 1 u 1 transport dominate over the dissipative actions of T 4 and T 1 in the middle of the flame brush, whereas the contributions of T 4 and T 1 dominate at the leading and trailing edges. As a result of this, ρu 1 u 1 decreases from the leading edge of the flame brush before increasing significantly at the middle of the flame brush, and it eventually decreases as the burned gas side is approached. By contrast, T 3 and T 4 roughly balance each other for the major part of the flame brush in the context of ρu 2 u 2 transport, but the dissipative action of T 4 dominates over the contribution of T 3 towards the leading and trailing edges of the flame brush. This gives rise to the decay of ρu 2 u 2 from the leading edge of the flame brush but it does not change significantly for the major part of the flame brush before decaying again towards the burned gas side. The smaller values of (u 1 u 1 ) R and (u 2 u 2 ) R in comparison to ρu 1 u 1 /ρ and ρu 2 u 2 /ρ, respectively, are associated with the influence of the second and third terms on the right-hand side of (10a). For instance, the third term on (10a) (i.e., −(ρu 1 c ) 2 /ρ 2 (1 − c) 2 ) is negative and its magnitude increases with c when (10a) is written for (u 1 u 1 ) R . However, the contribution of the third term on the right-hand side (i.e. −(ρu 2 c ) 2 /ρ 2 (1 − c) 2 ) vanishes when (10a) is written for (u 2 u 2 ) R and the difference between (u 2 u 2 ) R and ρu 2 u 2 /ρ in Figure 3(b) is significantly smaller than the difference between (u 1 u 1 ) R and ρu 1 u 1 /ρ in Figure 3(a). It is surprising that poor agreement (at c > 0.6) between the DNS data and (10a) is observed in this case that is, under the approximate conditions for that the equation was derived. It is evident from Figure 3(a) that (10a) in this case significantly overpredicts the magnitude of (u 1 u 1 ) R /S 2 L towards the burned gas side of the flame brush. Although the pdf of c remains roughly bimodal (see Figure 1(c)), the departure from bimodal distribution can be substantiated from Figure 1(b) which shows that c 2 remains slightly smaller than c(1 − c). The departure from bimodal distribution seems to have major influences on the predictive capability of (10a) because the reactive contribution to (u 1 u 1 ) R is not accurately represented by −ρu 1 2 when P(c) departs from the bimodal distribution.
The variations of (u 1 u 1 ) P /S 2 L and ρu 1 u 1 /ρS 2 L with c are shown in Figure 4(a) along with the prediction of (10b). It can be seen from Figure 4

Equation 10b
(b) Figure 4: (a) Variations of (u 1 u 1 ) P /S 2 L and ρu 1 u 1 /ρS 2 L with c across the flame brush along with the predictions of (10b). (b) Variations of (u 2 u 2 ) P /S 2 L and ρu 2 u 2 /ρS 2 L with c across the flame brush along with the predictions of (10b).
the magnitude of (u 1 u 1 ) P /S 2 L for the major portion of the flame brush. Figure 4(a) also shows that the variation of (u 1 u 1 ) P /S 2 L remains markedly different from ρu 1 u 1 /ρS 2 L for the case considered here due to the contribution of The variations of (u 2 u 2 ) P /S 2 L and ρu 2 u 2 /ρS 2 L with c are shown in Figure 4(b). The variations of (u 3 u 3 ) P /S 2 L and ρu 3 u 3 /ρS 2 L with c are both qualitatively and quantitatively similar to (u 2 u 2 ) P /S 2 L and ρu 2 u 2 /ρS 2 L , respectively, and thus are not shown here for the sake of brevity.
Figures 3(b) and 4(b) indicate that the predictions of both (10a) and (10b) agree reasonably well with the corresponding quantities obtained from DNS data for the major part of the flame brush. Although the differences (u 2 u 2 ) R − ρu 2 u 2 /ρ and (u 2 u 2 ) P − ρu 2 u 2 /ρ are of the same order, the relative difference [(u 2 u 2 ) R − ρu 2 u 2 /ρ]/(u 2 u 2 ) R is significantly greater than the relative difference [(u 2 u 2 ) P − ρu 2 u 2 /ρ]/(u 2 u 2 ) P . In an actual RANS simulation the quantities ρu i c and ρu i u j c are modelled and the accuracy of their modelling is also likely to affect the performance of the models given by (10a) and (10b). Modelling of ρu i c and ρu i u j c has been discussed elsewhere [1,24,27,30,33,34] and thus will not be repeated here for the sake of brevity.
The variations of (u 1 ) Rs /S L and (u 1 ) R /S L with c along with the predictions of (11c) according to c = c and c = c are shown in Figure 5(a). It is worth stressing that these velocities were evaluated in the coordinate framework selected so that the mean velocity of unburned gas upstream the flame brush is equal to S L . It is evident from Figure 5(a) that the simplest model given by (11b) [10] captures the behaviour of (u 1 ) Rs /S L for the major portion of the flame brush, which is in agreement with earlier DNS results [6,7]. Figure 5(a) show that the performance of (11c) is comparable for c = c and c = c and the model given by (11c) underpredicts (u 1 ) Rs /S L for the case considered here. It is worth reminding that (11c) is a linear interpolation between two limit cases, that is, (u 1 ) P /σ at c = 0 (or c = 0) to (u 1 ) R at c = 1 (or c = 1). Figure 5(a) indicates that the linear interpolation does not capture the behaviour of (u 1 ) Rs /S L in the major part of the flame brush, but, at the leading edge of the flame brush (0 < c 1), (11c) agrees with the DNS data in a better manner than the simple model given by (11c).
Lee and Huh [7] proposed a model for (u 1 ) Rs in the following manner: where K 1 (∂c/∂x 1 ) is a tuning parameter that is proportional to but significantly larger than the kinematic eddy viscosity C μ ( k 2 / ε) R conditional in unburned gas [7] where C μ = 0.09 is a model constant. Lee and Huh [7] found that the ratio K 1 (∂c/∂x 1 )/[C μ ( k 2 / ε) R ] was equal to 6.7 and 20 in the two flames studied by them. The prediction of (17) for K 1 = 2C μ ( k 2 / ε) R /(∂c/∂x 1 ) is also plotted in Figure 5(a), which shows that K 1 ∂ 2 c/∂x 2 1 remains much smaller than (u 1 ) R and (u 1 ) Rs and the prediction of (15) remains comparable to (11c) for the major part of turbulent flame brush with the exception of the leading edge. 2C μ ( k 2 / ε) R (∂ 2 c/∂x 2 1 )/[(∂c/∂x 1 )S L ] with c are plotted in Figure 5(b) which indicates that K 1 (∂ 2 c/∂x 2 1 )/S L = 2C μ ( k 2 / ε) R (∂ 2 c/∂x 2 1 )/[(∂c/∂x 1 )S L ] satisfactorily captures the qualitative behaviour of [(u 1 ) R − (u 1 ) Rs ]/S L . The quantitative agreement between K 1 (∂ 2 c/∂x 2 1 satisfactory towards the unburned gas side of the flame brush but the quantitative agreement is not satisfactory towards the burned gas side of the flame brush. The variations of (u 1 ) Ps /S L and (u 1 ) P /S L with c across the flame brush are shown in Figure 5(c) which demonstrates that (u 1 ) Ps remains close to but smaller than (u 1 ) P , which is also consistent with DNS results of Im et al. [6]. A physical mechanism that makes (u 1 ) P greater than (u 1 ) Ps is associated with the acceleration of burned gas under the action of the mean pressure gradient from the flame surface to the burned gas side of the turbulent flame brush. This physical mechanism has been discussed in detail elsewhere [12] and will therefore not be elaborated here for the sake of brevity.
The variations of (u 1 u 1 ) Rs /S 2 L and (u 1 u 1 ) R /S 2 L = [(u 1 ) R (u 1 ) R + (u 1 u 1 ) R ]/S 2 L with c are shown in Figure 6(a) along with the predictions of (12a) and (12b) (for c = c and c = c). The corresponding variations for (u 2 u 2 ) Rs /S 2 L and (u 2 u 2 ) R /S 2 L = [(u 2 ) R (u 2 ) R +(u 2 u 2 ) R ]/S 2 L with c are shown in Figure 6(b). It is evident from both Figures 6(a) and 6(b) that (u 1 u 1 ) Rs and (u 2 u 2 ) Rs are appropriately captured by (u 1 u 1 ) R and (u 2 u 2 ) R , respectively. A comparison between Figure 2(a) and Figure 3(a) shows that the magnitude of (u 1 ) R (u 1 ) R is much larger than the magnitude of (u 1 u 1 ) R . Therefore, the behaviour of (u 1 u 1 ) Rs appears to be well captured by (u 1 ) R (u 1 ) R .
It is evident from Figure 6(a) that (12b) for both c = c and c = c overpredicts the magnitude of (u 1 u 1 ) Rs for the major part of the flame brush. According to (12b) the quantity (u 1 u 1 ) Rs is bound between (u 1 u 1 ) P at c = 0 (or c = 0) to (u 1 u 1 ) R at c = 1 (or c = 1), and, as (u 1 u 1 ) P is greater than (u 1 u 1 ) R in this case because of greater values of (u 1 u 1 ) P and (u 1 ) P than (u 1 u 1 ) R and (u 1 ) R , respectively (see Figure 2 and compare Figures 3 and 4), the model given by (12b) overpredicts the magnitude of (u 1 u 1 ) Rs for the major portion of the flame brush. Moreover, in this case the value of (u 2 u 2 ) P remains greater than (u 2 u 2 ) R , which leads to an overprediction for the model given by (12b) whereas (12a) satisfactorily predicts (u 2 u 2 ) Rs throughout the flame brush.
The variations of (u 1 u 1 ) Ps /S 2 L and (u 2 u 2 ) Ps /S 2 L with c are shown in Figures 7(a) and 7(b), respectively, along with the variations of (u 1 u 1 ) P /S 2 L and (u 2 u 2 ) P /S 2 L . In the present case (u 1 u 1 ) P overpredicts (u 1 u 1 ) Ps whereas (u 2 u 2 ) P underpredicts (u 2 u 2 ) Ps . The former trend is consistent with the trend demonstrated in Figure 5(c), that is, (u 1 ) P > (u 1 ) Ps , and is explained when discussing Figure 5(c).
It is worth noting that the above analysis is carried out using a simplified chemistry-based DNS database representing the corrugated flamelets regime combustion for a moderate value of turbulent Reynolds number following several previous studies [16,24,25]. More analysis is needed based on detailed chemistry-based DNS and experimental data for higher values of turbulent Reynolds number Re t to develop more efficient models for (u i u j ) R , (u i u j ) P , and (u i u j ) Ps .

Conclusions
The statistics of mean fluid velocities conditional in reactants and products in the context of Reynolds Averaged Navier Stokes simulations are analysed in detail using a DNS database of statistically planar turbulent premixed flame representing the corrugated flamelets regime of premixed turbulent combustion. It has been found that the conditional mean velocities (u i ) R and (u i ) P and the velocity fluctuation correlations (u i u j ) R and (u i u j ) P differ substantially from the corresponding Favre averaged quantities.
While the conditional velocities (u i ) R and (u i ) P evaluated from DNS data agree well with the predictions obtained from the BML model, the conditional Reynolds stresses (u i u j ) R and (u i u j ) P are not appropriately predicted by the expressions derived based on a presumed bi-modal distribution although the segregation factor g = c 2 / c(1 − c) remains close to unity. It is also demonstrated that the conditional surface averaged velocities (u i ) Rs and the velocity correlations (u i u j ) Rs can be effectively modelled by (u i ) R and (u i u j ) R , respectively, for the major part of the flame brush with the exception of the leading edge. However, (u i ) Ps and (u i u j ) Ps cannot be modelled by (u i ) P and (u i u j ) P , respectively.
As the present analysis has been carried out for moderate value of turbulent Reynolds number Re t and simplified chemistry, detailed chemistry-based DNS and experimental data for higher values of Re t will be necessary for more comprehensive physical understanding and the development of more efficient models. Moreover, it remains to be seen how these models perform for nonunity Lewis number and low Damköhler number conditions where the pdf of c is unlikely to be bimodal in nature. Some of these issues will form the basis of future investigations.
The influence of reactive scalar mixing physics on turbulent premixed flame propagation is studied, within the framework of turbulent flame speed modelling, by comparing predictive ability of two algebraic flame speed models: one that includes all relevant physics and the other ignoring dilatation effects on reactive scalar mixing. This study is an extension of a previous work analysing and validating the former model. The latter is obtained by neglecting modelling terms that include dilatation effects: a direct effect because of density change across the flame front and an indirect effect due to dilatation on turbulence-scalar interaction. An analysis of the limiting behaviour shows that neglecting the indirect effect alters the flame speed scaling considerably when u /s o L is small and the scaling remains unaffected when u /s o L is large. This is evident from comparisons of the two models with experimental data which show that the quantitative difference between the two models is as high as 66% at u /s o L = 0.

Introduction
The propagation of a deflagration wave in turbulent medium is a classical problem of turbulent premixed combustion and is of great practical significance in the current climate of energy and environment. Expressions for the mean propagation velocity, referred to as "turbulent flame speed" or "turbulent burning velocity", are often required in practical computational fluid dynamics (CFDs) codes employed in the design of combustion systems. Furthermore, from a theoretical standpoint, turbulent flame speed is a useful analytical tool to assess the general validity of turbulent combustion models [1] for premixed flames. Analytical expressions of turbulent flame speed have been studied in many previous works [2][3][4][5][6][7][8]. It is well known from the classical theories [9] that in regimes of practical interest, the turbulent burning rate, and hence turbulent flame speed, is dictated by the rate of turbulent mixing at scales relevant to sustain combustion on the flame surface. A direct relationship between the average burning rate and the small scale mixing rate, known as scalar dissipation rate c , in large Damköhler number flames can be written as [9] where C m is a model constant of order unity and ρ is the average density. The Damköhler number is the ratio of the integral time scale of turbulence to the chemical time scale, Da = (Λ/δ o L )/(u /s o L ), where s o L and δ o L are, respectively, the flame speed and thermal thickness of an unstrained planar laminar flame, and the root-mean-square value of turbulent velocity fluctuations with integral length scale Λ is u . The subscript c denotes a reaction progress variable, which can be defined using scalar mass fractions or temperature or sensible enthalpy [10]. Here, it is defined using temperature, and it varies from zero in the unburnt mixture to one in the fully burnt mixture. The scalar dissipation rate is defined as c = ρD(∇c · ∇c )/ρ, where D is the diffusivity of progress variable c and c is the Favre fluctuation of c. The scalar dissipation rate, which denotes the fine scale mixing rate of hot and cold fluid parcels, embodies the scalar field dynamics, and in turn, influences the rate of chemical reactions, since the mixing of the fluid parcels is what ensures sustained combustion on the flame surface. Equation (1) suggests that the modelling of mean scalar dissipation rate yields a model for the mean burning rate. This approach has been used in [11,12] in conjunction with the Kolmogorov-Petrovskii-Piskunov (KPP) analysis to obtain an expression for the turbulent flame speed.
The KPP theorem [13], which is based on eigenvalue analysis of the Favre averaged transport equation for the progress variable, gives the propagation speed of the leading edge of the flame brush as where ν t is the turbulent kinematic viscosity and Sc c is the turbulent Schmidt number for the scalar c. Strictly, this analysis applies for statistically planar flames when the density and turbulence diffusivity are taken to be constant. Furthermore, it requires gradient flux approximation. However, Lipatnikov and Chomiak [1] have shown that the KPP analysis is applicable equally when the density and the diffusivity are not constant. The analysis [11] of statistically nonplanar flames using Taylor's series expansion of Hakberg and Gosman [2] results in an expression similar to (2).
In the presence of counter gradient flux, which is known to occur inside the flame brush when the thermochemical effects are stronger than the turbulence effects, Corvellec et al. [14] have noted that the solution to the KPP analysis is limited by the condition at the burnt side ( c → 1) rather than by the condition at the unburnt side ( c → 0). However, it is well known that the turbulent scalar flux is gradient at the leading edge of the flame brush where (2) applies. Furthermore, Bray [8] has pointed out that other premixed flame theories including nongradient transport of turbulent scalar flux give flame speed expressions which are equivalent to the KPP expression in (2). The KPP analysis has been used in many past turbulent flame speed studies [1-4, 8, 14]. Substituting (1) into (2), one obtains Now, the modelling of turbulent flame speed depends on the scalar dissipation rate modelling. Also, the behaviour of scalar dissipation rate at the leading edge of the flame brush controls the propagation speed, S T as per (3).
Early theories for turbulent flows have taken the scalar dissipation rate to be purely a function of the turbulence parameters. While this is true for turbulent flows with passive scalars, recent studies [15][16][17][18][19][20] have shown that the scalar dissipation rate is strongly influenced by chemical reactions in premixed flames. The heat release directly affects the local density, which induces alterations in the local dynamics of turbulence, scalar fields, and their interaction and thus, one can envisage a two-way coupling between heat release and the turbulence. This two-way coupling has been investigated in recent studies, and both of these are observed to be leading order [15] effects. A general review of thermal expansion effects is provided by Lipatnikov and Chomiak [21]. Also, algebraic models for the mean scalar dissipation rate and turbulent flame speed that appropriately include these couplings have been proposed recently [12] using the models developed [22] for various terms in the scalar dissipation rate transport equation. These models were validated using turbulent flame speed data covering a wide range of flame configurations and conditions [11].
The aim of the present work is to illustrate the influence of the two-way coupling of heat release effects and the associated physics on turbulent premixed flame propagation. This is done, within the turbulent flame speed framework, by comparing the predictive capability of the flame speed model with and without the terms signifying the physics behind the two-way coupling. In principle, one could also do this within the scalar dissipation rate framework, since the two models have the same terms with the same model parameters [12] as one will find in Section 3. However, the dissipation rate is a very difficult quantity to measure, and experimental data of this quantity are scarce. On the other hand, turbulent flame speed has been widely studied, and there is a wealth of experimental data available for this quantity.
The outline of this paper is as follows. The two-way coupling and the physics associated with the turbulent mixing of a reactive scalar in premixed flames are discussed in Section 2. A brief background on the turbulent flame speed model and the physical significance of various terms is presented in Section 3. The limiting behaviour of the flame speed model which excludes one of the physical effects of dilatation is discussed in Section 4. The analysis of the flame speed expression including the effects of twoway coupling has been reported earlier [11], and thus, it not repeated here. However, relevant results are quoted for convenience. The revised and original flame speed expressions are compared with experimental data in Section 5, and the conclusions of this study are summarised in Section 6.

Physics of Reactive Scalar Mixing
The physics behind the two-way coupling is explained best using the scalar dissipation rate transport equation, which also helps us to identify the mathematical terms signifying these physical processes. This transport equation has been derived in earlier studies for unity [15] and nonunity Lewis Journal of Combustion 3 number [23] flames. This equation can be written [15,16,22] as for unity Lewis number, where ρ is the density, u is the velocity vector, andω c is the reaction rate of c. The left hand side of the above equation represents the unsteady, convective and diffusive flux of the scalar dissipation rate in a control volume. On the right hand side, the first term denotes the molecular dissipation effects, and the direct influence of density change across the flame front is denoted by T 2 . These two are leading order terms [15]. The indirect influence of density change through the interaction of turbulence and scalar fields is denoted by the three T 3 terms. Out of these three terms, T 32 is identified [15] to be a leading order term in premixed flames. Using Eigen decomposition, this term can be written as T 32 = 2ρ c (e α cos 2 θ 1 + e β cos 2 θ 2 + e γ cos 2 θ 3 ), where e α > e β > e γ are the principal components of the turbulence strain tensor ∇u and θ i is their orientation angle with ∇c . The contribution of chemical reactions is denoted by T 4 , which is also a leading order term [15], and it will prevail even if the reaction is passive (without heat release). This transport equation has been studied in the past [15,17] to develop physical understanding and to develop models [22,24] for its various terms. Details can be found in these studies. For a passive chemical reaction, T 2 is zero. Furthermore, in turbulent flows T 32 is dictated by the alignment of scalar gradient with the principal turbulent strain rates as noted above. The alignment characteristics are opposite when the chemical reaction is passive or active; the scalar gradient aligns preferentially with the most compressive strain when the chemical reactions are passive, and the scalar gradient aligns preferentially with the most extensive strain when there is strong heat release from the reactions. These situations are schematically shown in Figure 1, and this change is because of the fact that the dilatation due to heat release is strong compared to the turbulent strain rate. Evidence for such a behaviour has been found in direct numerical simulation (DNS) studies of statistically planar [16][17][18][19][20]22], spherically symmetric [25], and Bunsen [26] flames and also experimental bluff body stabilised flames [27]. As a consequence, the isoscalar surfaces are brought together by the compressive strain resulting in an increase of scalar gradient when the reaction is passive. Because of the change in the alignment, the opposite is true for premixed flames with high heat release. In the present study, we are interested in elucidating the effect of this change in the physics on the propagation speed of the flame brush leading edge. As noted earlier, modelling of these physical process has been done in previous studies [22], and these models have been used [12] to obtain a model for the turbulent flame speed. We start from this model for this study.

Turbulent Flame Speed Model
A model for the scalar dissipation rate, c , accounting for the direct and indirect effects of heat release rate is written as [12] Using this model in (3), the algebraic model for the flame speed, S T , presented in [12] is obtained, and it is written as follows: for turbulent premixed flames with large Reynolds, Re, and Damköhler numbers. The constant C μ in (6) is the standard k − ε turbulence modelling constant with a value of 0.09 and the constant C m = 0.7 for hydrocarbon-air flames [9]. The physical meaning of other constants is as follows. The symbol β = 6.7 represents the contributions from processes related to flame front curvature, and it comes from the modelling of (T 4 − D 2 ) in (4). The direct contribution of dilatation, T 2 , at leading order is signified through the term containing K * c . The numerical value of K * c is 0.85τ for hydrocarbon-air mixtures and 0.65τ for hydrogen-air mixtures, where τ ≡ (T ad − T u )/T u is the heat release parameter, where T u and T ad are the unburnt mixture and adiabatic flame temperatures, respectively. The leading order effects of indirect influence of dilatation through the scalar gradient alignments, T 32 , is signified by the term with C 4 and C 3 represents the contribution of turbulence. The values of these two constants are, respectively, where Ka is the Karlovitz number, which is defined as the ratio of chemical to Kolmogorov time scales, and can be evaluated using Ka = (D/s 0 L η) 2 , where η is the Kolmogorov length scale. The ratio of the integral turbulent velocity  time scale to scalar time scale is denoted by C 3 , and its dependence on Ka is introduced to capture the variation in the value of the time scale ratio [12]. The Ka dependence for C 4 is introduced in [24] to recover the classical alignment behaviour of the scalar gradient with the principal turbulent strain rates when Ka becomes large.
It is evident that the terms pertaining to the two-way coupling are K * c and C 4 , and both are consequences of dilatation because of chemical heat release. It is worth pointing out that these terms are relevant for compressible fluids and may not be important for the case of, say, chemical reactions occurring in turbulent liquids. The limiting behaviour of (6) without these terms will be studied in the next section.

Influence of Dilatation Effects on Limiting Behaviour
The case of passive chemical reaction is trivial for this study and has been discussed already in [11]. Since there is no heat release, both K * c and C 4 are zero for this scenario. Hence, , as has been noted in earlier studies [11,28]. This classical expression, without Ka dependence, has been suggested by many turbulent premixed flame theories and is not particularly interesting for the present study.
If one were to retain only the direct part of the two-way coupling, then K * c is nonzero and C 4 = 0. This physically implies that the leading order contribution of dilatation due to density change across the flame front is important, but the induced alterations in the dynamics of turbulence and scalar mixing to be ignored. This situation is an improvement from the classical cases noted above, which assume that the alignment characteristics of reactive scalar gradients are the same as those of passive scalars. Now, from (6), one obtains It is worth noting that this feature would appear in a flame speed expression that could be obtained using the mean scalar dissipation rate, c , model proposed in the earlier work of Swaminathan and Bray [15] in (3).
To study the limiting behaviour of (8) in the limits of low and high u /s o L , following [11], we assume Λ/δ o L to be constant and recast (8) (i) In the limit of small u /s o L , Ka is small, and hence, The value of τ for most fuels lies typically in the range 5 to 10, and hence, one would expect A B, which yields the scaling Da, which is different from Da 1/4 discussed in [1]. The scaling in (9) suggests a square root dependence on (u /s o L ) for the normalised turbulent flame speed, S T /s o L , and also the contribution of thermochemistry overwhelms the contribution of turbulence for small u /s o L values. On the other hand, the opposite is true for (6) which yields a linear scaling for S T /s o L with u /s o L . This is because the contribution of direct effect of dilatation, signified by K * c term, is balanced by the contribution of C 4 making A to be of order zero for (6) as has been noted in [11]. The linear scaling is consistent with the scaling arguments proposed by Damköhler [29].
for both (6) and (8), since C 4 0 for large Ka. The above scaling can also be written as S T It is interesting to note that neglecting C 4 -related term affects the scaling of turbulent flame speed in the limit of small u /s o L but not in the limit of large u /s o L . This is not surprising, since C 3 and C 4 represent the competing effects of turbulent straining and thermochemistry respectively on the interaction between turbulence and scalar fields [16,22]. When u /s o L is very large, the turbulence is expected to overwhelm the thermochemical process, and the model parameters reflect this behaviour. Hence, neglecting C 4 makes no difference in the large u /s o L limit. (It is also important to note that these scaling does not say how large is large.) This does not imply that the direct effect of dilatation through K * c is also negligible, and its influence appears as the first term in (10). Also, the above two scaling expressions for the turbulent flame speed suggests a qualitative difference when compared to the corresponding classical scaling relations. This is somewhat different from the observation in [30], suggesting only a quantitative influence of density change for the propagation speed of premixed flames in large-scale and low-intensity turbulence.
Finally, it is worth commenting on why we have not considered the case where K * c is neglected but C 4 is retained. Such a model is not physically correct, since it ignores the leading order contribution of dilatation, and it would not satisfy the realisability condition for the mean scalar 6 Journal of Combustion . . dissipation rate discussed in [12]. Physically, setting K * c = 0 implies that there is no density change across the flame front, and thus, it is meaningless to consider the indirect influence of density change through the dynamics of turbulencescalar interaction. Furthermore, for low values of Ka, such a model might yield negative value of mean scalar dissipation rate which is unphysical. In the next section, the predictive capability of (8) is compared to that of (6) by making comparisons with experimental data.

Comparison to Measurements
The original flame speed model, (6), was comprehensively validated in [11] with experimental data from a wide range of turbulent flame configurations and conditions spanning all regimes of practical interest. Particularly, due rigour was paid to the definition of turbulent flame speed given by the KPP analysis-the propagation speed of flame brush leading edge-and the experimental data were chosen accordingly. Other definitions of turbulent flame speed are possible and the reader is referred to the review paper by Driscoll [32] or [11] for a discussion on this topic. For comparisons in the present study, we choose the following subset of the experimental data considered in [11]: (1) the planar flame data from the Taylor-Couette apparatus of Aldredge et al. [31], (2) the high pressure Bunsen flame data of Kobayashi et al. [33], (3) the very high turbulence intensity data of Il'yashenko and Talantov [34,35]. This is to emphasis the relative roles of direct and indirect effects of density change on the propagation speed of the flame brush leading edge for wide conditions of turbulent flames. As noted earlier, ignoring the effects of density change completely will give the classical result analysed in many earlier studies. Comparisons of (6) and (8) with the above three data sets are shown in Figures 2, 3 and 4, respectively. The values of s o L , δ o L , and τ required for the comparisons are calculated for each case in the same fashion as described in [11]. In general, the predictions of (8) are worse compared to the predictions of (6), and this is expected, since (8) does not include all the effects of density change across the flame front. The quantitative difference in predictions is more pronounced in Figure 2 compared to the other cases. This is because all the data in Figure 2 correspond to low u /s o L values (≤4). However, the trend in Figure 3 clearly illustrates the difference between the two expressions. The quantitative difference between (6) and (8) varies from 66% at u /s o L = 0.3 in Figure 3(a) to 12% at u /s o L = 10 in Figure 3(d). This is further evident in Figure 4, where the quantitative difference decreases from 62% to 4% as u /s o L varies from 1.1 to 52.4. These trends are consistent with the analysis in the previous section, which suggested that neglecting the C 4 term in (6) affects the scaling appreciably in the low u /s o L limit, while it is negligible in the large u /s o L limit. This observation underscores the effects of density change on the flame speed and the importance of validating turbulent flame speed models with data over a wide range of conditions. If one were to consider experimental data selectively, for instance, data at only large u /s o L in the present case, then the comparisons can be misleading. Furthermore, the dominant effect of gas expansion or density change across the flame front in weak turbulence limit is consistent with the observation by Peters et al. [36] for the corrugated flamelets regime combustion.
Nevertheless, the direct effect of density change prevail for large u /s o L values relevant for thin reaction zones regime of combustion with Da > 1. This is becomes clear when the two dashed lines in Figure 4 are compared. The short dashed line in this figure is for (8) with K * c = 0, which implies that both the direct and indirect effects of heat release are ignored. However, the indirect effect of dilatation becomes negligible only when u /s o L is larger than 50 and Da is of order unity. This is also clear from Figure 4 (compare the solid and long dashed lines).

Summary and Conclusions
The present study illustrates the influence of dilatation effects on the propagation speed of turbulent premixed flames, by assessing their contribution to turbulent flame speed calculation. An algebraic model for turbulent flame speed, (6), proposed and validated in earlier studies [11,12], is a suitable benchmark, since it incorporates the two-way coupling of the effects of density change across the flame front. The direct effect of density change on the local average mixing rate of hot products and cold reactants is through the term involving K * c in (6), and the indirect influence is through the alignment of reactive scalar gradient with the principal components of turbulent strain tensor. The term involving C 4 in (6) signify the indirect influence.
If both of these terms are ignored, then the original model proposed in [11,12] yields a classical form S T /s o , when the molecular viscosity is included in the analysis, that has been proposed by many early theories. Neglecting only the K * c term but retaining the C 4 term is unphysical, since it implies the influence of indirect effect when the density change across the flame front is not allowed. Furthermore, this situation leads to negative values for the scalar dissipation rate, c , violating its realisability condition [12]. On the other hand retaining K * c but neglecting C 4 yields a flame speed model as in (8), which has a different limiting behaviour compared to the model in (6). The scaling of (8) is different from (6) in the small u /s o L limit but not in the large u /s o L limit. Such a behaviour is also evident in comparisons with experimental data.
The turbulent flame speed values obtained using these two equations differ largely in the limit of weak turbulence (about 66% for u /s o L ≈ 0.3) or in the corrugated flamelets regime. This is consistent with the observations of Peters et al. [36]. However, the direct effect of dilatation, signified through K * c , prevails for large u /s o L values, but the indirect influence of dilatation thorough the dynamics of turbulencescalar interaction weakens as u /s o L increases. These two effects of heat release is observed to prevail for practically relevant values of u /s o L (∼20).

Introduction
The determination of the burning velocity of turbulent premixed flames S T has been the subject of many experimental studies. As far as a given modelling of such reacting flows is concerned, the a priori value of S T for statistically stationary one-dimensional turbulent flame is often estimated through the direct application of the results of the Kolmogorov-Petrovskii-Piskounov (KPP) theory [1] which showed that, for a one-dimensional reaction zone without heat release described through the evolution of a single progress variable c and with a constant diffusion coefficient D t , the steady propagation of the reaction zone is controlled by the behavior of the reaction rate at the reactants side and by the value of the diffusion coefficient. In many situations, experiments (Moss [2], Cheng and Shepherd [3], Troiani et al. [4], Zimmer et al. [5]) as well as direct numerical simulations (DNS) (Veynante et al. [6], Nishiki [7], Hauguel [8], Lee and Huh [9]) show that an expression of the turbulent transports via a gradient expression is not appropriate as some mechanisms may promote a transport with a sign identical to that of the mean scalar gradient. An extensive review of relevant experimental and DNS data can be found in the paper of Lipatnikov and Chomiak [10]. In such situations, the introduction of the Bray number proposed by Veynante et al. [6] proved to be helpful to discriminate between flow configurations featuring gradient turbulent flux (GTF) only (N B ≤ 1) or counter-gradient turbulent flux (CGTF) (N B ≥ 1).The Bray number was defined as N B = τS L /2α v u where τ designates the heat release parameter and α v some order unity efficiency function which depends mainly on (and is an increasing function of) the ratio between the turbulence integral length scale and the flame thermal thickness. For N B ≥ 1, one of the most prominent mechanism leading to CGTF is directly related to the differential effect of the pressure gradient on the pockets of reactants (heavy) and products (lighter). But it should, nevertheless, be stressed, that even when such a mechanism is predominant, there still exists a gradient turbulent flux acting at scales much smaller than those of the pockets and for which the turbulent transport coefficient cannot be estimated through a turbulence model which integrates the whole turbulence scales. Accordingly, and following in that respect Veynante et al. [6], Zimont and Biagioli [11] or Lipatnikov and Chomiak [12], this suggests to consider the turbulent scalar flux as resulting from the competition between gradient and counter-gradient mechanisms. Along these lines, we propose to study the propagation properties of a turbulent premixed flame by considering such a flux decomposition. First, an extended KPP analysis is developed. It combines all the possible scenarios with the change of N B and provides a general analysis whose results are not qualitatively dependent anymore on the modelling closures chosen but permit a clear and comprehensive discrimination between all the various propagation regimes.
We note that, in the framework of an eddy-breakup model, an analytical attempt aimed at explaining the influence of the presence of CGTF on the turbulent flame propagation properties has been done by Corvellec et al. [13] who considered an idealized one-dimensional premixed flame that propagates through high-Reynoldsnumber frozen turbulence. These authors dealt only with a flame brush formed by a GTF zone followed by a CGTF one. The scalar flux in each zone was expressed by employing a classical gradient formulation with the peculiarity of using a strictly negative diffusion coefficient for the CGTF zone. Thanks to this mathematical "trick", it was then possible to employ the KPP technique to study the characteristics of the governing equation at both sides of the flame brush and at the intermediate point at which the flux is zero. This approach did not put into evidence a qualitative change in the nature of the burning velocity domain which appears to be still a semi-infinite interval. The lower bound though appeared to be controlled by both sides in such a situation. It drew also the attention on the fact that the KPP results cannot be extrapolated directly in such situations (in accordance to some extent with the numerical simulations of Bradley et al. [14]). The results of Corvellec et al. [13] cannot be considered as being of general implication since they are contingent to the recourse to an effective diffusion coefficient and are limited to the analysis of only one flame structure, that is, a GTF zone followed by a CGTF one. The Bray number is absent in the theory developed in [13] whereas DNS [6,9] and recent experiments on turbulent premixed flames [4,5] show that a different mean flame structure may be encountered that is, a CGTF zone followed by a GTF one. Thus, the present approach is more general than those followed previously since it permits to cover all the two-zone mean flame structures that can be thought of. It shows in particular that as far as the burning velocity domain is concerned, qualitative differences with the KPP and Corvellec et al. [13] results are obtained.

Governing
Equations. An unsteady 1-D freely developing turbulent isenthalpic premixed flame is considered here and is described through the evolution of a single progress variable c = (T − T r )/(T b − T r ), where subscripts r and b denote the states of the reactant mixture and the fully burnt products, respectively, and τ = (T b − T r )/T r defines the heat release parameter. Thus, in a fixed coordinate system (O, x) the Favre averaged equation for c and the continuity equation to be used later on, are given by ρu c is the turbulent flux, and w is the mean reaction rate whose exact expression is not required at this stage. The flame is propagating from right to left, that is, ∂ c/∂x 0. We consider here the case for which counter-gradient mechanism is largely present in the flame as it is the case when the turbulence intensity is sufficiently low and/or the heat release parameter sufficiently large (see [6] for further details). Accordingly, the turbulent flux is expressed as F GTF < 0 represents the gradient contribution while F CGTF > 0 corresponds to the counter-gradient part. Such a decomposition is compatible with the findings of Veynante et al. [6] who derived an algebraic expression for the flux as a function of the laminar flame velocity S L , the heat release parameter τ and the rms velocity u , namely, ρu c = ρ c(1− c)(τS L − 2α v u ). In this particular case, the counter-gradient contribution is equal to and the gradient one reads as The ratio F CGTF /|F GTF | is just equivalent to the Bray number N B . If we introduce D t ≈ l t u , where l t is the turbulence integral length scale and suppose that ∂ c/∂x ≈ c(1 − c)/l t , then the gradient contribution can be classically written as F GTF = −ρD t ∂ c/∂x. Here, the dependency of α v on the ratio between the flame thermal thickness and the turbulence integral length scale has been neglected. In the following analysis, we shall consider F GTF as given by (4) combined with a more general expression for F CGTF , namely, where f is a positive continuous function of c such that f (0) = f (1) = 0. In the case of Veynante et al. [6], f ( c) is a quadratic symmetric function, namely, If we now rewrite the c-equation combined with the continuity equation, one obtains It is worth noting that the analysis could be done also in the case, not considered here, where F would be expressed without the introduction of a turbulent diffusion coefficient (e.g., using directly the decomposition of Veynante et al. [6]).
We rewrite (7) in a nondimensional form by defining the following reference parameters: (1) a velocity scale u = u , (4) a reference turbulent diffusion coefficient D t0 = l t u (l t is chosen in such a way that the proportionality coefficient between D t0 and l t u is unity), Thus, one has In the following and whenever unambiguous, we shall drop superscript * . So, with that convention, (8) reads as and the nondimensional form of the continuity equation (1) is given by At this stage, it is important to underline that the form of (7) allows us to give a physical interpretation of the influence of CGTF on the flame structure (following mainly the reasoning adopted in [6]). The term of CGTF is a nonlinear convection term which tends to move both sides toward each other because of the change of sign of ∂ρ f ( c)/∂ c (positive at the fresh reactants side and then negative at the burnt products side). Such a "thinning" effect has to be counter-balanced by the gradient transport term in order to obtain a steady state flame propagation regime.

Steady-State Regime of Flame
Propagation. Let us consider a steady state regime of flame propagation. Accordingly, in a coordinate system attached to the mean flame brush and introducing the flame mass consumptionṁ and the turbulent flame speed S t such asṁ = ρ u S t , the preceding continuity and progress-variable nondimensional equations are written as where Λ = S t /u . We can rewrite (11) as a first-order differential equation Introducing N B = N B 2α v and Ω = wRD, (12) is expressed as where () denotes the derivation with respect to c.

Analysis of the P -Equation in the Vicinity of the Singular Points
Now, our primary objective is to determine the domain of existence of Λ. In principle, we need for that to analyse the trajectories in the (P, c) phase plane in the KPP-like manner. It is clear that a priori, due to the presence of the CGTF derivative f , it is difficult to perform the KPP analysis in such a case. Thus, leaving aside the question of determining the trajectories themselves, we focus on the determination of the associated Λ domain by examining the characteristic equations of (13) at both sides of the brush. The implication of the relative position of P and F = ( N B /Λ)R f for the P-trajectories when they are tangent to the characteristics lines is also analyzed. We denote (a) Reactants Side. In the limit c → 0 + and assuming that P = s c and F = α c (s > 0 and α > 0), the characteristic equation associated with (13) is given by The required positiveness of the discriminant Δ 0 + = (1 + α) 2 − (4/Λ 2 )Ω (0) is assured as soon as Λ 2 > f 1 (α)Ω (0). The corresponding roots s 1 = (1 + α)/2 + Δ 0 + /4 and s 2 = (1+α)/2− Δ 0 + /4 are both positive with s 1 > s 2 . Accordingly, the reactants side is an improper unstable node and all the trajectories but one (the characteristic line s 1 c itself) are tangent to the s 2 c line. Two cases are to be considered in the vicinity of c = 0: (1) F GTF + F CGTF < 0 (the gradient contribution overcomes the counter-gradient one) and (2) the reverse situation prevails, that is, F GTF + F CGTF > 0.
(i) α > s 1 . Solving this inequality yields α > 1 and In both cases (1, 2), we must satisfy also Thus, for case 1 we have and for case 2 Using the definitions of f 1 (α) = 4/(1 + α) 2 , f 2 (α) = 1/α and For case 2, we have Introducing the variables we can rewrite for case 1 and for case 2 It is important to note that relations in (23) are the sum of two sets, but not the product. Figure 1 illustrates the domain of existence of GTF and CGTF in the plane ( Λ, N B ) that followed from the analysis of (13) (for P) in the vicinity c → 0.
In the limit β → 0, the original KPP behavior is recovered, that is, the absence of constraint. (2) CGTF prevails, that is,

Journal of Combustion
Using the relation we transform the above inequalities into the following form: Table 1 Configuration for case 2 (CGTF) where There exists a maximum of four possible configurations, listed in Table 1, which correspond to the total turbulent flux behavior (gradient (G) or counter-gradient (CG)) in the vicinity of the singular points c = 0 and c = 1.
(iv) CG/CG configuration: The different regions of the ( Λ, N B ) plane corresponding to these four configurations are presented in Figure 4.Ñ We note that the case f B < 1 is realised for the classical eddybreak-up model where case [6]. In that case, one has Ω( c) = c(1 − c)/(1 + τ c) 2 and consequently Ω (0) = 1 and Ω (1) = −1/(1 + τ) 2 . Thus, The case f B > 1 is obtained if where κ is some positive contant. So we have , For Veynante et al. [6], f ( c) = c(1 − c) and so It should be noted that experimental data [3,15] suggest that w * is shifted towards the burnt gas side. Such a situation can be modelled by relation (36) with κ > 0, and, in this case, we can have f B > 1. As was shown previously, in such a situation, the configuration G/CG is absent. Such a result is compatible with DNS results for 1-D freely propagating turbulent premixed flame [6,9]. It is also in correspondance with the recent experimental data of Troiani et al. [4] who investigated bluff-body stabilized turbulent premixed flames. These authors showed that the radial component of the scalar flux (which is the counterpart of F) could exhibit a transition between a counter-gradient to a gradient behavior when moving along the normal to the mean front from reactants toward products.

A Particular Configuration:
The Murray's Equation The Murray's equation [16] is a good example of the convection-diffusion equation, of the type considered here, to demonstrate that the common viewpoint that at the leading edge it is the diffusion term which always drives the flame propagation is not universal. Such an equation is obtained as a particular case of (9) corresponding to R = 1, We note that formally though, R = 1 corresponds to τ = 0 leading to N B = 0. But we will consider here N B like a parameter not linked to τ. In such a situation, by denoting K = 2 N B and c w = Λ + N B , (11) reduced to that analyzed by Murray [16] to determine the wave speed c w , namely, Murray [16] showed that travelling wave solutions exist for all c w ≥ c w (K) where or equivalently with our notations Thus, it appears that in that particular case, for N B 1 and in the zone c → 0, it is the counter-gradient transport that prevails (see the domain CG/G in Figures 3 and 4).

Concluding Remark
A methodology that permits to study the domain of the tubulent flame velocities for a steady regime of a turbulent premixed flame propagation has been presented. It has been shown how the introduction of the parameter f B (defined by 32) that combines both the counter-gradient flux and the mean reaction rate behavior at both edges of the flame can be used to determine the domain where a steady regime of propagation may exist for any combination of scalar fluxes that dominate at the flame edges. It is worth noticing though that the present analysis of the sole velocity domain does not prove the existence/unicity of the steady solutions. Consequently, future work will concentrate on the use of numerical simulations to study these questions.

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 t ), while employing at the same time complex physicochemical models to describe turbulence, molecular transport, and chemistry [2][3][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][7][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 t 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 multicomponent 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).

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, u j the components of the hydrodynamic velocity, p the pressure, τ i j the stress tensor, N s the total number of species, V k j the component of the diffusion velocity of species k, in the direction j,ω k the chemical production rate of species k and q j the jthcomponent of the heat flux vector. The volume viscosity coefficient κ is a function of the fluid local properties and appears explicitly within τ i j 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 m 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 m, the more expensive the underlying algorithm and the more accurate the expression for κ, hence our choice of m = 6 (EGSK6()) for the computations in the present study. The simplest and cheapest of the models (m = 2, EGSK2()) considers only the transport system matrix associated with the internal energy, whereas m = 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 fourthorder 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 t 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 for writing/reading all files by a factor of at least 3) and single (restart/solution) data files.

Flame Configuration and Initialization
A stoichiometric spherical premixed methane-air flame is considered in all computations, within a cubic computational domain of sides L = 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.  [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 Φ ≤ 1.0 for the present study.
The initial system is a stationary hot (T b ≈ 2 200 K) 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 atmospheric premixed mixture of methane and air at an unburnt temperature T u = 300 K. The initial profile of any primitive variable Φ is prescribed 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 of the fresh/burnt gas interface. The initial mass fraction values are Y CH4 = 0.055 and Y O2 = 0.220 at T u outside the kernel, and Y CO2 = 0.120 and Y H2O = 0.137 at T b within. An appropriate nitrogen complement is added everywhere.
To investigate the influence of the turbulent Reynolds number Re t 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 u , the turbulent Reynolds number Re t , the eddy turn-over time τ used to quantify flame/turbulence interaction for time-decaying turbulence [35], and the Karlovitz number Ka ≈ [(u /s L ) 3 (l t /δ th ) −1 ] 1/2 are given in Table 1 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 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 t/τ = 1. The longest computation requires about two months. For the results presented below, a nondimensional time t = 1.2τ, needed to obtain appropriate equilibrium conditions between turbulence and chemistry [35], is systematically retained for the analysis.

Results and Discussion
The obtained DNS datasets are processed using an inhouse 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 (u ) 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.

Laminar Flame Structure.
First, a laminar case (referred hereafter as case 0) was computed in a smaller computational box (L = 3.0 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 (t = 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: Corresponding temporal evolution of the fuel consumption (S C ) and integrated heat release (H r ) 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 t = 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 2 O 2 and HO 2 as well as global flame quantities like S C and H r . 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. 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. 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 (t = 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 t > 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.

Turbulent Flame
The iso-contours of the OH radical closely follow those of the temperature and are shown in Figure 7 at the same nondimensional time t = 1.2τ for cases 1-4 when increasing turbulence intensity (Re t = 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 t 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. 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 t = 1.2τ. As in the laminar comparisons, the differences in both the structure and peak values are negligible.

Temporal Flame Evolution.
Two particularly important global flame quantities-the burning rate S c and the volumeintegrated heat release rate H r (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. S c 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ω c , and when suitably scaled, serves as a measure of the burning velocity of the flame [47,48]: 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. 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, H r is approximately constant up to about t ≈ 0.2τ after which it increases steadily for about 0.1τ. At t ≥ 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 t ). 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 t/τ = 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.
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 κ.

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 |∇C| 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 |∇C| and H 2 O 2 mass fraction show a noticeable dependency on turbulence intensity, with progressively lowered peaks. The decrease of the conditional mean |∇C| reveals that flame thickening is predominant, especially around 0.1 ≤ C ≤ 0.7 (the active flame region) when increasing Re t . 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 |∇C| 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.

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 t = 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).

Introduction
Although turbulent flame speed S t , by analogy with laminar burning, is often considered to be a basic characteristic of premixed combustion and was the main subject of numerous studies reviewed elsewhere [1][2][3][4], a precise definition of this quantity is difficult due to the following two key differences between laminar and turbulent premixed flames.
First, while a typical laminar flame is thin as compared with both the flame curvature and scales of spatial nonuniformities of the incoming flow of fresh reactants, the mean thickness Δ t of a turbulent flame is comparable with length scales that characterize mean flow nonuniformities or the flame curvature. Due to the significant thickness of a turbulent flame brush, both the mean convective flux of the deficient reactant toward the flame and the mean flow velocity vary substantially across the flame brush, thus making a choice of a reference flux or a reference flow velocity (with respect to that flame speed should be defined) difficult [5,6]. For instance, Figure 2 in a paper by Gouldin [6] indicates more than threefold decrease in the mean normal convective flux from the leading to the trailing edges of a V-shaped flame brush. Other experimental data that show the same effect were recently reviewed by Lipatnikov and Chomiak [7] (see Section 5.2 in the cited paper). It is also worth remembering that the difference in the speeds of planar and curved laminar flames scales as a ratio δ L /R f of the flame thickness to the radius of its curvature [8], and this ratio is much less than unity in a typical laminar premixed flame. If a similar scaling is assumed to hold for a premixed turbulent flame, then the difference between the speeds of statistically planar and curved flames could be large [9], because a ratio of Δ t /R f is of unity order in a typical curved premixed turbulent flame.
Second, while a typical laminar flame has a fully developed structure and speed, because a time scale of laminar flame development is very short in comparison with time scales characterizing the flow, a typical turbulent flame is a developing wave, with the growth of Δ t being the most striking manifestation of premixed turbulent flame development, as stressed elsewhere [3]. Due to the growth of a turbulent flame brush thickness, the observed speeds of different isotemperature (or isoconcentration) surfaces differ from one another.
In the statistically planar one-dimensional case, the latter problem may be resolved by considering turbulent burning ]}|, that is, the normalized total mass rate of consumption (or production) of the deficient reactant (or a main product) Y per unit mean flame area. However, in the case of a curved flame brush, the definition of the area is ambiguous due to a large thickness Δ t . For instance, a ratio of the areas of surfaces associated with the leading and trailing edges of a statistically spherical turbulent flame brush may be as large as (R f + Δ t ) 2 /R 2 f ≈ 4 if the mean flame radius and thickness are of the same order.

Journal of Combustion
Several methods have been proposed to resolve the above issues and to determine turbulent flame speed S t and/or burning velocity U t in statistically stationary [5,6] or expanding [10][11][12] flames. In particular, in [10], an expanding statistically spherical premixed turbulent flame with a self-similar mean structure (i.e., the mean density ρ(ξ) depends on a single normalized distance ξ = [r − R f (t)]/Δ t (t), rather than on two independent variables, time t and radial distance r) was theoretically studied and the following two radii: were introduced. Here, c is the well-known combustion progress variable [13] and overbar designates Reynolds average. It was analytically shown that (i) if turbulent burning velocity is defined as follows U t ≡ ∞ 0 Wr 2 dr/(ρ u R 2 f ,u ), then it is not affected straightforwardly by the rate of the growth of the mean flame brush thickness, (ii) if turbulent flame speed is defined as follows S t = (ρ b /ρ u )dR f ,s /dt, then S t = U t , and (iii) in order for the difference in the observed flame speed dR f ,s /dt and a mean velocity U u of the unburned mixture to be equal to S t or U t , the velocity U u should be evaluated by extrapolating the mean flow velocity distribution u(r) in the unburned mixture ahead of the flame brush to the burning velocity surface characterized by r = R f ,u . Here, W is the mean mass rate of product creation, that is, the source term in the well-known balance equation [13] for the Favreaveraged combustion progress variable.
It is of interest to note that, as shown elsewhere [12], the above purely theoretical results are very close to the definition of the burning velocity of expanding statistically spherical premixed turbulent flames, proposed by Bradley et al. [11] based on different reasoning.
Although the theoretical analysis in [10] was solely performed for expanding flames, either statistically planar or statistically spherical ones, it was also hypothesized that the proposed methods for determining U t , S t , and U u may be useful in other cases, for example, in statistically stationary premixed turbulent flames stabilized in divergent flows such as impinging jets. In particular, turbulent burning velocity in such a flame was hypothesized to be approximately equal to a velocity obtained by extrapolating the velocity distribution in the unburned gas ahead of the flame to a burning velocity surface characterized by the following reference value of the Favre-averaged combustion progress variable: where and ξ = (x − x f )/Δ t is normalized distance in the direction normal to the mean flame brush. In a subsequent paper [14], this hypothesis was indirectly supported by simulating a premixed turbulent flame stabilized in an impinging jet and by invoking two different RANS models of premixed turbulent combustion. However, in the cited paper, only numerical results were reported, and the discussed hypothesis has never been tested against experimental data, to the best of the author's knowledge.
The present work is aimed at filling this gap and testing the hypothesis against experimental data obtained by four research groups [15][16][17][18] from seven premixed turbulent flames each stabilized in an impinging jet.
In the next section, the hypothesis is tested. Obtained results are discussed in the third section followed by conclusions.

Test
Because the studied flames (see Figure 1 and Table 1) look statistically planar in the vicinity of the axis of an impinging jet, turbulent burning velocity may be evaluated by integrating the mean mass rate of product creation along the axis, that is, where the burning velocity u t = U t /U, burning rate Ω = dW/ρ u U, and axial distance z = x/d (z = 0 at the wall) are normalized using the distance d between the burner exit and the wall, the mean flow velocity U at the burner exit, and the density ρ u of the unburned gas.
In the present work, the burning rate Ω is evaluated using the following well-known balance equation [19][20][21][22][23][24]: which holds in the vicinity of the jet axis. Note that the radial turbulent scalar flux, that is, the last term on the left hand Journal of Combustion 3 side (LHS) of (6), was taken into account in [23,24], but was neglected in the other cited papers [19][20][21][22]. Here, we set ρ u = 1 for simplicity, that is, symbol ρ designates the density normalized using the density of unburned gas, w = u/U is the normalized Favre-averaged axial velocity, which depend only on z on the axis (∂c/∂r = ∂ρ/∂r = ∂ w/∂r = 0 at r = 0 for symmetry reasons), and g = (d/U)(∂υ/∂r) r=0 is the normalized radial gradient of the radial velocity, calculated at the axis, with g depending solely on x. The following wellknown BML relation [25]: was invoked to close the turbulent scalar flux ρu c in the balance equation for c, with u u and u b designating velocity vectors conditioned to unburned and burned mixture, respectively. In order to evaluate the normalized burning rate by solving (6), the axial distributions of c(z), ρ(z), w(z), w u (z), w b (z), g u (z), and g b (z) should be known. In the present work, the following well-known BML state equation [13]: and the experimental data reported by Cho et al. [15], Cheng and Shepherd [16], Li et al. [17], and Stevens et al. [18] were used for these purposes. Here, σ = ρ −1 b is the density ratio. However, because the straightforward use of the experimental data could result in significant numerical errors when differentiating the data measured in a few points separated by substantial distances, smooth numerical approximations of these experimental data were obtained using the following method [24].

to compute c[c(z)].
Here and in the following, ξ = (z − z f )/δ t , with the values of the normalized flame coordinate z f and thickness δ t = Δ t /d being reported in Table 2. As shown elsewhere [3,27], this approximation works well in various premixed turbulent flames. The readers interested in a theoretical justification of (9) are referred to [28,29]. Second, the profiles of the mean axial velocity and either axial conditioned velocities or axial turbulent scalar flux were computed by solving the following well-known [19][20][21][22][23] continuity: and radial impulse equations supplemented with (i) the following boundary conditions: (ii) the well-known BML relation [13] and (iii) balance equations for g u and the difference in w b and w u , derived by the present author [30] and discussed in detail elsewhere [7,31]. Here, Q = −d 2 /U 2 (∂ 2 p/∂r 2 ) r=0 and g 1 are tuning parameters reported in Table 2.
The readers interested in discussion and validation of this numerical model are referred to a recent paper [24]. For the goal of the present study, the validity of this model is a side issue, but the most important point consists of the fact that curves plotted in Figures 3, 4, and 5 well approximate experimental data shown in symbols. Therefore, these curves could be considered to be smooth approximations of the data, while further details of the approximation techniques are of minor importance as far as the evaluation of Ω by numerically integrating (6) is concerned.
Thus, in the present work, the reference method of evaluating turbulent burning velocity consists of numerically integrating (5) and (6) by invoking the smooth approximations of experimental data shown in lines in Figures 2, 3,  4, and 5. The values of u t obtained using this method are reported in the first row of Table 3.
These reference values are used in order to test the method put forward in [10]. The latter method consists of the following.
First, the normalized velocity profile w(z) obtained in the unburned mixture is extrapolated by invoking the following analytical solution: to (10) and (11) written in the constant density case (ρ = 1) and supplemented with boundary conditions given by (12).
In the present work, the parameters Q and g 1 were taken from Table 2. If the method is applied to process raw experimental data, then these parameters may be tuned by comparing the profile of w(z) obtained from the oncoming flow of unburned mixture (e.g., see thin solid line in Figure 6) with profiles calculated using (14) with various Q and g 1 (e.g., see dashed line in Figure 6). Second, the reference value c 0 associated with the burning velocity surface is calculated using (3), (4), and (9). It can easily be shown that the reference surface is located at z = z f if (8) holds. Indeed, using (8) and (9), we have  [18], set 3 a d is the distance between the jet exit and the wall (see Figure 1). b U is the mean axial flow velocity in the jet exit. c Φ is the equivalence ratio. d σ is the density ratio. e u is the rms turbulent velocity. f For CH 4 -air mixtures, the values of the laminar flame speeds S L , reported in Table 1, were estimated using recent experimental data [26]. g This value was reported by Bray et al. [19] (see footnote on page 645 in the cited paper). h This value was reported by Cheng and Shepherd [16]. i The main difference between the two flames studied by Li et al. [17] was the diameters (h = 4 and 6 mm) of holes in grids used to generate turbulence. Subsequently, (4) and (15) read Therefore, and z = z f . In the present work, the normalized distance z = z f was taken from Table 2.
If the method is applied to process raw experimental data, then, the coordinate z 0 of the burning velocity surface may be calculated from the following equality: using the measured profile of the Reynolds-averaged combustion progress variable. If the Favre-averaged combustion progress variable is measured, then the proper reference value may be calculated substituting c = 0.5 into (8). Under typical conditions, the reference value of c ≈ 0.1 could be used (e.g., see bold solid line in Figure 6). Third, turbulent burning velocity is considered to be equal to w e (z f ) calculated from (14); see filled circle in Figure 6.
Results of the test, shown in Table 3, indicate that the values of w e (z f ) (the second row) are very close to the normalized turbulent burning velocity u t calculated using (5) in all the studied flames, thus validating the considered method.

Discussion
From the physical standpoint, the above results could be interpreted as follows. Although the turbulent burning velocity U t is controlled by processes within (or, maybe, at the leading edge of) turbulent flame brush, U t can be evaluated by properly analyzing the velocity field in front of the flame brush, because the changes in this field, induced by the flame, are controlled by the burning velocity and the density ratio. For instance, in the simplest hypothetical case of an expanding spherical infinitely thin laminar flame, the velocity in unburned mixture is equal to (σ − 1)S L (r/R f ) 2 , and the burning velocity can easily be evaluated by analyzing the radial velocity profile measured ahead of the flame.
The approximate equality of u t and |w e (z f )| may further be supported by invoking the so-called Zimont [32][33][34]  Stevens et al. [18], set 1 Stevens et al. [18], set 2 Stevens et al. [18], set 3 (c) Figure 2: Axial profiles of (a) the Reynolds-averaged combustion progress variable c(x), measured by Cho et al. [15] and by Cheng and Shepherd [16], and the Favre-averaged combustion progress variable c(z), measured (b) by Li et al. [17] and (c) by Stevens et al. [18]. Symbols show experimental data. Curves have been obtained using (9). and validated by the present author [35][36][37]. As reviewed elsewhere [3,27], the model is successfully used by a number of different groups to simulate various turbulent flames.
The model deals with the following balance equation: for the Favre-averaged combustion progress variable, where D t is turbulent diffusivity. Note that a basically similar balance equation was first introduced into the combustion literature by Prudnikov [38]. Along the axis of a statistically stationary flame stabilized in an impinging jet, this equation reads where d t = D t /(Ud) is the normalized turbulent diffusivity. Substitution of (8) and (9) into (20) yields Therefore, within the flame,    Differentiating this equation and using (8), (10), and (15), we obtain Let us consider the leading edge of a turbulent flame brush and match (22) and (23), which are valid within the flame, with (14), which is valid in the unburned mixture.
Because ρ → 1, c → 0, and ξ 1 at the leading edge z → z * , we obtain Simple algebraic manipulations with these two equations and (14) result in According to the present simulations (see the third row in Table 3 and curves 1 in Figure 7), the last term in (25) is small at the leading edges of all the seven flames addressed in the present paper. Therefore, the model predicts that u t ≈ −w e (z f ) in a premixed turbulent flame stabilized in a circular impinging jet. Note that a similar analysis of a flame stabilized in a planar impinging jet yields exactly u t = −w u (z f ). Certainly, the above substantiation of the studied method (i.e., u t = −w e (z f )) by invoking the extended Zimont model may be put into question by referring to a paper by Bray et al. [21] who claimed that the model failed in predicting the burning rate within flames stabilized in impinging jets, particularly, in flames nos. 2-5 studied in the present paper. More specifically, Bray et al. [21] have found that the following simple expression: yields too strong dependence of W on c if the turbulent burning velocity U t is evaluated using local Favre-averaged  (25) (1) and mean axial velocities calculated using (2) equation (14) and (3) equation (22) versus the Reynolds-averaged combustion progress variable. Results shown in bold and thin lines have been obtained for flames h4 and h6, respectively, investigated by Li et al. [17].
turbulence characteristics (turbulent kinetic energy k and its dissipation rate ε), which depend substantially on c.
As pointed out elsewhere [39], the model tested by Bray et al. [21] differed substantially from the model developed and validated by Zimont and Lipatnikov [35][36][37]. In particular, the latter model provides a joint closure of the sum of the transport and reaction terms in the following well-known [13] balance equation: that is, within the framework of the extended Zimont model. If applied to a flame stabilized in an impinging jet, (28) reads Different terms in (29) are shown in Figure 8 for the flame h4 investigated by Li et al. [17], and similar results were obtained for the other six flames. In these calculations, the burning velocity used by the extended Zimont model was evaluated by integrating (29), that is,  (29), normalized using d and U, versus the Favre-averaged combustion progress variable. Curves have been computed in the flame h4 investigated by Li et al. [17]. Term numbers are specified in legends. Curve 6 shows the difference between the LHS and RHS of (29). and the normalized turbulent diffusivity was set in order to satisfy (29) at a point characterized by the maximum magnitude of the diffusion term 4. Thus, neither u t nor d t depend on z. As discussed earlier [39], the use of a constant turbulent burning velocity within such flames appears to be more consistent with the underlying physics than substitution of k( c) into an expression for U t = U t (S 2 L / k, . . .), because the Favre-averaged k strongly overestimates the turbulent kinetic energy due to the substantial contribution of the unburnedburned intermittency into k, for example, see Figure 1 in [39]. Figure 8 shows that although the source terms on the LHS and RHS of (29) depend differently on the mean combustion progress variable (compare with curves 3 and 5), (29) is sufficiently accurate in all the studied flames (see curve 6 which shows the difference between the LHS and RHS).
Moreover, (23), which results from the extended Zimont model, yields a constant radial gradient g of the Favreaveraged radial velocity across a premixed turbulent flame stabilized in an impinging jet provided that u t = const and d t = const. Numerical results plotted in Figure 9 show that g is approximately constant in the range c < 0.8 of turbulent flame brush, thus further validating the model in the largest part of turbulent flame brush.
Furthermore, curves 2 and 3 in Figure 7 show velocities calculated using (14) and (22), respectively, at the leading edges of the two flames investigated by Li et al. [17]. The latter curves have been calculated by neglecting term proportional to e −πξ 2 , as it vanishes at the leading edge. The fact that curves 2 and 3 are very close to one another further supports the extended Zimont model. Similar results were computed for other studied flames.  Table 1.
Thus, the extended Zimont model is consistent with the experimental data obtained from premixed turbulent flames stabilized in impinging jets and, therefore, may be used to substantiate the tested method for evaluating turbulent burning velocity in such flames, as done above (Although the model was also put into question by Peters [40], his arguments were rebutted in earlier papers [28,29].) .
Finally, it is worth raising the following issue. To the the best of the present author's knowledge, the key equation of the extended Zimont model, that is, (19), has yet been substantiated by basic arguments [3,28,29,32,34,38] only in the case of a statistically planar one-dimensional premixed turbulent flame. Therefore, strictly speaking, the extended Zimont model yields the following joint closure: where the ξ-axis is locally normal to the mean flame brush and u n is the flow velocity in this direction. Accordingly, the model may be considered to avoid closing the turbulent scalar flux in the tangential direction, and the following balance equation: may be applied to a flame stabilized in an impinging jet by referring to the same extended Zimont model. While (29) invokes this model and a gradient closure of the transverse turbulent scalar flux (the closure terms vanishes along the jet axis for symmetry reasons), (32) invokes the same model and the BML (7) in order to close the transverse flux (the corresponding terms are not shown on the LHS and RHS of (32), because they cancel one another). Terms I and III on the LHS of (32) are equal to terms 1 and 3, respectively, on the LHS of (29), but the former equation does not involve a counterpart to term 2, which is much smaller than terms 1 and 5 (see Figure 8). Our simulations show that terms IV and V on the RHS of (32) are close to terms 4 and 5, respectively, on the RHS of (29), with the difference between the two sets of terms resulting from different d t and u t substituted into the two equations. Indeed, integration of (32) leads us to (5), with the difference between it and (30) consisting of the second term on the RHS of the latter equation. Comparison of the first and fourth rows in Table 3 shows that the integrated transverse scalar flux is much less than the integrated heat release rate in all the seven flames addressed in the present paper, with the largest difference (10%) being computed in flame no. 6 (set 2 investigated by Stevens et al. [18]). Thus, Figure 8 qualitatively describes the behavior of both terms 1-5 in (29) and terms I-V in (32) in all the flames studied here, but term II vanishes in the latter equation.
If (29) is substituted with (32), then (21) should also be substituted with the following equation: that is, the second term is introduced into the LHS. Accordingly, (22)- (25), which result from (21), do not result from (33). However, because the introduced term is much less than the mean convection term at the leading edges of all the seven flames simulated here (see the fifth row in Table 3), (25) appears to be sufficiently accurate even if (32) and (33) are invoked. At the same time, the variations in g within the studied flames (see Figure 9) may in part be caused by the transverse flux introduced into (33), as well as by the eventual dependencies of u t and d t on z due to variations in turbulence characteristics within the flames.

Conclusions
A method for evaluating burning velocity in premixed turbulent flames stabilized in divergent mean flows, proposed earlier by Lipatnikov and Chomiak [10], was quantitatively validated using numerical approximations of measured axial profiles of the mean combustion progress variable, mean and conditioned axial velocities, and axial turbulent scalar flux, obtained by Cho et al. [15], Cheng and Shepherd [16], Li et al. [17], and Stevens et al. [18] from seven different flames each stabilized in an impinging jet. The method is further substantiated by analyzing the combustion progress variable balance equation that is yielded by the extended Zimont model of premixed turbulent combustion. The compatibility of the model with the aforementioned experimental data is also shown.

Introduction
Displacement speed S d is an important quantity for the modelling of turbulent premixed flames using both levelset [1] and flame surface density (FSD) [2] approaches. The displacement speed is defined as the speed at which the flame surface moves normal to itself with respect to an initially coincident material surface. Statistics of displacement speed for turbulent premixed flames have been studied extensively in previous work based on two-dimensional direct numerical simulations (DNS) with detailed chemistry [3][4][5][6][7][8][9][10][11][12][13] and on three-dimensional DNS with simplified chemistry [14][15][16][17][18][19][20][21][22]. These studies have looked at the effects of local strain rate and curvature effect on S d [13][14][15][16][17][18][19][20][21][22][23][24] and have included the influence of root-mean-square (rms) turbulent velocity fluctuation magnitude u and Lewis number Le [9,15,17,22] on these effects. Peters [1] showed using an order of magnitude analysis that the curvature dependence of displacement speed becomes the leading order contributor in the thin reaction zones regime, whereas the curvature dependence remains relatively weak in the corrugated flamelets regime. This observation was supported using DNS data by Chakraborty [21] who also studied the individual response of the reaction, normal displacement, and curvature components of displacement speed. However, the effects of turbulent Reynolds number on displacement speed have yet to be addressed in detail. Hence, in this paper, the influence of Re t on the statistics of S d has been analysed based on a three-dimensional compressible DNS database of statistically planar turbulent 2 Journal of Combustion premixed flames. The main objectives of this study are as follows: (1) to demonstrate the influence of Re t on the statistical distribution of S d (2) to demonstrate the influence of Re t on the local curvature and strain rate dependencies of S d (3) to indicate the implications of the dependence of S d on Re t in the context of flame surface density-based turbulent combustion modelling.
In order to address the above objectives, a series of DNS runs has been carried out for turbulent Reynolds numbers ranging from 20 to 100. In one group of simulations, the change in Re t is obtained by modifying the Karlovitz number Ka for a given value of Damköhler number Da, whereas in the other group of simulations the change in turbulent Reynolds number is induced by modifying the Damköhler number for a given value of Karlovitz number. The rest of the paper is organised as follows. The mathematical background and the numerical details relevant to this work are discussed in the next two sections of this paper. Following this, results will be presented and subsequently discussed. In the final section, the main findings are summarised and conclusions are drawn.

Mathematical Background
Direct numerical simulations (DNS) of turbulent reacting flows ideally should be carried out in three-dimensions with detailed chemistry. However, such simulations remain extremely expensive and are not yet feasible for parametric studies. Previous work has indicated that displacement speed statistics obtained from three-dimensional DNS with simplified chemistry [14][15][16][17][18][19][20][21][22] are in agreement with those obtained from two-dimensional DNS with complex chemistry [3][4][5][6][7][8][9][10][11][12][13]. In the present study, three-dimensional DNS with simplified chemistry is adopted in order to allow for an extensive parametric study without excessive computational cost. The chemical reaction is represented in terms of single step Arrhenius type chemistry, and the species field is uniquely represented by a reaction progress variable c, which may be defined based on a suitable product mass fraction Y P according to where subscripts 0 and ∞ are used to refer to quantities in unburned reactants and fully burned products, respectively. The transport equation for reaction progress variable is given by whereẇ is the chemical reaction rate, ρ is the fluid density, and D is the diffusivity of reaction progress variable. The transport equation for reaction progress variable can be written in kinematic form for a given isosurface defined by c = c * [6,7] ρ ∂c ∂t + ρu j ∂c ∂x j = ρS d |∇c| c=c * .
Comparing (2a) and (2b), the expression for the displacement speed of a given c = c * isosurface can be obtained as [1,[3][4][5][6][7] It is often useful to decompose the displacement speed S d into three separate components [3][4][5][6][7] where S r , S n , and S t are the reaction, normal diffusion, and tangential diffusion components which are given as where κ m is the arithmetic mean of the two principal curvatures, which is defined as According to (6), flame surface that is curved convex to the unburned reactants has a positive curvature and vice versa. The definition of displacement speed and its components in (4)- (6) suggests that the statistics of the surface density function (SDF = |∇c|) [15,[25][26][27] and curvature κ m are each likely to have a significant effect on the statistical behaviour of displacement speed. The transport equation for the SDF for an isosurface of c is given by [15,[25][26][27] ∂|∇c| ∂t where a T is the tangential strain rate given by It is evident from (7) and (8) that tangential strain rate is likely to affect the statistical behaviour of |∇c|, which in turn can affect the local curvature and strain rate dependence of S r and S n . It is worth noting that in some previous studies [3][4][5][6][7][8][9][10][11][12], the statistics of the density-weighted displacement speed S * d = ρS d /ρ 0 were analysed instead of the statistics of S d itself. For low Mach number unity Lewis number flames, S * d and S d are related by S * d = S d (1 + τc * ) for a given isosurface of c, and thus the statistics of S * d are similar to those of S d . For this reason, and also since S d appears explicitly in the SDF and FSD transport equations [2,15,[25][26][27], the statistics of S d are addressed in the present work. The local strain rate and curvature dependence of S d and its components for different values of turbulent Reynolds number will be discussed in Section 4 of this paper.
It is essential for the purposes of the present study to note the relationship between the turbulent Reynolds number and both the Damköhler and Karlovitz numbers. The turbulent Reynolds number is defined as Re t = ρ 0 u l/μ, where ρ 0 is the density and μ is the dynamic viscosity of the unburned reactants, and l is the integral length scale of the turbulence. For unity Lewis number flames, the turbulent Reynolds number can be scaled as [1] Re where Da = lS L /u δ th is the Damköhler number and Ka = (u /S L ) 3/2 /(l/δ th ) 1/2 is the Karlovitz number. Here, δ th = (T ad − T 0 )/ Max |∇ T| L is the thermal flame thickness with T ad , T 0 and T being the adiabatic flame temperature, unburned gas temperature, and instantaneous dimensional temperature, respectively, while the subscript "L" refers to the unstrained planar laminar flame quantities. Equation (9) indicates that a change of either Da or Ka will lead to a change in the turbulent Reynolds number Re t . Thus, the influence of Re t on the statistics of S d is essentially induced by changes in Da and/or Ka.

Numerical Implementation
Three-dimensional compressible DNS runs for statistically planar turbulent premixed flames have been carried out using a DNS code called SENGA [28] in which the conservation equations of mass, momentum, energy, and species are solved in nondimensional form. The simulations have been carried out in a rectangular domain of size 36.6δ th × 24.1δ th × 24.1δ th . The simulation domain has been discretised using a Cartesian mesh of size 345 × 230 × 230 with uniform grid spacing in each direction. The mean direction of flame propagation is assumed to align with the x 1 -direction and the transverse directions are assumed to be periodic. The boundaries in the direction of mean flame propagation are taken to be partially nonreflecting and the boundary conditions are specified using the Navier-Stokes characteristic boundary conditions (NSCBC) technique [29]. The spatial discretisation is carried out using a 10th order central difference scheme for the internal points and the order of discretisation drops gradually to a one-sided 2nd order scheme at nonperiodic boundaries. The time advancement is carried out using a low storage 3rd order Runge-Kutta scheme [30]. The turbulent velocity fluctuations are initialised using an initially homogeneous isotropic field using a pseudospectral method proposed by Rogallo [31]. The combustion is initailised by a steady unstrained planar laminar premixed flame solution. For the present study, the thermophysical properties such as dynamic viscosity μ, thermal conductivity λ, specific heats C P and C V and the density-weighted mass diffusivity ρD are taken to be constant and independent of temperature following several previous studies [11][12][13][14][15][16][17][18][32][33][34]. For the present study, standard values are taken for Prandtl number (Pr = 0.7) and the ratio of specific heats (γ = C P /C  Table 1. From Table 1, it is evident that Ka remains greater than unity for all cases, which indicates that the combustion belongs to the thin reaction zones regime according to the regime diagram by Peters [1]. In cases B, C, and D the values of u /S L and l/δ th were chosen in such a manner as to vary the turbulent Reynolds number Re t by changing the Damköhler number Da while keeping the Karlovitz number Ka constant (see (9)). On the other hand, in cases A, C, and E, u /S L and l/δ th were chosen to vary Re t by changing Ka, while Da was kept constant. The range of turbulent Reynolds number values considered in this study remains modest although several previous studies [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22] with comparable values of Re t have made valuable contributions to the fundamental understanding of displacement speed in turbulent premixed flames. It is also worth noting that the range of Re t considered here is comparable to that of previous laboratory-scale experiments (e.g., Kobayashi et al. [35]). Despite the limited range, a number of important Re t effects on the displacement speed statistics have been identified in the present parametric study and are discussed in detail in Section 4 of this paper. In all cases, the flame-turbulence interaction takes place under conditions of decaying turbulence, in which simulations should be carried out for a minimum time t sim = Max(t f , t c ), where t f = l/u is the initial eddy turn over time and t c = δ th /S L is the chemical time scale. All the cases considered here were run for one chemical time scale t c , which corresponds to a time equal to 2.0t f in case D, 3.0t f in cases A, C, and E, and 4.34t f for case B. The present simulation time is comparable to the simulation times used for previous DNS studies which focused on the analysis of displacement speed statistics [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][23][24][25][26][27]. By the time statistics were extracted, the global turbulent kinetic energy and volume-averaged burning rate were no longer changing rapidly with time. The temporal variations of these quantities were presented in [36], and thus are not repeated here for the sake of conciseness. Analysis of the volume-averaged burning rate indicates that the effects of unsteadiness have become weak by the time when statistics were extracted. At this time,  the global level of turbulent velocity fluctuation had decayed by 52.66%, 61.11%, 45%, 24%, and 34% in comparison to the initial values for cases A to E, respectively. By contrast, the integral length scale had increased by factors of between 1.5 to 2.25, ensuring that sufficient numbers of turbulent eddies were retained in each direction to obtain useful statistics. The values for u /S L , l/δ th , and δ th /η at the time when statistics were extracted have been presented in Table 2 of [36]. For all cases, the thermal flame thickness δ th remains greater than the Kolmogorov length scale η at that time, which indicates the combustion in all cases still belongs to the thin reaction zones regime. Moreover, the turbulent Reynolds number Re t for all cases remains greater than unity (i.e., Re t > 1) for the values of u /S L and l/δ th at the time when statistics were extracted, which suggests that all the flames remain turbulent according to the regime diagram by Peters [1].

Flame-Turbulence
Interaction. The contours of the reaction progress variable c in the x 1 -x 2 mid-plane of the domain at time t sim = 1.0 δ th /S L are shown in Figures 1(a)-1(e) for cases A to E, respectively. It is evident that the extent of the wrinkling of the c isosurfaces increases with increasing u , that is, in going from case A through to case E (see Table 1). Moreover, the contours of c representing the preheat zone (i.e., c < 0.5) are much more distorted than those representing the reaction zone (i.e., 0.7 < c < 0.9). This tendency is more prevalent for the case with the highest Karlovitz number (case E), since the scale separation between δ th and η increases with increasing Ka, allowing more energetic eddies to enter into the preheat zone and causing greater distortion of the flame. This situation arises when the negative contribution of the molecular diffusion rate dominates over the positive semidefinite chemical reaction rate, as discussed in detail by Gran et al. [4]. Moreover, Figure 2(e) indicates that the probability of finding negative values of S d /S L is greatest for the highest Karlovitz number (case E). This behaviour can be explained in terms of the scaling analysis of Peters [1] for unity Lewis number flames, which suggested that the different components of displacement speed scale as

Probability Density Functions of Displacement Speed
where v η is the Kolmogorov velocity scale. Here, in accordance with Peters [1], (S r + S n ) is taken to scale with S L and the curvature κ m is taken to scale with the Kolmogorov length scale 1/η. These scalings of (S r + S n ) and κ m were confirmed by DNS data presented in [21]. Equation (10) clearly suggests that in the thin reaction zones regime (i.e., Ka > 1) the effects of (S r + S n ) are likely to weaken progressively in comparison to the contribution of S t with increasing Karlovitz number. This suggests that for large values of Karlovitz number the negative contribution of S t more readily overcomes the predominantly positive contribution of (S r + S n ) to yield a negative value of S d . This eventually gives rise to an increased probability of finding negative values of S d for flames with increasing value of Karlovitz number, as suggested by Figure 2. In order to demonstrate that the observed negative values of S d originates principally due to S t , the pdfs of (S r + S n )/S L are shown in Figure 3, respectively, for all cases. Comparing Figures 2 and 3 reveals that the probability of finding negative values of (S r + S n ) is significantly smaller than that of finding negative values of S d .

Comparing Figures 2(a)-2(e) and Figures 3(a)-3(e)
, it can be seen that the most probable values of S d /S L and (S r +S n )/S L remain of the order of ρ 0 /ρ and that the widths of the displacement speed pdf for a given c isosurface increases with increasing u (i.e., going from case A to case E). In order to explain this behavior, it is useful to examine the pdfs of the normalised tangential component of displacement speed S t /S L which are depicted in Figure 4 and show that these pdfs remain almost symmetrical with a peak close to a value of zero. This behaviour is consistent with the statistically planar nature of these flames. However, the greater extent of flame wrinkling observed in Figure 1 for larger values of u results in a greater spread of flame curvature, which in turn increases the width of the S t /S L pdfs according to (5). The negligible mean value of S t essentially leads to the mean values of S d and (S r + S n ) becoming almost equal to each other, but the statistical variation of S d induced by S t leads to broader pdfs of S d than those of (S r +S n ). These results clearly suggest that the distributions of S d and its components S r , S n and S t generally broaden with u /S L , which essentially leads to broadening of the distributions with increasing Re t ∼ (u /S L ) 2 Da ∼ (u /S L ) 4 /Ka 2 ∼ (l/δ th ) 4/3 Ka 2/3 when either Da or Ka is held constant.

Strain Rate and Curvature Dependences of |∇c|. For low
Mach number unity Lewis number flames, both density ρ and nondimensional temperature T remain uniform on a given isosurface of c and thus the statistical behaviours of S r and S n (see (5)) are mainly affected by the dependence of the SDF |∇c| on strain rate and curvature. The values of |∇c| × δ th conditional on normalised tangential strain rate a T × δ th /S L are shown in Figure 5(a) for the isosurface at c = 0.8 (close to the location of maximum reaction rate), which suggests clearly that |∇c| and a T are positively correlated for all cases. The correlation coefficients between |∇c| and a T are shown in Figure 5(b) for values of c between 0.1 and 0.9, and these indicate that the positive correlation is maintained throughout the flame brush. This correlation can be explained by noting that the dilatation rate ∇ · − → u can be expressed as For unity Lewis number flames, the dilatation rate can be scaled as ∇ · − → u ∼ τ S L /δ th [21], whereas a T can be scaled as a T ∼ u /l [37], which suggests that Equation (11b) essentially shows that for low values of the Damköhler number, a T remains much greater than ∇ · − → u . This essentially suggests that an increase in tangential strain rate is associated with a decrease in the normal strain rate, since the normal strain rate scales as a n ∼ −a T when dilatation rate effects are much weaker than those of turbulent straining. Under decreasing normal straining, the isoscalar lines are brought closer to each other and hence give rise to an increase in |∇c|. The decreasing trend of normal strain rate with increasing tangential strain rate is reflected in the positive correlation between |∇c| and a T throughout the flame brush. It can be seen from Figure 5 that the correlation between |∇c| and a T does not show any monotonic trend with either Karlovitz number (Damköhler number) or turbulent Reynolds number variations when Damköhler number (Karlovitz number) is held constant.
The variation of |∇c| × δ th conditional on normalised curvature κ m × δ th is shown in Figure 6(a) for the isosurface at c = 0.8. Both positive and negative correlation trends can be observed, which is consistent with several previous studies [14-22, 26, 27]. As a result of these opposing trends, the net correlation between |∇c| × δ th and κ m × δ th turns out be weak throughout the flame brush, as evident from the correlation coefficients between |∇c| and κ m presented for different values of c across the flame brush in Figure 6(b). Both the positively and negatively correlating branches can be explained using (11a). In turbulent premixed flames the dilatation rate ∇ · − → u is negatively correlated with curvature κ m , because of focusing (defocusing) of heat in the negatively (positively) curved regions of the flame surface [18,21]. The tangential strain rate a T is also found to be negatively correlated with curvature [18,21]. The correlation coefficients for the ∇ · − → u − κ m and a T − κ m correlations for different c values across the flame brush are shown in Figures 7(a) and 7(b), respectively. Comparing the correlation coefficients of the ∇ · − → u − κ m and a T − κ m correlations for different cases reveals that the strength of both of these negative correlations decreases with increasing u /S L and Re t ∼ (l/δ th ) 4/3 Ka 2/3 when either Da or Ka is held constant.
The combination of the positive correlation between |∇c| and a T and the negative correlation between a T and κ m leads to small values of |∇c| in regions of high positive curvature. This gives rise to the negative correlation trend between |∇c| and κ m , as observed in Figure 6. However, the dilatation rate ∇ · − → u may locally attain large values in the highly negatively curved regions of the flame due to strong focussing of heat. Hence, ∇ · − → u may locally exceed a T and induce an increase in a n according to (11a). As isoscalar lines move apart from each other under increasing normal strain rate, the increase in a n with increasing negative curvature leads to a positive correlation between |∇c| and κ m in negatively curved regions, as also observed in Figure 6.

Curvature
Dependence of S d . The mean normalised displacement speed S d /S L conditional on normalised curvature κ m × δ th is shown in Figure 8(a) for c = 0.8, which indicates that S d /S L is negatively correlated with curvature and that the correlation is nonlinear in nature. These observations are consistent with several previous DNS studies with both detailed and simplified chemistry [3-11, 13-18, 20-22]. The variation of the correlation coefficient between S d and κ m for different c isosurfaces across the flame brush is shown in Figure 8(b), which reveals that strength of the negative correlation decreases with increasing Re t ∼ (u /S L ) 2 Da ∼ (u /S L ) 4 /Ka 2 ∼ (l/δ th ) 4/3 Ka 2/3 when either Da (see cases A, C, and E) or Ka (see cases B, C, and D) is held constant (see Table 1). The tangential diffusion component of displacement speed S t is deterministically negatively correlated with curvature with correlation coefficient equal to unity for unity Lewis number flames (see (5)). Thus, the nonlinearity of the correlation between S d and κ m originates due to the curvature dependence of S r and S n .
Mean values of the reaction component S r /S L conditional on normalised curvature κ m × δ th are shown in Figure 9(a) for the isosurface at c = 0.8. The plot shows both positive and negative correlation trends. For low Mach number unity Lewis number flames, the curvature dependence of S r is principally determined by the curvature dependence of |∇c|. The negative (positive) correlation between |∇c| and κ m ultimately gives rise to the positive (negative) correlation between S r and κ m as observed in Figure 9. As a result of the weak net correlation between |∇c| and κ m , the net correlation between S r and κ m also turns out to be    weak, as indicated by the correlation coefficients plotted in Figure 9(b) Note that Figure 9(b) includes isosurfaces of c from 0.5 to 0.9 only, since the reaction rateẇ is negligible in the preheat zone (see also (5)).
Corresponding mean values of the normal diffusion component S n /S L conditional on normalised curvature κ m × δ th are shown in Figure 10(a) for c = 0.8, while the correlation coefficients between S n and κ m are shown in Figure 10(b) for values of c from 0.1 to 0.9. The magnitude of the normal molecular diffusion rate in the thin reaction zones regime can be scaled as | − → N · ∇(ρD − → N · ∇c)| ∼ ρD/δ 2 where δ is the characteristic flame thickness which in turn can be scaled as |∇c| ∼ 1/δ. This suggests that S n scales as D/δ towards the reactant side of the flame brush and as −D/δ towards the product side. As a result of this, the positively (negatively) correlating branch of the |∇c| − κ m correlation gives rise to a positive (negative) correlation between S n and κ m towards the reactant side, where S n assumes predominantly positive values. By contrast, the positive (negative) correlation between |∇c| and κ m results in a negative (positive) correlation between S n and κ m towards the product side where S n assumes predominantly negative values. Both positive and negative correlation trends can be discerned for all the c isosurfaces shown in Figure 10(b). As a result of this, the net correlation between S n and κ m turns out to be weak throughout the flame brush. The positive and negative correlating trends with κ m observed for both S r and S n with κ m give rise to positive and negative correlating branches in the variation of (S r + S n )conditional on κ m × δ th as shown in Figure 10(c) for c = 0.8. This nonlinear correlation between (S r + S n ) and κ m ultimately induces a nonlinear curvature dependence of displacement speed S d (see Figure 8(a)), noting that the correlation coefficient arising from the deterministic relation between S t and κ m remains close to −1.0 throughout the flame brush for all cases. The correlation coefficients between (S r + S n ) and κ m are shown in Figure 10(d) for values of c from 0.1 to 0.9, and demonstrate that the (S r + S n ) − κ m correlation is much weaker than the S t − κ m correlation in all cases. Hence, the net correlation between S d and κ m turns out to be negative throughout the flame brush for all the cases considered here (see Figure 8(a)).

Tangential Strain Rate Dependence of S d .
The mean values of normalised displacement speed S d /S L conditional on normalised tangential strain rate a T × δ th /S L are shown in Figure 11(a) for c = 0.8 which suggests that S d and a T are positively correlated for all cases. The observed strain rate dependence of S d is found to be qualitatively consistent with previous DNS studies for flames with global Lewis number close to unity [3, 5, 12, 14, 15, 17-19, 21, 22]. The variations of the correlation coefficients between S d and a T for different c isosurfaces are shown in Figure 11(b), which reveals that the strength of the correlation between S d and a T decreases with increasing u /S L when either Da or Ka is held constant. This is consistent with the analytical treatment by Joulin [38] who indicated that the strain rate dependence of S d weakens with increasing frequency of straining. In order to explain the observed behaviour it is useful to look into the dependence of S r , S n , and S t on a T . For low Mach number unity, Lewis number flames the strain rate dependence of S r is essentially determined by the strain rate dependence of |∇c|. It has been shown already that |∇c| and a T are positively correlated (see Figure 5), which gives rise to a negative correlation between S r and a T according to (5). This is substantiated by the correlation coefficients between S r and a T as shown in Figure 11(c) for the isosurfaces from c = 0.5 to 0.9. The correlation between S r and a T does not show any monotonic trend with either Karlovitz number (Damköhler number) or turbulent Reynolds number variations when the Damköhler number (Karlovitz number) is held constant, due to the lack of any such trend in the |∇c| − a T correlation. It has been discussed already that S n scales as S n ∼ D/δ on the reactant side of the flame brush, and as S n ∼ −D/δ on the product side. Since there is a positive correlation between |∇c| ∼ 1/δ and a T , this gives rise to a positive (negative) correlation between S n and a T on the reactant (product) side. This behaviour can be seen in the correlation coefficients plotted in Figure 11(d), which also indicates that the correlation coefficient between S n and a T does not show any monotonic trend with changes in either Da or Ka. However, it has been found that the correlation between S n and a T in general remains stronger for smaller values of u when either Da or Ka is held constant.
Since S r remains negligible towards the reactant side of the flame brush and S n is positively correlated with a T , the correlation coefficient between (S r + S n ) and a T remains positive towards the reactant side (see Figure 11(e)). By contrast, both S r and S n are negatively correlated with a T towards the product side, which leads to a locally negative correlation between (S r + S n ) and a T .
Since a T and κ m are negatively correlated, the tangential diffusion component S t is positively correlated with a T because S t is proportional to the negative of curvature (i.e., S t = −2Dκ m ). This positive correlation is evident from Figure 11(f), in which the correlation coefficient between S t and a T is shown for isosurfaces of c across the flame brush. The negative correlation between a T and κ m (see Figure 7(a)) arises principally due to heat release, and thus, the strength of this correlation weakens with increasing Karlovitz number, since the effects of turbulent velocity fluctuations are likely to mask the heat release effects as the combustion tends more towards the broken reaction zones regime [1]. Moreover, with decreasing Damköhler number the effects of heat release weaken in comparison to those of turbulent straining, as suggested by (11b). As a result of this, the strength of the negative correlation between a T and κ m weakens with decreasing (increasing) Damköhler number (Karlovitz number) when Ka (Da) is held constant (see Figure 7(a)), and this gives rise to a corresponding change in the correlation between S t and a T (see Figure 11(f)). This essentially suggests that the strength of the positive correlation between S t and a T decreases with increasing turbulent Reynolds number Re t when either Ka or Da is held constant. The effects of the positive S t − a T and negative a T − κ m correlations remain strongest in the reaction zone and the strengths of these correlations decrease towards both the reactant and product sides of the flame brush (see Figures 11(f) and 7(b), resp.).
It is evident from Figures 11(e) and 11(f) that the positive correlations between (S r + S n ) and a T , and between S t and a T , aid each other towards the reactant side of the flame brush. By contrast, towards the product side the positive correlation between S t and a T overcomes the negative correlation between (S r +S n ) and a T to result in a net positive correlation between S d and a T in all the cases considered here (see Figure 11(a)). The weakening of the correlation between S t and a T with increasing Re t ultimately leads to decrease in the strength of the correlation between S d and a T when either Da or Ka is held constant (see Figures 11(a) and 11(f)).

Implications for FSD Modelling.
In the context of Reynolds Averaged Navier Stokes (RANS) modelling, the generalised FSD Σ gen is defined as [39] Σ gen = |∇c| , where Q indicates the Reynolds averaged value of a general quantity Q. It has been demonstrated in [15,23,26,27] that the statistical behaviour of |∇c| and the terms of its transport (7) in the thin reaction zones regime depend significantly on the choice of progress variable isosurface c = c * . Hence, it is more appropriate to use the generalised FSD (see(12a)) which does not depend on the choice of isosurface. Reynolds averaging of (7) yields the transport equation for the generalised FSD [2,23,24] ∂Σ gen ∂t where Q s = Q|∇c| / |∇c| indicates a surface averaging operation for a general quantity Q [39]. The first term on the left hand side is the transient term while the second term is the advection term. The three terms on the right hand side represent the effects of strain rate, curvature, and propagation, respectively. It is evident from (12b) that the statistical behaviour of the curvature and propagation terms depends on the statistics of displacement speed. These terms can be rewritten in terms of the displacement speed components as The variations of all of the terms appearing in these equations with Favre averaged reaction progress variable c = ρc / ρ are shown in Figure 12. The Reynolds averages are evaluated by ensemble averaging the relevant quantities in the directions normal to the direction of mean flame propagation.
It can be seen from Figure 12 that the FSD curvature term S d ∇ · − → N s Σ gen (line with tick symbols) assumes mostly negative value throughout the flame brush for all cases although small positive contributions are found towards the reactant side of the flame brush for small Reynolds numbers (i.e., cases A and B). The magnitude of the term 2 (S r + S n )κ m s Σ gen (line with filled triangles) remains comparable to the contribution of the term −4 Dκ 2 m s Σ gen (line with open squares) for small values of Re t (i.e., cases A and B). However, the contribution of 2 (S r + S n )κ m s Σ gen remains smaller than that of −4 Dκ 2 m s Σ gen for larger values of Re t and the statistical behaviour of the FSD curvature term is then principally governed by the latter contribution. It is also evident from Figure 12 that the magnitude of (−4 Dκ 2 m s Σ gen ) becomes increasingly more important than the magnitude of 2 (S r + S n )κ m s Σ gen with increasing Karlovitz number (Damköhler number)  In some previous studies, the FSD curvature and propagation terms have been modelled by replacing S d with S L to yield the modelled expressions 2 S L κ m s Σ gen and −∇ · S L − → N s Σ gen , respectively; [40,41]. Figure 12 further indicates that there are significant differences between the original terms and these modelled expressions, which are shown, respectively, by a line with filled circles and a line with open triangles. This essentially indicates that S d cannot be approximated simply by S L in the modelling of the FSD curvature and propagation terms. Moreover, Figure 12 suggests that modelling of κ 2 m s is necessary for the closure of the FSD curvature term especially in the thin reaction zones regime since the term −4 Dκ 2 m s Σ gen remains close to 2 S d κ m s Σ gen throughout the flame brush for high values of the turbulent Reynolds number. The effects of turbulent Reynolds number on the modelling of the FSD curvature and propagation terms are beyond of the scope of this study and will be addressed in future work.

Conclusions
The turbulent Reynolds number dependence of the statistics of displacement speed in turbulent premixed flames has been studied using a DNS database of statistically planar flames, in which the variation of turbulent Reynolds number Re t from 20 to 100 has been brought about by modifying either Damköhler number Da or Karlovitz number Ka independently of each other. As both Da and Ka are functions of u /S L and l/δ th , the variation of Re t for the present database has been achieved by modifying the initial values of u /S L and l/δ th simultaneously.
The mean value of displacement speed remains positive throughout the flame brush, but there is nonzero probability of finding negative values of displacement speed, in accordance with the findings of previous work. It has been shown here that the probability of finding negative displacement speed increases with increasing turbulent Reynolds number when the Damköhler number is held constant.
The dependences of tangential strain rate and dilatation rate on flame curvature are shown to have a significant influence on the strain rate and curvature dependences of the SDF |∇c|, which in turn affects the statistical behaviour of displacement speed in response to strain rate and curvature. It has been found that the variation of turbulent Reynolds number does not alter the qualitative nature of the correlations between tangential strain rate and dilatation rate with curvature, but the strength of these correlations is found to weaken with increasing turbulent Reynolds number when either Damköhler or Karlovitz number is held constant. Similarly, the strain rate and curvature dependences of displacement speed are shown to weaken with increasing turbulent Reynolds number when either Damköhler or Karlovitz number is kept unaltered.
Detailed physical explanations are provided to explain the influences of turbulent Reynolds number on the strain rate and curvature dependences of displacement speed in terms of the individual response of the reaction, normal diffusion and tangential diffusion components of displacement speed.
The implications of the turbulent Reynolds number dependence of displacement speed for the modelling of the FSD transport equation has been explored. The statistical behaviour of the FSD curvature term contributions arising due to the combined reaction and normal diffusion components of displacement speed and to the tangential diffusion component of displacement speed have been analysed in detail. At low turbulence Reynolds numbers, the magnitudes of these contributions remain comparable. As the turbulence Reynolds number increases, the curvature dependence of the tangential component of displacement speed ensures that this contribution becomes increasingly important and is the predominant influence on the FSD curvature term for larger values of turbulent Reynolds number. Hence, the curvature dependence of the tangential diffusion component of displacement speed cannot be ignored for accurate modelling of the FSD curvature term for flames within the thin reaction zones regime.
It is worth noting that the effects of differential diffusion of heat and mass are not taken into account in the present study. However, some previous DNS studies with detailed chemistry [3][4][5][6][7]9] suggested that intermediate species may play an important role even for the flames with global Lewis number close to unity. Moreover, only a moderate range of turbulent Reynolds number has been considered in this study, so three-dimensional DNS studies with detailed chemistry and at higher values of turbulent Reynolds number will be needed for deeper understanding and for the purpose of quantitative predictions.

Introduction
Rising oil prices, national security concerns, and protection of the environment have resulted in research and development of alternative fuels. Biofuels form one category of alternative fuels that can be used with minor modifications to current automobiles. Biofuels are renewable and available in the US, and they have energy content similar to petroleumbased gasoline and diesel. One category of biofuels is derived from the transesterification of vegetable oils, such as canola and soy oil; examples of such biofuels are soy methyl ester (SME) and canola methyl ester (CME). These fuels are low in sulfur and can be produced locally; also, they result in reduced net carbon emissions.
There have been a number of studies on the use of these biofuels in diesel engines. Durbin et al. [1] measured the pollutant emissions and particulate matter (PM) from four commercially available diesel vehicles. The engines were operated on diesel, B100 (100%) biofuel, B80 (80% biofuel by volume) biofuel, and synthetic diesel fuels. In three of the four vehicles tested, there were no significant differences in the particulate matter (PM) emitted from engines using B100 and diesel fuels. Scholl and Sorenson [2] tested soy methyl ester in a four-cylinder, four-stroke cycle, normally aspirated direct-injection diesel engine and compared the results to those obtained with conventional diesel fuel. The authors found that the results varied with engine operating conditions, such as ignition delay and cylinder pressure rise. McCormick et al. [3] tested various biofuel and diesel fuels in a six-cylinder, direct-injected, turbocharged, four-stroke cycle engine. Particulate matter per unit power was measured and was found to be significantly less (by 300%) than that obtained with the use of petroleum-based diesel. Wang et al. [4] tested B35 (35% biofuel by volume) soy methyl ester biofuel blends in unaltered Cummins tractor truck engines and measured the resulting PM and pollutant emissions. It was shown that 25% less PM was produced when the engines were operated using the blended mixture than when operated on petroleum diesel.
Dorado et al. [5] used transesterified waste olive oil in a 2500 cc, three-cylinder, four-stroke cycle, water-cooled, direct-injection diesel engine. The authors showed that by using biofuels, CO emission was reduced by up to 58.9%, and NO emission was increased up to 37.5% when compared to the emissions obtained with the use of petroleum diesel. Also, CO emissions were shown to correlate with the production of PM. Labeckas and Slavinskas [6] conducted studies on naturally aspirated, four-stroke, four-cylinder, water-cooled direct-injection diesel engine with diesel fuel and cold pressed rapeseed oil (RO). An increase in the NOx emissions was observed with the use of rapeseed oil. The increase in NOx formation was attributed to the double bonds present in the fuel and the longer ignition delay caused due to the presence of rapeseed fatty acid. Agarwal [7] reviewed the production, characterization, and current statuses of vegetable oil, biofuel, and alcohols as well as the experimental research work carried out in various countries. The low CO emissions from methyl esters in comparison to raw vegetable oil were attributed to better spray quality. Mueller et al. [8] investigated the increase in NOx emission in a heavyduty compression ignition engine with soy biofuel. Results suggested that the NOx increase observed with the use of biofuels was not quantitatively determined by a change in a single fuel property, but rather was the result of a number of coupled mechanisms whose effects may tend to reinforce or cancel one another under different conditions, depending on the specific combustion and fuel characteristics. Due to the oxygenated nature of biofuel, the soot production was suppressed, which reduced the total radiative heat losses in the cylinder, thereby increasing the temperatures and thermal NOx production. Kousolidou et al. [9] studied the performance of B10 biodiesel fuel of palm oil in a lightduty common-rail Euro 3 engine. It was found that the use of biodiesel blends resulted in a reduction in PM emissions whereas the effects on NOx emissions were only marginally different from those measured with petroleum diesel. It was concluded that blends of up to 10% biodiesel could be used in current diesel engines without any significant effects on the performance.
A limited number of studies have been conducted to document the combustion properties of biofuels in an openflame environment. Love et al. [10] measured the radiative heat flux and the emissions from laminar partially premixed flames of prevaporized mixtures of various fuels with air. The emission indices of NOx and CO of both petroleum-derived and biofuels agreed well with those measured in engines, as reported in the literature. In a companion study, Love et al. [11] used the same experimental configuration to study the effect of iodine number on NOx formation in laminar partially premixed flames of canola methyl ester, soy methyl ester, and methyl stearate. It was found that the peak NOx concentration increased significantly with the fuel iodine number, indicating a strong correlation between the chemical structure of the fuel and the NOx emission. It was also concluded that the Fenimore mechanism (not the Zeldovich mechanism) was the dominant route in NOx production in laminar fuel rich biofuel flames. Erazo et al. [12] investigated the combustion properties of canola methyl ester and petroleum diesel spray flames at atmospheric pressure. The drop size and velocity distribution, in-flame temperature, and emissions were measured. Whereas the diesel spray flame burned in a mostly heterogeneous combustion mode, the canola methyl ester underwent homogeneous gas-phase combustion. The NO concentration was lower in the CME flame at all flame heights than that of the diesel spray flame. The CME flame was also cooler (by 200 K) than the diesel spray flame. The reduction in NOx emissions observed was in contrast to the NOx increase recorded in previous engine studies and was attributed to the continuous combustion in the spray flame as opposed to the intermittent combustion that occurred in engines.

Objective
The present experimental investigation is an extension of the studies by Love et al. [10,11] and is aimed at understanding the effects of turbulence on the combustion characteristics of vaporized blends of CME/diesel flames. The specific objectives of this research project are to measure the global radiation, global emissions, in-flame emissions, and temperature profiles in flames of blended CME and No. 2 diesel fuel under turbulent conditions.

Experimental Setup and Test Conditions.
A schematic diagram of the setup is presented in Figure 1. Hightemperature heating tape was used to heat the flow lines carrying air to desired temperatures in order to completely vaporize the liquid fuel. The temperature was regulated using a dual zone temperature controller connected to two K-type thermocouples. One thermocouple was placed just prior to the fuel injection and the other was positioned downstream right before the burner. The air flow was monitored using a digital mass flow meter. The liquid fuel was stored in a tank which was pressurized using a nitrogen cylinder and was metered using a Rotameter. After the air reached the desired temperature, the fuel was injected into the air. The air temperature was set to 460 • C for all the experiments; this temperature was considerably above the highest boiling point, 420 • C, but below the pyrolysis temperature of CME to enable the complete vaporization of the fuel without coking. The burner had an inner diameter of 9.5 mm and was fed symmetrically through four bypass ports. All experiments were conducted in a vertical steel test chamber of cross-section 76 cm × 76 cm and 150 cm height. The top of the combustion chamber was connected to the atmosphere through an exhaust duct. The ambient pressure of the laboratory was maintained slightly above (20 Pa) the atmospheric pressure to provide a positive draft inside the test chamber so that the combustion products did not leak into the main laboratory. Additional details of the setup are provided by Dhamale [13].
Three blends of CME and No. 2 diesel fuel in different proportions, B25, B50, and B75 (e.g., B25 is 25% CME by volume splash blended with 75% diesel) were tested. The properties of the blends are presented in Table 1. The oxygen content in the fuel increases with increase in SME concentration in the blend. Also, the heating value decreases as the CME concentration is increased because CME has a lower heating value (by about 10%) than diesel. Measurements were made for three Reynolds numbers: 2700, 3600, and 4500. The Reynolds numbers (Re) were  computed based on the bulk fluid velocity, burner inner diameter, and the mixture viscosity. The mixture viscosity was computed using Wilke's equation [14] at the burner exit temperature. The initial equivalence ratio was maintained at 7, thus, the fuel/air mixture was fuel rich. This condition was chosen to represent the local equivalence ratio in the turbulent diffusion-controlled combustion zone away from the injector in a diesel engine at full-load conditions. Since the mass flow rate was changed to vary the Reynolds number and the equivalence ratio remained constant, the change in Reynolds number was proportional to the change in the bulk velocity of the fuel-air mixture. The carbon input rate increased with the increase in Reynolds number for each fuel, as seen in Table 1. The measurements described below were repeated four to six times, and uncertainties were computed using standard procedures; details of the calculations are presented by Dhamale [13]. The uncertainties are shown as error bars in the relevant figures.

Radiation Measurement.
Radiation measurements were made with a high-sensitivity pyrheliometer. The pyrheliometer (with a view angle of 150 • and absorptivity of 0.96) was placed at a distance of 150 cm from the burner so that its view angle covered the entire flame length and the inverse square law was satisfied. A data-acquisition board along with suitable software was used to sample the measured radiative heat flux. Each test was run for a duration of 1 minute with a sampling rate of 10 Hz. The background radiative heat flux was sampled for 30 seconds at 10 Hz. The background radiation was subtracted from the measured heat flux to obtain the modified radiative heat flux of the flame. The modified radiative heat flux, R, was used to compute the radiative fraction of heat release, F, as given below: Here, L is the distance from the flame centerline to the pyrheliometer,ṁ is the mass flow rate of the liquid fuel, and LHV is the lower heating value of the fuel. The radiative fraction of heat released is the ratio of the total energy radiated by the flame to the net energy content of the burnt fuel.

Global and In-Flame Emission Measurement.
The volumetric concentrations of CO 2 , CO, and NOx in the exhaust were measured using a gas emissions analyzer. A conicalshaped flue was placed above the flame to collect the exhaust gases. At the top of the flue, a quartz probe with a 1 mm inner diameter was mounted to sample the exhaust gases. The exhaust gases were then passed through an ice bath to remove moisture and then through a filter to eliminate particulate matter. The concentration measurements of CO and NOx were then converted to emission indices on a mass basis (g of species/kg of fuel) [15] to compensate for dilution: Here, X i represents the mole fraction, N is the number of carbon atoms in the fuel, and MW is the molecular weight. It is assumed that all the carbon in the fuel is converted into CO and CO 2 and that the soot concentration is small.

3.4.
In-Flame Temperature Measurement. The radial temperature profiles were measured with an R type thermocouple with a bead size of 0.3 mm. The thermocouple was mounted on a traverse mechanism and was manually traversed through the flame. The thermocouple output was read and recorded with the help of a LabVIEW program at a sampling rate of 1 Hz, and the temperatures were averaged over a period of 40 seconds. The temperature measurements were corrected for radiation, conduction, and convection losses [13]. was obtained by measuring the laser beam intensity with no flame and through the flame. The flame intensity was also measured and subtracted from the total intensity measurements to get only the laser attenuation through the flame. The attenuation was then converted into soot volume fraction, f v : Here, I o is the attenuated laser intensity, I s is the incident laser intensity, λ is the laser beam wavelength, and δ is the flame thickness along the laser beam length. The spectral extinction coefficient, k λ , was assumed to be constant corresponding to that of petroleum diesel soot [16] in the calculations.
3.6. Exit Velocity Measurement. The exit vertical velocity of the fuel air mixture was measured with a laser Doppler velocimeter (LDV). The system was based on the 514.5 nm wavelength from an argon ion laser. Frequency shifting was used with a Bragg cell which shifted the frequency of one beam by 40 MHz; the signal was then downshifted by 38 MHz, yielding a net shift of 2 MHz. The scattered light was collected in the off-axis forward-scatter mode. Magnesium oxide (MgO) particles of nominal diameter of 4 microns were used as seed particles. A fraction (20%) of the air supply was diverted into the cylinder containing the MgO powder. The air intake port in the cylinder was positioned for tangential entry; the fine particles were carried through the exhaust port located at the center of the cylinder. The air carrying the seed particles was introduced 20 cm below the burner exit where it was mixed with the fuel-air mixture. numbers. The flames were brushy in appearance due to the turbulent conditions at the burner exit; this can be seen in the pictures in Figure 3 taken with an exposure time of 1/4000 s. The presence of significant soot content in the flame provides the yellow hue. The exit equivalence ratio was 7; therefore, enough air was not supplied at the exit for complete combustion of the fuel and the air needed to be entrained from the surroundings. The flame heights for all the fuels were comparable. Near-burner exit images of pure CME and diesel flames are shown in Figure 4. These images were captured with an exposure time of 1 s. The nearburner regions appear blue due to the dominant gas-phase reactions in this region. The blue region extended axially to larger heights as the biofuel content in the fuel was increased. The presence of oxygen in the fuel molecule of CME helps in the oxidation of the fuel, resulting in an increase in the homogeneous reaction zone near the burner.

Exit Velocity Profiles.
The mean velocities of the fuelair mixture at the exit of the burner for the pure diesel and CME flames are given in Figure 5 for the three Reynolds numbers. In order to maintain an initial equivalence ratio of 7, the amount of air supplied for the biofuel was smaller for each Reynolds number (due to the presence of oxygen in the fuel molecule); thus, the bulk velocity was lower for the biofuel flame than the diesel flame. This is observed in the mean velocity profiles. The peak axial velocity was attained at the centerline; the profiles are not flat because of the transitional nature of the flows at these Reynolds numbers. The mean centerline axial velocity was 22% higher than the bulk velocity. The corresponding turbulence intensities are presented in Figure 6. The turbulence intensities were higher near the edge than at the center (typical of a pipe flow). The maximum turbulence intensity (near the edge) increased with Reynolds number.

Flame Radiation.
The radiative heat fraction for the various flames at the three Reynolds numbers is plotted in Figure 7. The flame radiation is emitted by the burning soot particles and the high-temperature gases present in the flame; it, thus, depends on the local temperature, gas emissivity, and soot emissivity. The radiative fraction of the pure CME flames was slightly lower than that of the diesel flames at all Reynolds numbers. In general, the soot content in the diesel flames was higher than that of the pure CME flames, but the maximum temperatures in the far-burner region were higher in the pure biofuel flames than in the diesel flames. Thus, the slightly lower radiative fraction in the pure CME flames could be attributed to the dominant effect of soot content in the flame. In the case of the flames of the blends, the radiative fraction of the B75 flame was higher than that of the diesel flame, whereas the radiative fractions of the B50 and B25 flames were comparable to the radiative fraction of the diesel flame. While the soot content in these flames was lower than that of the pure fuel flames, the flame temperatures differed from those measured in the pure fuel flames. The in-flame temperatures were highest for the B75 flame at 25 and 50% flame heights, compared to the flames of the B25 and B50 blends [13]; the high temperatures resulted in the B75 flame having the highest radiation fraction. In general, the radiative fraction increased as the Reynolds number was increased from 2700 to 3600 and then decreased with a further increase in Reynolds number to 4500 (except for the B50 flame). As noted earlier, the Reynolds number increase was achieved by increasing the bulk velocity, thereby increasing the air entrainment as the Reynolds number was increased. Thus, the highest radiative fraction at the Reynolds number of 3600 could be attributed to the tipping of the balance between the increased soot formation due to the larger carbon input rate and the higher soot oxidation rate with increased air entrainment. Figure 8 shows the global emission index of CO for all the flames. The CO emission index was lowest for the pure CME flame (one-third of the value of the diesel flame at a Reynolds number of 2700 and onehalf the value of the diesel flame at Reynolds numbers of 3600 and 4500). The B25 and B50 flames had a CO emission index comparable to that of the diesel flame whereas the CO emission index of the B75 flame was slightly lower than that of the diesel flame. As the biofuel content in the fuel is increased, the oxygen content in the fuel is increased, and the increased oxygen would be expected to help in the oxidation of CO to CO 2 . However, the global emission index of CO of the B25 flame was slightly higher than that of the diesel flame. This is a result of the coupling effect between CME and diesel. The carbon input rate in the diesel and B25 flames was comparable and higher than that of the B75 and pure CME flame (Table 1); the flame temperatures were lowest in the B25 flame at midflame and at 75% flame heights. The soot volume fraction in the B25 flame was lower than that of the diesel flame (as discussed later). Thus, while less soot was formed in the B25 flame, the comparable carbon input rate resulted in higher CO formation. Also, the CO emission index decreased with an increase in Reynolds number for all the fuels due to the higher air entrainment that occurred as the Reynolds number was increased.

Global Emissions.
The global NOx emission index measured in the various flames is presented in Figure 9. It was highest for the pure CME flame; in the case of the flames of the blends, the CME B25 flame had a higher NOx emission index than that of the diesel flame. NOx production in flames is governed by the thermal, Fenimore, and N 2 O mechanisms [15]. Thus, the local temperatures and the fuel-air mixture ratios play a significant role in the generation of NOx. In general, the maximum temperature in the flames of CME and CME blends were higher than that of the diesel flame, as described below. Therefore, NOx generation by the thermal mechanism could be significant at these conditions, in contrast to laminar flames at the same initial equivalence ratio [11]. The NOx emission index was found to increase with an increase in Reynolds number for most cases, except for the B50 and B75 flames, in which the NOx emission  index dropped at the Reynolds number of 3600 and increased at a Reynolds number of 4500. The larger air entrainment at higher Reynolds numbers leading to leaning of the local mixtures, and the shorter residence times (due to higher velocities) affect the thermally formed NOx in opposite directions, and thus could be a reason for the minimum observed at the Reynolds number of 3600 for the B50 flame. Figure 10 shows the radial temperature profiles in the flames at three-quarter flame height. The effect of Reynolds number (Re = 2700, 3600, and 4500) can be seen by comparing the three figures. In each figure, the effect of blend ratio at a given Reynolds number is presented. At this axial location, the temperature reached a maximum near the centerline and decreased as the radial distance was increased for all flames except for the B75 Turbulence intensity (%)

Temperature Profiles.
Diesel Re 4500 CME Re 4500 SME Re 4500 (c) Figure 6: Axial velocity turbulence intensity profiles at the burner exit.
flame, which shows an off-axis peak that is particularly dominant at Re = 4500. This trend is characteristic of sootdominated heterogeneous combustion in the far-burner regions of high initial equivalence ratio premixed flames. At the highest Re (4500), the combined effect of soot reduction and cumulative air entrainment could make homogeneous gas-phase reactions significant, resulting in off-axis peaks in radial temperature profiles. That seems to be the case for the flame of the CME B75 blend.
Since the local equivalence ratio is determined by the soot concentration and cumulative air entrainment up to this region, it can approach stoichiometric or even lean condition. Hence, the maximum temperature that could occur is limited by the adiabatic flame temperature (2303 K for pure diesel, 2300 K for CME B25, 2297 K for CME B50, 2294 K for CME B75, and 2291 K for CME). The temperature in the pure CME flame was the highest at this location. At a Reynolds number of 2700, the temperatures in the diesel and B75 flame were comparable and higher than those of the B25 and B50 flames; the temperatures in the B25 and diesel flames were the lowest at the Reynolds numbers of 3600 and 4500. The peak temperature increased with Reynolds number due to better mixing between the surrounding air and the partially premixed fuel/air mixture.
4.6. Soot Volume Fraction. The axial variation of soot volume fraction in the flames is presented in Figure 11. In general, the soot volume fraction increased with downstream distance. In the flames of the CME blends, the soot volume fraction decreased near the flame tip, particularly at the CME CME B75 CME B25 Diesel Reynolds number of 4500. The diesel flame had the highest soot volume fraction, followed by the pure biofuel flame. The presence of oxygen in the fuel molecule of CME would be expected to reduce the formation of soot in the biofuel flames. It is interesting that the soot volume fraction in the flames of the CME blends was lower than that of the pure CME flames. The fuel with high carbon/hydrogen ratio (diesel) and oxygen in the fuel molecule that can enhance soot pyrolysis rate (CME) can yield high soot concentrations overall. Thus, the pure diesel and CME flames have higher soot concentrations than their blends. In the near-burner region (x/L < 0.35), the residence time of the fuel at Re = 2700 is much larger than that at Re = 4500. Fuel pyrolysis, soot nucleation, and particle growth are dominant in this region compared to oxidation reactions. Hence, at Re = 2700, the differences in the fuel molecular structure manifest themselves as significant differences in soot volume fraction. Due to higher velocities and entrainment of additional air that promote oxidation reactions also, these difference are mitigated at Re = 4500. The synergistic coupling effect between CME and diesel fuels appears to be dominant in these flames. The insignificant change in soot volume fraction with the increase in Reynolds number was a result of the cumulative overall effect of higher carbon input rate and higher dilution due to increased air entrainment. Figure 12 shows the inflame concentration measurement of NOx at three-quarter flame height. At this flame height, the peak NOx concentration occurred at the centerline and decreased towards the flame edge, although slight off-axis hump is seen for the B75 flame at Re = 2700. The peak NOx concentrations were recorded for the pure CME flames, which also had the highest in-flame temperatures [13]. Further, the Zeldovich thermal NOx mechanism shows that NOx formation rate also depends upon the local oxygen concentration The NOx concentrations were lowest for the pure diesel flames. Soot volume fraction CME Re 4500 CME B75 Re 4500 CME B50 Re 4500 CME B25 Re 4500 Diesel Re 4500 (c) Figure 11: Axial profiles of soot volume fraction. the length of the homogeneous gas-phase reactions (indicated by blue hue) was longer for the CME flames. The presence of oxygen in the fuel molecule of CME promoted the oxidation of the fuel in the nearburner region.

In-Flame NOx Concentration.
(2) The soot volume concentration in the flames of the CME blends was lower than that of the flames of the pure fuels, highlighting the combined effect of the reduction in soot precursor oxidation (due to the presence of oxygen in the biofuels) and the higher degree of pyrolysis of diesel. As the Reynolds number was increased, the soot volume concentration did not change significantly due to the cumulative overall effect of higher carbon input and higher air entrainment.  (3) The radiative heat fraction of the B75 flame was higher than that of the other flames. Whereas the soot content in the flames of the blends was lower than that of the pure fuel flames, the temperatures at mid flame and three-quarter flame heights were higher for the flames of the blends. The lower amount of soot appeared to burn at a higher rate, thereby increasing the energy release rate with the same amount of air entrained at a given Reynolds number.
(4) The CO emission index decreased with an increase in Reynolds number due to the increase in air entrainment.
(5) The global NOx emission index was highest for the pure CME flame. At three-quarter flame height, the NOx concentrations in the flames were lowest for the pure diesel flame, and increased with Reynolds number. These trends corresponded to those of