Effects of Turbulent Reynolds Number on the Performance of Algebraic Flame Surface Density Models for Large Eddy Simulation in the Thin Reaction Zones Regime: A Direct Numerical Simulation Analysis

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 (Σ_gen) is defined as [5] Σ_gen=|∇𝑐|, (1) 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 Σ_gen as 𝑤̇+∇⋅(𝜌𝐷_𝑐∇𝑐)=(𝜌𝑆_𝑑)_𝑠Σ_gen, 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, Σ_gen is expressed in terms of the wrinkling factor Ξ, which is defined as Ξ=|∇𝑐||∇𝑐|=Σ_gen|∇𝑐|. (2) Hence the prediction of Σ_gen 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 Re_𝑡 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 Re_𝑡. The variation of Re_𝑡 is brought about by varying the Damköhler number Da and the Karlovitz number Ka independently of each other, using the following relation [17]: Re_𝑡∼Da²Ka²∼𝑢^′𝑆_𝐿²𝑢Da∼^′𝑆_𝐿⁴Ka⁻². (3) Here, Re_𝑡=𝜌₀𝑢^′𝑙/𝜇₀ is the turbulent Reynolds number, Da=𝑙𝑆_𝐿/𝑢^′𝛿_th is the Damköhler number and Ka=(𝑢^′/𝑆_𝐿)^3/2(𝑙𝑆_𝐿/𝛼_𝑇0)^−1/2 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 𝛿_th 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 Σ_gen models in the TRZ regime,(2)to identify the influence of Re_𝑡 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] Ξ=(𝜂₀/𝜂_𝑖)^𝐷−2, where 𝜂₀ and 𝜂_𝑖 are the outer and inner cut-off scales and 𝐷 is the fractal dimension of the flame surface. For LES, 𝜂_o is taken to be equal to the filter width Δ. According to Peters [17], 𝜂_𝑖 scales with the Gibson scale 𝐿_𝐺=𝑆³_𝐿/𝜀 (the Obukhov-Corrsin scale 𝜂_OC=(𝐷_𝑐³/𝜀)^1/4) 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 𝛿_𝑧=𝛼_𝑇0/𝑆_𝐿. A recent DNS analysis [8] found that 𝜂_𝑖 does indeed scale with 𝐿_𝐺 and 𝜂_OC for the CF and TRZ regimes, respectively, but also scales with 𝛿_𝑧 for both regimes.

