Stretch Rate and Displacement Speed Correlations for Increasingly-Turbulent Premixed Flames

Probability density functions of the components of stretch rate are investigated using a previously-published Direct Numerical Simulation dataset spanning a range of turbulence intensities in the Thin Reaction Zones (TRZ) regime. The dataset was generated by varying the turbulence intensity across five different simulations while maintaining fixed the remaining physico-chemical input parameters such as integral length scale and laminar flame thickness and speed. Across the entire dataset, the joint probability density function of stretch rate and displacement speed displays a distinctive shape with two branches consistent with previous studies at low turbulence intensities. This joint probability density function is analysed further by extracting individual contributions of stretch rate components to determine their relative importance across the branches. The curvature dependence of displacement speed appears to play an important role in shaping these branches. Implications of this result with regard to evaluation of the components of stretch rate in the TRZ regime are discussed.


Flamelet modelling: stretch rate and the displacement speed
In flamelet modelling of turbulent premixed flames, the rate at which reactants are transformed into products is a function of the mean stretch rate ṡ of the flame surface [1]. The instantaneous stretch rateṡ of the flame surface is defined as the fractional rate of change of its area A [2] as followsṡ with the flame surface itself represented by an isosurface of the reaction progress variable c which is defined as where Y is the local mass fraction of a suitable reacting species, the subscripts R and P denoting its value in the reactant and product mixtures, respectively. At any point on a progress variable isosurface c = c * , the local instantaneous stretch rateṡ is written as [3,4] a linear combinationṡ = a t + 2s d h m , of its components a t and s d h m that represent individually: 1. the tangential strain rate a t imposed by the flow field u on the local tangent plane to the flame surface as given by where n ≡ −∇c/|∇c| is the local unit normal vector defined to be positive in the direction of the reactant mixture, ∇ is the gradient operator and | · | is the Euclidean norm, and 2. the propagation rate s d h m of the flame surface which operates in regions of non-zero mean curvature h m ≡ −(1/2)∇ · n, (5) by virtue of its local displacement speed whereω is the reaction rate, ρ is the local mixture density and D is the molecular diffusivity of the same chemical species that defines the progress variable c.
Given that the mean stretch rate ṡ governs the overall transformation rate of reactants as mentioned at the outset, the key task of flamelet modelling is the estimation of ṡ using approximations for the individual stretch-rate components a t , h m , and s d . These components, however, are strongly coupled and the underlying inter-dependences are highly non-linear. For example, instantaneous as well as persistent strain a t can affect both displacement speed s d [5] and flame curvature h m [6]. In turn, the flame curvature has a strong influence on the displacement speed [7,8]. Subsequently, the resulting displacement speed may enhance or inhibit the flame curvature further [9]. Finally, the propagation of the flame surface changes the strain field experienced by it [10]. These effects highlight the recursive nature of the relationship between instantaneous stretch rateṡ and its components, particularly that betweenṡ and s d . This renders flamelet modelling and the associated task of estimating the mean stretch rate ṡ challenging.
In order to simplify the estimation procedure for ṡ , the expected dependences betweeṅ s and its components are usually delineated into regimes. These regimes are based on the ratio u /s L of turbulence intensity u and the unstrained laminar flame speed s L , and the ratio 0 /δ L of integral length scale 0 and the unstrained laminar flame thickness δ L [2]. In the low turbulence intensity regime where u /s L ∼ 1, approximate expressions have been developed for s d as a function of the stretch rateṡ [11]. However, such expressions are generally not applicable to higher intensity regimes where u /s L 1 [12]. Consequently, the question remains: what relationships hold between the stretch rate components for increasing values of u /s L ?

