Statistics of Conditional Fluid Velocity in the Corrugated Flamelets Regime of Turbulent Premixed Combustion: A Direct Numerical Simulation Study

The statistics of mean fluid velocity components conditional in unburned reactants and fully burned products in the context of Reynolds Averaged Navier Stokes (RANS) simulations have been studied using a Direct Numerical Simulation database of statistically planar turbulent premixed flame representing the corrugated flamelets regime combustion. Expressions for conditional mean velocity and conditional velocity correlations which are derived based on a presumed bimodal probability density function of reaction progress variable for unity Lewis number flames are assessed in this study with respect to the corresponding quantities extracted from DNS data. In particular, conditional surface averaged velocities (ui)Rs and the velocity correlations (uiuj)Rs in the unburned reactants are demonstrated to be effectively modelled by the unconditional velocities (ui)R and velocity correlations (uiuj)R, respectively, for the major part of turbulent flame brush with the exception of the leading edge. By contrast, conditional surface averaged velocities (ui)Ps and the velocity correlations (uiuj)Ps in fully burned products are shown to be markedly different from the unconditional velocities (ui)P and velocity correlations (uiuj)P , respectively.


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) 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 u, c; x = α c x P R u; 0; x δ(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 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 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] ( 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 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(a) reveals that the contours 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), 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 c 2 remains smaller than c(1 − c), which demonstrates that O(γ c ) contribution is small but not negligible in this case.This behaviour is expected as the reaction progress variable contours representing the leading edge of the preheat zone are affected by the turbulent eddies and introduce nonnegligible O(γ c ) contributions.
The variation of (u 1 ) R /S L with c is shown in Figure 2(a), which demonstrates that (u 1 ) R /S L increases from unburned gas side to the burned gas side of the flame brush.
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 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 of L .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 term 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 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 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(a) that (10b) overpredicts 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 )/[(∂c/∂x 1 )S L ] and [(u 1 ) R − (u 1 ) Rs ]/S L remain 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 L with c across the flame brush along with the predictions of (12b) according to c = c (shown as Model 1 in the legend) and c = c (shown as Model 2 in the legend).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 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 ) 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.

Figure 1 :
Figure 1: (a) Contours of reaction progress variable c for c = 0.1-0.9(left to right) in steps of 0.1 at the x 1 -x 2 mid-plane.(b) Variations of c 2 with c across the flame brush.(c) Pdfs of c at the location corresponding to c = 0.5.

Figure 2 :
Figure 2: (a) Variations (u 1 ) R /S L and u 1 /S L with c across the flame brush along with the predictions of (8a).(b) Variations (u 1 ) P /S L and u 1 /S L with c across the flame brush along with the predictions of (8b).Variations of (c) [(u 1 ) P − (u 1 ) R ]/S L and (d) ρu 1 c × ∂ c/∂x 1 × δ th /ρ 0 S L with c across the flame brush.

Figure 3 :
Figure 3: (a) Variations of (u 1 u 1 ) R /S 2L 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).

2 (u 2 u 2 ) 2 L(u 2 u 2 ) R /S 2 LFigure 6 :
Figure 6: (a) 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 across the flame brush along with the predictions of (12b) according to c = c (shown as Model 1 in the legend) and c = c (shown as Model 2 in the legend).(b) Variations of (u 2 u 2 ) Rs /S 2 L and (u 2u 2 ) R /S 2 L = [(u 2 ) R (u 2 ) R + (u 2 u 2 ) R ]/S 2L with c across the flame brush along with the predictions of (12b) according to c = c (shown as Model 1 in the legend) and c = c (shown as Model 2 in the legend).

Table 1 :
List of initial values of turbulence and combustion parameters for the present DNS database.