A Priori Assessment of Algebraic Flame Surface Density Models in the Context of Large Eddy Simulation for Nonunity Lewis Number Flames in the Thin Reaction Zones Regime

Abstract

The performance of algebraic flame surface density (FSD) models has been assessed for flames with nonunity Lewis number (Le) in the thin reaction zones regime, using a direct numerical simulation (DNS) database of freely propagating turbulent premixed flames with Le ranging from 0.34 to 1.2. The focus is on algebraic FSD models based on a power-law approach, and the effects of Lewis number on the fractal dimension D and inner cut-off scale 𝜂_𝑖 have been studied in detail. It has been found that D is strongly affected by Lewis number and increases significantly with decreasing Le. By contrast, 𝜂_𝑖 remains close to the laminar flame thermal thickness for all values of Le considered here. A parameterisation of D is proposed such that the effects of Lewis number are explicitly accounted for. The new parameterisation is used to propose a new algebraic model for FSD. The performance of the new model is assessed with respect to results for the generalised FSD obtained from explicitly LES-filtered DNS data. It has been found that the performance of the most existing models deteriorates with decreasing Lewis number, while the newly proposed model is found to perform as well or better than the most existing algebraic models for FSD.

1. Introduction

Reaction rate closure based on flame surface density (FSD) is one of the most popular approaches to combustion modelling in turbulent premixed flames [1–11]. In the context of LES the generalised FSD (Σ_gen) is defined as follows [3–11]: Σ_gen=||∇𝑐||.(1) where the overbar denotes the LES filtering operation. The reaction progress variable 𝑐 may be defined in terms of a reactant mass fraction 𝑌_𝑅, for example, 𝑐=(𝑌_𝑅0−𝑌_𝑅)/(𝑌_𝑅0−𝑌_𝑅∝) such that 𝑐 rises monotonically from zero in fresh reactants (subscript 0) to unity in fully burned products (subscript ∝).

In the context of LES, several models have been proposed for the wrinkling factor Ξ [12–16], which is often used in the context of thickened flame modelling [13, 14]. The wrinkling factor Ξ is closely related to Σ_gen according to [12–16]:Ξ=Σ^gen||∇𝑐||.(2a)

Often, Ξ is expressed in terms of a power-law expression [7, 9, 13, 14] Ξ=(𝜂₀/𝜂_𝑖)^𝐷−2 in which 𝜂₀ and 𝜂_𝑖 are the outer and inner cut-off scales and 𝐷 is the fractal dimension. This leads to a power-law expression for Σ_gen as: Σ_gen=Ξ||∇𝑐||=Δ𝜂_𝑖^𝐷−2||∇𝑐||,(2b)where, for LES, the outer cut-off scale 𝜂_𝑜 is taken to be equal to the filter width Δ. According to Peters [17], 𝜂_𝑖 scales with the Gibson length scale 𝐿_𝐺=𝑆^3𝐿/𝜀 in the corrugated flamelets (CF) regime, and with the Kolmogorov length scale 𝜂=(𝜈³/𝜀)^1/4 in the thin reaction zones (TRZ) regime. Here, 𝑆_𝐿 is the unstrained laminar burning velocity, 𝑣 is the kinematic viscosity in the unburned gas, and 𝜀 is the dissipation rate of turbulent kinetic energy. Experimental analyses by Knikker et al. [7] and Roberts et al. [18] indicated that 𝜂_𝑖 scales with the Zel’dovich flame thickness 𝛿_𝑍=𝛼_𝑇0/𝑆_𝐿, where 𝛼_𝑇0 is the thermal diffusivity in unburned gases. A recent a priori DNS analysis [9] demonstrated that 𝜂_𝑖 scales with 𝐿_𝐺 and 𝜂 for the CF and TRZ regimes, respectively, as suggested by Peters [17]. However, 𝜂_𝑖 is also found to scale with thermal flame thickness 𝛿_th in both the CF and TRZ regimes [9]. North and Santavicca [19] parameterised 𝐷 in terms of the root-mean-square (rms) turbulent velocity fluctuation 𝑢′ as: 𝐷=2.05/(𝑢′/𝑆_𝐿+1)+2.35/(𝑆_𝐿/𝑢′+1), whereas Kerstein [20] suggested that 𝐷 increases from 2 to 7/3 for increasing values of 𝑢′/𝑆_𝐿, where 𝐷=7/3 is associated with the material surface.

Since combustion is set to remain a major practical means of energy conversion for the foreseeable future, it has become necessary to find novel ways to reduce carbon emissions from relatively conventional combustion systems. One such approach is the use of hydrogen-blended hydrocarbon fuels in IC engines, aeroengines, and furnaces. Increased abundance of fast diffusing species such as H and H₂ leads to significant effects of differential diffusion of heat and mass in hydrogen-blended flames [21, 22], whereas these effects are relatively weaker in conventional hydrocarbon flames [22, 23]. The differential rates of thermal and mass diffusion in premixed flames are often characterised by the Lewis number Le which is defined as the ratio of the thermal diffusivity to mass diffusivity (i.e., Le=𝛼_𝑇/𝐷_𝑐). Assigning a global characteristic value of Le is not straightforward since many species with different individual values of Le are involved in actual combustion. Often the Lewis number of the deficient reactant species is used as the characteristic Le [21, 24–28] and this approach has been adopted here. It is worth noting that, to date, most FSD-based modelling has been carried out for unity Lewis number flames (e.g., [1–11]) and the effects of differential diffusion of heat and mass on the statistical behaviour of FSD have rarely been addressed [28]. More specifically the effects of Le on 𝐷 and 𝜂_𝑖 have not yet been analysed in detail, or in the context of power-law FSD reaction rate models. Moreover, most algebraic models for Σ_gen have been proposed for the CF regime where the effects of Le are not accounted for. Thus, it is important to assess the performance of existing models for combustion in the TRZ regime with nonunity Lewis number.