Thin reaction zones regime: scope of present work
In recent years, Direct Numerical Simulation (DNS) has provided evidence of correlations between stretch rate components for moderate-and high-intensity turbulence [13][14][15][16]. In these investigations, the Karlovitz number, a dimensional group defined as based on the two ratios mentioned above, has been used for delineating the various possible relationships betweenṡ and its components. The regime addressed by these recent simulations, known as the Thin Reaction Zones (TRZ) regime, corresponds to a Karlovitz number 1 < Ka < 100. In this regime, the dependence of displacement speed s d on curvature h m has been observed to be significant towards estimation of the overall stretch rate [13,15]. Furthermore, the effects of h m on s d have been found to be regulated by the Lewis number Le, the ratio of the thermal diffusivity to the molecular diffusivity of the reacting mixture [17]. In multi-species mixtures, these effects may be exaggerated due to differential diffusion [18]. As a result, the overall relationship between s d andṡ in the TRZ regime is necessarily complicated and differs significantly from the expressions [19] derived for low-intensity regimes.
In this context, a relatively early 2D DNS study of methane-air combustion with detailed chemistry [9] has elaborated upon a distinctive branched shape of the distribution betweeṅ s and s d in the low to moderate intensity region of the TRZ regime. This shape was also seen, but not discussed, in a more recent 3D DNS study with single-step chemistry at similar turbulence intensities [20]. Most recently, the branched shape has also been observed in 3D DNS of planar turbulent methane-air flames using detailed chemistry at moderate turbulence intensities [21]. The shape of this distribution has implications with regard to estimation of the propagation term in the mean stretch rate ṡ . However, since these simulations have considered a relatively narrow range of Ka values, it remains to be confirmed if the branched distribution holds for a broader range of u /s L values for a given value of 0 /δ L . If it holds, then what factors determine such a branched shape? The present study takes up these open questions, the objective being to discern the relationship between stretch rate components as observed for increasing values of u /s L in the TRZ regime. Correlations between stretch rate components are investigated within a 3D DNS dataset [22,23] generated by varying u systematically from 1.5s L ≤ u ≤ 30s L to span vertically the entire TRZ regime (1 < Ka < 100) by maintaining fixed the ratio 0 /δ L as shown in Fig. 1.
In the following section, the DNS dataset investigated here is described briefly. In Section 3, stretch rate component correlations are examined, and these correlations are discussed in the context of modelling in Section 4.1.

Direct Numerical Simulation Dataset
The dataset investigated here consists of five separate DNS cases obtained by solving the Navier-Stokes equations in 3D compressible form along with equations for the transport of a reacting species. In each simulation, homogeneous isotropic turbulence of desired intensity is initialized in an inflow-outflow configuration. An identical copy of the turbulent flow field is convected through the inlet at a mean axial velocity equal in magnitude to the unstrained laminar flame speed so as to maintain the turbulence ahead of the flame. A statisticallyplanar flame brush is allowed to propagate freely towards the inlet. At any given snapshot in time, the instantaneous flame brush defined as the spatial region corresponding to 0.05 < c < 0.95 remains far from the inlet and outlet boundaries.
Flame chemistry is represented using a single-step Arrhenius reaction mechanism with the Lewis number of the reacting species set equal to unity. The laminar flame speed s L = 0.39 ms −1 is computed as the consumption speed of an unstrained 1D laminar premixed flame, and is intended to be representative of stoichiometric methane-air combustion. Similarly, the laminar flame thickness δ L = 0.36 mm based on the maximum temperature gradient for the unstrained 1D laminar premixed flame. Further details are provided in the references [22,23]. The key feature of the dataset is that all parameters pertaining to the flame thermochemistry are maintained fixed as well as maintaining fixed turbulent flow parameters such as the integral length scale while varying only the turbulence intensity u . The value of u at the inlet ranges from as low as 1.5s L to as high as 30s L as reported in Table 1. The turbulence decays in intensity slightly from the inlet up to the leading edge of the flame brush. The turbulence intensity u | le /s L computed at the leading edge is presented alongside in Table 1. In accordance with earlier work [22,23], it is the inlet turbulence intensity u that is employed throughout the analysis and presentation of results. The integral length scale 0 remains nearly constant from the inlet up to the leading edge of the flame brush. The values for the integral length scale 0 and the Taylor length scale λ computed for the (initial) Batchelor-Townsend energy spectrum are reported in Table 1 [22,23]. The Karlovitz number Ka defined according to Eq. 7 based on the parameters u and 0 varies across the dataset as shown in Fig. 1.

