Research Article  Open Access
Mohit Katragadda, Nilanjan Chakraborty, "A Priori Direct Numerical Simulation Modelling of the Curvature Term of the Flame Surface Density Transport Equation for Nonunity Lewis Number Flames in the Context of Large Eddy Simulations", International Journal of Chemical Engineering, vol. 2012, Article ID 103727, 14 pages, 2012. https://doi.org/10.1155/2012/103727
A Priori Direct Numerical Simulation Modelling of the Curvature Term of the Flame Surface Density Transport Equation for Nonunity Lewis Number Flames in the Context of Large Eddy Simulations
Abstract
A Direct Numerical Simulation (DNS) database of freely propagating statistically planar turbulent premixed flames with Lewis numbers ranging from 0.34 to 1.2 has been used to analyse the statistical behaviours of the curvature term of the generalised Flame surface Density (FSD) transport equation, in the context of the Large Eddy Simulation (LES). Lewis number is shown to have significant influences on the statistical behaviours of the resolved and subgrid parts of the FSD curvature term. It has been found that the existing models for the subgrid curvature term do not capture the qualitative behaviour of this term extracted from the DNS database for flames with . The existing models of only predict negative values, whereas the subgrid curvature term is shown to assume positive values within the flame brush for the and 0.6 flames. Here the subgrid curvature terms arising from combined reaction and normal diffusion and tangential diffusion components of displacement speed are individually modelled, and the new model of the subgrid curvature term has been found to capture extracted from DNS data satisfactorily for all the different Lewis number flames considered here for a wide range of filter widths.
1. Introduction
Flame Surface Density (FSD) based reaction rate closure is well established in the context of Reynolds Averaged NavierStokes (RANS) simulations of turbulent premixed flames [1, 2]. The increased affordability of high performance computing has made Large Eddy Simulation (LES) an alternative simulation tool, where the largescale physical processes are resolved, but modelling is still required for the subgrid quantities. The FSDbased reaction rate closure has recently been successfully extended for the purpose of LES [3–14]. In LES simulation of premixed combustion, a Favrefiltered reaction progress variable transport equation is solved alongside other filtered conservation equations. The reaction progress variable is defined as , where is the mass fraction of a suitable reactant and the subscripts 0 and ∞ denote the values in the fully unburned and burned gases, respectively. The generalised FSD is defined as [3–14], where the overbar indicates the LES filtering process. The Favrefiltered reaction progress variable transport equation takes the following form: where indicates the Favre filtered value of a general variable Q, is the velocity component in the jth direction, ρ is the density, D is the molecular diffusivity, and is the filtered reaction rate. The first two terms on right hand side of (1) denote the filtered molecular diffusion and reaction rates, respectively, and their combined contribution can be modelled using in the following manner: where indicates the surfaceweighted filtered value of a general quantity Q and is the displacement speed, which denotes the speed at which a given c isosurface moves normal to itself with respect to an initially coincident material surface. The generalised FSD is an unclosed quantity and is closed either by using an algebraic expression or by solving a modelled transport equation alongside other conservation equations. The algebraic closure is valid when the generation rate of flame surface area remains in equilibrium with its destruction rate, but this assumption is rendered invalid under unsteady conditions (e.g., combustion instabilities). Under unsteady conditions, it is often advantageous to solve a modelled transport equation of. The exact transport equation for the generalised FSD is given as [1, 4–7, 9, 10, 12]: where is the ith component of flame normal vector. The terms on the left hand side of (3) denote transient and mean advection effects, respectively. The first three terms on the right hand side of (3) denote the effects of subgrid convection, flame surface area generation due to fluiddynamic straining, and flame normal propagation, respectively. The last term of (3) describes the production/destruction of due to flame curvature and thus referred to as the FSD curvature term [4–7, 9, 10, 12]. It has been found in several previous studies [5–7, 9, 14] that the FSD curvature term remains a leading order contributor to the FSD transport for both unity and nonunity Lewis number turbulent premixed combustion. As the curvature term remains a leading order contributor to the FSD transport, the modelling of is crucial for the transport equationbased FSD closure. The statistical behaviour of is significantly affected by curvature dependence on [9, 10, 12]. Earlier a priori Direct Numerical Simulation (DNS) analyses [9, 10, 12] showed that existing models for the subgrid curvature term do not adequately capture the qualitative behaviour of this term obtained from DNS data. Moreover, the model parameters for the existing subgrid curvature term models are found to be strong functions of the LES filter width Δ [9, 10, 12].
To date, most existing FSDbased models have been proposed for unity Lewis number flames where the differential diffusion of heat and mass has been ignored. The Lewis number is defined as the ratio of thermal diffusivity to mass diffusivity (i.e., ). The effects of on the statistical behaviour of the FSD curvature term are yet to be analysed in detail, and this paper aims to bridge this gap in the existing literature. It is worth noting that, in a premixed flame, different species have different values of Lewis number. Thus, specifying a global Lewis number characterising the whole combustion process is not straightforward. The Lewis number of the deficient reactant is often considered to be the characteristic of the combustion process in question [15, 16]. Moreover, several previous studies [16–29] analysed the effects of differential diffusion of heat and mass by modifying the characteristic Lewis number in isolation, and the same procedure has been adopted here. In the present study, a simplified chemistrybased DNS database of statistically planar turbulent premixed flames with global Lewis numbers ranging from 0.34 to 1.2 has been considered to analyse the statistical behaviour of the FSD curvature term in the context of LES. In this context, the main objectives of this study are as follows: (1)to analyse the statistical behaviours of the subgrid FSD curvature term in the context of LES, for flames with different values of Lewis number; (2)to propose models for different components of the subgrid FSD curvature terms and assess their performances in comparison to the corresponding quantities extracted from DNS data.
The rest of the paper will be organised as follows. The necessary mathematical background will be provided in the next section. This will be followed by a brief description of the numerical implementation related to the DNS database. Following this, results will be presented and subsequently discussed. The main findings will be summarised, and conclusions will be drawn in the final section.
2. Mathematical Background
The curvature term of the FSD is often decomposed in the following manner [4–7, 9, 10, 12]: where and are the resolved and subgrid components of the FSD curvature term, respectively. The resolved curvature term can be expressed in three different manners [5, 9, 10, 12]: where is the ith component of the resolved flame normal vector. It was demonstrated by Chakraborty and Cant [10, 12] that provides the best option for the resolved curvature term , as it gives rise to the smallest magnitude of among all the possibilities shown in –. Equation was found to perform the best among the three possibilities shown in , , for this database. This is advantageous from the perspective of efficient modelling of the FSD curvature term as most of the modelling uncertainty is associated with . Moreover, has also been used for the modelling of in previous LES simulations [5–7, 13]. For the present analysis , (i.e., ) will be considered for the resolved curvature term .
It is often useful to decompose the flame displacement speed in the following manner for the purpose of modelling the FSD curvature term [9–12, 30, 31]: where and are the reaction and normal diffusion components of displacement speed and is the tangential diffusion component of displacement speed. The following expression for can be obtained using – and (i.e., ): where Equation indicates that curvature () dependences of and significantly influence the statistical behaviour of . Equation suggests that is expected to assume negative values throughout the flame brush.
Hawkes and Cant [6, 7] modified a version of the Coherent Flamelet Model (CFM) by Candel et al. [2] for the purpose of LES as: where is a resolution parameter which vanishes when the flow is fully resolved and is a model parameter. Hawkes [5] discussed a possibility of modifying a RANS model proposed by Cant et al. [1] for the purpose of LES as: where , , is the subgrid turbulent velocity fluctuation, is the subgrid kinetic energy, and is a model parameter. Another model of was proposed by Charlette et al. [4]: where is a model parameter. The models given by – (henceforth will be referred to as CSGCFM, CSGCPB, and CSGCHAR, resp.) ensure that vanishes when the flow is fully resolved (i.e., and ). A priori DNS assessment of the CSGCFM, CSGCPB, and CSGCHAR models and the modelling of and will be addressed in Section 4 of this paper.
3. Numerical Implementation
In principle combustion, DNS should account for both three dimensionality of turbulence and detailed chemical mechanism. However, until recently, most combustion DNS studies were carried out either in two dimensions with detailed chemistry or in three dimensions with simplified chemistry due to the limitation of computer storage capacity. Although it is now possible to carry out threedimensional DNS with detailed chemistry, they remain extremely expensive (e.g., millions of CPU hours and thousands of processors [32]) and the cost of an extensive parametric analysis based on threedimensional detailed chemistrybased DNS often becomes prohibitive. As the present analysis concentrates on an extensive parametric variation in terms of Lewis number, the chemical mechanism is simplified here by an Arrheniustype irreversible singlestep chemical reaction (i.e., Reactants→Products) following several previous studies [1–12, 14]. It has been found that the strain rate and curvature dependences of and obtained from threedimensional simplified chemistry DNS [25–27, 33, 34] are found to be qualitatively similar to the corresponding behaviours obtained from detailed chemistrybased DNS simulations [16, 30, 31, 35]. As the statistical behaviours of the FSD curvature term are strongly dependent on the curvature dependences of and , the results for this analysis are expected to be valid even for detailed chemistry based simulations at least in a qualitative sense without much loss of generality. Several studies [3–7, 9–12] have concentrated on a priori DNS modelling of FSD based on simplified chemistry in the past and the same approach has been adopted here.
A compressible threedimensional DNS code SENGA [36] was used for the simulations where the conservation equations of mass, momentum, energy, and species are solved in nondimensional form. A cubic domain of each side equal to is considered for the present DNS database where is the thermal flame thickness, which is defined as , and the subscript refers to quantities in an unstrained planar laminar flame with , , and being the adiabatic flame, unburned gas, and instantaneous gas temperatures, respectively. The computational domain was discretised using a Cartesian grid of with equal grid spacing in each direction. The grid spacing is determined based on the flame resolution, and about 10 grid points are kept within the thermal flame thickness for all the cases considered here. This grid spacing corresponds to , where is the Kolmogorov length scale. The boundaries in the mean flame propagation were taken to be partially nonreflecting and were implemented using the NavierStokes Characteristic Boundary Conditions (NSCBC) technique [37]. The boundary conditions in the transverse direction were taken to be periodic. The spatial derivatives for the internal grid points were evaluated using a tenthorder central differencing scheme, and the order of differentiation gradually decreases to a onesided 2nd order scheme at the partially nonreflecting boundaries. The time advancement was carried out using an explicit low storage thirdorder RungeKutta scheme [38].
For the current DNS database, the turbulent velocity field was initialised using a pseudospectral method [39] following the BatchelorTownsend turbulent kinetic energy spectrum [40]. The flame is initialised using a steady planar unstrained laminar flame solution. The initial values of normalised root mean square (rms) turbulent velocity fluctuation , integral length scale to thermal flame thickness ratio , heat release parameter , Damköhler number , and Karlovitz number are listed in Table 1. According to Peters [41], all the cases considered here can be taken to represent the thin reaction zone regime combustion, as Ka remains greater that unity. Standard values are considered for Prandtl number and the Zel’dovich number (i.e., and , where is the activation temperature).

