Abstract

The effects of turbulent Reynolds number on the statistical behaviour of the displacement speed have been studied using three-dimensional Direct Numerical Simulation of statistically planar turbulent premixed flames. The probability of finding negative values of the displacement speed is found to increase with increasing turbulent Reynolds number when the Damköhler number is held constant. It has been shown that the statistical behaviour of the Surface Density Function, and its strain rate and curvature dependence, plays a key role in determining the response of the different components of displacement speed. Increasing the turbulent Reynolds number is shown to reduce the strength of the correlations between tangential strain rate and dilatation rate with curvature, although the qualitative nature of the correlations remains unaffected. The dependence of displacement speed on strain rate and curvature is found to weaken with increasing turbulent Reynolds number when either Damköhler or Karlovitz number is held constant, but the qualitative nature of the correlation remains unaltered. The implications of turbulent Reynolds number effects in the context of Flame Surface Density (FSD) modelling have also been addressed, with emphasis on the influence of displacement speed on the curvature and propagation terms in the FSD balance equation.

1. Introduction

Displacement speed 𝑆𝑑 is an important quantity for the modelling of turbulent premixed flames using both level-set [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 [313] and on three-dimensional DNS with simplified chemistry [1422]. These studies have looked at the effects of local strain rate and curvature effect on 𝑆𝑑 [1324] and have included the influence of root-mean-square (rms) turbulent velocity fluctuation magnitude 𝑢 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𝑡 on the statistics of 𝑆𝑑 has been analysed based on a three-dimensional compressible DNS database of statistically planar turbulent premixed flames. The main objectives of this study are as follows:(1)to demonstrate the influence of Re𝑡 on the statistical distribution of 𝑆𝑑(2)to demonstrate the influence of Re𝑡 on the local curvature and strain rate dependencies of 𝑆𝑑(3)to indicate the implications of the dependence of 𝑆𝑑 on Re𝑡 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𝑡 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.

2. 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 [1422] are in agreement with those obtained from two-dimensional DNS with complex chemistry [313]. 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𝑐, which may be defined based on a suitable product mass fraction 𝑌𝑃 according to𝑌𝑐=𝑃𝑌𝑃0𝑌𝑃𝑌𝑃0,(1) 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𝜌𝜕𝑐𝜕𝑡+𝜌𝑢𝑗𝜕𝑐𝜕𝑥𝑗=̇𝜕𝑤+𝜕𝑥𝑗𝜌𝐷𝜕𝑐𝜕𝑥𝑗,(2a) where ̇𝑤 is the chemical reaction rate, 𝜌 is the fluid density, and 𝐷 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 𝑐=𝑐 [6, 7] 𝜌𝜕𝑐𝜕𝑡+𝜌𝑢𝑗𝜕𝑐𝜕𝑥𝑗=𝜌𝑆𝑑||||𝑐𝑐=𝑐.(2b)

Comparing (2a) and (2b), the expression for the displacement speed of a given 𝑐=𝑐 isosurface can be obtained as [1, 37]𝑆𝑑=̇𝑤+(𝜌𝐷𝑐)𝜌||||||||𝑐𝑐=𝑐.(3) It is often useful to decompose the displacement speed 𝑆𝑑 into three separate components [37]𝑆𝑑=𝑆𝑟+𝑆𝑛+𝑆𝑡,(4) where 𝑆𝑟, 𝑆𝑛, and 𝑆𝑡 are the reaction, normal diffusion, and tangential diffusion components which are given as 𝑆𝑟=̇𝑤𝜌||||||||𝑐𝑐=𝑐,𝑆𝑛=𝑁(𝜌𝐷𝑁𝑐)𝜌||||||||𝑐𝑐=𝑐,𝑆𝑡=2𝐷𝜅𝑚,(5) where 𝜅𝑚 is the arithmetic mean of the two principal curvatures, which is defined as 𝜅𝑚=12𝑁.(6) 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=|𝑐|) [15, 2527] and curvature 𝜅𝑚 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 𝑐 is given by [15, 2527] 𝜕||||𝑐+𝜕𝑢𝜕𝑡𝑗||||𝑐𝜕𝑥𝑗=𝑎𝑇||||𝑐+2𝑆𝑑𝜅𝑚||||𝜕𝑆𝑐𝑑𝑁𝑖||||𝑐𝜕𝑥𝑖,(7) where 𝑎𝑇 is the tangential strain rate given by𝑎𝑇=𝛿𝑖𝑗𝑁𝑖𝑁𝑗𝜕𝑢𝑖𝜕𝑥𝑗𝑐=𝑐.(8)

