Analysis of the Closures of Sub-grid Scale Variance of Reaction Progress Variable for Turbulent Bunsen Burner Flames at Different Pressure Levels

The statistical behaviour and modelling of the sub-grid variance of reaction progress variable have been analysed based on a priori analysis of direct numerical simulation (DNS) data of turbulent premixed Bunsen burner flames at different pressure levels. An algebraic expression for sub-grid variance, which can be derived based on a presumed bi-modal sub-grid distribution of reaction progress variable with impulses at unburned reactants and fully burned products, has been found to be inadequate for the purpose of prediction of sub-grid variance even for the flames in the wrinkled flamelets/corrugated flamelets regime. Moreover, an algebraic model, which is often used for modelling sub-grid variance of passive scalars, has been found to significantly overpredict the sub-grid variance of reaction progress variable for all the cases considered here. The modelling of the unclosed terms of the sub-grid variance transport equation has been analysed in detail. Suitable model expressions have been identified for the sub-grid flux of variance, reaction rate contribution and scalar dissipation rate based on a priori analysis of DNS data. It has been found that the alternation of pressure does not have any significant impact on the closures of sub-grid flux of variance but a model parameter for the Favre-filtered scalar dissipation rate needs to be modified to account for the variation of pressure.


Introduction
In gas turbines and internal combustion engines, premixed turbulent combustion often takes place under elevated thermodynamic pressures. However, most combustion models have been proposed and validated at atmospheric conditions because of relative ease of conducting experiments and simulations for atmospheric pressure. The flame thickness decreases and turbulent Reynolds number increases (for a given set of values of root-mean-square velocity fluctuation and integral length scale) with increasing pressure and thus it often becomes challenging to analyse high-pressure premixed flames. The advances in high performance computing have enabled detailed analyses of high pressure premixed turbulent combustion and its associated modelling using direct numerical simulations (DNS) data (Savard et al. 2017;Klein et al. 2018a, b;Chakraborty et al. 2019;Alquallaf et al. 2019). The present paper focuses on the assessment of the closure of sub-grid scale (SGS) reaction progress variable variance under elevated pressure conditions for turbulent Bunsen burner flames. The SGS reaction progress variable variance 2 v is defined as (Pierce and Moin 1998;Jimenez et al. 2001;Ranjan et al. 2016) where c is the reaction progress variable, and Q = Q∕̄ and Q denote Favre-filtered and filtered value of a general variable Q , respectively with being the gas density.
The closure of sub-grid scalar variance plays a crucial role in the modelling of micro-mixing in the context of Large Eddy Simulations (LES) (Pierce and Moin 1998;Jimenez et al. 2001). Moreover, 2 v is a quantity of pivotal importance in turbulent premixed combustion modelling Ranjan et al. 2016). The present analysis utilises an existing DNS database of turbulent Bunsen flames at different thermodynamic pressures representing the flamelet regime of combustion [flames with low Karlovitz number (i.e. Ka < 1)], which has been extended for this study by simulating a Bunsen flame with 25 times the baseline pressure. A recent study (Nilsson et al. 2019) focussed on the closure of the SGS variance 2 v and its transport for high Karlovitz number flames but such an analysis is yet to be carried out for the low Karlovitz number conditions, located at the boundary of the wrinkled and the corrugated flamelets regimes, for a range of different thermodynamic pressures. This deficit in the existing literature has been addressed in the present work with a particular focus on the modelling of pressure effects on the unclosed terms in the transport equation of 2 v for a range of different filter widths.

