#### Abstract

A direct numerical simulation (DNS) database of freely propagating statistically planar turbulent premixed flames with a range of different turbulent Reynolds numbers has been used to assess the performance of algebraic flame surface density (FSD) models based on a fractal representation of the flame wrinkling factor. The turbulent Reynolds number Re_{t} has been varied by modifying the Karlovitz number Ka and the Damköhler number Da independently of each other in such a way that the flames remain within the thin reaction zones regime. It has been found that the turbulent Reynolds number and the Karlovitz number both have a significant influence on the fractal dimension, which is found to increase with increasing Re_{t} and Ka before reaching an asymptotic value for large values of Re_{t} and Ka. A parameterisation of the fractal dimension is presented in which the effects of the Reynolds and the Karlovitz numbers are explicitly taken into account. By contrast, the inner cut-off scale normalised by the Zel’dovich flame thickness does not exhibit any significant dependence on Re_{t} for the cases considered here. The performance of several algebraic FSD models has been assessed based on various criteria. Most of the algebraic models show a deterioration in performance with increasing the LES filter width.

#### 1. Introduction

Large eddy simulation (LES) is becoming increasingly popular for computational fluid dynamics (CFD) analysis of turbulent reacting flows due to the advancement and increased affordability of high-performance computing. The exponential temperature dependence of the chemical reaction rate poses one of the major challenges in LES modelling of turbulent reacting flows [1, 2]. Reaction rate closure models based on the concept of flame surface density (FSD) are well established in the context of the Reynolds averaged Navier Stokes simulations [3, 4] of turbulent premixed flames. However, the application of FSD-based modelling in LES is relatively recent [5–10]. The generalised FSD () is defined as [5] where is the reaction progress variable. The overbar indicates the LES filtering operation in which the filtered value of a general quantity is evaluated as , where is a suitable filter function [5]. The combined contribution of the filtered reaction and molecular diffusion rates can be modelled using as , where is the fluid density, is the progress variable diffusivity, is the density-weighted surface-filtered displacement speed , and is the surface-filtered value of a general quantity . Often, is expressed in terms of the wrinkling factor , which is defined as Hence the prediction of depends on the accuracy of the modelling of [8]. Several models for have been proposed in the context of LES [11–16], many of which make use of a power-law scaling involving a fractal representation of the flame surface.

To date, most of the algebraic models for FSD based on the wrinkling factor have been developed for the corrugated flamelet (CF) regime [17] in which the flame thickness remains smaller than the Kolmogorov length scale. However, there has been no detailed assessment of the performance of these models in the thin reaction zone (TRZ) regime [17], in which the Kolmogorov length scale remains smaller than the flame thickness. Moreover, the effects of turbulent Reynolds number on these models also remain to be investigated in detail. In order to address these gaps in the existing literature, the performance of several algebraic FSD models has been assessed here for TRZ regime combustion based on a direct numerical simulation (DNS) database of freely propagating statistically planar turbulent premixed flames with different values of . The variation of is brought about by varying the Damköhler number and the Karlovitz number independently of each other, using the following relation [17]: Here, is the turbulent Reynolds number, is the Damköhler number and is the Karlovitz number [17]. Subscript 0 denotes quantities evaluated in the unburned reactants, is the turbulent velocity fluctuation magnitude, is the turbulence integral length scale, is the dynamic viscosity, is the thermal diffusivity, is the laminar burning velocity, and is the thermal thickness of the laminar flame.

The main objectives of the present study are(1)to assess the performance of a range of wrinkling factor based algebraic models in the TRZ regime,(2)to identify the influence of on the performance of these models in the TRZ regime.

A summary of the existing wrinkling factor models will be provided in the next section. This will be followed by a discussion on the numerical implementation. Following this, the results will be presented and discussed. Finally, the main findings will be summarised and conclusions will be drawn.

#### 2. Mathematical Background

The wrinkling factor can be expressed in the form of a power law [11–16, 18] , where and are the outer and inner cut-off scales and is the fractal dimension of the flame surface. For LES, is taken to be equal to the filter width . According to Peters [17], scales with the Gibson scale (the Obukhov-Corrsin scale ) in the CF (TRZ) regime, with being the dissipation rate of turbulent kinetic energy (TKE). Knikker et al. [16] found by experiment that scales with the Zel’dovich flame thickness . A recent DNS analysis [8] found that does indeed scale with and for the CF and TRZ regimes, respectively, but also scales with for both regimes.