Results
All five simulations of the dataset [22] are analysed in the current work, but only three cases representative of the key trends are presented. These three cases correspond to low (1.5s L ), moderate (10s L ), and high (30s L ) levels of u , respectively. The cases are each run for a time t = 4τ 0 eddy turn over times corresponding to each u level. At this point, the turbulent burning velocity s T reaches a statistically-stationary value [22,23]. The eddy turnover time τ 0 and the Damköhler number Da = τ 0 /τ c are reported in Table 1. All statistics presented have been examined at this instant of time over the progress variable isosurface c * = 0.8, which is close to the peak reaction rate isosurface in the corresponding unstrained 1D laminar flame. The statistics have been extracted by first finding the locations corresponding to the c = c * contour by Lagrangian interpolation followed by the interpolation of required quantities to the locations found on the contour. Choice of the instant of time as well as the progress variable isosurface in the neighbourhood t ∈ (3τ 0 , 4τ 0 ) and c * ∈ (0.6, 0.9) has little influence on the statistics obtained. The particular instant of time and progress variable serve as a representative choice for the relevant statistics captured by the dataset.

Stretch rate components
In Fig. 2, the joint probability density functions (PDFs) of stretch rate components: a t , s d and h m sampled on the c * = 0.8 progress variable isosurface are shown for the three representative cases (u increases from left to right columns). The components have been normalized as in previous studies: strain rate a t is normalized using the Taylor scale strain rate [24], mean curvature h m using the laminar flame thickness δ L , the overall stretch rateṡ is scaled based on the chemical time scale δ L /s L [9], and the displacement speed is weighted by the local mixture density as s * d = ρs d /ρ R (where ρ R is the density of the unburned reactant mixture) and normalized by the laminar flame speed s L [7].
The joint PDF of a t and h m in Fig. 2a shows that the majority of points (dark concentrated regions) correspond to zero mean curvature h m ≈ 0 (horizontal line) and an extensive strain rate a t > 0 (vertical line). This conforms with the long-established picture of statisticallyplanar turbulent flames [24]. The correlation between a t and h m is negative for all u values examined (left to right in Fig. 2a). The correlation coefficient χ , defined as the covariance of the two quantities normalized by the product of their individual standard deviations, denotes the strength of the correlation. The value of χ is high at low u levels (left panel in Fig. 2a), similar to previous studies [7,14], and low at high u levels. In conclusion, the correlation strength diminishes as u increases (panels to the right in Fig. 2a).
The joint PDF of a t and s * d for the low u case in Fig. 2b shows a positive correlation: s * d tends to be high in regions of extensive tangential strain rate a t > 0, whereas s * d is low under compressive strain a t < 0. The correlation strength diminishes for successively higher u values (panels to the right in Fig. 2b). In the highest u case (right panel in Fig. 2b), negative displacement speeds are noticeable (points below the horizontal line). The majority of points correspond to near-unity displacement speeds s * d /s L (horizontal line in all panels of Fig. 2b).
The joint PDF of h m and s * d in Fig. 2c shows a negative correlation. High displacement speeds (s * d > s L ) are found exclusively in negatively curved regions (h m < 0), whereas low (or negative) displacement speeds (s * d ≤ 0) are found in positively curved regions of the flame. The correlation strength does not diminish significantly with increasing intensities in contrast to the distributions involving a t (Fig. 2a and b). Moreover, on close examination, it is seen that the joint PDF of h m and s * d has two branches (with slightly different negative slopes) that originate close to h m ≈ 0 with s * d ≈ s L . A similar distribution was computed previously for a lean methane-air flame at u = 10s L modelled using a detailed chemical-kinetic model [9]. The present result shows that this branched distribution is maintained at higher levels of u (left to right panels in Fig. 2c). Furthermore, it suggests that the distribution is maintained even for single-step chemistry flames.