It is evident from (7) and (8) that tangential strain rate is likely to affect the statistical behaviour of |𝑐|, which in turn can affect the local curvature and strain rate dependence of 𝑆𝑟 and 𝑆𝑛. It is worth noting that in some previous studies [312], the statistics of the density-weighted displacement speed 𝑆𝑑=𝜌𝑆𝑑/𝜌0 were analysed instead of the statistics of 𝑆𝑑 itself. For low Mach number unity Lewis number flames, 𝑆𝑑 and 𝑆𝑑 are related by 𝑆𝑑=𝑆𝑑(1+𝜏𝑐) for a given isosurface of c, and thus the statistics of 𝑆𝑑 are similar to those of 𝑆𝑑. For this reason, and also since 𝑆𝑑 appears explicitly in the SDF and FSD transport equations [2, 15, 2527], the statistics of 𝑆𝑑 are addressed in the present work. The local strain rate and curvature dependence of 𝑆𝑑 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𝑡=𝜌0𝑢𝑙/𝜇, where 𝜌0 is the density and 𝜇 is the dynamic viscosity of the unburned reactants, and 𝑙 is the integral length scale of the turbulence. For unity Lewis number flames, the turbulent Reynolds number can be scaled as [1]Re𝑡Da2Ka2𝑢𝑆𝐿2𝑢Da/𝑆𝐿4Ka2𝑙𝛿th4/3Ka2/3,(9) where Da=𝑙𝑆𝐿/𝑢𝛿th is the Damköhler number and Ka=(𝑢/𝑆𝐿)3/2/(𝑙/𝛿th)1/2is the Karlovitz number. Here, 𝛿th=(𝑇ad𝑇0)/Max|𝑇|𝐿 is the thermal flame thickness with 𝑇ad, 𝑇0 and 𝑇 being the adiabatic flame temperature, unburned gas temperature, and instantaneous dimensional temperature, respectively, while the subscript “𝐿” 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𝑡. Thus, the influence of Re𝑡 on the statistics of 𝑆𝑑 is essentially induced by changes in Da and/or Ka.