Under decaying turbulence, DNS simulations should be carried out for a simulation time [42], where is the initial eddy turn over time and is the chemical time scale. For this database, the statistics were extracted after about three eddy turn over times (i.e., ), which corresponded to one chemical time scale (i.e., ). This simulation time remains small but comparable to several studies [3, 24, 28, 43–47] which contributed significantly to the FSDbased modelling in the past. The statistics presented in this paper did not change significantly since halfway through the simulation (i.e., ). The value of in the unburned gas ahead of the flame had decayed by 50% of its initial value when the statistics were extracted. By this time, the normalised integral value had increased to around 1.7 times of its initial value. The values of and at the time statistics were extracted are also representative of the thin reaction zones regime combustion [41]. This DNS database was used extensively earlier for the purpose of RANS modelling [27, 28, 48, 49], and the interested readers are referred to these papers for further details.
The DNS data was explicitly LES filtered using a Gaussian filter kernel in physical space for the purpose of a priori analysis. The filtered quantity , is given by where is the Gaussian filter kernel, which is defined in the following manner:
The filtered quantities of interest were extracted for filter widths ranging from to in steps of . These filter sizes are comparable to the range of used in several previous studies [3, 4, 9–12, 14] for a priori DNS analysis and span a useful range of length scales (i.e., from , where the flame is partially resolved, up to , where the flame becomes fully unresolved and is comparable to the integral length scale).
4. Results and Discussion
The instantaneous isosurfaces of ranging from 0.01 to 0.99 at are shown in Figure 1, which indicates that the flame wrinkling increases with decreasing Lewis number and this tendency is particularly prevalent for the flames due to thermodiffusive instabilities [17–29]. The unburned reactants diffuse into the reaction zone at a faster rate than the rate at which heat diffuses out in the flames. This gives rise to simultaneous presence of high temperature and reactant concentration in the reaction zone for the flames, which in turn leads to greater burning rate and flame surface area generation in comparison to the unity Lewis number flame. By contrast, heat diffuses faster than the diffusion rate of reactants into the reaction zone in the case of , which reduces the burning rate and the rate of flame area generation in comparison to the unity Lewis number flame. The increase in burning rate and flame area generation with decreasing Lewis number can be substantiated by the values of normalised turbulent flame speed and normalised flame surface area which are presented in Table 2. The values of have been evaluated by volume integrating the reaction rate using the expression , where is the projected area of the flame in the direction of mean flame propagation, while the values of have been evaluated by volume integrating (i.e., ) under both turbulent and laminar conditions. Table 2 shows that both and increase with decreasing Lewis number, and this effect is particularly prevalent in the flames with Le < 1 due to the presence of thermodiffusive instabilities [17–29]. The increase in flame wrinkling with decreasing Lewis number is also visually evident from the isosurfaces presented in Figure 1.