Stretch rate and displacement speed
The joint PDF of normalized stretch rateṡδ L /s L and displacement speed s * d /s L is shown in Fig. 3 for the same three cases with low, moderate, and high turbulence intensities. In all three cases, the majority of the points correspond to a region with low positive stretchṡ. The overall shape of the distribution does not change between the three cases: a boomerang shape with two distinct branches is observed. A similar branched distribution was also observed in 2D [9] and 3D [21] DNS of planar methane-air flames using detailed chemistry, and across a range of flame kernel sizes in 3D DNS of spherically-propagating flames [20].
In the upper branch, s * d is negatively correlated withṡ (all panels in Fig. 3). The values of s * d in this branch are high and positive. The lower branch shows a positive correlation, with relatively-lower values of s * d andṡ. In the low intensity case (left panel in Fig. 3), most of the distribution belongs to this branch, supporting the use of a linear relation s d andṡ as consistent with relevant theory [19]. The magnitude ofṡ in the lower branch increases with u (left to right panels Fig. 3). The stretch rateṡ peaks at a low positive value in each case that increases with u . This peak value, however, remains much smaller in magnitude than the negative stretch values attained. As noted previously for both single-step chemistry [20] and detailed chemistry flames [9,21], the occurrence of high negative stretch rates is rare through the range of u /s L values considered. The joint PDF ofṡ and s * d (shown in Fig. 3) is examined further in Fig. 4 by colouring the distribution based on strain rate a t and mean curvature h m . In all three cases, the upper branch corresponds to negative curvature h m < 0 and positive strain rate a t > 0 (red squares in Fig. 4). The lower branch corresponds instead to regions with positive curvature h m > 0 (blue and yellow squares in Fig. 4). In the low intensity case (left panel in Fig. 4), the lower branch predominantly is characterised with compressive strain rate a t < 0 (yellow squares). In higher intensity cases, positive strain rate a t > 0 dominates regions with low stretch (blue squares) and negative strain rate a t > 0 dominates regions with relatively higher stretch (yellow squares). The fourth configuration, negative curvature and compressive strain rate, has the least frequent contribution to the entire distribution. The shift in contributions from different regions of the flame surface points to the underlying competition between strain and curvature in determining the overall mean stretch rate ṡ of the flame surface (here, the averaging · implies a surface mean).
The overall distribution in any of the cases may be interpreted by following the distribution curve from its upper left corner downwards towards the lower branch. The upper left corner of the distribution corresponds to negative curvature regions h m < 0 where s * d takes on large positive values and the stretch rateṡ is high but negative. Asṡ decreases in magnitude along the upper branch, the normalized displacement speed s * d /s L tends to unity. Wheṅ s is close to zero and s * d /s L is near unity (lower right corner), the contributions from positive curvature h m > 0 begin to gain prominence. As stretch rateṡ increases further along the lower branch, h m > 0 and s * d > 0 so that propagation rate s d h m alone cannot account for the increasingly negative stretch rate values. This indicates that compressive strain rate is becoming important.

Displacement speed components
The local instantaneous displacement speed s d can be decomposed into separate components [25] as follows which can be expressed as s d = s r +s n +s t in terms of the reaction s r , normal (diffusion) s n and tangential (curvature) s t components. In the following, each of these components s r,n,t have been weighted by the local mixture density (s * i = ρs i /ρ R for i = {r, n, t} where ρ R is the density of the unburned reactant mixture) and normalized using the unstrained laminar flame speed s L . The profiles of the different s * d components across the instantaneous flame brush (the region corresponding to 0.05 < c < 0.95) are presented in Fig. 5a alongside the corresponding laminar case. In the laminar case (first panel in Fig. 5a), the reaction component s * r (blue line) is non-zero downstream of the preheat zone 0.6 < c < 0.95 where heat release occurs. The normal diffusion component s * n (green line) is positive in the preheat zone 0.05 < c < 0.6 and negative in the reaction zone 0.6 < c < 0.95. The tangential component s * t (red line) remains zero throughout. In the turbulent cases presented alongside, the range of the curvature component s * t is widened (red squares in the second panel onwards to the right in Fig. 5a) as compared to the laminar case (first panel in Fig. 5a). The range of s * t , however, remains within a band of fairly constant width across most of the instantaneous flame brush, indicating that isosurfaces of progress variable remain parallel to each other in this region. The width of this band increases with u , particularly near the leading edge 0.05 < c < 0.2. In the low-intensity case (second panel in Fig. 5a) the curvature component s * t /s L at any progress variable value ranges between slightly negative to unity, the majority exhibiting a negative value, especially towards the trailing edge 0.8 < c < 0.95. The s * t distribution tends to favour negative values as the turbulence intensity u is increased (second and last panel in Fig. 5a).
The normal diffusion component s * n stays close to the laminar profile in the low intensity case (green squares in the second panel in Fig. 5a), with a minor spread near the leading edge. The spread near the leading edge 0.05 < c < 0.2 increases as u is increased (green squares from second panel to last panel). In the highest intensity case, the normal diffusion  Fig. 5a).
The normal and reaction components taken together s * r + s * n are compared with the contribution of the curvature component s * t in Fig. 5b. In the low intensity case (second panel), reaction and normal components (blue squares) constitute together nearly the entire displacement speed. The sum s * r +s * n tends to be close to unity, i.e. s * r +s * n ∼ s L in line with the argument of Peters [26] for the TRZ regime. With increasing turbulence intensity (second panel to the last panel), the curvature component s * t contributes increasingly to the overall displacement speed supporting the inferences of 2D DNS conducted more recently [13]. In the highest intensity cases (right panels in Fig. 5a), the sum (s * r + s * n )/s L tends to be less than unity. This is true especially near the leading edge, highlighting the role of negative diffusion component s * n values in this region of the instantaneous turbulent flame brush, i.e. the space corresponding to 0.01 < c < 0.99.