The present study aims to bridge this gap in the existing literature. In this respect the main objectives of the work are as the following.(i)To understand the effects of Lewis number on 𝐷 and 𝜂_𝑖 in the context of LES modelling.(ii)To assess the performance of existing wrinkling factor-based algebraic models of FSD in the context of LES for flames with nonunity global Lewis number based on a priori DNS analysis.(iii) To identify or develop a power-law-based algebraic model for FSD in the context of LES which is capable of predicting the correct behaviour of FSD even for nonunity Lewis number flames.

The rest of the paper is organised as follows. An overview of the different algebraic FSD models considered here are presented in the next section. This will be followed by a brief discussion of the numerical implementation. Following this, results will be presented and subsequently discussed. Finally the main findings will be summarised and conclusions will be drawn.

2. Overview of Power-Law-Based FSD Models

A model for Ξ suggested by Angelberger et al. [4] (FSDA model) can be written in terms of Σ_gen as follows:Σ_gen=1+𝑎Γ𝑢^Δ𝑆_𝐿||∇𝑐||,(3a) where 𝑎=1.0 is a model parameter, 𝑢^Δ=2𝑘̃_Δ/3 is the subgrid turbulent velocity fluctuation, 𝑘̃_Δ=(𝑢_𝑖𝑢_𝑖−𝑢_𝑖̃𝑢_𝑖̃)/2 is the subgrid turbulent kinetic energy and 𝑄=𝜌𝑄/𝜌 denotes the Favre-filtered value of a general quantity 𝑄. In (3a), Γ is an efficiency function which is given by: Γ=0.75exp−1.2𝑢^Δ𝑆_𝐿^−0.3⋅𝛿Δ_𝑧^2/3(3b)Weller et al. [12] also presented an algebraic model for Ξ, which can be recast in the form (FSDW model): Σ_gen=[1+2𝑐̃(Θ−1)]||∇𝑐||,(4) where Θ=1+0.62𝑢^Δ/𝑆_𝐿Re_𝜂 and Re_𝜂=𝑢^Δ⋅𝜂/𝜈 with 𝜂 and 𝜌₀ denoting the Kolmogorov length scale and unburned gas density respectively. Colin et al. [13] proposed an algebraic model for Ξ, which can be expressed in terms of FSD (FSDC model) as: Σ_gen=1+𝛼Γ𝑢^Δ𝑆_𝐿||∇𝑐||,(5) where Γ is given by (3b), 𝛼=𝛽×2ln(2)/[3𝑐_𝑚𝑠(Re_𝑡1/2−1)] with Re_𝑡=𝜌₀𝑢′𝑙/𝜇₀, where 𝜇₀ is the unburned gas viscosity and 𝑙 is the integral length scale, 𝛽=1.0 and 𝑐_𝑚𝑠=0.28. The FSDC model requires three input parameters, namely 𝑢^Δ/𝑆_𝐿, Δ/𝛿_𝑧, and Re_𝑡. Charlette et al. [14] reduced the input parameters to only 𝑢^Δ/𝑆_𝐿 and Δ/𝛿_𝑧 by using (FSDCH model): Σ_gen=1+minΔ𝛿_𝑧,Γ_Δ𝑆𝑢^Δ_𝐿^𝛽1||∇𝑐||,(6) with the efficiency functionΓ_Δ=(𝑓^−𝑎1_𝑢+𝑓^−𝑎1_Δ)^−1/𝑎1^−𝑏1+𝑓^−𝑏1_Re^−1/𝑏1,(7a) where Re_Δ=𝑢^ΔΔ/𝜈 and with model constants 𝑏₁=1.4, 𝛽₁=0.5, 𝐶_𝑘=1.5, and functions 𝑎₁, 𝑓_𝑢, 𝑓_Δ, and 𝑓_Re are defined by: 𝑎₁=0.60+0.20exp−0.1𝑢^Δ𝑓_𝑢𝑆_𝐿−0.20exp−0.01𝛿Δ_𝑧,=427110𝐶^𝑘^1/25185𝐶^𝑘𝑆𝑢^Δ_𝐿²𝑓_Δ,=27110𝐶^𝑘𝜋^4/3Δ𝛿_𝑧^4/3𝑓_Re−1^1/2,=955exp(−1.5𝐶^𝑘𝜋^4/3Re^Δ−1)^1/2Re^Δ1/2.(7b) Knikker et al. [7] proposed a model for Σ_gen (FSDK model) as: Σ_gen=Δ𝜂_𝑖^𝛽𝑘||∇𝑐||,(8) where the inner cut-off scale 𝜂_𝑖 is taken to be 𝜂_𝑖=3𝛿_𝑧 and 𝛽_𝑘 is estimated based on a dynamic formulation as 𝛽_𝑘=[log⟨|∇𝑐|⟩−log⟨|∇𝑐̂|⟩]/log𝛾, where 𝑐̂‌ denotes the reaction progress variable at the test filter level γΔ. Fureby [16] proposed a model for Ξ which can be written in terms of Σ_gen (FSDF model) as: Σ_gen=Γ𝑢^Δ𝑆_𝐿^𝐷−2⋅||∇𝑐||,(9) where Γ is given by (3b), and 𝐷 is specified according to the parameterisation 𝐷=2.05/(𝑢^Δ/𝑆_𝐿+1)+2.35/(𝑆_𝐿/𝑢^Δ+1) [19].