North and Santavicca [19] parameterised as , which suggests that increases with . Kerstein [20] indicated that increases from 2 to 7/3 for increasing values of , where is associated with a nonpropagating material surface. The *a priori* DNS analysis by Chakraborty and Klein [8] demonstrated that indeed increases from a value slightly greater than 2.0 in the CF regime to 7/3 for high unity Lewis number flames in the TRZ regime.

Weller et al. [11] proposed a model for , (denoted here as FSDW), which can be recast in terms of as where and with being the Kolmogorov length scale. The subgrid scale velocity fluctuation magnitude is defined as . A model for proposed by Angelberger et al. [12] (denoted as FSDA) can be written as where is a model parameter of the order of unity and the efficiency function is defined as Colin et al. [13] proposed a slightly different model (denoted as FSDC): where is given by and with . Charlette et al. [14] reduced the input parameters to only and using the expression (denoted as FSDCH) with the efficiency function where , and with model constants , , and functions , , , and defined by Fureby [15] proposed a model (denoted as FSDF) which can be written as: where is given by and [19]. Knikker et al. [16] proposed a model (denoted as FSDK) as where and is dynamically evaluated as where denotes the filtered at the test filter level . The angled bracket indicates a volume averaging operation, as often used in dynamic models [16, 21].

In Section 4, the performance of these models is assessed with respect to obtained from DNS based on the following criteria.

*Criterion 1. *As FSD represents the flame surface area to volume ratio [3], the volume-averaged value of the generalised FSD over the DNS domain represents the total flame surface area within the domain and therefore should be independent of . Thus the model predictions of should not change with .

*Criterion 2. *The models for should be able to capture the correct variation of the mean values of conditional on across the flame brush.

*Criterion 3. *The correlation coefficient between the modelled and actual values of should be as close as possible to unity, in order to capture correctly the local strain rate and the curvature effects on in the context of LES.

#### 3. Numerical Implementation

Three-dimensional DNS of freely propagating statistically planar turbulent premixed flames has been carried out using the DNS code SENGA [22]. The domain of size was discretised using a Cartesian mesh of size with uniform mesh spacing in each direction. The grid spacing was determined by the flame resolution, and in all cases, about 10 grid points are kept within . The boundaries in the direction of the mean flame propagation (i.e., -direction) were taken to be partially nonreflecting whereas the transverse directions were taken to be periodic. Higher order finite-difference and Runge-Kutta schemes were employed for spatial discretisation and time advancement, respectively. The scalar fields were initialised using a steady unstrained planar laminar premixed flame solution. For the present study, standard values were taken for the Prandtl number (), the Zel’dovich number (), and the ratio of specific heats (), where , and are the activation temperature and the specific heats at constant pressure and constant volume, respectively. The initial values of , , , , and are listed in Table 1. In cases A, C, and E (B, C, and D), the values of and were chosen to vary by changing () while keeping () constant (see (3)). The heat release parameter and the Lewis number are taken to be equal to 4.5 and 1.0, respectively. Table 1 shows that remains greater than unity for all cases indicating the TRZ regime combustion [17]. In all cases the flame-turbulence interaction takes place under decaying turbulence, and the simulation time corresponds to one chemical time scale . This time is equal to in case D, in cases A, C, and E, and for case B where is the initial eddy turn-over time. The present simulation time remains comparable with several previous DNS studies which focused on the modelling of [5, 7–9, 14, 21, 23]. By the time statistics were extracted, the global TKE and volume-averaged burning rate were no longer changing rapidly with time [24]. The global level of turbulent velocity fluctuation had decayed by 53%, 61%, 45%, 24%, and 34% in comparison to the initial values for cases A–E, respectively. By contrast, the integral length scale increased by factors of 1.5 to 2.25, ensuring that a sufficient number of turbulent eddies was retained in each direction. Values for , , and at the time when statistics were extracted were presented in [24, Table 2] and are not repeated here (note that the Karlovitz number was defined in [24] in terms of the thermal flame thickness as , whereas in this paper is defined in terms of the Zel’dovich flame thickness as ).

In the present analysis, the DNS data is explicitly filtered using [5] to obtain the relevant filtered quantities. The results will be presented for ranging from to , where is the DNS mesh size . These filter sizes are comparable to the range of used in *a priori* DNS analysis in several previous studies [5, 7, 9, 14, 18] and span a useful range of length scales (i.e., from comparable to , where the flame is partially resolved, up to where the flame becomes fully unresolved and is comparable to the integral length scale). For this range of filter widths, the underlying combustion process varies from the “laminar flamelets-*G* DNS” [25] combustion regime (for ) to well within the TRZ regime (for ) on the regime diagrams by Pitsch and Duchamp de Lageneste [25] and Düsing et al. [26].