North and Santavicca [19] parameterised 𝐷 as 𝐷=2.05/(𝑢^′/𝑆_𝐿+1)+2.35/(𝑆_𝐿/𝑢^′+1), which suggests that 𝐷 increases with 𝑢^′/𝑆_𝐿. Kerstein [20] indicated that 𝐷 increases from 2 to 7/3 for increasing values of 𝑢^′/𝑆_𝐿, where 𝐷=7/3 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 Ka 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 Σ_gen as Σ_gen=[1+2𝑐̃(Θ−1)]|∇𝑐|, (4) where Θ=1+0.62𝑢^′_Δ/𝑆_𝐿Re_𝜂 and Re_𝜂=𝜌₀𝑢^′_Δ𝜂/𝜇₀ with 𝜂 being the Kolmogorov length scale. The subgrid scale velocity fluctuation magnitude is defined as 𝑢^′_Δ=(𝑢_𝑖𝑢_𝑖−𝑢_𝑖𝑢_𝑖)/3. A model for Ξ proposed by Angelberger et al. [12] (denoted as FSDA) can be written as Σ_gen=𝑢1+𝑎Γ^′Δ𝑆_𝐿|∇𝑐|, (5i) where 𝑎 is a model parameter of the order of unity and the efficiency function Γ is defined as 𝑢Γ=0.75exp−1.2^′Δ𝑆_𝐿^−0.3Δ𝛿_𝑧^2/3. (5ii)Colin et al. [13] proposed a slightly different model (denoted as FSDC): Σ_gen=𝑢1+𝛼Γ^′Δ𝑆_𝐿|∇𝑐|, (6) where Γ is given by (5ii) and 𝛼=2ln(2)/[3𝑐_ms(Re_𝑡^1/2−1)] with 𝑐_ms=0.28. Charlette et al. [14] reduced the input parameters to only 𝑢^′_Δ/𝑆_𝐿 and Δ/𝛿_𝑐 using the expression (denoted as FSDCH)Σ_gen=Δ1+min𝛿_𝑐,Γ_Δ𝑢^′_Δ𝑆_𝐿^𝛽1|∇𝑐| (7i) with the efficiency function Γ_Δ=𝑓^−𝑎1_𝑢+𝑓^−𝑎1_Δ^{−1/𝑎1−𝑏1}+𝑓^−𝑏1_Re^−1/𝑏1, (7ii) where 𝛿_𝑐=4.0𝜇₀/𝜌₀𝑆_𝐿, Re_Δ=4(𝑢^′_Δ/𝑆_𝐿)⋅(Δ/𝛿_𝑐) and with model constants 𝑏₁=1.4, 𝛽₁=0.5, 𝐶_𝑘=1.5 and functions 𝑎₁, 𝑓_𝑢, 𝑓_Δ, and 𝑓_Re defined by 𝑎₁𝑢=0.60+0.20exp−0.1^′Δ𝑆_𝐿Δ−0.20exp−0.01𝛿_𝑐,𝑓_𝑢=427𝐶110^𝑘1/218𝐶55^𝑘𝑢^′_Δ𝑆_𝐿²,𝑓_Δ=27𝐶110^𝑘𝜋^4/3Δ𝛿_𝑐^4/3−1^1/2,𝑓_Re=955exp−1.5𝐶_𝑘𝜋^4/3Re_Δ^−11/2Re_Δ^1/2. (7iii)Fureby [15] proposed a model (denoted as FSDF) which can be written as: Σ_gen=Γ𝑢^′_Δ𝑆_𝐿^𝐷−2|∇𝑐|, (8) where Γ is given by (5ii) and 𝐷=2.05/(𝑢^′_Δ/𝑆_𝐿+1)+2.35/(𝑆_𝐿/𝑢^′_Δ+1) [19]. Knikker et al. [16] proposed a model (denoted as FSDK) as Σ_gen=Δ𝜂_𝑖^𝛽𝑘|∇𝑐|, (9) where 𝜂_𝑖=3𝛿_𝑧 and 𝛽_𝑘 is dynamically evaluated as 𝛽_𝑘=[log⟨|∇𝑐|⟩−log⟨|∇c|̂⟩]/log𝛾 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 Σ_gen 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 ⟨Σ_gen⟩ represents the total flame surface area within the domain and therefore should be independent of Δ. Thus the model predictions of ⟨Σ_gen⟩ should not change with Δ.

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

Criterion 3. The correlation coefficient between the modelled and actual values of Σ_gen should be as close as possible to unity, in order to capture correctly the local strain rate and the curvature effects on Σ_gen 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 36.6𝛿_th×24.1𝛿_th×24.1𝛿_th was discretised using a Cartesian mesh of size 345×230×230 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 𝛿_th. 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 (Pr=0.7), the Zel’dovich number (𝛽_𝑍=𝑇_ac(𝑇_ad−𝑇₀)/𝑇²_ad=6.0), and the ratio of specific heats (𝛾_𝑔=𝐶_𝑃/𝐶_𝑉=1.4), where 𝑇_ac,𝐶_𝑃, and 𝐶_𝑉 are the activation temperature and the specific heats at constant pressure and constant volume, respectively. The initial values of 𝑢^′/𝑆_𝐿, 𝑙/𝛿_th, Da, Ka, and Re_𝑡 are listed in Table 1. In cases  A, C, and E (B, C, and D), the values of 𝑢^′/𝑆_𝐿 and 𝑙/𝛿_th were chosen to vary Re_𝑡 by changing Ka (Da) while keeping Da(Ka) 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 Ka 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 𝑡_𝑐=𝛿_th/𝑆_𝐿. This time is equal to 2.0𝑡_𝑓 in case D, 3.0𝑡_𝑓 in cases A, C, and E, and 4.34𝑡_𝑓 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 Σ_gen [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 𝑢^′/𝑆_𝐿, 𝑙/𝛿_th, and 𝛿_th/𝜂 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 𝛿_th as Ka=(𝑢^′/𝑆_𝐿)^1.5(𝑙/𝛿_th)^−1/2, whereas Ka in this paper is defined in terms of the Zel’dovich flame thickness 𝛿_𝑧=𝛼_𝑇0/𝑆_𝐿 as Ka=(𝑢^′/𝑆_𝐿)^1.5(𝑙𝑆_𝐿/𝛼_𝑇0)^−1/2).

Table 1 

List of initial simulation parameters and nondimensional numbers.