In the present study, the performance of each algebraic model described above is assessed with respect to Σ_gen obtained from DNS. There are three requirements for each model. Firstly, the volume-averaged value of Σ_gen represents the total flame surface area, and therefore this quantity should not change with Δ. Secondly, the model should be able to capture the correct variation of the averaged value of Σ_gen conditional on 𝑐 across the flame brush. Thirdly, the correlation coefficient between the modelled and actual values of Σ_gen should be as close to unity as possible in order to capture the effects of local strain rate and curvature on Σ_gen.

3. Numerical Implementation

For the purposes of the analysis, a DNS database of three-dimensional turbulent premixed flames has been generated using the compressible DNS code SENGA [29]. Until recently most combustion DNS was carried out either in three dimensions with simplified chemistry or in two dimensions with detailed chemistry due to the limitations of available computational power. Although it is now possible to carry out three-dimensional DNS with detailed chemistry, such computations remain extremely expensive [30] and are not practical for a parametric study as in the present case. Thus three-dimensional DNS with single-step Arrhenius type chemistry has been used in the present study in which the effects of Lewis number are to be investigated in isolation.

For the present DNS database, the computational domain is considered to be a cube of size 24.1𝛿_th×24.1𝛿_th×24.1𝛿_th, which is discretised using a uniform grid of 230×230×230. The grid spacing is determined by the flame resolution, and in all cases, about 10 grid points are kept within 𝛿_th=(𝑇_ad−𝑇₀)/max|∇𝑇|_𝐿, where 𝑇_ad,𝑇₀ and 𝑇 are the adiabatic flame, unburned reactant and instantaneous dimensional temperatures respectively, and the subscript 𝐿 is used to refer to unstrained planar laminar flame quantities. The boundaries in the direction of mean flame propagation are taken to be partially nonreflecting and are specified using the Navier Stokes Characteristic Boundary Conditions formulation [31], while boundaries in the transverse direction were taken to be periodic. A 10th order central difference scheme was used for spatial discretisation for internal grid points and the order of differentiation gradually decreases to a one-sided second-order scheme at non-periodic boundaries [29]. A low storage 3rd-order Runge-Kutta scheme [32] is used for time advancement. The turbulent velocity field is initialised by using a standard pseudo-spectral method [33], and the flame is initialised using an unstrained planar steady laminar flame solution.

The initial values of 𝑢′/𝑆_𝐿 and 𝑙/𝛿_th for all the flames considered here are shown in Table 1 along with the values of heat release parameter 𝜏=(𝑇_ad−𝑇₀)/𝑇₀, Damköhler number Da=𝑙𝑆_𝐿/𝑢′𝛿_th, Karlovitz number Ka=(𝑢′/𝑆_𝐿)^3/2(𝑙𝑆_𝐿/𝛼_𝑇0)^−1/2 and turbulent Reynolds number Re_𝑡=𝜌₀𝑢′𝑙/𝜇₀. For all cases Ka remains greater than unity, which indicates that combustion is taking place in the TRZ regime [17]. Standard values are taken for Prandtl number (Pr=0.7), ratio of specific heats (𝛾_𝐺=𝐶_𝑃/𝐶_𝑉=1.4), and the Zel’dovich number (𝛽_𝑍=𝑇_ac(𝑇_ad−𝑇₀)/𝑇²_ad=6.0), where 𝑇_ac is the activation temperature.

Table 1 

In all cases, statistics were collected after three eddy turn-over times (i.e., 3𝑡_𝑓=3𝑙/𝑢′), which corresponds to one chemical time scale (i.e., 𝑡_𝑐=𝛿_th/𝑆_𝐿). The turbulent kinetic energy and its dissipation rate in the unburned reactants ahead of the flame were slowly varying at 𝑡_sim=3.0𝑙/𝑢′ and the qualitative nature of the statistics was found to have remained unchanged since 𝑡=2.0𝑙/𝑢′ for all cases. By the time the statistics were extracted, the value of 𝑢′/𝑆_𝐿 in the unburned reactants ahead of the flame had decayed by about 50%, while the value of 𝑙/𝛿_th had increased by about 1.7 times, relative to their initial values. Further details on the flame-turbulence interaction of this DNS database may be found in [27, 28]. The present simulation time is short, but remains comparable to several studies [3, 8–10, 14, 34–37] which have contributed significantly to the fundamental understanding and modelling of turbulent premixed combustion in the past. The DNS data was explicitly filtered according to the integral 𝑄(^→𝑥)=∫𝑄(^→𝑥−^→𝑟)𝐺(^→𝑟)𝑑^→𝑟 using a Gaussian kernel given by the expression 𝐺(^→𝑟)=(6/𝜋Δ²)^3/2exp(−6^→𝑟⋅^→𝑟/Δ²). The results will be presented for Δ ranging from Δ=4Δ_𝑚≈0.4𝛿_th to Δ=24Δ_𝑚≈2.4𝛿_th, where Δ_𝑚 is the DNS grid spacing (Δ_𝑚≈0.1𝛿_th). These filter sizes are comparable to the range of Δ used in a priori DNS analysis in several previous studies [3, 8–10, 14], and span a useful range of length scales from Δ comparable to 0.4𝛿_th≈0.8𝛿_𝑧, where the flame is partially resolved, up to 2.4𝛿_th≈4.8𝛿_𝑧, where the flame becomes fully unresolved and Δ is comparable to the integral length scale. For these filter widths, the underlying combustion process ranges from the “laminar flamelets-G DNS” [38] combustion regime (for Δ=0.4𝛿_th≈0.8𝛿_𝑧) to well within the TRZ regime (for Δ≥0.5𝛿_th≈𝛿_𝑧) on the regime diagrams by Pitsch and Duchamp de Lageneste [38] and Düsing et al. [39]. However, these regime diagrams have been proposed based on scaling arguments for unity Lewis number flames and the likely effects of nonunity Lewis number on these regime diagrams have yet to be ascertained. This topic is the subject of a separate investigation and will not be taken up in this paper.