#### 4. Results and Discussion

##### 4.1. Flame Turbulence Interaction

Contours of in the mid-plane at time are shown in Figures 1(a)–1(e) for cases A–E, respectively. Figure 1 shows that the extent of flame wrinkling increases with increasing (see Table 1). Moreover, the contours of representing the preheat zone (i.e., , see colour scale) are much more distorted than those representing the reaction zone (i.e., ). This tendency becomes more prevalent with the increase in the Karlovitz number , since the scale separation between and increases with increasing Ka, allowing more energetic eddies to enter into the preheat zone causing greater distortion of the flame.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

##### 4.2. Behaviour of and in Response to

The power-law produces the expression where the angled brackets indicate the volume-averaging operation. The quantity decreases with increasing due to an increase in the proportion of flame wrinkling that occurs at the sub-grid scales with increasing . By contrast, indicates the total flame surface area in the computational domain, thus remaining independent of . As a result, increases with increasing , which can be substantiated from Figures 2(a)–2(e). The variation of with is linear when but becomes nonlinear for . The slope of the best-fit straight line representing the greatest slope of the variation of with gives while the intersection point of this straight line with the abscissa gives . It has been found that remains about 2.0 for all cases (i.e., ) and that for the present thermochemistry. This is qualitatively consistent with previous experimental [16, 27] and computational [8] findings.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

Figures 2(a)–2(e) demonstrate that is greater for flames with higher and that attains an asymptotic value of 7/3 for flames with (e.g., cases C, D, and E). The extent of the flame wrinkling and the flame surface area generation increases with increasing . This can be substantiated from the values of normalised flame surface area obtained by volume integrating (i.e., , where is an infinitesimal volume element), which yields 1.15, 1.33, 1.87, 3.63, and 3.70 for cases A–E, respectively, at the time when statistics were extracted. This behaviour is consistent with Figure 1 which demonstrates that the wrinkling of isosurfaces increases with increasing . Figure 2 suggests that is expected to increase with increasing before assuming an asymptotic value when either Da or is held constant. This behaviour is also qualitatively similar to the trend predicted by the parameterisation of North and Santavicca [19]. The parameterisation by Chakraborty and Klein [8], that is, , predicts for all the cases considered here because this parameterisation accounts for the dependence of on Ka only. Based on the above findings, is parameterised here as where is a model parameter. The prediction of with obtained from DNS and obtained from (11) is also shown in Figures 2(a)–2(e), which indicates that (11) satisfactorily captures the best-fit straight line corresponding to the power law. According to (11), increases from a value close to 2 for small values of (e.g., cases A and B) to an asymptotic value of 7/3 for large values of and (e.g., cases D and E), according to Kerstein [20].

It has been found in previous analyses [8, 19, 28] that approaches for unity Lewis number flames within the flamelet regime when approaches a value of about and the Karlovitz number remains greater than unity (i.e., ). The fractal dimension for the present case C is found to be , and this small deviation from is within the statistical noise and should not be over-interpreted. The fractal dimensions for cases D and E are equal to which indicates that for high Karlovitz numbers (i.e., ) within the flamelet regime, approaches an asymptotic value of for according to the present analysis.

Note that and in (11) were evaluated here based on and in the unburned reactants. However, in actual LES simulations, needs to be evaluated based on local velocity and length scale ratios (i.e., and ). Here is estimated from the sub-grid TKE (i.e., ) following previous studies [6, 8, 10–12, 15, 18]. The local Karlovitz number can be evaluated as where is a model parameter.

A power-law-based expression for is proposed here based on the observed behaviours of and in Figure 2: where is a bridging function which increases monotonically from 0 for small (i.e., or ) to 1 for large (i.e., or ). The expression given by (12) will be denoted as FSDNEW. Equation (12) ensures that approaches (i.e., ) for large (i.e., or ) and at the same time approaches (i.e., ) for small (i.e., ). It has been found from Figure 2 that provides better agreement with obtained from the DNS data for , whereas the power-law starts to predict more accurately for . Based on this observation, the bridging function is taken to be which ensures a smooth transition . As is found to scale with (i.e., according to the present thermochemistry), in (12) is taken to be the thermal flame thickness . In the context of LES, needs to be evaluated based on local quantities, which is achieved by replacing Ka and in (11) with their local values and , respectively. The choice of model constants and ensures an accurate prediction of for and yields the value of obtained based on the global quantities according to (11).

##### 4.3. Assessment of Model Performance