In the present analysis, the DNS data is explicitly filtered using 𝐺(^→𝑟)=(6/𝜋Δ²)^3/2exp(−6^→𝑟⋅^→𝑟/Δ²) [5] to obtain the relevant filtered quantities. The results will be presented for Δ ranging from Δ=4Δ_𝑚≈0.4𝛿_th to Δ=24Δ_𝑚≈2.4𝛿_th, where Δ_𝑚 is the DNS mesh size (Δ_𝑚≈0.1𝛿_th). 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 0.4𝛿_th, where the flame is partially resolved, up to 2.4𝛿_th 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 Δ=0.4𝛿_th≈0.7𝛿_𝑧) to well within the TRZ regime (for Δ≥0.5𝛿_th≈𝛿_𝑧) 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 𝑡_sim=1.0𝛿_th/𝑆_𝐿 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 𝑢^′/𝑆_𝐿∼Re_𝑡^1/2/Da^1/2∼Re_𝑡^1/4Ka^1/2 (see Table 1). Moreover, the contours of 𝑐 representing the preheat zone (i.e., 𝑐<0.5, see colour scale) are much more distorted than those representing the reaction zone (i.e., 0.7<𝑐<0.9). This tendency becomes more prevalent with the increase in the Karlovitz number Ka, since the scale separation between 𝛿_th 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)

Figure 1 

Contours of 𝑐 in the 𝑥1-𝑥2 midplane at time 𝑡=𝛿th/𝑆𝐿 for: (a–e) cases A–E.

4.2. Behaviour of 𝐷 and 𝜂_𝑖 in Response to Re_𝑡

The power-law Ξ=Σ_gen/|∇𝑐|=(𝜂₀/𝜂_𝑖)^𝐷−2 produces the expression Σlog^gen⟨|∇𝑐=()Δ|⟩𝐷−2log𝜂_𝑖, (10) 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, ⟨Σ_gen⟩ indicates the total flame surface area in the computational domain, thus remaining independent of Δ. As a result, log(⟨Σ_gen⟩/⟨|∇𝑐|⟩) increases with increasing Δ, which can be substantiated from Figures 2(a)–2(e). The variation of log(⟨Σ_gen⟩/⟨|∇𝑐|⟩) with log(Δ/𝛿_𝑧) is linear when Δ≫𝛿_𝑧 but becomes nonlinear for Δ≪𝛿_𝑧. The slope of the best-fit straight line representing the greatest slope of the variation of log(⟨Σ_gen⟩/⟨|∇𝑐|⟩) with log(Δ/𝛿_𝑧) gives 𝐷 while the intersection point of this straight line with the abscissa gives log𝜂_𝑖. It has been found that 𝜂_𝑖/𝛿_𝑧 remains about 2.0 for all cases (i.e., 𝜂_𝑖/𝛿_𝑧≈1.785) and that 𝛿_th=1.785𝛿_𝑧 for the present thermochemistry. This is qualitatively consistent with previous experimental [16, 27] and computational [8] findings.

(a)

(b)

(c)

(d)

(e)

Figure 2 

Variation of ⟨Σgen⟩/⟨|∇𝑐|⟩ () with Δ/𝛿𝑧 on a log-log plot for (a–e) cases A–E. The prediction of ⟨Σgen⟩/⟨|∇𝑐|⟩=(Δ/𝜂𝑖)𝐷−2 () with 𝜂𝑖 obtained from DNS and 𝐷 according to (11) () is also shown.