4. Results and Discussion

4.1. Effects of Le on 𝐷 and 𝜂_𝑖

The power law expression (2b) for Σ_gen may be rewritten as: logΣ^gen||∇𝑐||=(𝐷−2)logΔ−(𝐷−2)log𝜂^𝑖,(10) where the angled brackets indicate a volume-averaging operation. The variation of ⟨Σ_gen⟩/⟨|∇𝑐|⟩ with the ratio (Δ/𝛿_𝑧) is shown in Figure 1 on a log-log plot for all the different Lewis number cases. The quantity ⟨Σ_gen⟩ denotes the total flame surface area which remains independent of filter size Δ. By contrast, the quantity ⟨|∇𝑐|⟩ denotes the resolved portion of the flame wrinkling, which decreases with increasing Δ. As a result, log[⟨Σ_gen⟩/⟨|∇𝑐|⟩] increases with increasing Δ. The variation of log[⟨Σ_gen⟩/⟨|∇𝑐|⟩] with log(Δ/𝛿_𝑧) is linear when Δ≫𝛿_𝑧 but becomes nonlinear for Δ≪𝛿_𝑧. The best-fit straight line representing the greatest slope of the linear variation has been used to obtain values of 𝐷 and 𝜂_𝑖. It has been found that 𝜂_𝑖/𝛿_𝑧 remains independent of Le, and for all cases 𝜂_𝑖 remains on the order of thermal flame thickness 𝛿_th (i.e., 𝜂_𝑖/𝛿_th≈1.0), which is about twice the Zel’dovich flame thickness 𝛿_𝑧 for the present thermochemistry (i.e., 𝜂_𝑖=1.79𝛿_𝑧≈𝛿_th). The scaling of the inner cut-off scale 𝜂_𝑖 with 𝛿_𝑧 is consistent with previous DNS [9] and experimental [7, 18] findings. Figure 1 shows that the slope of the linear region decreases with increasing Lewis number (i.e., in moving from case a to case e), which suggests that the fractal dimension 𝐷 decreases with increasing Le.

(a)  Le = 0.34

(b) Le = 0.6

(c) Le = 0.8

(d) Le = 1.0

(e) Le = 1.2

Figure 1 

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

Contours of reaction progress variable 𝑐 in the 𝑥₁−𝑥₂ midplane are shown in Figure 2 for all cases and show that the extent of flame wrinkling is significantly greater at lower Lewis number. The rate of flame area generation increases with decreasing Le, and this behaviour is particularly noticeable for the cases with Le = 0.34 and Le = 0.6 because of the occurrence of thermo-diffusive instabilities [21, 24–28]. This can be substantiated from values of the ratio of turbulent to laminar flame surface area 𝐴_𝑇/𝐴_𝐿 obtained by volume integration of |∇𝑐| (i.e., 𝐴=∫_𝜗|∇𝑐|𝑑𝜗). This produces the values 𝐴_𝑇/𝐴_𝐿=3.93, 2.66, 2.11, 1.84, and 1.76 for the cases with Le = 0.34, 0.6, 0.8, 1.0, and 1.2, respectively, at the time when statistics were extracted. The experimental findings of North and Santavicca [19] suggested that 𝐷 increases with increasing 𝑢′/𝑆_𝐿∼Re^𝑡1/4Ka^1/2, which indicates that 𝐷 is expected to have a dependence on both Re_𝑡 and Ka. Moreover, the analysis of Kerstein [20] suggested that 𝐷 is expected to assume an asymptotic value of 7/3 for large values of Re_𝑡 and Ka. The present findings indicate that Le also has an influence on 𝐷 in addition to Re_𝑡 and Ka, and that 𝐷 can assume values greater than 7/3 for flames with Le≪1.0 (see Figure 1). The Karlovitz number Ka dependence of 𝐷 for unity Lewis number flames has been analysed in detail by Chakraborty and Klein [9] and they parameterised 𝐷 as: 𝐷=2+(1/3)erf(2Ka), which does not account for the effects of Re_𝑡 and Le. The parameterisation proposed by Chakraborty and Klein [9] has been extended here by accounting for the effects of Karlovitz number, turbulent Reynolds number, and global Lewis number (i.e., Ka, Re_𝑡, and Le) according to the following: 𝐷=2+13erf(3.0Ka)1−exp−0.1R𝐴e^𝑡_𝑚^1.6Le^−0.45,(11) where 𝐴_𝑚≈7.5 is a model parameter. Further details on the basis of this parameterisation are given in Appendix A.