Mathematical Background and Numerical Implementation
An algebraic expression is often used for the evaluation of SGS variance of passive scalars in the following manner (Pierce and Moin 1998;Ranjan et al. 2016): where C = 0.5 is the model constant (Ranjan et al. 2016) and Δ is the LES filter width. It is expected that 2 v is affected by both chemical and turbulent processes Lai and Chakraborty 2016), which is not accounted for in Eq. 2. It is worth noting that 2 v assumes the following expression: (1) 2 v = CΔ 2 |∇c| 2 1 3 subject to the assumption of a bimodal sub-grid probability density function (PDF) of c with impulses at c = 0.0 and c = 1.0 . Here Da Δ = ΔS L ∕u � Δ th is the sub-grid Damköhler number with u � Δ = √ (� u i u i −ũ iũi )∕3 = � 2k sgs ∕3 and k sgs = (� u i u i −ũ iũi )∕2 being the SGS velocity fluctuation and the SGS turbulent kinetic energy, respectively. The term O Da −1 Δ can be ignored only for Da Δ ≫ 1.0 , and under that condition, the maximum possible value of the SGS variance (i.e. 2 v =c(1 −c) ) is obtained. However, the sub-grid PDF of c may not be bi-modal (Bray et al. 1985) and it might be necessary to solve a modelled transport equation of 2 v in the cases where Eq. 2 and 2 v =c(1 −c) do not provide satisfactory performance. The exact transport equation of 2 v is given by: where u j ,ẇ and D c are the jth component of velocity, reaction rate of progress variable and reaction progress variable diffusivity, respectively. In Eq. 4, � c = D c ∇c ⋅ ∇c∕̄− � D c ∇c ⋅ ∇c is the sub-grid scalar dissipation rate (SDR), whereas � N c = D c ∇c ⋅ ∇c∕̄ will henceforth be referred to as the Favre-filtered SDR in this paper. The term T 1 accounts for turbulent transport of SGS variance, whereas the term T 2 denotes generation/destruction of SGS variance due to ∇c . The term T 3 accounts for generation/ destruction of SGS variance due to chemical reaction rate, whereas T 4 depicts molecular diffusion of 2 v . The term D v = −2̄� c originates due to molecular dissipation of SGS variance under the action of SDR. In LES, the SGS flux of reaction progress variable u j c −̄ũ jc is modelled in order to solve the transport equation of c , and thus the term T 2 can be considered to be closed. The modelling of u j c −̄ũ jc has been discussed elsewhere (Gao et al. 2015a, b;Klein et al. 2016;2018c) and thus is not repeated here. The molecular diffusion term T 4 is also a closed term and thus the closure of the transport equation of 2 v depends on the modelling of T 1 , T 3 and D v . In the present analysis, the terms of Eq. 4 have been extracted to analyse their statistical behaviours by explicitly filtering the DNS data and T 1 , T 3 and D v extracted from DNS data have been utilised for assessment of their respective models.
For the purpose of this analysis, 6 different turbulent premixed Bunsen burner flames have been considered and these flames are taken from a database (which has been extended for this study) consisting of 16 different cases (Klein et al. 2018b). The chemical mechanism is simplified here using a single step irreversible reaction for the sake of computational economy so that a detailed parametric analysis involving different thermodynamic pressures can be conducted. The Reynolds number Re = U B d n ∕ u (based on the bulk inlet velocity U B , nozzle diameter d n , and the unburned-gas kinematic viscosity u ); turbulent Reynolds number Re t = u � l∕ u ; integral length scale to thermal flame thickness ratio l∕ th ; Damköhler number Da = lS L ∕ th u � ; and Karlovitz number Ka = u � ∕S L 3∕2 l∕ th −1∕2 are listed in Table 1. Here th = T ad − T 0 ∕max|∇T| L is the thermal flame thickness with T ad and T 0 being the adiabatic flame temperature and the unburned reactant temperature, respectively, and the subscript 'L' refers to the unstrained laminar flame quantities. The normalised mean inlet velocity is set to U B ∕S L = 6.0 and the normalised turbulent rootmean-square (rms) velocity fluctuation equals u � ∕S L = 1.0 , except for case D where u � ∕S L at the inlet is 3.10. It is worth noting that there was a special reason for considering moderate turbulence intensities: due to an increasing ratio of the hydrodynamic to critical length scale, the flames become increasingly hydrodynamically unstable (e.g. Darrieus-Landau (DL) instability) with increasing pressure, which is an important point for modelling high pressure turbulent premixed flames. As the effects of DL instability are masked by turbulence effects for high turbulence intensities, a small value of normalised turbulent rootmean-square (rms) velocity fluctuation u � ∕S L has been used in this work in order to be able to analyse the effects of DL instability at elevated pressure levels. The integral length scale to Bunsen burner nozzle diameter ratio is given by l∕d n = 1∕5 except for case E where l∕d n is 3∕5 . For case F this corresponds to a ratio of l∕ th = 25 and in this sense the present database has much more realistic scale separation in terms of integral scale to thermal flame thickness ( l∕ th ) than most existing DNS databases where l is of the same order as th . The heat release parameter = T ad − T 0 ∕T 0 and the Zel'dovich number = T ac T ad − T 0 ∕T 2 ad are considered to be 4.5 and 6.0 respectively, and T ac is the activation temperature. Standard values of Prandtl number ( Pr = 0.7) and ratio of specific heats ( g = 1.4 ) have been used. All non-dimensional numbers in Table 1 are to be understood as the inlet values.
It can be noted from Table 1 that the turbulent Reynolds number Re t increases from case A to case C to F. Furthermore, it can be seen from Table 1 that cases C, D and E have same values of Re t but cases D and E have one tenth of the pressure of that of case C. The value of u � ∕S L is higher in case D than in case C and E, whereas u � ∕S L and l∕ th values are exactly the same for cases C and E and thus they fall on the same point on the regime diagram (but behave very differently in terms of flame morphology). Although turbulent Reynolds number remains moderate for the cases considered here, it has been shown in the past that the model parameters for SDR reach an asymptotic value for Re t ≈ 50.0 (Chakraborty and Swaminathan 2013). Moreover, SDR models proposed using a priori analysis of moderate values of DNS data (Dunstan et al. 2013;Gao et al. 2014a) have been found to be perform well for flow configurations with much larger turbulent Reynolds number Re t based on a posteriori analysis using LES (Ma et al. 2014a).
The simulations of high pressure Bunsen flames have been conducted by adjusting the viscosity and the Arrhenius parameters in such a way that S L ∼ P −0.5 , and u ∼ P −1 are satisfied, as that of methane-air flames (Turns 2011). Thus, the Zeldovich flame thickness Z scales as Z ∼ P −0.5 (where Z = T0 ∕S L with T0 being the thermal diffusivity in the unburned gas) and the numerical resolution must be adjusted accordingly. The  (Jenkins and Cant 1999) in which the governing equations are solved using high order finite difference and Runge-Kutta time-advancement. Inflow data has been generated using a parallelised and modified version of the digital filter-based inflow methodology proposed by Klein et al. (2003) where the Gaussian filter in the temporal space has been replaced by an autoregressive AR1 process in order to avoid excessive filter length in this direction caused by the small-time step in the compressible flow solver. Interested readers are referred to Klein et al. (2018a) for further information regarding the DNS database considered here and its numerical implementation. The reacting scalar and flow fields have been initialised using an unstrained premixed laminar flame solution, which is specified as a function of radius from the nozzle centre. All the domain faces except for the inlet are specified using the partially non-reflecting NSCBC formalism (Poinsot and Lele 1992). The simulation time, when statistics were first taken, is chosen to be larger than two flow-through which corresponds to nearly four eddy-turnover times.
In this analysis, the DNS data has been explicitly LES filtered using a Gaussian filter kernel G(r) so that the LES filtered values of a general quantity Q can be calculated as follows (Peters 2000;Nilsson et al. 2019;Lai and Chakraborty 2016;Dunstan et al. 2013;Gao et al. 2014a;Ma et al. 2014a;Boger et al. 1998;Reddy and Abraham 2012;Cant 2007, 2009): Results will be presented from Δ ≈ 0.03d n (corresponding for case A to the smallest filter width to flame thickness ratio of Δ∕ th = 0.8 ) where the flame is partially resolved, up to Δ ≈ 0.15d n (corresponding for case F to the largest filter width to flame thickness ratio Δ∕ th = 19.0 ) where the flame becomes fully unresolved and Δ becomes comparable to the integral length scale l.

