Algebraic Flame Surface Density Modelling of High Pressure Turbulent Premixed Bunsen Flames

Performance of representative algebraic flame surface density (FSD) models have been a-priori assessed based on a direct numerical simulation database consisting of four turbulent premixed Bunsen flames at four different pressure levels. Results indicate that for a given resolution of the flame front, the considered algebraic FSD closures perform in a qualitatively similar manner irrespective of pressure variation. However, for a given computational mesh, the performance deteriorates as flame thickness decreases with increasing pressure. Additionally, the contribution of the subgrid scale surface-filtered density-weighted tangential diffusion component of the displacement speed tends to increase with pressure increment, further complicating the modelling of FSD. Further analysis also indicate that the inner cut-off scale normalized by the thermal flame thickness increases with increasing pressure, possibly due to the presence of Darrieus–Landau instability which is more likely to occur at higher pressure due to a larger ratio of hydrodynamic length scale to thermal flame thickness.


Introduction
In many engineering applications like industrial gas turbines and spark ignition engines, premixed combustion takes place at elevated pressures under turbulent conditions. Increased pressure affects both, flame and turbulence characteristics alongside chemical kinetics. As an example, for hydrocarbon fuels, laminar burning velocities and flame thicknesses typically decrease with a pressure rise (Fragner et al. 2015). This in turn affects kinematic viscosity, the Reynolds, Damköhler and Karlovitz numbers, reaction rates and flame speed, extinction, ignition delay time and pollutant formation (Bougrine et al. 2011).
The reduction of the flame thickness promotes the hydrodynamic instability which has an impact on flame morphology, flame area and in particular turbulent flame speed (Creta et al. 2016;Klein et al. 2018) which is a quantity of fundamental interest. However, a large part of the existing literature on numerical and theoretical analyses, as well as model development in the context of premixed combustion has been carried out at ambient pressure conditions. Further, in high-pressure test-rigs optical access is limited, laser and signal beams must pass through several windows often resulting in a signal deterioration (Stopper et al. 2013) making the acquisition of reliable results more difficult and more expensive (Bougrine et al. 2011;Stopper et al. 2013). In this regard the present analysis, together with previous investigations by these authors (Klein et al. 2018a, b), aims at understanding and modelling the effects of elevated pressure on flame-turbulence interaction based on Direct Numerical Simulations (DNS) data. Under several simplifying assumptions, the complexity of chemistry can be reduced by introducing a reaction progress variable c assuming the values c = 0 on the reactant side and c = 1 in the fully burned products (Peters 2000). In such a framework, the transport equation for the filtered reaction progress can be written as: Here , u, D c and ̇c denote density, velocity, progress variable diffusivity and reaction rate respectively. Further, Q and Q = Q∕ refer to filtering and Favre filtered values of a general quantity Q , repectively. Both terms on the right hand side of Eq. 1 are unclosed. For the treatment of the turbulent scalar flux u c − ũc in the context of LES, the reader is referred to Gao et al. (2015), Klein et al. (2016Klein et al. ( , 2018, Allauddin et al. (2017), and references therein, whereas the combined modelling of subgrid transport and filtered flame propagation is discussed in Klein et al. (2016). Several methodologies exist for turbulent premixed combustion modelling such as G-equation (Pitsch and Duchamp de Lageneste 2002), artificially thickened flame (Colin et al. 2000) and FSD (Boger et al. 1998) based closures. In the FSD framework, the sum of filtered molecular diffusion and filtered reaction rate are expressed as: where The quantities Σ gen and Ξ are referred to as the generalized FSD (Boger et al. 1998) and the wrinkling factor, respectively and (Q) s = Q |∇c|∕Σ gen denotes the surface-weighted filtered value of the general quantity Q (Boger et al. 1998). The displacement speed S d in Eq. 2 can be decomposed into reaction, normal and tangential diffusion components (Echekki and Chen 1999;Peters et al. 1998) with the mean curvature of the respective isosurface m : This gives rise to the expression S d s = S r + S n s + S t s in Eq. 2. Often S r + S n s is approximated as 0 S L , with the unburned gas density 0 and the unstrained (2) ∇ ⋅ D c ∇c +̇c = S d s Σ gen = S d s Ξ |∇c| , (3) Ξ = |∇c|∕|∇c| and Σ gen = |∇c| .
laminar burning velocity S L . Assuming D c is constant (Peters 2000) on a given c-isosurface, S d s can be rewritten as: It has been shown recently (Chakraborty et al. 2018) that the curvature contribution to the flame displacement speed potentially has a non-negligible influence on the propagation of turbulent premixed Bunsen flames, however the quantity m s is unclosed. It is possible to split the surface-filtered curvature into a resolved contribution (superscript res ) given by the divergence of the surface-filtered flame normal vector and a subgrid scale contribution (superscript sgs ) as follows (Chakraborty and Cant 2007;Chakraborty and Stewart Cant 2009;Ma et al. 2013): where the normal and the surface-filtered normal vectors are defined as Here, it is important to note that the flame curvature can be expressed as the divergence of the flame normal vector i.e. m = ∇ ⋅ N∕2 . The surface-filtered flame normal vector is itself unknown and depends on the modelling of c and Σ gen according to Eq. 8. Alternatively, the resolved contribution of curvature could be approximated in terms of the filtered reaction progress variable c as follows: It is worth noting that the expression for c depends only on c . Hence, the resolved normal can alternatively be expressed as (− c∕ x i )∕|∇c| . The modelling of the sum of filtered molecular diffusion and filtered reaction rate (see Eqs. 2, 3) requires modelling of both S d S and Σ gen or Ξ . These aspects will be discussed in Sects. 3-5. A large variety of algebraic models for Σ gen closure is available in the existing literature and has been analyzed in Chakraborty and Klein (2008). The focus of this work is to (1) analyze the effect of pressure on the performance of (a subset of) these models; (2) to quantify the magnitude of the resolved and subgrid scale curvature contribution (see Eqs. 6,9) to the surface-filtered density-weighted displacement speed; (3) to contribute to the understanding of flame wrinkling in high pressure flames, possibly under the influence of the Darrieus-Landau instability.

Numerical Method and Database
A compressible DNS code, known as SENGA (Jenkins and Stewart Cant 1999), has been used for the present analysis, which employs high order finite difference and Runge-Kutta schemes for spatial discretization and explicit time advancement, respectively. The chemical mechanism is simplified by a single-step irreversible reaction (i.e. Reactants → Products ) for the purpose of computational economy. Four turbulent premixed Bunsen flames at different pressure levels are considered for this analysis. Cases A, B, C and D correspond to pressure levels of P 0 , 5P 0 , 10P 0 and 25P 0 , respectively and for simplicity they are henceforth referred to 1, 5, 10 and 25 bar cases in the following. The effect of pressure has been taken into account by adjusting viscosity and Arrhenius parameters in such a way that S L ∼ p −0.5 , u ∼ p −1 , similar to the behavior of methane-air flames, where u is the kinematic viscosity in the unburned gas. Consequently, the Zel'dovich flame thickness L scales as L ∼ p −0.5 and the grid resolution is adjusted accordingly. The high computational cost of simulating elevated pressure flames justifies the simplified treatment of chemistry used in this work and it has been shown in the past that the qualitative nature of the flame-turbulence interaction does not change if detailed chemistry and transport are used instead (see details mentioned in Klein et al. 2018a 1272 3 ) . The nozzle diameter D corresponds roughly to half the domain length. Inflow data has been generated using a modified version of the method suggested by Klein et al. (2003), where the Gaussian filter in axial direction has been replaced by an autoregressive AR1 process to avoid excessive filter length in this direction caused by the small time step in the compressible flow solver. All other boundaries are taken to be partially non-reflecting. An unstrained premixed laminar flame solution with the geometry of a semi-sphere located at the inflow is used as initial data. The simulation time, when first statistics are taken, is chosen to be larger than the maximum of two flow through and two eddy turnover times. The global Lewis number has been set to Le = 1.0 for the present analysis and standard values have also been used for the Prandtl number (Pr = 0.7) and ratio of specific heats ( g = 1.4) . The heat release parameter and the Zel'dovich number have been set as = (T ad − T 0 )∕T 0 = 4.5 and = 6 respectively, where T ad and T 0 denote the adiabatic flame and fresh gas temperature, respectively. The Reynolds number Re = U B D∕ u based on the bulk inlet velocity U B = 6 S L , nozzle diameter D and u ranges from 399 for case A to 1995 for case D. The values of Damköhler number Da th = l S L ∕ th u � , and Karlovitz number Ka th = (u � ∕S L ) 3∕2 (l∕ th ) −1∕2 are given in Table 1. Here, the velocity fluctuation intensity has been set to u � = S L , the integral length scale to nozzle diameter is l∕D = 1∕5 and th = (T ad − T 0 )∕max|∇T| L denotes the thermal flame thickness, where subscript L refers to the unstrained laminar flame quantities. All non-dimensional numbers mentioned before have to be understood as inlet values. All the considered cases fall on the boundary of the wrinkled and the corrugated flamelets regimes according to the regime diagram by Peters (2000) but due to decreasing flame thickness, l∕ th increases from 5 for case A to almost 26 for case D, offering a considerable range of scale separation for high pressure flames in particular. Moderate turbulence intensities have been considered to effectively assess the influence of Darrieus-Landau instability which is masked by strong velocity fluctuations (Klein et al. 2018). Figure 1a shows instantaneous views of iso-c surfaces for the considered cases, while the curvature PDFs for the four Bunsen flames are shown in Fig. 1b. The relatively strong negative skewness of the curvature PDF in cases C and D, together with the shown flame surface morphology provide a clear indication of Darrieus-Landau instability for these cases, as discussed in detail in Klein et al. (2018). Due to excessive postprocessing overheads associated with very large filter-width evaluations, results for case D are presented exemplarily. . 1 a Instantaneous iso-c surfaces for cases A, B, C and D (left to right). The value of c increases from 0.1 (dark orange) to 0.9 (transparent yellow). Cases A, B, C and D represent a 1 Bar, 5 bar, 10 bar and 25 bar Bunsen flame. b Curvature PDFs for the considered cases. The curvature m is normalized with the thermal flame thickness th of the corresponding flame

Performance of Existing Algebraic FSD Models
A set of various algebraic FSD models were assesed for statistically planar flames at ambient pressure conditions in Chakraborty and Klein (2008) and Klein and Chakraborty (2019). In this paper, we assess a subset of three algebraic FSD models, as mentioned in Table 2, for the Bunsen flames given in Table 1. It is remarked that the FZ model, was not precisely suggested in this form for the generalized FSD by Zimont and Lipatnikov (1995), which was adopted for LES analysis by Fureby (2005) for estimating the wrinkling factor from the turbulent flame speed expression provided by the RANS model proposed in Zimont and Lipatnikov (1995). This will henceforth be referred to as the Fureby extended Zimont model (i.e. FZ model) (Fureby 2005).
Beside these models, three more models (Angelberger et al. 1998;Charlette et al. 2002;Gülder 1990) have been tested for the present database. The performance of the model by Gülder (1990) is qualitatively similar to the FZ model (Fureby 2005). The model by Charlette et al. (2002) behaves similar to the model by Fureby (2005), whereas the performance of the model by Angelberger et al. (1998) is similar to the model by Pocheau (1994). The subset of these 6 models has been analyzed in Klein and Chakraborty (2019) and for the sake of brevity only results for the Pocheau (1994), Fureby (2005) and FZ model (Fureby 2005) models are presented here. Their performances have been tested for the considered cases at four different pressure levels. For this purpose, the DNS database has been explicitly filtered using a Gaussian filter function with a filter size ranging from ∕ th = 0.8 where the flame is still partially resolved, up to ∕ th = 19.2 (for case D) where the flame is fully unresolved. For this analysis all gradients of filtered variables have been calculated on the DNS grid. Figure 2 exemplarily shows correlation coefficients between model predictions and Σ gen evaluated from DNS for cases A and D. The difference in correlation coefficients is small between the different models and for this reason only average values are shown. It is well known and also evident from Fig. 2 that the correlation coefficient decreases with increasing filter size. Comparing the 1 bar and 25 bar cases (i.e. cases A and D), it becomes clear that correlation coefficients are comparable, for similar values of ∕ th . However, it is important to understand that the same LES grid for a 1 bar and 25 bar flame results in a five times larger value of ∕ th for case D compared to case A as flame thickness decreases with increasing pressure (i.e. th ∼ p −0.5 ) . Therefore, assuming the same computational grid for cases A and D, the same model is likely to exhibit worse performance for higher values of pressure.
It is important for a FSD model to predict the flame surface area accurately. Figure 3 shows the volume averaged value of modelled FSD (volume-averaging is represented by angled brackets) normalized with the corresponding DNS value for all cases, with three Table 2 List of considered algebraic FSD models, where L = T 0 ∕S L is the Zel'dovich flame thickness with T 0 being the thermal diffusivity in the unburned gas Fureby (2005) Σ  Fig. 3. The three considered FSD models reflect a qualitatively similar behavior in relation to filter size variation, i.e. as filter size increases, the flame surface area predicted by a particular model increasingly deviates from the DNS value. For the same value of ∕ th , a particular model reflects similar level of accuracy across pressure variations. However, it is important to note that the thermal flame thickness decreases with increasing pressure (i.e. th ∼ p −0.5 for the current database). As computing resources are often limited, higher pressure cases are likely to have higher values of ∕ th , as compared to low pressure cases. Therefore, for the same grid resolution, the models will have higher uncertainty and inaccuracy for increasingly higher pressure. The local performance of the models can be assessed by analyzing the behavior of the FSD conditionally averaged over the span of the reaction progress variable, i.e. for 0 ≤ c ≤ 1 . The models considered in this work depend on a wrinkling factor multiplied by the resolved FSD (i.e. |∇c| ). Hence, the profile of the modelled FSD (i.e. Σ mod gen ) depends largely on the shape of |∇c| , which is identical for all the three considered models. For ambient pressure flames, the models discussed in this paper (see Table 1) and their local behaviour (i.e. Σ mod gen versus c ) have been shown in references Chakraborty and Klein (2008) and Klein and Chakraborty (2019) and varying pressure does not change the general behaviour discussed there. Therefore, the profiles of Σ gen conditioned upon c are not explicitly shown in this paper and interested readers are referred to Chakraborty and Klein (2008) and Klein and Chakraborty (2019) for further information.

Resolved and SGS Curvature Contribution to Displacement Speed
In the context of FSD modelling, surface-filtered density-weighted displacement speed is often modelled as S d s = 0 S L (Boger et al. 1998;Cant 2000, 2001). However, this approximation ignores the direct effect of the curvature on the displacement speed, via Eq. 4. High fidellity modelling of ∇ ⋅ D c ∇c +̇c requires also the modelling of S d s (Chakraborty and Klein 2009;Klein and Chakraborty 2006). It is preferable to approximate S r + S n s = 0 S L and accordingly S d s can be modelled as: The influence of mean curvature in turbulent premixed Bunsen flames on S d s has been discussed in Chakraborty et al. (2018) and it was found for the present database that the negative mean curvature may yield a value of S d s which is up to 10% greater than 0 S L . As pointed out earlier m s is an unknown quantity in LES and has to be approximated with the help of either Eq. 6 or 9. Figure 4 shows the mean of the normalized magnitude of different decompositions of m s for cases A, B and C, i.e. mean( |C ∕ m s | ) . Results for case D follow the same trend but are not shown here explicitly. In Fig. 4, C res s (C sgs s ) refers to the resolved (SGS) part in the decomposition given by Eq. 6, where Σ gen has been determined based on the explicitly filtered DNS data. Ultimately, the value of Σ gen needs to be modelled in the expression C res s and the resulting expression using the Fureby model is called C res s,m . Finally, the first part in Eq. 9, where resolved curvature is calculated based on filtered reaction progress variable c , is henceforth denoted as C res c . Figure 4 shows that both expressions of resolved curvature represent reasonably good approximations to m s for small filter widths, but all models tend to deviate for larger filter widths. The expressions C res s and C res c are more or less equivalent in performance and C res s,m obviously depends on the underlying model for FSD. The SGS contribution to curvature |C sgs s ∕ m s | is zero for small filter sizes but can (in a mean magnitude sense) reach about 50% of m s . It is noted that the range of filter sizes is again increased from case A to C to reflect the fact that the same resolution of the flame front will require considerably smaller grid spacing for high pressure flames.
Beside quantifying the magnitude of the curvature contributions in Fig. 4, it is instructive to look at the average of the true values of curvature components normalized by the laminar flame speed, i.e. mean( − D c C res/sgs s ∕ 0 S L ) . In this manner the net effect of resolved and SGS contributions to S t ∕ 0 is quantified. Figure 5a shows that the resolved contribution of flame curvature (i.e. − D c C res/sgs s ∕ 0 ) amounts to 5-10% of S L for the flames considered here, and this value is more or less independent of the filter width. Although the magnitude of the subgrid contribution can attain considerable values as shown in Fig. 4, it can happen that positive and negative contributions cancel each other in an average sense resulting in a small mean contribution. However, Fig. 5b demonstrates that this is not the case entirely and that there is a net mean negative contribution of (the signed) C sgs s of about 5% the laminar flame speed for large filter width. A possible explanation for this behavior is the asymmetry of the curvature PDF (see Fig. 1b). The values of C sgs s seem to be comparable for different pressure levels at the same normalized filter size, but again, for the high pressure case it is likely that one has to deal with higher values of ∕ th . It is also worth noting that the Darrieus-Landau instability observed in case C (see Klein et al. 2018 for a discussion in this regard) with its localized regions of sharp negative curvature does not seem to have an additional effect on the magnitude of C sgs s . Figure 4 revealed that both definitions of resolved flame curvature, i.e. C res s (respectively C res s,m ) and C res c yield reasonable and comparable representations of mean m s . In order to assess their local behavior, correlation coefficients between these quantities and m s were computed over a range of filter sizes for cases A, B and C, where it is evident that both decompositions of m s lead to comparable results. Figure 6 shows correlation coefficients between these quantities and m s for a range of filter sizes for cases A, B and C. It is evident from Fig. 6 that both decompositions of m s work comparably. It is important to note that the FSD transport equation modelling typically makes use of the surfacefiltered normal vectors rather than the normalized gradient of resolved reaction progress variable. Hence, the latter quantity cannot directly be used in this context.

Power Law Based Wrinkling Factor Modelling
In power-law based modeling (e.g. Fureby 2005) the wrinkling factor Ξ is often expressed as Ξ = 0 ∕ i D f −2 , where i and 0 are inner and outer cut-off scales, respectively, and D f is the fractal dimension of the premixed flame. In the context of LES, the outer cut-off scale is taken to be the LES filter width . It will be interesting to analyze the fractal dimension and the inner cut-off scale for high pressure flames and in particular the possible influence of the Darrieus-Landau instability. Figure 7 shows the variation of ⟨Σ gen ⟩∕⟨�∇c�⟩ as a function of ∕ th on a log-log scale. It is evident from the figure that the SGS flame wrinkling is not only a property of the physical flame front, but also depends on the filter-width . For the largest filter-width corresponding to the case D, a value of 1.7 is obtained for the SGS flame wrinkling. It can be seen from Fig. 7 that the variation of log � ⟨Σ gen ⟩∕⟨�∇c�⟩ � with log ∕ th becomes almost linear for relatively large filter sizes, and the best fit straight line to this range can be used to estimate the inner cut-off scale and fractal dimension following a procedure explained in Chakraborty and Klein (2008). Table 3 shows the fractal dimension and inner cut-off scale for cases A, B, C and D, and it becomes clear that although the fractal dimension increases to a very small extent, the inner cut-off scale, normalized with thermal flame thickness, increases significantly with increasing pressure. Statistical averaging and curve fitting give rise to some uncertainty in these numbers but the qualitative trends remain unaffected. It becomes clear that i ∕ th follows a scaling of the form Ka −b with a positive number b as discussed in Chakraborty and Klein (2008). The observed increment in i ∕ th can be possibly attributed to the effects of Darrieus-Landau instability, which give rise to sharply negative localized curved cusps with large wavelengths positively curved bulges in between. Cases A-D have the same value of u � ∕S L = 1 , and thus the value of u � ∕S L tends to unity at larger filter width for these cases. As a consequence, the parameterization of D f given in Table 2  Correlation strength between m s and different terms representing resolved mean curvature (see Eqs. 6, 9) for a case A; b case B; c case C. The data presented here is calculated conditionally over the range 0.7 < c < 0.9 independent of pressure. However, the present analysis suggests a moderate increase in D f with increasing pressure.

Conclusions
Algebraic flame surface density based modelling of high pressure turbulent premixed Bunsen flames has been analyzed in this work using a-priori analysis of DNS data. The main findings can be summarized as follows: (1) Table 3  pressure. Results indicate that for a fixed value of ∕ th , existing algebraic FSD models show qualitatively similar performances. However, achieving the flame front resolution corresponding to ambient pressure is either extremely demanding or impossible for high pressure flames. As an example, the 25 bar case has a flame thickness which is 25 −0.5 ≈ 0.2 times smaller than that of the 1 bar case, which necessitates roughly a five times finer mesh in each coordinate direction in the 25 bar case than in the 1 bar case. Consequently, FSD models are expected to play a much more important role under high pressure conditions, and in this sense possible modelling errors are amplified in high pressure premixed turbulent combustion; (2) Surface-filtered, density-weighted displacement speed should be corrected by the resolved mean curvature contribution unless either the Markstein diffusivity or the mean flame curvature vanishes (for a detailed discussion of this the reader is referred to Chakraborty et al. 2018). Due to the asymmetry of the curvature PDF, the mean value of the subgrid scale contribution of the surface-filtered density-weighted tangential component of displacement speed does not vanish. However, although it can assume locally large values, the mean (true) values remain of the order of 5% of the laminar flame speed for the present configuration independent of pressure. This means that the Darrieus-Landau instability effects are implicitly present in the modelling of surface filtered mean curvature; (3) Assuming a power-law based modelling, it has been shown that, for a given set of values of normalized turbulent velocity fluctuations and integral length scale in cases A, B, C and D, the fractal dimension increases by a small extent, while the inner cut-off scale significantly increases with increasing pressure. This effect possibly can be attributed to an evolving Darrieus-Landau instability which is clearly present in case C and case D but absent in case A (see Klein et al. 2018 for more details), but a more detailed analysis is required in order to confidently validate this conjecture. Further, it will be necessary to consider a much wider parameter range, before these influences can be parameterised for modelling purposes. Development of model modifications in order to account for flame instabilities and for contributions of the subgrid scale surface-filtered curvature will be addressed in the future. Moreover, FSD transport based modelling for high pressure flames will be considered in the future analysis.