(a)

(b)

(c)

(d)

(e)

Figure 2 

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

The prediction of ⟨Σ_gen⟩/⟨|∇𝑐|⟩=(Δ/𝜂_𝑖)^𝐷−2 with 𝜂_𝑖 obtained from DNS and 𝐷 obtained from (11) is also shown in Figure 1, which indicates that (11) satisfactorily captures the best-fit straight line corresponding to the power law. It is worth noting that Re_𝑡 and Ka in (11) were evaluated for this purpose 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 subgrid turbulent kinetic energy as 𝑢^Δ=2𝑘̃_Δ/3 following previous studies [12, 15, 16]. The local Karlovitz number Ka_Δ can be evaluated as Ka_Δ=𝐶_Ka(√𝑘_Δ/𝑆_𝐿)^3/2(𝛿_𝑧/Δ)^1/2, where 𝐶_Ka is a model parameter. Similarly, the local turbulent Reynolds number Re_𝑡Δ can be evaluated using Re_𝑡Δ=𝐶_Re(𝜌₀𝑢^ΔΔ/𝜇₀). 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).

Based on the observed behaviour of 𝐷 and 𝜂_𝑖, a power-law expression for Σ_gen is proposed here (model FSDNEW): Σ_gen=||∇𝑐||(1−𝑓)+𝑓Δ𝜂_𝑖^𝐷−2,(12) where 𝑓 is a bridging function which increases monotonically from zero for small Δ (i.e., Δ/𝛿_th→0 or Δ≪𝛿_th) to unity for large Δ (i.e., Δ≫𝜂_𝑖 or Δ≫𝛿_th). Equation (12) ensures that Σ_gen approaches |∇𝑐|(Δ/𝜂_𝑖)^𝐷−2 for large Δ and at the same time Σ_gen approaches |∇𝑐| (i.e., lim_Δ→0Σ_gen=lim_Δ→0|∇𝑐|=|∇𝑐|) for small Δ. It has been found that Σ_gen≈|∇𝑐| provides better agreement with Σ_gen obtained from DNS data for Δ≤0.8𝜂_𝑖, whereas the power-law Σ_gen=|∇𝑐|(Δ/𝜂_𝑖)^𝐷−2 starts to predict Σ_gen more accurately for Δ≥1.2𝜂_𝑖 (see Figure 1). Based on this observation, the bridging function 𝑓 is taken to be 𝑓=1/[1+exp{−60(Δ/𝜂_𝑖−1.0)}], which ensures a smooth transition between 0.8𝜂_𝑖<Δ<1.2𝜂_𝑖. As 𝜂_𝑖 is found to scale with 𝛿_𝑧 (i.e., 𝜂_𝑖≈1.79𝛿_𝑧≈𝛿_th according to the present thermochemistry), 𝜂_𝑖 in (12) is taken to be the thermal flame thickness 𝛿_th.

The performance of the various algebraic models for Σ_gen will be assessed next, using the model requirements stated earlier.

4.2. Performance of Models for the Volume-Averaged FSD ⟨Σ_gen⟩

The inaccuracy in the model predictions of ⟨Σ_gen⟩ can be characterised using a percentage error (PE): PE=Σ^gen_model−Σ_genΣ_gen×100,(13) where ⟨Σ_gen⟩_model is the volume-averaged value of the model prediction of ⟨Σ_gen⟩. Results for the PE for a range of filter size Δ are shown Figure 3. These demonstrate that the models denoted by FSDA (see (3a) and (3b)) and FSDC (see (5)) overpredict ⟨Σ_gen⟩ for all the Lewis number cases, and that the level of overprediction increases with increasing Δ. The FSDW model (see (4)) also overpredicts ⟨Σ_gen⟩, although the level of overprediction decreases for Δ≫𝛿_th, especially for cases with Le≥0.6 (i.e., cases B–E). The FSDC model has greater PE than both the FSDA and FSDW models for all Δ in the same cases. However, the FSDW model has the highest PE relative to both the FSDA and FSDC models for all Δ in the Le=0.34 case.

(a)

(b)

(c)

(d)

(e)

Figure 3 

Percentage error (13) 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; Δ=24Δ𝑚=2.4𝛿th for (a–e) cases A–E.

The FSDCH (6), FSDA, and FSDC models provide accurate predictions of ⟨Σ_gen⟩ at small values of Δ (i.e., Δ≪𝛿_𝑧) but they overpredict ⟨Σ_gen⟩ for large values of Δ (i.e., Δ≫𝛿_𝑧). The FSDF model (9) predicts accurately for small Δ, and marginally underpredicts for larger Δ, for cases with Le≥0.6. However, the FSDF model remains better than the FSDA, FSDC, FSDCH, and FSDW models. The FSDNEW model (12) provides an accurate prediction of ⟨Σ_gen⟩ for all filter sizes because this model is designed to do so for all values of Le. The PE for the FSDCH model remains small for cases with Le≈1.0 (i.e., cases C–E), although the FSDCH model overpredicts ⟨Σ_gen⟩ for Δ≫𝛿_th for cases with Le≪1 (i.e., cases A and B). The FSDK model (see (8)) underpredicts the value of ⟨Σ_gen⟩ for all Δ for all cases. However, the level of underprediction of the FSDK model decreases for larger Δ.