Results and Discussion
The isosurfaces of c = 0.5 for cases A, C, E and F for Δ∕d n = 0.06 and 0.15 are exemplarily shown in Fig. 1 along with instantaneous views of c = 0.5 isosurfaces of the corresponding cases. It can be seen from Fig. 1 that the flame morphology changes significantly with increasing pressure. The extent of wrinkling in cases C and F is significantly greater than in cases A and E due to DL instability (Klein et al. 2018a, b;Alquallaf et al. 2019) and interested readers are referred to Klein et al. (2018a, b) and Alquallaf et al. (2019) for further discussion on this aspect. Figure 1 further shows that the extent of wrinkling of the isosurface decreases with increasing Δ because of the smearing of local information due to the convolution operation associated with the filtering. For the sake of brevity, results will be shown for cases A, C and F only but results for cases B, D and E are qualitatively similar. The variations of mean values of 2 v conditional upon bins of c for cases A, C, F for different values of Δ∕d n are shown in Fig. 2. It can be seen from Fig. 2 that 2 v increases with increasing Δ before approaching an asymptotic value which is considerably smaller than the maximum possible value given by c(1 −c) . The variation of mean values of Da Δ conditional upon bins of c for cases A, C, F for different normalised filter widths Δ∕d n is shown in Fig. 3, which shows that Da Δ increases with increasing Δ but Da Δ does not assume large (i.e. Da Δ ≫ 1 ) values in all cases even for the largest value of Δ considered here. This is consistent with the analysis in Keil et al. (2019). Thus, the O Da −1 Δ contribution can only be ignored for Da Δ ≫ 1 and under that condition the maximum possible value of the SGS variance (i.e. 2 v =c(1 −c) ) is obtained. Differences between 2 v and c(1 −c) indicate the extent of the departure of the sub-grid PDF of c from the bi-modal distribution with impulses at c = 0.0 and 1.0, except for the largest filter width for case F which corresponds to Δ∕ th = 19 . This is consistent with previous findings (Dunstan et al. 2013;Gao et al. 2014a;Ma et al. 2014a;Cant 2007, 2009).
In fact, Eq. 2 significantly overpredicts 2 v and provides an unphysical prediction (i.e. 2 v >c(1 −c) ) for all filter widths for the model parameter C = 0.5 . Moreover, the ratio 2 v ∕Δ 2 |∇c| 2 has been found to increase with increasing Δ and these values are different for each of the cases considered here (not shown). Thus, an alteration of C from 0.5 to a different value will not be sufficient for algebraic closure of 2 v for all cases because the optimum value of C will not be a priori known.
The variations of {T 1 , T 2 , T 3 , T 4 and D v } × th ∕ 0 S L (with 0 being the unburned gas density) conditional upon bins of c for Δ∕d n = 0.03 and 0.15 are shown in Fig. 4 for cases A, C, F. It can be seen from Fig. 4 that D v acts as a leading order sink term for all cases irrespective of the filter width. By contrast, the reaction rate contribution T 3 acts as a leading order source term apart from the negative contribution towards the burned gas side of the flame brush for all cases for all filter widths considered here. For small values of Δ , � c = D c ∇c ⋅ ∇c∕̄−D∇c ⋅ ∇c tends to disappear as both D c ∇c ⋅ ∇c∕̄ and D c ∇c ⋅ ∇c approach D c ∇c ⋅ ∇c (i.e. lim Δ→0 D c ∇c ⋅ ∇c∕̄= D c ∇c ⋅ ∇c and lim Δ→0Dc ∇c ⋅ ∇c = D c ∇c ⋅ ∇c ) and as a result of this the magnitude of the mean contribution of D v remains small for small filter widths but the relative magnitude of D v in comparison to T 3 increases with increasing Δ as the sub-grid SDR ̃ c increases with an increase in filter width.
The mean contribution of the resolved scalar gradient term T 2 remains negative for all cases for all filter widths which is indicative of counter-gradient transport (i.e. u j c −̄ũ jc c∕ x j > 0 ). The cases considered here are subjected to weak turbulence and thus the value of Bray number N B = S L ∕u � remains greater than unity (i.e. N B > 1.0 ) for all cases. This suggests that the effects of flame normal acceleration dominate over those of turbulent velocity fluctuations leading to a dominant counter-gradient transport (Veynante et al. 1997). The magnitude of the mean contribution of T 2 remains smaller than that of T 3 for all cases considered here. The mean contribution of the molecular diffusion term T 4 remains positive on both unburned and burned gas sides of the flame brush with a negative dip in the middle of the flame brush in all cases irrespective of Δ , whereas the qualitative behaviour of the mean values of the turbulent transport term T 1 remains just the opposite to that of T 4 . This behaviour originates due to counter-gradient behaviour of the SGS flux of scalar variance, which will be discussed later in this paper. The magnitudes of mean contributions of T 4 in comparison to the magnitudes of D v decrease with increasing Δ for all cases. It can further be seen from Fig. 4 that T 3 and D v remain of the similar order of magnitude for large Δ (e.g. Δ∕d n = 0.15 where Δ ≫ th ) but this equilibrium is not maintained for small filter sizes (e.g. Δ∕d n = 0.03 where Δ ≤ th ).
The observed mean behaviours of T 1 , T 2 , T 3 , T 4 and D v have been found to be consistent with scaling estimates presented in Nilsson et al. (2019) and thus are not presented here. Moreover, the qualitative behaviours of the mean behaviours of T 1 , T 2 , T 3 , T 4 and D v have been qualitatively similar to those previously reported for high Karlovitz number detailed chemistry DNS (Nilsson et al. 2019) and simple chemistry DNS (Keil et al. 2019) data for the flames representing the thin reaction zones regime combustion. The mean behaviors of {T 1 , T 2 , T 3 , T 4 and D v } × th ∕ 0 S L conditional upon bins of c for Δ∕ th = 0.8 and 4.0 in cases A, C, F are shown in the "Appendix". It becomes obvious that the behavior of the terms in the transport equation of 2 v is very similar for different pressures when the same filter width to thermal flame thickness ratio is chosen. However, for a real LES simulation that will mean that the computational time increases drastically because the mesh size has to follow the decreasing flame thickness with increasing pressure.
Equation 4 indicates that the closure of T 1 depends on the appropriate closure of the SGS flux of variance F j = u j c 2 − 2 u j c −̄� u jc c −̄� u j � c 2 . According to the gradient hypothe- where C s = 0.18 is the Smagorinsky constant, S ij = 0.5 ũ i ∕ x j +ũ j ∕ x i is the resolved strain rate and Sc t is the turbulent Schmidt number. The mean values of F n ∕ 0 S L (where F n is the projection of the Flux F j in resolved flame normal direction −∇c∕|∇c| ) conditional  Figure 5 shows that F GHM j does not capture the qualitative behaviour of F n and it predicts the opposite sign to the value obtained from DNS, which is indicative of the counter-gradient behaviour of the SGS variance flux. These findings are expected for these flames because the effects of flame normal acceleration are likely to overcome the influences of turbulent fluctuation to yield counter-gradient transport for small values of u � ∕S L (Gao et al. 2015a, b;Klein et al. 2016Klein et al. , 2018cVeynante et al. 1997). This suggests that it will be desirable to have a model which is capable of predicting both gradient and counter-gradient behaviours of F j . It has been found in previous analyses that Clark's gradient model (CGM) (Clark et al. 1979) for SGS scalar flux (Gao et al. 2015a, b;Klein et al. 2016Klein et al. , 2018c provides satisfactory performance, which has been utilised here to propose an alternative model as: The predictions of F CGM j are also shown in Fig. 5 which indicates that this model is capable of capturing both gradient and counter-gradient behaviours of SGS flux of variance. However, the quantitative agreement between DNS data and model prediction deteriorates with increasing Δ. A model expression, which was originally proposed by Chakraborty and Swaminathan (CSM) (2011) for the Reynolds flux of variance in the context of RANS, has been extended for LES in the present work in the following manner. According to Bray Moss Libby modelling (Bray et al. 1985) the Reynolds flux of variance can be expressed based on the presumed bi-modal PDF of c with impulses at c = 0 and c = 1.0 as: where ⟨q⟩ and are the Reynolds averaged and Favre-averaged values of a general quantity q , respectively and ⟨q⟩ P and ⟨q⟩ R are conditional mean values of the quantity q in products and reactants, respectively. However, Eq. 8 does not perform well for Da < 1 where the PDF of c does not remain bi-modal (Bray et al. 1985). The departure from bimodal PDF can be quantified with the help of a segregation factor so that g increasingly deviates from unity with decreasing Damköhler number Da . In order to extend Eq. 8 for Da < 1 combustion, Chakraborty and Swaminathan (2011) proposed the following model expression and validated it with the help of a priori analysis: It can be seen from Eq. 9 that this expression reduces to Eq. 8 when g = 1.0 which is likely to be realised for Da ≫ 1 where the PDF of c is likely to be bi-modal. In the context of LES, the model can be written as:  (Bray 1980), with f (c) being the burning mode PDF, which can be taken as any appropriate continuous function according to Bray (Bray 1980). This thermo-chemical parameter c m has been found to be 0.84 for cases A-F. Accordingly, T 3 can be expressed as: T 3 = 2̄ẇ c m −c . The variations of mean values of T 3 conditional upon bins of c in cases A, C, F are shown in Fig. 6 for Δ∕d n = 0.03 and 0.15, which shows that 2̄ẇ c m −c (with ̄ẇ extracted from DNS data) satisfactorily captures both qualitative and quantitative behaviours of T 3 for all Δ in all cases considered here. According to Bray (Bray 1980), the mean reaction rate in the context of RANS for Da ≫ 1 can be modelled by . It has been shown in Dunstan et al. (2013), Gao et al. (2014a) and Ma et al. (2014a) based on a priori analysis of DNS that the expression ̄ẇ = {2̄Ñ c ∕(2c m − 1)} performs well for Δ∕ th ≫ 1 and the agreement between ̄ẇ and {2̄Ñ c ∕ 2c m − 1 } improves with increasing Δ , as sub-grid Damköhler number Da Δ = ΔS L ∕u � Δ th increases with an increase in LES filter width (see Fig. 3). However, {2̄Ñ c ∕ 2c m − 1 } does not adequately predict ̄ẇ for small filter widths (i.e. Δ∕ th < 1 ) when Da Δ assumes small values. In order to obtain a single expression which predicts ̄ẇ accurately for both small and large filter widths the following model expression was suggested in Dunstan et al. (2013), Gao et al. (2014a) and Ma et al. (2014a): where T is the temperature, and = 0.56 th S L ∕ T0 is a model parameter with T0 being the thermal diffusivity in the unburned gas. The exponential terms exp − Δ∕ th and 1 − exp − Δ∕ th act as bridging function such that one obtains ̄ẇ = {2̄Ñ c ∕ 2c m − 1 } for Δ∕ th ≫ 1 , and for Δ → 0 , Eq. 11 reduces to lim c, T) . It has been shown elsewhere Ma et al. 2014a) that the prediction of Eq. 11 shows better agreement with ̄ẇ extracted from DNS data only for Δ > th and the quantitative agreement improves with increasing filter width. Using Eq. 11 in T 3 = 2 ẇc −̄ẇc = 2̄ẇ c m −c yields: Figure 6 shows that T 3 can be reasonably predicted by this expression and the quantitative agreement improves with increasing Δ . This behaviour arises due to improved prediction of ̄ẇ by Eq. 11 for large values of Δ (e.g. Δ = 0.15d n ) (Dunstan et al. 2013;Gao et al. 2014a;Ma et al. 2014a).
The transport equation of Ñ c takes the following form (Gao et al. , 2015c where u j is the jth component of velocity vector and the terms on the left-hand side denote the transient effects and resolved advection of Ñ c respectively. The term D 1 = ∇ ⋅ ( D∇N c ) represents the molecular diffusion of Ñ c and the terms T � 1 , T � 2 , T � 3 , T � 4 , −D 2 and f (D) are unclosed and given by: The term T ′ 1 represents the effects of sub-grid convection, whereas T ′ 2 denotes the effects of density-variation due to heat release. The term T ′ 3 is determined by the alignment of ∇c with local strain rates e ij = 0.5 u i ∕ x j + u j ∕ x i and this term is commonly referred to as the scalar-turbulence interaction term. The term T ′ 4 arises due to the correlation between ∇ẇ and ∇c , whereas −D 2 denotes the molecular dissipation of scalar dissipation rate and these terms will henceforth be referred to as the reaction rate term and dissipation term respectively. The term f (D) denotes the effects of D variation. It has been demonstrated based on scaling arguments by Gao et al. (2014b) that an equilibrium is maintained between the T � 2 + T � 3 + T � 4 + f (D) and −D 2 for large filter widths (i.e. Δ∕ th ≫ 1 ). A similar conclusion can be reached if the transport equation of is analysed and the terms equivalent to T � 2 + T � 3 + T � 4 + f (D) and −D 2 are considered (Chakraborty et al. 2008). Kolla et al. (2009) and  considered the modelled expressions for the terms equivalent to T � 2 , T � 3 , T � 4 − D 2 + f (D) in the context of RANS and utilised their balance (i.e. T � 2 + T � 3 + T � 4 + f (D) − D 2 ≈ 0 ) to propose a model expression: where C * 3 , C * 4 and ′ c are the model parameters and N c L f (c)dc is a thermo-chemical parameter . Equation 15 has been extended to LES based on the equilibrium between the T � 2 + T � 3 + T � 4 + f (D) and −D 2 for large filter widths (i.e. Δ∕ th ≫ 1 ) in the following manner (Dunstan et al. 2013;Gao et al. 2014a;Ma et al. 2014a): The model parameters C * 3 , C * 4 and c are given by the following expressions Ma et al. 2014a): where This suggests that the model parameter c (and ′ c ) is likely to be dependent on flame curvature description in the underlying computational methodology. As flame wrinkling is partially resolved in LES in contrast to the total sub-grid curvature contribution in RANS, the modification from ′ c to c in the model expression of scalar dissipation rate is not completely unexpected. It is important to note that the model parameters C * 3 , C * 4 and c have been calibrated using a large number of DNS cases including variations of , Le, Da and Ka (see Gao et al. 2014a) in the past, and thus, the dependence of Ka is explicitly addressed.
In Eq. 16i (Dunstan et al. 2013;Gao et al. 2014a;Ma et al. 2014a), N c L f (c)dc is a thermo-chemical parameter , which is found to be K * c = 0.77 for cases A-F for the present thermo-chemistry and f b = exp −0.7 Δ∕ th 1.7 is a bridging function, that ensures that Ñ c approaches D∇c ⋅ ∇c , when the flame is fully resolved (i.e. lim Δ→0Ñc = D∇c ⋅ ∇c ). The model parameters C * 3 , C * 4 and c are given by Eq. 16ii Ma et al. 2014a The predictions of Eq. 16i are compared to the mean values of Ñ c extracted from DNS conditional upon the bins of c in Fig. 7 for Δ∕d n = 0.03 and 0.15, which shows that this model satisfactorily predicts Ñ c for small filter widths (e.g. Δ∕d n = 0.03 ) where the resolved component (i.e. D ∇c ⋅ ∇c ) plays the dominant role but over predictions can be observed for large filter widths (e.g. Δ∕d n = 0.15 ). This is consistent with previous findings Ma et al. 2014a), which also showed that Eq. 16i overpredicts the mean values of Ñ c conditional upon c for small values of u � ∕S L but this overprediction disappears for large values of u � ∕S L . The extent of this overprediction is particularly severe for large filter widths (e.g. Δ∕d n = 0.15 ) for high pressure cases B, C and F and the extent of overprediction increases with increasing pressure.
A given value of Δ∕d n implies larger value of Δ∕ th for higher pressures as the flame thickness th decreases with increasing pressure. The overprediction of the mean values of (16i) 1 3 N c conditional upon c according to Eq. 16i in cases B, C and F has been found to be comparable to that in cases A, D and E for a given value of Δ∕ th (not shown here). It has been shown elsewhere Chakraborty et al. 2008) that the model parameter c originates due to the curvature m contribution to the SDR transport and m dependence of S d plays a key role in this contribution. It has been shown and discussed elsewhere (Klein et al. 2018b;Alquallaf et al. 2019) that the flame curvature statistics are significantly affected by pressure. The likelihood of obtaining DL instability increases with increasing pressure (Klein et al. 2018a, b;Alquallaf et al. 2019) and this is reflected in the increased skewness of negative curvature values in cases B, C and F (not shown here but refer to Klein et al. 2018a, b;Alquallaf et al. 2019). Thus, some pressure dependence of c is not unexpected. It has been found that the value of c , which ensures that the volume-integrated ̄Ñ c obtained from DNS data can be satisfactorily predicted, needs to increase with increasing pressure. This has been empirically parameterised as: It can be seen from Fig. 7 that the use of c expression given by Eq. 18 in Eq. 16 provides better agreement with DNS data. Furthermore, Fig. 7 shows that using Eq. 16i with c given by Eq. 18 in Eq. 12 provides satisfactory prediction of T 3 in all cases. It is worth noting that T 3 is overpredicted if Eq. 16i is used without the pressure dependence of c and thus is not explicitly shown here. Finally, as D ∇c ⋅ ∇c is a resolved quantity, Eq. 16i with c according to Eq. 18 can be used to close the molecular dissipation term In the present analysis all statistics are taken from the whole flame. This is based on the finding that there is little variation between different axial locations for a variety of quantities (Klein et al. 2018d). Finally, it is worth noting that the models discussed in this paper need to be implemented in actual LES simulations because numerical and modelling errors interact in a complex manner so the applicability of these models cannot entirely be assessed based on a priori DNS analysis. Interested authors are referred to (Vervisch et al. 2008) for discussion on the coupling between the numerical discretization of scalar field transport and the modelling of unresolved sub-grid scale fluctuations of chemical species. The novelty of the current analysis lies in the analysis of transport equation-based closure of SGS scalar variance for turbulent premixed flames under elevated pressures, which is yet to be carried out in the existing literature. Therefore, model capabilities are a priori assessed in this analysis without the interference (18) c = max 2∕ 2c m − 1 , p∕p 0 0.37 1.05 ∕( + 1) + 0.51 4.6 Fig. 7 Variations of mean values of Ñ c × th ∕S L (red) and the predictions of Eq. 16i with c given by Eq. 16ii (green) and c according to Eq. 18 (blue) conditional upon bins of c for Δ∕d n = 0.03 and 0.15 in cases A, C, F of numerical errors because a posteriori assessment of models using LES is intrinsically code dependent as a result of the interaction of numerical and modelling errors (Vervisch et al. 2008). It is worth noting that the same approach has been adopted in several previous studies (Ranjan et al. 2016;Nilsson et al. 2019;Reddy and Abraham 2012;Veynante et al. 1997;Charlette et al. 2002a;Selle and Bellan 2007;Lignell et al. 2009;Chatakonda et al. 2013;Lapointe and Blanquart 2017;Aspden et al. 2019;Devaud et al. 2019;Bushe et al. 2020) by other authors and there are several examples where the models proposed based on a priori DNS analysis (e.g. Gao et al. 2014a;Charlette et al. 2002a) have been demonstrated to perform well based on a posteriori assessment (e.g. Charlette et al. 2002b;Ma et al. 2013Ma et al. , 2014aButz et al. 2015). A recent analysis by different authors (Nilsson et al. 2019) on the modelling of SGS variance transport has also been conducted entirely by a priori assessment and the same approach has been adopted in the current study.
It is also worth mentioning that SGS variance transport has been adopted in several LES simulations Langella et al. , 2018Massey et al. 2018;Chen et al. 2019a;b, 2020;Massey et al. 2019) in the past for atmospheric flames, and these simulations successfully used several model expressions considered in this analysis (e.g. Eqs. 6, 16i and 16ii) for closures of SGS flux of variance and SDR. The modelling of the chemical reaction contribution T 3 is reliant upon the SDR based filtered reaction rate closure given by Eq. 11. This SDR based reaction rate closure was also previously successfully implemented in LES albeit for atmospheric premixed turbulent combustion (Ma et al. 2014a;Butz et al. 2015) and the SDR closure given by Eqs. 16i and 16ii has been used in several previous LES studies (Ma et al. 2014a;Butz et al. 2015;Langella et al. 2018;Massey et al. 2018;Chen et al. 2019aChen et al. , b, 2020Massey et al. 2019) for atmospheric premixed combustion. The current analysis justifies the choices of the model expressions given by Eqs. 11, 16i and 16ii based on a priori DNS analysis. The only two expressions which are yet to be a posteriori analysed by implementing them in LES are the SGS variance flux model given by Eq. 10 and the pressure correction in the SDR modelling given by Eq. 18. A close inspection of Eq. 10 reveals that the success of this model depends on the successful modelling of SGS flux of reaction progress variable u j c −̄� u jc and the modelling of u j c −̄� u jc has been addressed elsewhere (Allauddin et al. 2017) by simultaneous a priori DNS and a posteriori LES analyses.
It is important to note that an experimental database, which offers information regarding SGS variance, and its sub-grid flux and dissipation rate for a range of elevated pressure conditions, is rare in the existing literature. Thus, a posteriori assessment of the models of the unclosed terms in the SGS variance transport equation is not a straightforward task because a discrepancy between LES results and experimental measurements may not originate from the inaccuracy in the SGS variance modelling and similarly a good agreement between LES and experimental results might be obtained due to fortuitous cancellation of errors. However, notwithstanding the above comments, the model expressions identified in this analysis need to be implemented in LES for an experimental configuration for which measurements are available for the assessment of the performances of the SGS variance closure.