(a)
(b)
(c)
(d)
(e)
The variations of , , and conditionally averaged in bins of isosurfaces for cases (a)–(e) are shown in Figure 2 for filter widths and , where Δ_{m} is the DNS grid size. It is evident from Figure 2 that Le significantly affects the statistical behaviours of the curvature terms. The filter widths , and span a useful range of length scales (i.e., from , where the flame is partially resolved, up to where the flame becomes fully unresolved and is comparable to the integral length scale). In the Le ≪ 1 flames (e.g., cases (a) and (b)), the FSD curvature term behaves as a source term for the major part of the flame brush before assuming negative values towards the burned gas side for . For , the FSD curvature term acts as a source (sink) term towards the unburned (burned) gas side of the flame brush in the flames. In the case of flames (i.e., cases (c)–(e)) the curvature term behaves as a sink type term throughout the flame brush for all filter widths. It can be seen from Figure 2 that acts as a source (sink) term for cases (a)(b) ((c)–(e)). The magnitude of () decreases (increases) with increasing in all cases, and for large filter widths is principally made up of . The LES filtering is a convolution process, and the weighted averaging involved in the filtering process leads to a decrease in the magnitude of with increasing filter width . The flow becomes increasingly unresolved with increasing filter width , and this is reflected in the rise in magnitude with increasing filter width .
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
The resolved curvature term can be seen to capture the behaviour of the curvature term , well at small filter widths (i.e., ) for flames with (i.e., cases (c)–(e)). However, the magnitude of decreases with increasing and it does not capture the behaviour of the FSD curvature term for the Le ≪ 1.0 flames (i.e., cases (a) and (b)). The subgrid curvature term, C_{sg} follows the qualitative behaviour of the FSD curvature term for all filter widths. The subgrid curvature term almost entirely makes up the FSD curvature term for , and this is especially true for the Le ≪ 1.0 cases (i.e., cases (a) and (b)). It can further be observed from Figure 2 that assumes positive values towards the unburned gas side of the flame brush in the Le ≪ 1 flames (e.g., cases (a) and (b)), whereas the existing models for allow for only negative values (see –). This suggests that new models for are warranted to account for the influences of nonunity Lewis number (i.e., ).
In order to be able to model the subgrid curvature term , the decomposition prescribed in  has been used here. The variations of , , and conditionally averaged in bins of isosurfaces for cases (a)–(e) are shown in Figure 3 for filter widths and . It is evident from Figure 3 that remains negative throughout the flame brush for all cases and follows the qualitative behaviour of . A comparison between and reveals that remains the major contributor to for all the flames at all values of , which is consistent with the expected behaviour in the thin reaction zones regime [41]. The contribution of remains significant for the cases (i.e., cases (a), (b) and (c)), but its contribution remains weak in comparison to the magnitude of in the and 1.2 flames (i.e., cases (d) and (e)). Figure 3 demonstrates that remains close to the magnitude of for all for the flame (i.e., case (d)), indicating that does not play a major role in capturing the behaviour of . However, there is a significant difference in magnitudes of and for small values of (i.e., ) in the nonunity Lewis number flames (i.e., cases (a)–(c) and (e)), which indicates that plays a key role for small values of filter width in these flames. For large values of filter width (i.e., ) remains the major contributor to for all cases considered here, indicating that plays progressively less important role for increasing values of .
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
Figure 3 shows that there is a significant difference between and for all cases for small values of , and the difference between these quantities decreases with increasing . As most of the contribution of remains unresolved for large values of , the subgrid curvature term remains the major contributor to , indicating that plays progressively less important role for increasing values of where the flame is fully unresolved. However, the contribution of remains significant for small values of , where the flame is partially resolved. Figure 3 further shows that the order of magnitudes of both and remains comparable and thus accurate modelling of and is necessary for precise predictions of .
As the range of values obtained on a flame surface increases with increasing flame wrinkling, the magnitude of increases with decreasing , which in turn leads to increasing magnitude of and (see Figure 3). The positive contribution of overcomes the negative contribution of towards the unburned gas side of the flame brush for the and 0.6 flames (i.e., cases (a) and (b)) and yields a net positive contribution of towards the reactant side of the flame brush (see Figure 2).
The statistical behaviours of and depend on the nature of the curvature dependences of and , and the variation of across the flame brush. The correlation coefficients for and for five different isosurfaces across the flame brush for all the cases are shown in Figures 4(a) and 4(b), respectively. For all cases, remains negatively correlated with with a correlation coefficient close to (−1.0). However, Figures 4(a) and 4(b) demonstrate that significantly affects the curvature dependences of and . It can be seen from Figures 4(a) and 4(b) that and remain positively (negatively) correlated with for the () flames, whereas both and show weak curvature dependences in the unity Lewis number flame. The positive correlations between and , and between and strengthen with decreasing for the flames. The physical explanations for the observed influences of Lewis number on the curvature dependence of and have been discussed elsewhere [25–27] and thus will not be discussed in this paper.
(a)
(b)
(c)
(d)
The variations of conditionally averaged in bins of isosurfaces for cases A–E are shown in Figures 4(c) and 4(d) for filter widths and , respectively. It is evident from Figures 4(c) and 4(d) that assumes positive (negative) values towards the unburned (burned) gas side of the flame brush. For small values of , the surfaceweighted filtered value of curvature approaches to (i.e., ) and thus the ensemble averaged value of remains small for small values of filter width as the ensemble averaged value of remains negligible for statistically planar flames. The difference between the ensemble averaged values of and increases with increasing filter width , as flame wrinkling increasingly takes place at the subgrid level. For the flame (i.e., case D), the combination of positive (negative) value of and weak and correlations gives rise to positive (negative) values of the ensemble averaged values of and towards the unburned (burned) gas side of the flame brush for all values of . The predominant positive and correlations give rise to positive values of the ensemble averaged values of and throughout the flame brush for small values of in the , 0.6, and 0.8 flames. By contrast, negative and correlations (see Figures 4(a) and 4(b)) give rise to negative values of the ensemble averaged values of and throughout the flame brush for small values of in the flame. These local dependences are progressively smeared with increasing because of the convolution operation associated with LES filtering process, and this leads to positive (negative) values of and towards the unburned (burned) gas side of the flame brush for all cases considered here, including the nonunity Lewis number flames where the curvature dependences of and are particularly strong.
The dependences of and on are likely to capture some of dependences of and at small values of filter widths (i.e., ), where the flame is partially resolved. This effect is particularly prevalent in the nonunity Lewis number flames where both and are strongly correlated with curvature even though the flames are statistically planar in nature. As a result of this, the contribution of remains close to that of for small filter widths (i.e., ) for the nonunity Lewis number flames, which is reflected in the small contribution of (see variations in Figures 3(a)–3(c) and Figure 3(e)). The correlation between the resolved quantities (e.g., dependences of and on ) weakens with increasing filter width due to smearing of local information. Moreover, physical processes take place increasingly at the subgrid level for , and thus does not capture the behaviour of for large filter widths in all cases considered here, including the nonunity Lewis number flames where the curvature dependences of and are particularly strong. This leads to for in all cases considered here (see variations in Figures 3(f)–3(j)). It can be seen from Figure 3 that the positive contribution of overcomes the negative contribution of towards the unburned gas side of the flame brush in the and 0.6 flames, which lead to positive value of towards the unburned gas side for all values of in these cases (see Figure 2). By contrast, negative values of overcome the positive contributions of towards the unburned gas side of the flame brush in the , 1.0, and 1.2 flames, which lead to negative values of throughout the flame brush in these cases (see Figure 2).
The subgrid fluctuations of the surfaceweighted contributions of and are scaled here using and , respectively, to propose the following model for : where , and are the model parameters. The function in is used to capture the correct qualitative behaviour of throughout the flame brush. In a compressible, LES simulation is readily available and needs to be extracted from .The methodology of extracting from in the context of LES was discussed elsewhere [9, 10, 12] and will not be discussed in detail in this paper. The model parameter ensures that the transition from positive to negative value of takes place at the correct location within the flame brush. The quantity vanishes when the flow is fully resolved (i.e., ), and thus becomes exactly equal to zero when the flow is fully resolved (i.e., ) according to . It has been found that enables to capture the qualitative behaviour of when the optimum values of and are chosen. The optimum values of () tend to increase with decreasing (increasing) . The curvature dependences of and are influenced by (see Figures 4(a) and 4(b)), and these local dependences also appear in the resolved scale but their strength diminishes with increasing due to filtering operation. As the resolved and subgrid curvature terms are closely related [9, 10, 12], the qualitative behaviour of is also affected by the curvature dependences of displacement speed components and scalar gradient at the resolved scale, which leads to the variation of the optimum values of , , and with and . The model parameter needs to be deceased for decreasing values of for satisfactory prediction of . The prediction of ensemble averaged on isosurfaces is compared with the ensemble averaged values of in Figure 5 for all cases considered here for the optimum values of , and for and when is taken to be . The optimum values of , and are estimated by calibrating the prediction of with respect to the ensemble averaged values of obtained from DNS data and the variation of the optimum values of , , and with for all cases are shown in Figure 6. The optimum values of , , and are parameterised here as where where
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(a)
(b)
(c)
(d)
(e)
Figure 5 shows that satisfactorily predicts when is taken to be , and the optimum values of , and are used. According to the parameterisation given by –, increases with decreasing Le, as the effects of chemical reaction strengthen with decreasing values of Lewis number (see Table 2). Moreover, , , and approach to asymptotic values for large values of and turbulent Reynolds number based on LES filter width , where and are the unburned gas density and viscosity, respectively.
Here, the contribution of is scaled with (i.e., ) where the subgrid fluctuations of are taken to scale with . The above relations are utilised here to propose a model for in the following manner: where is the wrinkling factor [8, 11, 43, 50, 51], and are the model parameters, and is used to capture the correct qualitative behaviour of . The subgrid curvature term vanishes when the flow is fully resolved according to , (i.e., ). It has been found that satisfactorily captures the behaviour of throughout the flame brush for in all cases when a suitable value of is used. The variation of the global mean optimum values of with is shown in Figure 6 for all cases considered here. The optimum values of have been parameterised here in the following manner: where
The predictions of ensemble averaged on isosurfaces are compared with ensemble averaged values of in Figure 5 for all cases at and , which show that satisfactorily predicts the statistical behaviour of when n is taken to be and the optimum value of is used. According to –, approaches to asymptotic values for large values of and .
Equations and – can be combined to propose a model for in the following manner:
The above model will henceforth be referred to CSGNEW model in this paper. Equation allows for a positive contribution of through the contribution of , which is absent in the CSGCAND, CSGCANT, and CSGCHAR models. The predictions of the CSGCAND, CSGCANT, CSGCHAR, and CSGNEW models for and are compared with