Strain rate and curvature response of displacement speed components
The components of displacement speed analysed above are further examined for their dependence on strain and curvature in Figs. 6 and 7, respectively.   Figure 6 shows the joint PDF between the the local strain rate a t and the components of displacement speed grouped as a) s * r + s * n and b) s * t . First, the terms a t and s * r + s * n (Fig. 6a) appear to be negatively correlated, especially at low intensities (left panel in Fig. 6a), with a diminished correlation strength at higher turbulence intensities (panels to the right in Fig. 6a). This means that the contribution to displacement speed s * d from the normal and reaction components taken together s * r +s * n reduces as an extensive strain rate is experienced by the flame surface. For compressive straining a t < 0, this contribution can be several times the unstrained laminar flame speed s L . The dependence of the tangential component s * t = −2Dh m on strain rate a t is shown in Fig. 6b. As expected, this distribution follows (with a negative sign) the dependence between curvature h m and strain rate a t as noted in Fig. 2a. Interestingly, however, Fig. 6b reveals another aspect this joint probability distribution which is the increasing tendency of s * t (i.e. h m in Fig. 2a) to attain near-zero values over a range of extensional strain a t > 0 as u is increased (panels from left to right in Fig. 6b). Figure 7 shows the joint PDF of the local mean curvature h m and the components of displacement speed grouped as a) s * r + s * n and b) s * t . First, the joint probability distribution of s * r + s * n and h m (Fig. 7a) shows a definite positive correlation, with high probabilities of near-zero positive curvatures where the contribution of s * r +s * n tends to be of the order of s L . At higher intensities (panels to the right in Fig. 7a) the dependence is not so straightforward even though signatures of positive correlation are evident. The dependence of the tangential component s * t = −2Dh m on curvature h m is shown for sake of completeness in Fig. 7b. A deterministic linear behaviour with a negative slope is recovered as expected. Together, these joint PDFs as shown in Figs. 6 and 7 highlight the relatively complex and varied dependence of the components of s * d on the remaining components -a t and h m -of the local stretch rate.