Conclusions
The modelling of SGS variance of reaction progress variable and its sensitivity to thermodynamic pressure have been analysed based on a priori analysis of DNS data of turbulent premixed Bunsen burner flames at different pressure levels. An algebraic expression of the SGS variance based on the presumed bi-modal distribution with impulses at c = 0.0 and 1.0 has been found to be inadequate even for the flames in the corrugated/ wrinkled flamelets regime. An alternative algebraic expression of SGS variance, which is usually used for passive scalar mixing, has been shown to significantly overpredict the corresponding quantity extracted from DNS data. The statistical behaviours of the unclosed terms of the SGS variance transport equation have been analysed. It has been found that the reaction rate contribution plays a leading order role for all filter widths and remains of the same order as that of the molecular dissipation contribution for large filter widths. The modelling of the unclosed terms of the SGS variance transport equation has been assessed based on a priori DNS analysis and suitable models have been identified for the SGS flux of variance, reaction rate contribution and scalar dissipation rate. One of the model parameters of a widely used algebraic closure of the Favrefiltered SDR closure (Dunstan et al. 2013;Gao et al. 2014a;Ma et al. 2014a) has been found to be dependent on pressure, whereas the model performance for the SGS flux of variance has been found to be independent of pressure variation. As the effects of Darrieus-Landau instability are likely to become weak for high values of Ka , the modelling methodology needs to be assessed for high pressure flames under high Karlovitz number. Moreover, these model expressions need to be implemented in actual LES for a posteriori assessment for comprehensive validation, which will form the basis of future analyses.