*Criterion 1. *The inaccuracy in the model predictions of can be characterised using a percentage error (PE)
where is the model prediction of . Results for the PE for a range of are shown in Figure 3, which demonstrate that the models denoted by FSDA (see and and FSDC (see (6)) overpredict for all the cases and that the level of overprediction increases with increasing . The FSDW model (see (4)) also overpredicts the value of , although the level of overprediction for is smaller for this model, especially for the cases with higher (i.e., cases D and E). The FSDC model remains inferior to both the FSDA and FSDW models for all in all cases. The PE for the FSDCH model (see , , and remains small for the small and moderate cases (i.e., cases A–C), although the FSDCH model slightly overpredicts the value of for for the higher cases (i.e., cases D and E). Replacing (i.e., for the present thermochemistry) by in – leads to a deterioration of the performance of the FSDCH model especially for higher cases where it starts to overpredict the value of for . However, the performance of the FSDCH model even with in – remains better than the FSDA and FSDC models for all cases considered here. The FSDK model (see (8) underpredicts the value of for all for all cases. However, the level of underprediction decreases for larger .

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

The PEs for the FSDCH, FSDF (see (9)), and FSDNEW (see (12)) models remain negligible in comparison to the PEs for all the other models considered. Figure 3 shows that the FSDF model underpredicts slightly for small and moderate values of (i.e., cases A–C), but the prediction of this model remains very close to the DNS result for high values of (cases D and E). The PEs for the FSDCH, FSDNEW and FSDF models remain comparable. Note that should approach (i.e., ) when vanishes because the flow tends to be fully resolved (i.e., and ). Although the FSDF model performs well for all for all the cases considered here, does not tend to as approaches zero but instead predicts a finite value close to zero. Hence the FSDF model underpredicts the value of for small values of for all the cases considered here (see Figure 3). This limitation of the FSDF model can be avoided using a modified form of (8): where is a bridging function as before, the efficiency function is given by , and [19]. Equation (14) ensures that becomes exactly equal to when the flow is fully resolved (i.e., or ) where also vanishes (i.e., ). The model given by (14) is denoted as the MFSDF model. Figure 3 shows that the modification given by (14) does not appreciably alter the performance of (8), but this modification ensures that the model given by (8) will approach the correct asymptotic limit (i.e., ) for very small (i.e., ). Note that the combination of parameterisation of and according to [19] and , respectively, is essential for the satisfactory performance of the FSDF model. Using (13), for in the FSDF model is found to lead to a deterioration in its performance. Similarly, using as given by [19] in (12) worsens the performance of the FSDNEW model.

The FSDK model is based on a power-law scheme which is strictly valid only for which are sufficiently greater than (i.e., ), as can be seen from Figure 2. Hence the predictive capability of the FSDK model improves when (see Figure 3). However, the FSDK model underpredicts Σ_{gen} because the inner cut-off scale is taken to be in this model whereas for all the cases considered here. An accurate estimation of in the framework of the FSDK model results in comparable performance to the FSDNEW model for large (i.e., ). Moreover, vanishes when according to the FSDK model, whereas should approach when (i.e., ). This limitation can be avoided by modifying the FSDK model in the same manner as shown in (14) for the FSDF model (not shown here for conciseness).

The stretch rate represents the fractional rate of change of flame surface area [3], where is the displacement speed, is the local flame normal vector, and is the tangential strain rate. It is possible to decompose into the reaction, normal diffusion, and tangential diffusion components (i.e., , , and ) as [17, 24, 29, 30] , , and . It has been shown previously [6, 7, 9] that remains positive and thus acts to generate flame area, whereas is primarily responsible for flame area destruction. The equilibrium of flame area generation and destruction yields , which leads to [8] . The stretch rate induced by becomes the leading order sink term in the TRZ regime [17]. However, most algebraic models (e.g., FSDA, FSDC, and FSDW) were proposed in the CF regime based on the equilibrium of the stretch rates induced by and , and the flame surface area destruction due to the term was ignored [10–14]. Hence these models underestimate the flame surface area destruction in the TRZ regime, which leads to overprediction of for the FSDA, FSDC, and FSDW models. Although the FSDCH model was proposed for the CF regime (where the term was ignored), the efficiency function was tuned to capture the expected behaviour of the turbulent flame speed for both (as in the CF regime) and (as in the TRZ regime). Hence this model is somewhat capable of predicting the behaviour of for the TRZ regime flames considered here. However, this model starts to overpredict due to underestimation of the destruction of flame surface area in the TRZ regime for higher cases where the effects of are significant. Moreover, the performance of this model is sensitive to the definition of the flame thickness used in the efficiency function.

*Criterion 2. *To assess model performance with respect to Criterion 2, the variations of mean conditionally averaged on are shown in Figure 4 for and in Figure 5 for . These values have been chosen since they correspond to and respectively. The following observations can be made from Figure 4.(i)The FSDW model tends to overpredict the value of for flames with higher (e.g., cases D and E). However, the extent of this overprediction is relatively lower in the low cases (e.g., cases A and B). (ii)The FSDW model tends to predict a peak value of at , whereas the peak value of obtained from DNS is attained close to for all the cases.(iii)The FSDK model tends to underpredict the value of for all cases. The physical explanations provided earlier for the underprediction of by the FSDK model (see Figure 3) are also valid here.(iv)Models FSDA, FSDC, FSDCH, FSDF, and FSDNEW all tend to capture the variation with obtained from DNS data. The prediction of the MFSDF model remains comparable to that of the FSDF model.

Comparing Figure 4 with Figure 5 reveals that there is significantly greater spread between the predictions of the various models for than for . The following observations can be made from Figure 5. (i)Models FSDW, FSDA, and FSDC all overpredict the value of , and the overprediction increases with increasing . The FSDCH model captures the behaviour of for small values of (e.g., cases A–C), but it overpredicts the value of for higher cases (e.g., cases D–E).(ii)Similar to Figure 4, the FSDW model predicts a peak at > 0.6, whereas the peak value of , from DNS occurs at for all cases.(iii)Models FSDF, FSDK, FSDNEW and MFSDF predict satisfactorily throughout the flame brush.(iv)The difference between the models MFSDF and FSDF seems to be very small, and both predict satisfactorily.

Unlike any of the other models, the prediction of the FSDK model improves with increasing , which is consistent with observations made in the context of Figure 3. Moreover, the prediction of the FSDW model remains skewed towards the burned products due to the dependence of (i.e., ). The models FSDW, FSDA, and FSDC underestimate the destruction of flame surface area in the TRZ regime by ignoring destruction arising due to which eventually leads to the overprediction of . As discussed above, the efficiency function in the FSDCH model is parameterised to capture the turbulent flame speed behaviour in both the CF and TRZ regimes, and hence this model performs satisfactorily with respect to criterion 2 for all cases considered here.

The use of local Karlovitz number and turbulent Reynolds number enables local variation of , and a satisfactory performance of the FSDNEW model with respect to Criteria 1 and 2 indicates that (11) satisfactorily captures the local Reynolds and Karlovitz number dependences of for the present definitions of these numbers (i.e., and ).

*Criterion 3. *The correlation coefficients between the DNS result and the model prediction for are shown in Figure 6 for . For and is small and thus the correlation coefficients in these regions of the flame brush do not have much physical significance. On comparing cases A–E, it is clear that the correlation coefficients are significantly affected by and by . With increasing , the correlation coefficients decrease significantly with increasing which indicates that the accuracy of the models deteriorates for large , as flame wrinkling increasingly takes place at the sub-grid scale with increasing . A modelled transport equation [6, 7, 9, 10, 18] for may be advantageous in terms of capturing the correct strain rate and the curvature dependences of , provided the unclosed terms of transport equation are accurately modelled.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

It has been shown in Figures 3–6 that the FSDCH, FSDF, MFSDF, and FSDNEW models perform better than the other models in the thin reaction zones regime flames considered here. Note that the MFSDF model is now considered instead of the FSDF model since it removes the limitations of the FSDF model when the flow is fully resolved (i.e., and ). Moreover, the performance of the FSDCH, MFSDF, and FSDNEW models remains comparable, and thus any of these models is likely to provide the most satisfactory prediction within the thin reaction zones regime.

#### 5. Conclusions

The performance of several wrinkling-factor-based algebraic models for for flames within the TRZ regime has been assessed based on DNS of turbulent premixed flames over a range of different values of . It has been found that the fractal dimension increases with increasing and before reaching an asymptotic value. By contrast, the inner cut-off scale is not significantly affected within the range of considered. The observed behaviours of and have been incorporated into a new power-law model for in the context of LES. Various criteria have been used to assess the performance of this model along with the other existing models. Most models show a deterioration of performance with increasing , and in general the performance of the models is better for lower . Based on this assessment, models have been identified which predict satisfactorily for all the cases considered here and for different values of . It is worth noting that the present study has been carried out for a range of moderate without the effects of detailed chemistry and transport. Thus, three-dimensional DNS with detailed chemistry along with experimental data for higher values of will be necessary for a more comprehensive analysis.

#### Acknowledgment

The authors are grateful to EPSRC, UK, for financial assistance.