Profiles of stretch rate components
The results above suggest that the non-linear contribution of displacement speed s d to the stretch rateṡ is significant. Analysis of s d components suggests further that the curvature component s t contributes primarily to this non-linearity. Using the expression for s * d in Eq. 8 in the stretch rate equation Eq. 3, the following expression resultṡ where the tangential component s t = −2Dh m is written explicitly in terms of curvature to highlight the non-linear contribution of the term 2s t h m = −4Dh 2 m towards the stretch rate. The normalized variation of overall stretch rateṡ and its components across the instantaneous turbulent flame brush is shown in Fig. 8. Figure 8a shows the variation of strain rate, while the components of propagation rate s d h m are shown grouped separately as: 1) the sum of normal and reaction propagation rates s n h m + s r h m (Fig. 8b) and the tangential propagation rate s t h m (Fig. 8c). The variation of the overall stretch rate is shown in Fig. 8d. The profiles are coloured based on the turbulence intensity (red squares: low-intensity, blue squares: moderate-intensity, and green squares: high intensity).
The variation of strain rate a t across the instantaneous flame brush is shown in Fig. 8a. Consistent with the vertical lines in Fig. 2a and b, the (normalized) strain rate has a (slightly) positive mean conforming with previous investigations [24]. The spread of the distribution about zero increases with turbulence intensity (red through blue to green squares). This spread is greater near the leading edge 0.05 < c < 0.2. The increased spread in the leading edge region suggests that the turbulent flow upstream has higher vorticity possibly as it is prior to suppression by dilatation across the flame [27]. This is consistent with the Thin Reaction Zones regime.
The combination s r h m + s n h m tends to be negative in the preheat zone 0.05 < c < 0.6 and slightly positive in the reaction zone 0.6 < c < 0.95 region (red squares in Fig. 8b). This trend is more apparent at higher intensities (blue and green squares in Fig. 8b). The tangential component s t h m , however, is always negative (−4Dh 2 m < 0∀h m ∈ R, D ∈ R + ), and the range of s t h m increases with turbulence intensity (see green squares in Fig. 8c).
The mean values of s t h m (horizontal lines in Fig. 8c) become increasingly negative with increasing u . It is evident, therefore, that the contribution of s t h m increases in magnitude with u and it does not average out overall (right panel). Rather, it is the contribution of s r h m + s n h m that tends to average out over the flame brush (horizontal lines in Fig. 8b). This provides a contrasting data point with respect to various approaches for modelling the stretch rate that consider the curvature contribution s t h m to cancel out in the mean [5], and provides support for more recent investigations that highlight its importance [13].
The variation of overall stretch rateṡ is shown in Fig. 8d. For all three cases (red: low u , blue: moderate u , and green: high u ), the overall stretch rate has both positive and negative contributions. At low turbulence intensities (red squares), the most frequent contribution is from low positive stretch. This contribution appears to increase gradually across the instantaneous flame brush (low to high values of c in Fig. 8d). Negative contributions arise from Stretch rate and its components obtained by decomposing the displacement speed for increasing values of turbulence intensity (red: u = 1.5s L , blue: u = 10s L , and green: u = 30s L fewer regions in the flame brush (red squares below horizontal line). At higher intensities (blue and green squares), this profile of stretch rates seems to follow with a much wider spread of stretch rates.

Implications for flamelet modelling
The results discussed above have some interesting implications with regard to the modelling of turbulent premixed flames in the context of the Flame Surface Density (FSD) approach [3,4,11] and related flamelet-based modelling approaches.
Firstly, there appears to be no qualitative change in the statistical behaviour of the displacement speed and its separate components as the turbulence intensity u is increased over the range considered here (as per Fig. 3). Even though the scatter in all quantities increases with u as might be expected, no new correlations are seen to emerge. This implies that a well-formulated FSD model should be applicable equally across the range of u spanning the TRZ regime.
Secondly, accurate modelling of the displacement speed appears to be a requirement in order to capture the non-linearities observed in the stretch rate and its components. Indeed, separate models may be needed for each of the three components of displacement speed, owing to their rather different response to strain and curvature (as per Figs. 6 and 7). Previous modelling [13,17] has gone some way towards achieving this goal, but more work is needed still.

Conclusions
Joint probability density functions of mean curvature and displacement speed as well as stretch rate and displacement speed have been extracted from a DNS database wherein turbulence intensity is increased successively over five separate simulations. The shapes of both of these joint PDFs show little variation across the range of turbulence intensities in the Thin Reaction Zones regime. The distributions show a branched shape akin to previous observations in lean [9] as well as stoichiometric [20] methane-air flames. Components of displacement speed and, in particular, the tangential component through its dependence on curvature appears to play a significant role in effecting the branched structure of the distribution. The inherent non-linearities involved in the expression for stretch rate indicate that the recursive relationship between stretch rate and displacement speed is a worthy topic of further investigation.