Figures 2(a)–2(e) demonstrate that 𝐷 is greater for flames with higher Re_𝑡 and that 𝐷 attains an asymptotic value of 7/3 for flames with Re_𝑡≥50 (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 𝑢^′/𝑆_𝐿∼Re_𝑡^1/4Ka^1/2∼Re_𝑡^1/2/Da^1/2. Figure 2 suggests that 𝐷 is expected to increase with increasing Re_𝑡∼Da²Ka² before assuming an asymptotic value when either Da or Ka 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, 𝐷=2+1/3erf(2Ka), predicts 𝐷=7/3 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 1𝐷=2+3erf(3Ka)1−exp−0.1Re_𝑡𝐴_𝑚^1.6, (11) where 𝐴_𝑚≈7.5 is a model parameter. The prediction of ⟨Σ_gen⟩/⟨|∇𝑐|⟩=(Δ/𝜂_𝑖)^𝐷−2 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 Re_𝑡 and Ka (e.g., cases  D and E), according to Kerstein [20].

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

Note that Re_𝑡 and Ka in (11) were evaluated here based on 𝑢^′/𝑆_𝐿 and 𝑙/𝛿_th 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 𝑘_Δ=(𝑢_𝑖𝑢_𝑖−𝑢_𝑖𝑢_𝑖)/2 (i.e., 𝑢^′_Δ=2𝑘_Δ/3) following previous studies [6, 8, 10–12, 15, 18]. The local Karlovitz number Ka_Δ can be evaluated as Ka_Δ=𝐶_Ka(𝑘_Δ/𝑆_𝐿)^3/2(𝛿_𝑍/Δ)^1/2 where 𝐶_Ka is a model parameter.

A power-law-based expression for Σ_gen is proposed here based on the observed behaviours of 𝐷 and 𝜂_𝑖 in Figure 2: Σ_gen=|∇𝑐Δ|1−𝑓+𝑓𝜂_𝑖^𝐷−2, (12) where 𝑓 is a bridging function which increases monotonically from 0 for small Δ (i.e., Δ/𝛿_th→0 or Δ≪𝛿_th) to 1 for large Δ (i.e., Δ≫𝜂_𝑖 or Δ≫𝛿_th). The expression given by (12) will be denoted as FSDNEW. Equation (12) ensures that Σ_gen approaches |∇𝑐|(Δ/𝜂_𝑖)^𝐷−2 (i.e., Σ_gen=|∇𝑐|(Δ/𝜂_𝑖)^𝐷−2) for large Δ (i.e., Δ≫𝜂_𝑖 or Δ≫𝛿_𝑡ℎ) and at the same time Σ_gen approaches |∇𝑐| (i.e., lim_Δ→0Σ_gen=lim_Δ→0|∇𝑐|=|∇𝑐|) for small Δ (i.e., Δ/𝛿_th→0). It has been found from Figure 2 that Σ_gen≈|∇𝑐| provides better agreement with Σ_gen obtained from the DNS data for Δ≤0.8𝜂_𝑖, whereas the power-law Σ_gen=|∇𝑐|(Δ/𝜂_𝑖)^𝐷−2 starts to predict Σ_gen more accurately for Δ≥1.2𝜂_𝑖. Based on this observation, the bridging function 𝑓 is taken to be 𝑓=1/[1+exp{−60(Δ/𝜂_𝑖−1.0)}] which ensures a smooth transition 0.8𝜂_𝑖<Δ<1.2𝜂_𝑖. As 𝜂_𝑖 is found to scale with 𝛿_𝑧 (i.e., 𝜂_𝑖=𝛿_th=1.785𝛿_𝑧 according to the present thermochemistry), 𝜂_𝑖 in (12) is taken to be the thermal flame thickness 𝛿_th. In the context of LES, 𝐷 needs to be evaluated based on local quantities, which is achieved by replacing Ka and Re_𝑡 in (11) with their local values Ka_Δ=𝐶_Ka(𝑘_Δ/𝑆_𝐿)^3/2(Δ/𝛿_𝑧)^−1/2 and Re_𝑡Δ=𝐶_Re(𝜌₀𝑢^′_ΔΔ/𝜇₀), respectively. The choice of model constants 𝐶_Ka=6.6 and 𝐶_Re≈4.0 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 ⟨Σ_gen⟩ can be characterised using a percentage error (PE) ΣPE=^genmodel−Σ_genΣ_gen×100, (13) where ⟨Σ_gen⟩_model is the model prediction of ⟨Σ_gen⟩. Results for the PE for a range of Δ are shown in Figure 3, which demonstrate that the models denoted by FSDA (see (5i) and (5ii) and FSDC (see (6)) overpredict ⟨Σ_gen⟩ for all the cases and that the level of overprediction increases with increasing Δ. The FSDW model (see (4)) also overpredicts the value of ⟨Σ_gen⟩, although the level of overprediction for Δ≫𝛿_th is smaller for this model, especially for the cases with higher Re_𝑡 (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 (7i), (7ii), and (7iii) remains small for the small and moderate Re_𝑡 cases (i.e., cases  A–C), although the FSDCH model slightly overpredicts the value of ⟨Σ_gen⟩ for Δ≫𝛿_th for the higher Re_𝑡 cases (i.e., cases  D and E). Replacing 𝛿_𝑐=4𝜇₀/𝜌₀𝑆_𝐿 (i.e., 𝛿_𝑐=2.8𝛿_𝑧 for the present thermochemistry) by 𝛿_𝑧=𝐷₀/𝑆_𝐿 in (7i)–(7iii) leads to a deterioration of the performance of the FSDCH model especially for higher Re_𝑡 cases where it starts to overpredict the value of ⟨Σ_gen⟩ for Δ≫𝛿_th. However, the performance of the FSDCH model even with 𝛿_𝑧 in (7i)–(7iii) remains better than the FSDA and FSDC models for all cases considered here. The FSDK model (see (8) underpredicts the value of ⟨Σ_gen⟩ for all Δ for all cases. However, the level of underprediction decreases for larger Δ.

(a)

(b)

(c)

(d)

(e)

Figure 3 

Percentage error of the model prediction from ⟨Σgen⟩ obtained from DNS for LES filter widths Δ=4Δ𝑚=0.4𝛿th; Δ=8Δ𝑚=0.8𝛿th; Δ=12Δ𝑚=1.2𝛿th; Δ=16Δ𝑚=1.6𝛿th; Δ=20Δ𝑚=2.0𝛿th and Δ=24Δ𝑚=2.4𝛿th for (a–e) cases A–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 ⟨Σ_gen⟩ slightly for small and moderate values of Re_𝑡 (i.e., cases  A–C), but the prediction of this model remains very close to the DNS result for high values of Re_𝑡 (cases  D and E). The PEs for the FSDCH, FSDNEW and FSDF models remain comparable. Note that Σ_gen should approach |∇𝑐| (i.e., lim_{𝑢′Δ→0}Σ_gen=lim_{𝑢′Δ→0}|∇𝑐|=|∇𝑐|) when 𝑢^′_Δ vanishes because the flow tends to be fully resolved (i.e., lim_Δ→0𝑢^′_Δ=0 and lim_Δ→0Σ_gen=|∇𝑐|). Although the FSDF model performs well for all Δ for all the cases considered here, Σ_gen does not tend to |∇𝑐| as 𝑢^′_Δ approaches zero but instead predicts a finite value close to zero. Hence the FSDF model underpredicts the value of ⟨Σ_gen⟩ 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): Σ_gen=|∇𝑐|Γ𝑢1−𝑓+𝑓^′Δ𝑆_𝐿^𝐷−2, (14) where 𝑓=1/[1+exp{−60(Δ/𝛿_th−1.0)}] is a bridging function as before, the efficiency function Γ is given by (5ii), and 𝐷=2.05/(𝑢^′_Δ/𝑆_𝐿+1)+2.35/(𝑆_𝐿/𝑢^′_Δ+1) [19]. Equation (14) ensures that Σ_gen becomes exactly equal to |∇𝑐| when the flow is fully resolved (i.e., Δ≪𝜂_𝑖 or Δ→0) where 𝑢^′_Δ also vanishes (i.e., lim_Δ→0𝑢^′_Δ=0). 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., lim_Δ→0Σ_gen=lim_Δ→0|∇𝑐|=|∇𝑐|) for very small Δ (i.e., Δ→0). Note that the combination of parameterisation of 𝐷 and Γ according to [19] and (5ii), 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 Ξ=(𝜂₀/𝜂_𝑖)^𝐷−2 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 3𝛿_𝑧 in this model whereas 𝜂_𝑖=1.785𝛿_𝑧 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, Σ_gen vanishes when Δ→0 according to the FSDK model, whereas Σ_gen should approach |∇𝑐| when Δ→0 (i.e., lim_{𝑢′Δ→0}Σ_gen=lim_{𝑢′Δ→0}|∇𝑐|=|∇𝑐|). 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 𝐾=(1/𝛿𝐴)𝑑(𝛿𝐴)/𝑑𝑡=𝑎_𝑇+𝑆_𝑑∇⋅𝑁^→ 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 (𝐾)_𝑠=0, 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 ⟨Σ_gen⟩ 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 ⟨Σ_gen⟩ 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 Re_𝑡 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 Σ_gen conditionally averaged on 𝑐 are shown in Figure 4 for Δ=8Δ_𝑚=0.8𝛿_th and in Figure 5 for Δ=24Δ_𝑚=2.4𝛿_th. 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 Σ_gen for flames with higher Re_𝑡 (e.g., cases  D and E). However, the extent of this overprediction is relatively lower in the low Re_𝑡 cases (e.g., cases  A and B). (ii)The FSDW model tends to predict a peak value of Σ_gen at 𝑐>0.6, whereas the peak value of Σ_gen obtained from DNS is attained close to 𝑐≈0.6 for all the cases.(iii)The FSDK model tends to underpredict the value of Σ_gen for all cases. The physical explanations provided earlier for the underprediction of ⟨Σ_gen⟩ by the FSDK model (see Figure 3) are also valid here.(iv)Models FSDA, FSDC, FSDCH, FSDF, and FSDNEW all tend to capture the Σ_gen 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 Σ_gen models for Δ=2.4𝛿_th than for Δ=0.8𝛿_th. The following observations can be made from Figure 5. (i)Models FSDW, FSDA, and FSDC all overpredict the value of Σ_gen, and the overprediction increases with increasing Re_𝑡. The FSDCH model captures the behaviour of Σ_gen for small values of Re_𝑡 (e.g., cases  A–C), but it overpredicts the value of Σ_gen for higher Re_𝑡 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 Σ_gen, from DNS occurs at 𝑐≈0.5 for all cases.(iii)Models FSDF, FSDK, FSDNEW and MFSDF predict Σ_gen satisfactorily throughout the flame brush.(iv)The difference between the models MFSDF and FSDF seems to be very small, and both predict Σ_gen satisfactorily.

Figure 4 

Variation of the mean values of Σgen×𝛿𝑧 conditional on 𝑐 across the flame brush for Δ=8Δ𝑚=0.8𝛿th according to the DNS () and model () for cases A (1st row), B (2nd row), C (3rd row), D (4th row), and E (5th row).

Figure 5 

Variation of the mean values of Σgen×𝛿𝑧 conditional on 𝑐 across the flame brush for Δ=24Δ𝑚=2.4𝛿th according to the DNS () and model () for cases A (1st row), B (2nd row), C (3rd row), D (4th row), and E (5th row).

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., Ξ=1+1.24𝑐̃𝑢^′_Δ/𝑆_𝐿Re_𝜂). The models FSDW, FSDA, and FSDC underestimate the destruction of flame surface area in the TRZ regime by ignoring Σ_gen destruction arising due to −[𝐷_𝑐(∇⋅𝑁^→)²]_𝑠 which eventually leads to the overprediction of Σ_gen. 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 Ka_Δ and turbulent Reynolds number Re_𝑡Δ 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., Ka_Δ=6.66(𝑘_Δ/𝑆_𝐿)^3/2(Δ/𝛿_𝑧)^−1/2 and Re_𝑡Δ=4.0(𝜌₀𝑢^′_ΔΔ/𝜇₀)).