The PEs for the FSDF and FSDNEW models remain negligible in comparison to the PEs for all the other models. 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. This limitation of the FSDF model can be avoided using a modified form of (8) (MSFDF model): Σ_gen=||∇𝑐||(1−𝑓)+𝑓Γ𝑢^Δ𝑆_𝐿^𝐷−2,(14) where 𝑓=1/[1+exp{−60(Δ/𝛿_th−1.0)}] is a bridging function as before, the efficiency function Γ is given by (3b) 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). Figure 3 shows that the modification given by (14) does not appreciably alter the performance of (8) while ensuring the correct asymptotic behaviour. Note that the parameterisation of 𝐷 and Γ according to [19] and (3b), 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 the power-law Ξ=(𝜂₀/𝜂_𝑖)^𝐷−2 which is strictly valid only for filter sizes Δ which are sufficiently greater than 𝜂_𝑖 (i.e., Δ≫𝜂_𝑖), as can be seen from Figure 1. 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.79𝛿_𝑧 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 K=(1/𝛿𝐴)𝑑(𝛿𝐴)/𝑑𝑡=𝑎_𝑇+𝑆_𝑑∇⋅𝑁^→ represents the fractional rate of change of flame surface area A [1], 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 𝑆_𝑡) [8–10, 40, 41]: 𝑆_𝑟=𝑤̇𝜌||∇𝑐||,𝑆^𝑛=^→𝑁⋅∇𝜌𝐷𝑁_𝑐→∇𝑐𝑆𝜌||∇𝑐_𝑡||,=−𝐷_𝑐∇⋅𝑁^→.(15) It has been shown in several previous studies [5, 6, 8, 10, 25] that (𝑎_𝑇)_𝑠 remains positive throughout the flame brush and thus acts to generate flame surface area, whereas the contribution of curvature to stretch (𝑆_𝑑∇⋅𝑁^→)_𝑠=[(𝑆_𝑟+𝑆_𝑛)∇⋅𝑁^→]_𝑠−[𝐷_𝑐(∇⋅𝑁)²]_𝑠 is primarily responsible for flame surface area destruction. The equilibrium of flame surface area generation and destruction yields (K)_𝑠=0, which gives rise to [9]: 𝑎_𝑇_𝑠=−𝑆^𝑟+𝑆_𝑛∇⋅𝑁^→_𝑠+[𝐷_𝑐∇⋅𝑁^→)²_𝑠.(16)