3. 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 𝑥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 𝐶𝑃 and 𝐶𝑉 and the density-weighted mass diffusivity 𝜌𝐷 are taken to be constant and independent of temperature following several previous studies [1118, 3234]. For the present study, standard values are taken for Prandtl number (Pr = 0.7) and the ratio of specific heats (𝛾=𝐶𝑃/𝐶𝑉=1.4). The flame Mach number is taken to be Ma=𝑆𝐿/𝛾𝑅𝑇0=0.014. The Lewis number Le is taken to be equal to unity for all the cases considered here. As Da and Ka are often used to characterise the combustion situation on combustion regime diagram, the variation of Re𝑡 in the present study is brought about by modifying Da and Ka independently of each other. In order to achieve this, the values of both 𝑢/𝑆𝐿 and 𝑙/𝛿th are modified as both Da=𝑙𝑆𝐿/𝛿th𝑢 and Ka=(𝑢/𝑆𝐿)3/2(𝑙/𝛿th)1/2 are dependent on the values of 𝑢/𝑆𝐿 and 𝑙/𝛿th. The initial values of normalised rms turbulent velocity fluctuation 𝑢/𝑆𝐿, integral length scale to flame thickness ratio 𝑙/𝛿th, heat release parameter 𝜏=(𝑇ad𝑇0)/𝑇0, Damköhler number Da, Karlovitz number Ka and turbulent Reynolds number Re𝑡 are listed in 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 𝑢/𝑆𝐿 and 𝑙/𝛿th were chosen in such a manner as to vary the turbulent Reynolds number Re𝑡 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, 𝑢/𝑆𝐿 and 𝑙/𝛿th were chosen to vary Re𝑡 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 [322] with comparable values of Re𝑡 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𝑡 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𝑡 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 𝑡sim=Max(𝑡𝑓,𝑡𝑐), where 𝑡𝑓=𝑙/𝑢 is the initial eddy turn over time and 𝑡𝑐=𝛿th/𝑆𝐿 is the chemical time scale. All the cases considered here were run for one chemical time scale 𝑡𝑐, which corresponds to a time equal to 2.0𝑡𝑓 in case D, 3.0𝑡𝑓 in cases A, C, and E, and 4.34𝑡𝑓 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 [321, 2327]. 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 𝑢/𝑆𝐿, 𝑙/𝛿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𝑡 for all cases remains greater than unity (i.e., Re𝑡>1) for the values of 𝑢/𝑆𝐿 and 𝑙/𝛿th at the time when statistics were extracted, which suggests that all the flames remain turbulent according to the regime diagram by Peters [1].

4. Results and Discussion

4.1. Flame-Turbulence Interaction