Criterion 3. The correlation coefficients between the DNS result and the model prediction for Σ_gen are shown in Figure 6 for 0.1≤𝑐≤0.9. For 𝑐<0.1 and 𝑐>0.9,Σ_gen 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 Re_𝑡. With increasing Re_𝑡, 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 Σ_gen may be advantageous in terms of capturing the correct strain rate and the curvature dependences of Σ_gen, provided the unclosed terms of Σ_gen transport equation are accurately modelled.

(a)

(b)

(c)

(d)

(e)

Figure 6 

Correlation coefficients between modelled and actual values of Σgen in the 𝑐 range 0.1≤𝑐≤0.9 for Δ=4Δ𝑚=0.4𝛿th; Δ=8Δ𝑚=0.8𝛿th; Δ=12Δ𝑚=1.2𝛿th; Δ=16Δ𝑚=1.6𝛿th; Δ=20Δ𝑚=2.0𝛿th and Δ=24Δ𝑚=2.4𝛿th for (a–e) cases A–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., Δ→0 and 𝑢^′_Δ→0). 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 Σ_gen for flames within the TRZ regime has been assessed based on DNS of turbulent premixed flames over a range of different values of Re_𝑡. It has been found that the fractal dimension 𝐷 increases with increasing Re_𝑡 and Ka before reaching an asymptotic value. By contrast, the inner cut-off scale 𝜂_𝑖 is not significantly affected within the range of Re_𝑡 considered. The observed behaviours of 𝐷 and 𝜂_𝑖 have been incorporated into a new power-law model for Σ_gen 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 Re_𝑡. Based on this assessment, models have been identified which predict Σ_gen 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 Re_𝑡 without the effects of detailed chemistry and transport. Thus, three-dimensional DNS with detailed chemistry along with experimental data for higher values of Re_𝑡 will be necessary for a more comprehensive analysis.

Acknowledgment

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