Abstract
A simplified chemistry based three-dimensional Direct Numerical Simulation (DNS) database of freely propagating statistically planar turbulent premixed flames with a range of different values of turbulent Reynolds number has been used for the a priori modelling of the curvature term of the generalised Flame Surface Density (FSD) transport equation in the context of Large Eddy Simulation (LES). The curvature term has been split into the contributions arising due to the reaction and normal diffusion components of displacement speed and the term originating from the tangential diffusion component of displacement speed. Subsequently, these contributions of the curvature term have been split into the resolved and subgrid contributions. New models have been proposed for the subgrid curvature terms arising from the combined reaction and normal diffusion components and the tangential diffusion component of displacement speed. The performances of the new model and the existing models for the subgrid curvature term have been compared with the corresponding quantity extracted from the explicitly filtered DNS data. The new model for the subgrid curvature term is shown to perform satisfactorily in all cases considered in the current study, accounting for wide variations in LES filter size.
1. Introduction
Flame Surface Density (FSD) based reaction rate closure is one of the popular methods of turbulent premixed combustion modelling in the context of Reynolds Averaged Navier Stokes (RANSs) simulations [1, 2]. The FSD based modelling has recently been extended to Large Eddy Simulations (LESs) [3–12]. The generalised FSD is defined as [3–10] where is the reaction progress variable and the overbar indicates a LES filtering operation. The transport equation of is given by [1, 4–7, 9, 11]: where is the th component of flame normal vector and is the displacement speed, and are the surface-weighted and Favre filtered values of a general quantity . The final term on the right hand side of (1) originates due to flame curvature and thus this term (i.e., ) is referred to as the curvature term [4–7, 9, 11]. It is evident from (1) that the curvature dependence of plays a key role in the statistical behaviours of and this was confirmed in previous a priori Direct Numerical Simulation (DNS) based analyses [9, 11]. It was previously demonstrated [9, 11] that the existing models for the subgrid curvature term often do not capture its correct qualitative and quantitative behaviours, particularly in the Thin Reaction Zones (TRZ) regime flames. Moreover, the model parameters for the existing models were found to be strong functions of LES filter width [9, 11]. The modelling of therefore remains one of the weakest points in the LES modelling of the transport equation. This gap in the existing literature is addressed in this paper by explicitly LES filtering a DNS database of freely propagating statistically planar turbulent premixed flames with different values of turbulent Reynolds number . In this regard, the main objectives of the present study are as follows:(1)to analyse the statistical behaviours of the subgrid FSD curvature term in the context of LES for the flames with different values of ;(2)to propose models for the subgrid FSD curvature term and assess their performances in comparison to the corresponding quantities extracted from DNS data.
The necessary mathematical background and numerical details will be provided in the next section. Following this, the results will be presented and subsequently discussed. Finally the main findings will be summarised and conclusions will be drawn.
2. Mathematical Background and Numerical Implementation
Although three-dimensional DNS with detailed chemical mechanism is currently possible, it remains extremely computationally intensive [13] and is often not suitable for a detailed parametric analysis. Thus the chemical mechanism is simplified here using a single step Arrhenius type chemical reaction in order to carry out a parametric variation in terms of . For the convenience of modelling, is often split as [4–7, 9, 11] where and are the resolved and subgrid curvature terms, respectively. Chakraborty and Cant [9, 11] analysed the possibility of using three different expressions of : where and are the th component of surface-weighted and resolved flame normal vector, respectively. Previous a priori DNS analyses suggested that is the most preferred expression for the resolved curvature term out of the three options presented in (3), as it allows for the smallest magnitude of , while satisfactorily capturing the qualitative behaviour of [9, 11]. Moreover, was used in previous LES simulations [5–7, 12]. It is useful to split in the following manner [9–11, 14, 15] for obtaining further insight into : where is the reaction rate and is the progress variable diffusivity. The following expression for can be obtained using (4) and : where Equation indicates that the curvature dependences of and are likely to influence the statistical behavior of . According to , remains deterministically negative throughout the flame brush. Hawkes and Cant [6, 7] modified a version of the coherent flamelet model by Candel et al. [2] 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 the RANS model proposed by Cant et al. [1] for the purpose of LES as: where and is the subgrid kinetic energy and is a model parameter. Charlette et al. [4] modelled as where is a model parameter. The models given by (7)–(9) (henceforth will be referred to as CSGCFM, CSGCPB, and CSGCHAR, respectively) ensure that vanishes when the flow is fully resolved (i.e., and ). Modelling of and using a priori analysis of DNS data, and the assessment of the models given by (7)–(9), will be addressed in Section 3 of this paper.
In the present study a compressible DNS database of freely propagating statistically planar turbulent premixed flames under decaying turbulence has been considered. The simulation domain of size , was discretised using a Cartesian mesh of size with uniform mesh spacing in each direction where is the thermal flame thickness with , , and being the adiabatic flame, unburned gas, and instantaneous gas temperatures, respectively, and the subscript refers to the unstrained planar laminar flame quantities. The domain boundaries in the direction of mean flame propagation (i.e., -direction) are taken to be partially nonreflecting, whereas the transverse boundaries are taken to be periodic. The partially nonreflecting boundaries are specified using the well-known Navier Stokes Characteristic Boundary Conditions (NSCBC) technique [16]. The simulations have been carried out using a three-dimensional compressible DNS code called SENGA [17] A 10th order central difference scheme is used for spatial differentiation for internal grid points, and the order of differentiation decreases gradually to a one-sided 2nd order scheme towards nonperiodic boundaries [17]. A third order Runge-Kutta scheme was used for the purpose of time advancement [17]. For all cases, the reacting flow field is initialised by a steady unstrained planar laminar flame solution, and the initial turbulent velocity fluctuations are specified using an initially homogeneous isotropic velocity field. About 10 grid points are kept within the thermal flame thickness for all cases considered here. The initial values for the root-mean-square turbulent velocity fluctuation normalised by unstrained planar laminar burning velocity and the integral length scale to flame thickness ratio are presented in Table 1 along with the values of Damköhler number , Karlovitz number , and turbulent Reynolds number , where and are the unburned gas density and viscosity, respectively. The heat release parameter and Lewis number are taken to be 4.5 and 1.0 for all cases considered here. Standard values are taken for Prandtl number , ratio of specific heats , and the Zel’dovich number (i.e., , , and ), where is the activation temperature. The turbulent Reynolds number scales as , and thus the variation of in cases A–E is brought about by modifying and independently from each other. In cases A, C and E, is held constant, while is held constant in cases B, C, and D. For all cases the Karlovitz number remains greater than unity indicating the TRZ regime combustion according to the regime diagram by Peters [18]. The range of values considered in this study remains modest, although several previous studies [3, 9, 11, 15, 19–23] with comparable values of have made valuable contributions to the fundamental understanding and the modelling of turbulent premixed combustion. Moreover, the range of considered here is comparable to that of previous laboratory-scale experiments [24].
In all cases flame-turbulence interaction takes place under decaying turbulence. The simulations were run for a time equal to one chemical time scale (i.e., ), which is equivalent to in case D; in cases A, C, and E; for case B. The aforementioned simulation times remain comparable to several studies [3, 9, 11, 15, 19–23], which contributed to the FSD based modelling in the past. The global turbulent kinetic energy and burning rate were not varying significantly with time when statistics were extracted (see Figure 1 of [23]) and the qualitative nature of the statistics was found to have remained unchanged since for all cases [23]. At time , the global level of had decayed from the initial values by about 45%, 55%, 40%, 25%, and 32% in cases A–E, respectively. The values of had increased from their initial values by a factor of about 1.5–2.25 at , but there were still enough turbulent eddies on each side of the computational domain [23]. Values for and at the time when statistics were extracted were presented in Table 2 of [23] and are not repeated here. The flame thickness remained greater than the Kolmogorov length scale for all cases when the statistics were extracted (see Table 2 of [23]), confirming the TRZ regime combustion.
For the purpose of a priori DNS analysis, the relevant quantities are explicitly filtered using a Gaussian kernel . The filtered value of a general quantity is evaluated using the following convolution operation: . The statistical behaviours of the FSD curvature term have been analysed here for ranging from to , where is the DNS mesh size . These filter sizes are comparable to the range of used in a priori DNS analysis in several previous studies [3, 4, 9–11] and span a useful range of length scales (i.e., from comparable to , where the flame is partially resolved, up to , where the flame becomes fully unresolved and is comparable to the integral length scale).
3. Results and Discussion
The isosurfaces of ranging from 0.01 to 0.99 at time for all cases are shown in Figure 1. A comparison between Figures 1(a)–1(e) reveals that the wrinkling of the flame surface with increasing . Figures 1(a)–1(e) further demonstrate that the flame surfaces in all cases show a range of different curvatures and this range increases with increasing . This indicates that the interrelations between and and between and may lead to nonnegligible value of even for statistically planar flames.
The variations of the ensemble averaged values of conditional on isosurfaces for all cases are shown in Figure 2 for and . The filter widths and correspond to two representative situations, where the flame is partially resolved and where the flame is fully unresolved, respectively. Figure 2 shows that assumes predominantly negative values throughout the flame brush. Although the resolved curvature term captures the qualitative behaviour of throughout the flame brush, the magnitude of remains smaller than the magnitude of for the major portion of the flame brush in all cases for all values of (see Figure 2). This leads to predominantly negative values of , although the ensemble averaged values of on isosurfaces exhibits positive values towards the unburned gas side of the flame brush for the flames with low and moderate values of turbulent Reynolds number (e.g., cases A–C). By contrast, the variation of ensemble averaged values of on isosurfaces exhibits only negative values throughout the flame brush for the flames with high values of turbulent Reynolds number (e.g., cases D and E). The models of given by (7)–(9) only predicts negative values of and thus will not be capable of predicting the positive values of in cases A–C.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
Comparing , , and magnitudes for and reveals that the magnitudes of and decrease with increasing , whereas the relative magnitude of in comparison to increases with increasing . This observation is consistent with previous findings [9, 11]. The smearing of local information as a result of the weighted-averaging process involved in LES filtering leads to the decrease in the magnitudes of and for increasing values of . The flow becomes increasingly unresolved with increasing and thus the flame curvature and its influence on the FSD evolution are increasingly felt at the subgrid scale, which is reflected in the relatively high magnitudes of in comparison to for large values of .
It is useful to examine the statistical behaviours of and in order to explain the differences in the behaviours of for flames with different . The variations of the ensemble averaged values of , , and conditional on isosurfaces are shown in Figure 3 for cases A–E for and , respectively. Figure 3 demonstrates that both and remain predominantly positive (negative) towards the unburned (burned) gas side of the flame brush for all values of considered here. The contribution of and remains deterministically negative throughout the flame brush (see Figure 3). It is evident from Figure 3 that remains a leading order contributor to for all the flames at all values of (see Figure 3), which is consistent with the expected behaviour in the TRZ regime, where is expected to play an important role [18]. Figure 3 further shows that remains close to the magnitude of for all for all cases considered here, indicating that does not play a major role in capturing the behaviour of . By contrast, 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 , remains the leading order contributor to , indicating that plays a progressively less important role for increasing values of , where the flame is fully unresolved. However, the contribution of remains significant for small values of when the flame is partially resolved. Figure 3 further shows that the order of magnitudes of both and remains comparable for large values of (i.e., ) and thus accurate modelling of and is necessary for accurate modelling of . As the range of values obtained on a flame surface increases with increasing flame wrinkling at higher values of , the magnitude of increases with increasing , which in turn leads to increasing magnitude of and with increasing for a given value of or (see Figure 3). The positive contribution of overcomes the negative contribution of towards the unburned gas side of the flame brush for the flames with small and moderate values of turbulent Reynolds number (i.e., cases A–C) and yields a net positive contribution of towards the reactant side of the flame brush (see Figure 2).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
The statistical behaviours of and depend on the nature of the correlations between and and between , , and the variation of across the flame brush. The correlation coefficients for the and dependences for five different isosurfaces across the flame brush for all cases are shown in Figures 4(a) and 4(b), respectively. For unity Lewis number flames is deterministically negatively correlated with with a correlation coefficient equal to −1.0. Figure 4(a) suggests that correlation is much weaker than the correlation in all cases. Moreover, Figure 4(b) demonstrates that and remain weakly correlated throughout the flame brush for all cases. Figures 4(a) and 4(b) demonstrate that the curvature dependences of and remain qualitatively similar for all the flames [25]. The physical explanations of the observed curvature dependences of and have been discussed elsewhere [25] and will not be repeated here.
(a)
(b)
(c)
(d)
The variation of the ensemble averaged values of conditional on isosurfaces for all cases are shown in Figures 4(c) and 4(d) for and , respectively, which demonstrates that predominantly assumes positive (negative) values towards the unburned (burned) gas side of the flame brush and the magnitude of increases with increasing . The quantity approaches to for small values of (i.e., ) and the mean value of remains negligible for all the isosurfaces due to the statistical planar nature of the flames. However, subgrid level curvature increases with increasing and thus the magnitude of increases with increasing values of . Relatively weak curvature dependences of and lead to positive (negative) values of and towards the unburned (burned) gas side of the flame brush due to positive (negative) value of . The contribution of resolved curvature term remains negligible in comparison to , due to relatively small values of in comparison to in statistically planar flames. Thus the contributions of and remain close to each other for all values of (see Figure 3).
The subgrid fluctuations of the surface-weighted contributions of and are taken to scale with and , respectively, to propose the following model for in this analysis: where , and are the model parameters. The function in (10) is used to capture the correct qualitative behaviour of across 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, 11] 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 (10). It has been found that enables (10) to capture the qualitative behaviour of when the optimum values of and are chosen. The optimum value of tends to increase with decreasing (increasing) . The dependences of and are reflected mostly in the resolved scale but these effects weaken with increasing values of [9, 11]. As the resolved and subgrid curvature terms are closely related [9, 11], the qualitative behaviour of is also affected by the dependences of and , which leads to the variation of the optimum values of and . The model parameter is found to decrease with decreasing values of for satisfactory quantitative prediction of , which is accounted for by expressing as . The prediction of (10) ensemble averaged on isosurfaces is compared with the ensemble averaged values of in Figure 5 for all cases for the optimum values of and , for and , when and are taken to be and . The optimum values of and are estimated by calibrating the prediction of (10) with respect to the values of obtained from DNS data and the variation of the global mean optimum values of and with for all cases are shown in Figure 6. The optimum values of and are parameterised here as where Figure 5 shows that (10) satisfactorily predicts when and , and , , and are used for and .
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(a)
(b)
(c)
(d)
(e)
Here the contribution of is scaled with (i.e., ), where the subgrid fluctuations of is taken to scale with (i.e., ). The above relations are utilised here to propose a model for (see ) in the following manner: where is the wrinkling factor [8, 10, 19], is a model parameter, and is used to capture the correct qualitative behaviour of . According to (12), vanishes when the flow is fully resolved (i.e., ). It has been found that (12) satisfactorily captures the behaviour of throughout the flame brush for in all cases considered here when a suitable value of is used. The variation of the global mean optimum values of with for all cases is shown in Figure 6. The optimum values of has been parameterised here in the following manner: where The predictions of (12) ensemble averaged on isosurfaces are compared with the ensemble averaged values of in Figure 5 all for cases for and , which show that (12) satisfactorily predicts the statistical behaviour of when and , , and are used for .
Equations (10) and (12) 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 (14) 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 obtained from DNS in Figure 7 for optimum values of , , and , where the optimum values are estimated by calibrating the models based on the ensemble averaged value of obtained from DNS data. The variations of the optimum values of , , and with for cases A–E are also shown in Figures 6(a)–6(e), respectively, which demonstrate that the model constants , , and remain greater than unity for all cases. This is found to be consistent with the realisability analysis by Hawkes and Cant [26]. Figures 6(a)–6(e) demonstrate that the optimum values of , , and change with respect to , which is also consistent with earlier findings [9]. Moreover, the optimum values of , , and for a given value of vary between cases considered here (see Figures 6(a)–6(e)). The optimum values of , , and can also be parameterised in the same manner in which is parameterised in , , and . However, this is not presented here as the models given by (7)–(9) fail to capture the positive contribution of for cases A–C. Moreover, the CSGCAND, CSGCANT, and CSGCHAR models do not capture the correct qualitative behaviour of even when the optimum values of , , and are used. The CSGCHAR model tends to overpredict the negative values of towards the unburned gas side and this behaviour becomes more prominent with increasing filter size. Figure 7 shows that for , the CSGCHAR model predicts the maximum magnitude of near the middle of the flame whereas the actual maximum magnitude of is attained slightly towards the burned gas side. The CSGCAND and CSGCANT models predict the correct magnitude of for optimum values of and , but they do not satisfactorily capture the qualitative behaviour of and underpredict (overpredict) its magnitude towards the burned gas side (middle) of the flame brush. Figure 7 demonstrates that the CSGNEW model captures the qualitative behaviour of in a better manner than the CSGCAND and CSGCANT models, and the quantitative agreement between and the CSGNEW model remains better than the CSGCAND, CSGCANT, and CSGCHAR models for the major part of the flame brush for all cases, when optimum values of , , , and are used.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
4. Conclusions
The LES modelling of the curvature term of the transport equation has been analysed using a simplified chemistry based DNS database of freely propagating statistically planar turbulent premixed flames with a range of different turbulent Reynolds numbers . The variation of is brought about by modifying and independently from each other. The statistical behaviours of the subgrid curvature term for a range of different values of have been analysed in terms of its contributions and , which arise from and , respectively. Detailed physical explanations have been provided for the observed filter size dependences of the different components of . Models have been identified for individual components of the subgrid curvature term (i.e., and ) and the performances of these models have been compared to the corresponding quantities extracted from DNS data. It has been found that the new models for and satisfactorily capture the statistical behaviours of the corresponding terms extracted from DNS data. The performance of the new model for has been found to be either better than or comparable to the performances of the existing models. It is worth noting that the present analysis has been carried out using a DNS database with moderate values of in the absence of the effects of detailed chemistry and differential diffusion. Thus, three-dimensional DNS data with detailed chemistry and experimental data at higher values of will be necessary for more comprehensive modelling of and in the context of LES.
Nomenclature
| : | Model parameter for the newly proposed model |
| : | Model parameter for Cant et al. [1] Model |
| : | Mean curvature term of flame surface density transport equation |
| : | Subgrid curvature term of flame surface density transport equation |
| : | Subgrid curvature term component due to the combined reaction and normal diffusion components of displacement speed |
| : | Subgrid curvature term due to the tangential diffusion component of displacement speed |
| : | Reaction progress variable |
| : | Model parameter for the newly proposed model |
| : | Progress variable diffusivity |
| : | Damköhler number |
| : | Gaussian Kernel used for filtering DNS data |
| : | Karlovitz number |
| : | Subgrid turbulent kinetic energy |
| , and : | Model parameter for the newly proposed model |
| : | Lewis number |
| : | Resolved flame normal vector and its th component |
| : | Model parameter for the newly proposed model |
| : | Local flame normal vector and its th component |
| : | Model parameter for the newly proposed model |
| : | Prandtl number |
| : | Turbulent Reynolds number |
| : | Subgrid turbulent Reynolds number |
| , and : | Model parameter for the newly proposed model |
| : | Displacement speed |
| : | Laminar flame speed |
| : | Normal diffusion component of displacement speed |
| : | Reaction diffusion component of displacement speed |
| : | Tangential diffusion component of displacement speed |
| : | Nondimensional temperature |
| : | Instantaneous gas temperature (dimensional) |
| : | Activation temperature |
| : | Adiabatic flame temperature |
| : | Reactant temperature |
| : | Chemical time scale |
| : | Initial eddy turnover time |
| : | Simulation time |
| : | th component of nondimensional fluid velocity |
| : | Initial root-mean-square velocity fluctuation |
| : | Subgrid scale turbulent velocity fluctuation |
| : | th Cartesian coordinate |
| : | Reaction rate of reaction progress variable. |
| : | Resolution parameter of Cant et al. [1] and Candel et al. [2] models |
| : | Model parameter for the newly proposed model |
| : | Zel’dovich number |
| : | Model parameter for Candel et al. [2] model |
| : | Model parameter for Cant et al. [1] model |
| : | Model parameter for Charlette et al. [4] model |
| : | Model parameter for the newly proposed model |
| : | Model parameter for the newly proposed model |
| : | Ratio of specific heats |
| : | Large eddy simulation filter width |
| : | DNS mesh size |
| : | Thermal flame thickness |
| : | Kolmogorov length scale |
| : | Flame curvature |
| : | Dynamic viscosity |
| : | Dynamic viscosity of unburned gas |
| : | Wrinkling factor |
| : | Density |
| : | Unburned gas density |
| : | Generalised flame surface density |
| : | Heat release parameter |
| : | Kolmogorov eddy turn-over time. |
| : | LES filtering operation |
| : | LES Favre filtering operation |
| : | LES surface averaging operation. |
| CSGCFM: | Subgrid curvature model proposed by Candel et al. [2] |
| CSGCPB: | Subgrid curvature model proposed by Cant et al. [1] |
| CSGCHAR: | Subgrid curvature model proposed by Charlette et al. [4] |
| CSGNEW: | Newly proposed subgrid curvature model |
| DNS: | Direct Numerical Simulation |
| FSD: | Flame Surface Density |
| LES: | Large Eddy Simulation |
| NSCBCs: | Navier Stokes Characteristic Boundary Conditions |
| TRZ: | Thin Reaction Zones. |
Acknowledgment
The financial support by EPSRC, UK is gratefully acknowledged.