The contours of the reaction progress variable 𝑐 in the x1-x2 mid-plane of the domain at time 𝑡sim=1.0𝛿th/𝑆𝐿 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 𝑐 isosurfaces increases with increasing 𝑢, that is, in going from case A through to case E (see Table 1). Moreover, the contours of 𝑐 representing the preheat zone (i.e., 𝑐<0.5) are much more distorted than those representing the reaction zone (i.e., 0.7<𝑐<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.

4.2. Probability Density Functions of Displacement Speed 𝑆𝑑

The probability density functions (pdfs) of normalised displacement speed 𝑆𝑑/𝑆𝐿 on the 𝑐=0.1,0.3,0.5,0.7, and 0.9 isosurfaces are shown in Figures 2(a)2(e), respectively. It is evident in all cases that there is a small but finite probability of finding negative displacement speed although the mean value of 𝑆𝑑/𝑆𝐿 remains strongly positive throughout the flame brush. Negative values of displacement speed suggest that the progress variable isosurface in question is moving locally in a direction opposite to that of mean flame propagation. This situation arises when the negative contribution of the molecular diffusion rate dominates over the positive semi-definite 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 𝑆𝑑/𝑆𝐿 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𝑆𝑟+𝑆𝑛𝑣𝜂𝑆𝐿𝑣𝜂1𝑂,𝑆Ka𝑡𝑣𝜂2𝐷𝜅𝑚𝑣𝜂𝑂(1),(10) where 𝑣𝜂 is the Kolmogorov velocity scale. Here, in accordance with Peters [1], (𝑆𝑟+𝑆𝑛) is taken to scale with 𝑆𝐿 and the curvature 𝜅𝑚 is taken to scale with the Kolmogorov length scale 1/𝜂. These scalings of (𝑆𝑟+𝑆𝑛) and 𝜅𝑚 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 (𝑆𝑟+𝑆𝑛) are likely to weaken progressively in comparison to the contribution of 𝑆𝑡 with increasing Karlovitz number. This suggests that for large values of Karlovitz number the negative contribution of 𝑆𝑡 more readily overcomes the predominantly positive contribution of (𝑆𝑟+𝑆𝑛) to yield a negative value of 𝑆𝑑. This eventually gives rise to an increased probability of finding negative values of 𝑆𝑑 for flames with increasing value of Karlovitz number, as suggested by Figure 2. In order to demonstrate that the observed negative values of 𝑆𝑑 originates principally due to 𝑆𝑡, the pdfs of (𝑆𝑟+𝑆𝑛)/𝑆𝐿 are shown in Figure 3, respectively, for all cases. Comparing Figures 2 and 3 reveals that the probability of finding negative values of (𝑆𝑟+𝑆𝑛) is significantly smaller than that of finding negative values of 𝑆𝑑.

Comparing Figures 2(a)2(e) and Figures 3(a)3(e), it can be seen that the most probable values of 𝑆𝑑/𝑆𝐿 and (𝑆𝑟+𝑆𝑛)/𝑆𝐿 remain of the order of 𝜌0/𝜌 and that the widths of the displacement speed pdf for a given 𝑐 isosurface increases with increasing 𝑢 (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 𝑆𝑡/𝑆𝐿 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 𝑢 results in a greater spread of flame curvature, which in turn increases the width of the 𝑆𝑡/𝑆𝐿 pdfs according to (5). The negligible mean value of 𝑆𝑡 essentially leads to the mean values of 𝑆𝑑 and (𝑆𝑟+𝑆𝑛) becoming almost equal to each other, but the statistical variation of 𝑆𝑑 induced by 𝑆𝑡 leads to broader pdfs of 𝑆𝑑 than those of (𝑆𝑟+𝑆𝑛). These results clearly suggest that the distributions of 𝑆𝑑 and its components 𝑆𝑟,𝑆𝑛 and 𝑆𝑡 generally broaden with 𝑢/𝑆𝐿, which essentially leads to broadening of the distributions with increasing Re𝑡(𝑢/𝑆𝐿)2Da(𝑢/𝑆𝐿)4/Ka2(𝑙/𝛿th)4/3Ka2/3 when either Da or Ka is held constant.

4.3. Strain Rate and Curvature Dependences of |𝑐|

For low Mach number unity Lewis number flames, both density 𝜌 and nondimensional temperature 𝑇 remain uniform on a given isosurface of 𝑐 and thus the statistical behaviours of 𝑆𝑟 and 𝑆𝑛 (see (5)) are mainly affected by the dependence of the SDF |𝑐| on strain rate and curvature. The values of |𝑐|×𝛿th conditional on normalised tangential strain rate 𝑎𝑇×𝛿th/𝑆𝐿 are shown in Figure 5(a) for the isosurface at 𝑐=0.8 (close to the location of maximum reaction rate), which suggests clearly that |𝑐| and 𝑎𝑇 are positively correlated for all cases. The correlation coefficients between |𝑐| and 𝑎𝑇 are shown in Figure 5(b) for values of 𝑐 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 𝑢 can be expressed as 𝑢=𝑎𝑇+𝑎𝑛.(11a) For unity Lewis number flames, the dilatation rate can be scaled as 𝑢𝜏𝑆𝐿/𝛿th [21], whereas 𝑎𝑇 can be scaled as 𝑎𝑇𝑢/𝑙 [37], which suggests that𝑢𝑎𝑇𝑆𝜏𝐿𝑙𝛿th𝑢𝜏Da.(11b)

Equation (11b) essentially shows that for low values of the Damköhler number, 𝑎𝑇 remains much greater than 𝑢. 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 𝑎𝑛𝑎𝑇 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 |𝑐|. The decreasing trend of normal strain rate with increasing tangential strain rate is reflected in the positive correlation between |𝑐| and 𝑎𝑇 throughout the flame brush. It can be seen from Figure 5 that the correlation between |𝑐| and 𝑎𝑇 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 |𝑐|×𝛿th conditional on normalised curvature 𝜅𝑚×𝛿th is shown in Figure 6(a) for the isosurface at 𝑐=0.8. Both positive and negative correlation trends can be observed, which is consistent with several previous studies [1422, 26, 27]. As a result of these opposing trends, the net correlation between |𝑐|×𝛿th and 𝜅𝑚×𝛿th turns out be weak throughout the flame brush, as evident from the correlation coefficients between |𝑐| and 𝜅𝑚 presented for different values of 𝑐 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 𝑢 is negatively correlated with curvature 𝜅𝑚, because of focusing (defocusing) of heat in the negatively (positively) curved regions of the flame surface [18, 21]. The tangential strain rate 𝑎𝑇 is also found to be negatively correlated with curvature [18, 21]. The correlation coefficients for the 𝑢𝜅𝑚 and 𝑎𝑇𝜅𝑚 correlations for different 𝑐 values across the flame brush are shown in Figures 7(a) and 7(b), respectively. Comparing the correlation coefficients of the 𝑢𝜅𝑚 and 𝑎𝑇𝜅𝑚 correlations for different cases reveals that the strength of both of these negative correlations decreases with increasing 𝑢/𝑆𝐿 and Re𝑡(𝑙/𝛿th)4/3Ka2/3 when either Da or Ka is held constant.

The combination of the positive correlation between |𝑐| and 𝑎𝑇 and the negative correlation between 𝑎𝑇 and 𝜅𝑚 leads to small values of |𝑐| in regions of high positive curvature. This gives rise to the negative correlation trend between |𝑐| and 𝜅𝑚, as observed in Figure 6. However, the dilatation rate 𝑢 may locally attain large values in the highly negatively curved regions of the flame due to strong focussing of heat. Hence, 𝑢 may locally exceed 𝑎𝑇 and induce an increase in 𝑎𝑛 according to (11a). As isoscalar lines move apart from each other under increasing normal strain rate, the increase in 𝑎𝑛 with increasing negative curvature leads to a positive correlation between |𝑐| and 𝜅𝑚 in negatively curved regions, as also observed in Figure 6.

4.4. Curvature Dependence of 𝑆𝑑

The mean normalised displacement speed 𝑆𝑑/𝑆𝐿 conditional on normalised curvature 𝜅𝑚×𝛿th is shown in Figure 8(a) for 𝑐=0.8, which indicates that 𝑆𝑑/𝑆𝐿 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 [311, 1318, 2022]. The variation of the correlation coefficient between 𝑆𝑑 and 𝜅𝑚 for different 𝑐 isosurfaces across the flame brush is shown in Figure 8(b), which reveals that strength of the negative correlation decreases with increasing Re𝑡(𝑢/𝑆𝐿)2Da(𝑢/𝑆𝐿)4/Ka2(𝑙/𝛿th)4/3Ka2/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 𝑆𝑡 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 𝑆𝑑 and 𝜅𝑚 originates due to the curvature dependence of 𝑆𝑟 and 𝑆𝑛.

Mean values of the reaction component 𝑆𝑟/𝑆𝐿 conditional on normalised curvature 𝜅𝑚×𝛿th are shown in Figure 9(a) for the isosurface at 𝑐=0.8. The plot shows both positive and negative correlation trends. For low Mach number unity Lewis number flames, the curvature dependence of 𝑆𝑟 is principally determined by the curvature dependence of |𝑐|. The negative (positive) correlation between |𝑐| and 𝜅𝑚 ultimately gives rise to the positive (negative) correlation between 𝑆𝑟 and 𝜅𝑚 as observed in Figure 9. As a result of the weak net correlation between |𝑐| and 𝜅𝑚, the net correlation between 𝑆𝑟 and 𝜅𝑚 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 𝑐 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 𝑆𝑛/𝑆𝐿 conditional on normalised curvature 𝜅𝑚×𝛿th are shown in Figure 10(a) for 𝑐=0.8, while the correlation coefficients between 𝑆𝑛 and 𝜅𝑚 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 |𝑁(𝜌𝐷𝑁𝑐)|𝜌𝐷/𝛿2 where 𝛿 is the characteristic flame thickness which in turn can be scaled as |𝑐|1/𝛿. This suggests that 𝑆𝑛 scales as 𝐷/𝛿 towards the reactant side of the flame brush and as 𝐷/𝛿 towards the product side. As a result of this, the positively (negatively) correlating branch of the |𝑐|𝜅𝑚 correlation gives rise to a positive (negative) correlation between 𝑆𝑛 and 𝜅𝑚 towards the reactant side, where 𝑆𝑛 assumes predominantly positive values. By contrast, the positive (negative) correlation between |𝑐| and 𝜅𝑚 results in a negative (positive) correlation between 𝑆𝑛 and 𝜅𝑚 towards the product side where 𝑆𝑛 assumes predominantly negative values. Both positive and negative correlation trends can be discerned for all the 𝑐 isosurfaces shown in Figure 10(b). As a result of this, the net correlation between 𝑆𝑛 and 𝜅𝑚 turns out to be weak throughout the flame brush. The positive and negative correlating trends with 𝜅𝑚 observed for both 𝑆𝑟 and 𝑆𝑛 with 𝜅𝑚 give rise to positive and negative correlating branches in the variation of (𝑆𝑟+𝑆𝑛)conditional on 𝜅𝑚×𝛿th as shown in Figure 10(c) for 𝑐=0.8. This nonlinear correlation between (𝑆𝑟+𝑆𝑛) and 𝜅𝑚 ultimately induces a nonlinear curvature dependence of displacement speed 𝑆𝑑 (see Figure 8(a)), noting that the correlation coefficient arising from the deterministic relation between 𝑆𝑡 and 𝜅𝑚 remains close to −1.0 throughout the flame brush for all cases. The correlation coefficients between (𝑆𝑟+𝑆𝑛) and 𝜅𝑚 are shown in Figure 10(d) for values of 𝑐 from 0.1 to 0.9, and demonstrate that the (𝑆𝑟+𝑆𝑛)𝜅𝑚 correlation is much weaker than the 𝑆𝑡𝜅𝑚 correlation in all cases. Hence, the net correlation between 𝑆𝑑 and 𝜅𝑚 turns out to be negative throughout the flame brush for all the cases considered here (see Figure 8(a)).

4.5. Tangential Strain Rate Dependence of 𝑆𝑑

The mean values of normalised displacement speed 𝑆𝑑/𝑆𝐿 conditional on normalised tangential strain rate 𝑎𝑇×𝛿th/𝑆𝐿 are shown in Figure 11(a) for 𝑐=0.8 which suggests that 𝑆𝑑 and 𝑎𝑇 are positively correlated for all cases. The observed strain rate dependence of 𝑆𝑑 is found to be qualitatively consistent with previous DNS studies for flames with global Lewis number close to unity [3, 5, 12, 14, 15, 1719, 21, 22]. The variations of the correlation coefficients between 𝑆𝑑 and 𝑎𝑇 for different 𝑐 isosurfaces are shown in Figure 11(b), which reveals that the strength of the correlation between 𝑆𝑑 and 𝑎𝑇 decreases with increasing 𝑢/𝑆𝐿 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 𝑆𝑑 weakens with increasing frequency of straining. In order to explain the observed behaviour it is useful to look into the dependence of 𝑆𝑟, 𝑆𝑛, and 𝑆𝑡 on 𝑎𝑇. For low Mach number unity, Lewis number flames the strain rate dependence of 𝑆𝑟 is essentially determined by the strain rate dependence of |𝑐|. It has been shown already that |𝑐| and 𝑎𝑇 are positively correlated (see Figure 5), which gives rise to a negative correlation between 𝑆𝑟 and 𝑎𝑇 according to (5). This is substantiated by the correlation coefficients between 𝑆𝑟 and 𝑎𝑇 as shown in Figure 11(c) for the isosurfaces from c = 0.5 to 0.9. The correlation between 𝑆𝑟 and 𝑎𝑇 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 |𝑐|𝑎𝑇 correlation.

It has been discussed already that 𝑆𝑛 scales as 𝑆𝑛𝐷/𝛿 on the reactant side of the flame brush, and as 𝑆𝑛𝐷/𝛿 on the product side. Since there is a positive correlation between |𝑐|1/𝛿 and 𝑎𝑇, this gives rise to a positive (negative) correlation between 𝑆𝑛 and 𝑎𝑇 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 𝑆𝑛 and 𝑎𝑇 does not show any monotonic trend with changes in either Da or Ka. However, it has been found that the correlation between 𝑆𝑛 and 𝑎𝑇 in general remains stronger for smaller values of 𝑢 when either Da or Ka is held constant.

Since 𝑆𝑟 remains negligible towards the reactant side of the flame brush and 𝑆𝑛 is positively correlated with 𝑎𝑇, the correlation coefficient between (𝑆𝑟+𝑆𝑛) and 𝑎𝑇 remains positive towards the reactant side (see Figure 11(e)). By contrast, both 𝑆𝑟 and 𝑆𝑛 are negatively correlated with 𝑎𝑇 towards the product side, which leads to a locally negative correlation between (𝑆𝑟+𝑆𝑛) and 𝑎𝑇.

Since 𝑎𝑇 and 𝜅𝑚 are negatively correlated, the tangential diffusion component 𝑆𝑡 is positively correlated with 𝑎𝑇 because 𝑆𝑡 is proportional to the negative of curvature (i.e., 𝑆𝑡=2𝐷𝜅𝑚). This positive correlation is evident from Figure 11(f), in which the correlation coefficient between 𝑆𝑡 and 𝑎𝑇 is shown for isosurfaces of 𝑐 across the flame brush. The negative correlation between 𝑎𝑇 and 𝜅𝑚 (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 𝑎𝑇 and 𝜅𝑚 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 𝑆𝑡 and 𝑎𝑇 (see Figure 11(f)). This essentially suggests that the strength of the positive correlation between 𝑆𝑡 and 𝑎𝑇 decreases with increasing turbulent Reynolds number Re𝑡 when either Ka or Da is held constant. The effects of the positive 𝑆𝑡𝑎𝑇 and negative 𝑎𝑇𝜅𝑚 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 (𝑆𝑟+𝑆𝑛) and 𝑎𝑇, and between 𝑆𝑡 and 𝑎𝑇, aid each other towards the reactant side of the flame brush. By contrast, towards the product side the positive correlation between 𝑆𝑡 and 𝑎𝑇 overcomes the negative correlation between (𝑆𝑟+𝑆𝑛) and 𝑎𝑇 to result in a net positive correlation between 𝑆𝑑 and 𝑎𝑇 in all the cases considered here (see Figure 11(a)). The weakening of the correlation between 𝑆𝑡 and 𝑎𝑇 with increasing Re𝑡 ultimately leads to decrease in the strength of the correlation between 𝑆𝑑 and 𝑎𝑇 when either Da or Ka is held constant (see Figures 11(a) and 11(f)).

4.6. Implications for FSD Modelling

In the context of Reynolds Averaged Navier Stokes (RANS) modelling, the generalised FSD Σgen is defined as [39] Σgen=||||𝑐,(12a) where 𝑄 indicates the Reynolds averaged value of a general quantity 𝑄. It has been demonstrated in [15, 23, 26, 27] that the statistical behaviour of |𝑐| and the terms of its transport (7) in the thin reaction zones regime depend significantly on the choice of progress variable isosurface 𝑐=𝑐. 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+𝜕𝑢𝜕𝑡𝑗𝑠Σgen𝜕𝑥𝑗=(𝛿𝑖𝑗𝑁𝑖𝑁𝑗)𝜕𝑢𝑖𝜕𝑥𝑗𝑠Σgen+2𝑆𝑑𝜅𝑚𝑠Σgen𝜕𝑆𝑑𝑁𝑖𝑠Σgen𝜕𝑥𝑖,(12b) where 𝑄𝑠=𝑄|𝑐|/|𝑐| indicates a surface averaging operation for a general quantity 𝑄 [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 as2𝑆𝑑𝜅𝑚𝑠Σgen𝑆=2𝑟+𝑆𝑛𝜅𝑚𝑠Σgen4𝐷𝜅2𝑚𝑠Σgen,𝑆(12c)𝑑𝑁𝑠Σgen𝑆=𝑟+𝑆𝑛𝑁𝑠Σgen+2𝐷𝜅𝑚𝑁𝑠Σgen.(12d)

The variations of all of the terms appearing in these equations with Favre averaged reaction progress variable ̃𝑐=𝜌𝑐/𝜌 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 𝑆𝑑𝑁𝑠Σ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(𝑆𝑟+𝑆𝑛)𝜅𝑚𝑠Σgen (line with filled triangles) remains comparable to the contribution of the term 4𝐷𝜅2𝑚𝑠Σgen (line with open squares) for small values of Re𝑡 (i.e., cases A and B). However, the contribution of 2(𝑆𝑟+𝑆𝑛)𝜅𝑚𝑠Σgen remains smaller than that of 4𝐷𝜅2𝑚𝑠Σgen for larger values of Re𝑡 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𝐷𝜅2𝑚𝑠Σgen) becomes increasingly more important than the magnitude of 2(𝑆𝑟+𝑆𝑛)𝜅𝑚𝑠Σgen with increasing Karlovitz number (Damköhler number) when Da (Ka) is held constant. The statistical behaviour of 2(𝑆𝑟+𝑆𝑛)𝜅𝑚𝑠Σgen is determined by the correlation between |𝑐| and 𝜅𝑚 and between (𝑆𝑟+𝑆𝑛) and 𝜅𝑚 and the strength of both of these correlations tends to weaken with increasing (increasing) Karlovitz number (Damköhler number) when Da (Ka) is held constant (see Figures 6(b) and 10(d)), causing a corresponding reduction in the relative magnitude of the term 2(𝑆𝑟+𝑆𝑛)𝜅𝑚𝑠Σgen.

The FSD propagation term 𝑆𝑑𝑁𝑠Σgen (line with crosses) remains positive towards the reactant side, but it eventually becomes negative towards the product side. It can be seen that the term (𝑆𝑟+𝑆𝑛)𝑁𝑠Σgen (line with filled squares) remains almost equal to the FSD propagation term itself, because the contribution of 2𝐷𝜅𝑚𝑁𝑠Σgen (line with open circles) remains negligible throughout the flame brush due to the statistically planar nature of the flames considered here. The contribution of 2(𝑆𝑟+𝑆𝑛)𝜅𝑚𝑠Σgen remains positive towards the reactant side and assumes negative values towards the product side.

In some previous studies, the FSD curvature and propagation terms have been modelled by replacing 𝑆𝑑 with 𝑆𝐿 to yield the modelled expressions 2𝑆𝐿𝜅𝑚𝑠Σgen and 𝑆𝐿𝑁𝑠Σ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 𝑆𝑑 cannot be approximated simply by 𝑆𝐿 in the modelling of the FSD curvature and propagation terms. Moreover, Figure 12 suggests that modelling of 𝜅2𝑚𝑠 is necessary for the closure of the FSD curvature term especially in the thin reaction zones regime since the term 4𝐷𝜅2𝑚𝑠Σgen remains close to 2𝑆𝑑𝜅𝑚𝑠Σ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.

5. 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𝑡 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 𝑢/𝑆𝐿 and 𝑙/𝛿th, the variation of Re𝑡 for the present database has been achieved by modifying the initial values of 𝑢/𝑆𝐿 and 𝑙/𝛿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 |𝑐|, 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 [37, 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.

Acknowledgment

N. Chakraborty gratefully acknowledges the financial assistance of the EPSRC, UK.