The stretch rate induced by −[𝐷_𝑐(∇⋅𝑁^→)²]_𝑠 becomes the leading order sink term in the thin reaction zones regime [8–10, 42]. However, most algebraic models (e.g., FSDA, FSDC, FSDCH, 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 −[𝐷_𝑐(∇⋅𝑁^→)²]_𝑠 was ignored [4, 12–14]. As a result, these models underestimate the flame surface area destruction rate in the thin reaction zones regime, which leads to overprediction of ⟨Σ_gen⟩ for the FSDA, FSDC, FSDCH, and FSDW models.

The disagreement between the FSDF model prediction and DNS data originates principally due to the inaccuracy in estimating Γ and 𝐷, while the difference between the FSDK prediction and DNS data arises from inaccurate estimation of 𝜂_𝑖. Hence a more accurate estimation of Γ, 𝐷, and 𝜂_𝑖 will result in better performance of both the FSDF and FSDK models.

4.3. Performance of Models for the Variation of Σ_gen

It is important to assess the models based on their ability to capture the correct variation of Σ_gen with 𝑐 across the flame brush. The variation of mean Σ_gen conditionally averaged on 𝑐 is shown in Figure 4 for Δ=8Δ_𝑚=0.8𝛿_th and Figure 5 for Δ=24Δ_𝑚=2.4𝛿_th, respectively. These filter widths have been chosen since they correspond to Δ<𝜂_𝑖 and Δ>𝜂_𝑖 respectively. The following observations can be made from Figure 4 about the model predictions at Δ=8Δ_𝑚=0.8𝛿_th.(i)The models FSDA, FSDC, FSDCH, FSDF, and FSDNEW tend to capture the variation of the conditional mean value of Σ_gen with 𝑐 obtained from DNS data. The prediction of the MFSDF model remains comparable to that of the FSDF model for Δ=8Δ_𝑚=0.8𝛿_th.(ii)The FSDW model consistently overpredicts the conditional mean value of Σ_gen for all cases. The FSDW model also predicts a skewed shape, which fails to capture the trend predicted by DNS.(iii)The model FSDK underpredicts the conditional mean value of Σ_gen in all cases. The physical explanations provided earlier for the underprediction of ⟨Σ_gen⟩ by the FSDK model is also responsible for the underprediction seen here.

(a)

(b)

(c)

(d)

(e)

Figure 4 

Variation of mean values of Σgen×𝛿𝑧 conditional on 𝑐 across the flame brush for Δ=8Δ𝑚=0.8𝛿th according to DNS, FSDA, FSDC, FSDW, FSDCH, FSDK, FSDF, FSDNEW, and MFSDF predictions for (a–e) cases A–E.

(a)

(b)

(c)

(d)

(e)

Figure 5 

Variation of mean values of Σgen×𝛿𝑧 conditional on 𝑐 across the flame brush for Δ=24Δ𝑚=2.4𝛿th according to DNS, FSDA, FSDC, FSDW, FSDCH, FSDK, FSDF, FSDNEW, and MFSDF predictions for (a–e) cases A–E.

A comparison between Figures 4 and 5 reveals that the predictions of the various algebraic FSD models exhibit greater spread for Δ=24Δ_𝑚=2.4𝛿_th than in the case of Δ=8Δ_𝑚=0.8𝛿_th. The following observations can be made from Figure 5 about the model predictions at Δ=24Δ_𝑚=2.4𝛿_th.(i)Similar to Δ=8Δ_𝑚, the FSDW model predicts a peak at 𝑐> 0.6, whereas the peak value of conditionally averaged Σ_gen from DNS occurs at 𝑐≈ 0.5 for all the cases.(ii)The models FSDW, FSDA, FSDC, and FSDCH tend to overpredict the conditionally averaged value of Σ_gen and the level of the overprediction increases with decreasing Lewis number.(iii)The models FSDF, FSDK, FSDNEW, and MFSDF tend to predict the conditionally averaged value of Σ_gen satisfactorily throughout the flame brush.(iv)The difference in the predictions of the models MFSDF, and FSDF seem to be very small for all the flames considered here.

The inaccuracy in the predictions of the mean value of Σ_gen conditional on 𝑐 can be characterised once again using a percentage error (PE₂): PE₂=Σ^MODELcond−Σ^DNScondΣ^max_cond×100,(17) where Σ^MODELcond and Σ^DNScond are the mean values of Σ_gen conditional on 𝑐 as obtained from model prediction and DNS respectively, and Σ^max_cond is the maximum value of conditionally averaged Σ_gen obtained from DNS. The error in the model prediction according to (16) is shown in Figure 6 for filter size Δ=8Δ_𝑚=0.8𝛿_th and in Figure 7 for filter size Δ=24Δ_𝑚=2.4𝛿_th. Note that the models predicting PE₂ outside a margin of ±15% have been discarded. In the case of Le = 0.34 (case A) the models FSDNEW, FSDF, MFSDF, and FSDC stay within the ±15% error limit for Δ=8Δ_𝑚 whereas only the models FSDF, MFSDF, FSDK and FSDNEW remain within the ±15% error limit for Δ=24Δ_𝑚. As Le increases to 0.6 (case B), the models FSDNEW, FSDCH, FSDF, MFSDF, FSDC, FSDA, and FSDK predict within the ±15% error margin and have been listed in terms of decreasing accuracy for Δ=8Δ_𝑚. For case B only the predictions of FSDNEW, FSDF, MFSDF and FSDK remain within the ±15% error margin for Δ=24Δ_𝑚. In the Le = 0.8 case (case C), the models FSDF, FSDNEW, MFSDF, FSDCH, FSDA, FSDC, FSDK and FSDW all provide predictions within ±15% for Δ=8Δ_𝑚, whereas the predictions of FSDNEW, FSDF, MFSDF, FSDCH, FSDK and FSDW remain within ±15% for Δ=24Δ_𝑚. For Le = 1.0 and 1.2 (cases D and E) the models FSDF, MFSDF, FSDNEW, FSDCH, FSDA, FSDC, FSDK and FSDW all predict within the ±15% error margin for Δ=8Δ_𝑚, while the models FSDF, MFSDF, FSDNEW, FSDK and FSDCH predict within ±15% for Δ=24Δ_𝑚. The model FSDW was found to predict within the ±15% error margin for Δ=24Δ_𝑚 in the Le = 1.0 flame but its prediction remains marginally beyond the ±15% error margin for Δ=24Δ_𝑚 for the Le = 1.2 flame considered here (The maximum magnitude of PE₂ for the FSDW model in the Le = 1.2 case is 15.2%, and the variation of PE₂ with 𝑐 in this case is qualitatively similar to the Le = 1.0 case considered here).

(a)

(b)

(c)

(d)

(e)

Figure 6 

Variation of percentage error (17) on 𝑐 across the flame brush for Δ=8Δ𝑚=0.8𝛿th according to FSDA, FSDC, FSDW, FSDCH, FSDK, FSDF, FSDNEW, and MFSDF predictions for (a–e) cases A–E.

(a)

(b)

(c)

(d)

(e)

Figure 7 

Variation of percentage error (17) on 𝑐 across the flame brush for Δ=24Δ𝑚=2.4𝛿th according to FSDA, FSDC, FSDW, FSDCH, FSDK, FSDF, FSDNEW, and MFSDF predictions for (a–e) cases A–E.

Comparing the performance of the models at Δ=8Δ_𝑚 and Δ=24Δ_𝑚, it can be seen that FSDA, FSDCH and FSDC predict Σ_gen satisfactorily at Δ=8Δ_𝑚 but the agreement with DNS deteriorates at Δ=24Δ_𝑚. By contrast, the FSDK prediction is closer to DNS data at Δ=24Δ_𝑚 than at Δ=8Δ_𝑚. The models FSDF, MFSDF, and FSDNEW fare well at both Δ=8Δ_𝑚 and Δ=24Δ_𝑚 for all the Lewis number values considered here. It is worth noting that the FSDNEW model was designed to predict the volume-averaged value of generalised FSD ⟨Σ_gen⟩, but judging from Figures 4–7, this model also performs satisfactorily with respect to predicting the correct variation of Σ_gen across the flame brush.

The prediction of the model FSDK improves with increasing filter width Δ, unlike the other models, which is consistent with observations made in the context of Figure 3. The prediction of the FSDW model remains skewed towards the product side of the flame brush due to the 𝑐̃ dependence of Ξ (i.e., Ξ=1+1.24𝑐̃𝑢^Δ/𝑆_𝐿Re_𝜂) proposed in [12]. The FSDW, FSDA, FSDC, and FSDCH models underestimate the destruction rate of flame surface area in the thin reaction zones regime due to the underestimation of FSD destruction arising due to the curvature stretch contribution −[𝐷_𝑐(∇⋅𝑁^→)²]_𝑠, which eventually leads to the overprediction of conditionally averaged value of Σ_gen.

4.4. Performance of Models for the Local Σ_gen Behaviour

The FSD predicted by the models should have the correct resolved strain rate and curvature dependence in the context of LES and thus the correlation coefficient between the FSD obtained from DNS and from the model prediction should remain as close to unity as possible. The variation of the correlation coefficients between the model prediction and generalised FSD Σ_gen obtained from DNS in the range of filtered reaction progress variable 0.1≤𝑐≤0.9 are shown in Figure 8 for different filter widths. The regions corresponding to 0.1<𝑐 and 𝑐>0.9 have been ignored since the correlation coefficients have little physical significance in these regions due to the small values of Σ_gen obtained from both DNS and model predictions. Figure 8 indicates that the correlation coefficients decrease with increasing Δ due to increased unresolved subgrid wrinkling, which makes the local variation of Σ_gen different from |∇𝑐|. The extent of the deviation of the correlation coefficients from unity increases with decreasing Le for a given value of Δ. Figure 8 indicates that the models FSDA, FSDC, FSDCH, FSDF, MFSDF, FSDK, FSDNEW, and FSDW have comparable correlation coefficients, which deviate considerably from unity for large values of Δ. This indicates that algebraic models may not be able to predict FSD such that its local strain rate and curvature dependencies can be appropriately captured, especially in the TRZ regime. Hence a transport equation for FSD might need to be solved to account for the local strain rate and curvature effects on Σ_gen [5, 6, 8, 10, 11].

(a)

(b)

(c)

(d)

(e)

Figure 8 

Correlation coefficients between the modelled and the actual values of Σgen in the 𝑐 range 0.1≤𝑐≤0.9 for filter widths Δ=4Δ𝑚=0.4𝛿th; Δ=8Δ𝑚=0.8𝛿th; Δ=12Δ𝑚=1.2𝛿th; Δ=16Δ𝑚=1.6𝛿th; Δ=20Δ𝑚=2.0𝛿th; Δ=24Δ𝑚=2.4𝛿th for (a–e) cases A–E.

5. Conclusions

The performance of several wrinkling factor based LES algebraic models for Σ_gen has been assessed for nonunity Lewis number flames in the TRZ regime based on a DNS database of freely propagating statistically planar turbulent premixed flames with Le ranging from 0.34 to 1.2. It has been found that the fractal dimension 𝐷 increases with decreasing Le, whereas Le does not have any significant influence on the value of the normalised inner cut-off scale 𝜂_𝑖/𝛿_𝑧. For all Lewis number cases the inner cut-off scale is found to be equal to the thermal flame thickness (i.e., 𝜂_𝑖≈𝛿_th). Based on the analysis of DNS data, a new parameterisation of 𝐷 is proposed, where the effects of Le are explicitly accounted for. This new parameterisation of 𝐷 has been used to propose a power-law based model for Σ_gen to account for nonunity Lewis number effects. The performance of this new model has been assessed with respect to Σ_gen obtained from DNS data alongside other existing models. The new model was found to be capable of predicting the behaviour of Σ_gen in the TRZ regime with greater or comparable accuracy in comparison to the existing models for all values of Le considered here. However, the present study has been carried out for moderate values of turbulent Reynolds number Re_𝑡 and the effects of detailed chemistry and transport are not accounted for. Thus, three-dimensional DNS with detailed chemistry will be necessary, together with experimental data, for a more comprehensive assessment of LES algebraic models for Σ_gen.

Appendix

A. Effects of Re_𝑡 on Fractal Dimension 𝐷

The effects of Re_𝑡 on 𝐷 have been analysed based on a simplified chemistry based DNS database [43, 44], in which the variation of Re_𝑡∼Da²Ka² is brought about by modifying Da and Ka independently of each other. The initial values of 𝑢′/𝑆_𝐿 and 𝑙/𝛿_th for all the flames in this DNS database are shown in Table 1(a) along with the values of heat release parameter 𝜏=(𝑇_ad−𝑇₀)/𝑇₀, Damköhler number Da=𝑙𝑆_𝐿/𝑢′𝛿_th, Karlovitz number Ka=(𝑢′/𝑆_𝐿)^3/2(𝑙𝑆_𝐿/𝛼_𝑇0)^−1/2, and turbulent Reynolds number Re_𝑡=𝜌₀𝑢′𝑙/𝜇₀.

The variations of log(⟨Σ_gen⟩/⟨|∇𝑐|⟩) with log(Δ/𝛿_𝑧) for cases A1–E1 are shown in Figure 9, which demonstrate that 𝐷 is greater for flames with higher Re_𝑡, and that 𝐷 attains an asymptotic value of 7/3 for unity Lewis number flames with high values of Re_𝑡 (e.g., cases D1 and E1). The prediction of ⟨Σ_gen⟩/⟨|∇𝑐|⟩=(Δ/𝜂_𝑖)^𝐷−2 with 𝜂_𝑖 obtained from DNS and 𝐷 obtained from (11) is also shown in Figure 9, which indicates that (11) satisfactorily captures the slope of the best-fit straight line.

(a)

(b)

(c)

(d)

(e)

Figure 9 

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

Acknowledgment

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