Large Eddy Simulations of a turbulent premixed swirl ﬂame using an algebraic scalar dissipation rate closure

Large Eddy Simulations (LES) of a swirl-stabilised turbulent premixed ﬂame in the well-known TECFLAM burner conﬁguration have been carried out by solving transport equations of Favre-ﬁltered reaction progress variable and mixture fraction. A recently proposed closure for the Scalar Dissipation Rate (SDR) is used for the modelling of the ﬁltered reaction rate of reaction progress variable, whereas the Favre ﬁl-tered mixture fraction is used to account for mixture stratiﬁcation due to entrainment. The computational results are utilised to analyse the nature of stratiﬁcation at representative locations in the swirl ﬂame to gain physical insight into the ﬂame structure. Additionally, two algebraic Flame Surface Density (FSD) closures, which were found to perform well in a previous analysis (Ma et al., 2013), are used for the modelling of the ﬁltered reaction rate of reaction progress variable. The predictions of SDR closure are compared to the corresponding results obtained from algebraic FSD closures. The predictions of SDR based simulations show reasonably good agreement with experimental ﬁndings; the level of accuracy is at least comparable to that achieved with algebraic FSD models and to the results reported in the literature. (cid:2) 2015 The Authors. Published by Elsevier Inc. on behalf of The Combustion Institute. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
The Scalar Dissipation Rate (SDR) of mixture fraction n plays a pivotal role in the fundamental understanding and modelling of the rate of micro-mixing in non-premixed combustion; interested readers are referred to Bilger [1] for a detailed discussion in this regard.In premixed combustion, the SDR of reaction progress variable c is not only useful for the closure of the variance of the progress variable, but also closely related to the mean reaction rate in the context of Reynolds Averaged Navier Stokes (RANS) simulations [2][3][4][5].It has been demonstrated recently [6][7][8] that the Favre-filtered SDR e N c of the reaction progress variable c also remains proportional to the filtered reaction rate _ w according to the following relation when the filter width D remains much greater than the thermal flame thickness d th (i.e.D > d th ): In Eq. ( 1), q is the filtered gas density, N c ¼ Drc Á rc is the instantaneous SDR with D being the progress variable diffusivity, and Q and e Q ¼ qQ= q are the Large Eddy Simulation (LES) filtered and Favre-filtered values of a general quantity Q respectively.The parameter c m is given by [2]: where f ðcÞ is the burning mode probability density function (pdf); it has been demonstrated by Bray [2] that the numerical value of c m remains insensitive to any presumed continuous smooth function f ðcÞ, and for typical hydrocarbon flames c m ranges from 0.7 to 0.9 [2].Eq. ( 1) was originally proposed by Bray [2] in the context of RANS and this expression can be analytically derived from the transport equation of reaction progress variable variance for high values of Damköhler number (i.e.fast chemistry) where the probability density function of c can be approximated by a bi-modal distribution with impulses at c ¼ 0 and c ¼ 1:0.It was shown in Refs.[3,4] that Eq. (1) holds even for low Damköhler number (i.e.Da < 1) combustion based on an order of magnitude analysis of the terms in the progress variable variance transport equation, which was subsequently confirmed based on a priori DNS analysis [3].Dunstan et al. [6] discussed about the possibility of extending an algebraic RANS-SDR closure for the purpose of LES based on a priori DNS analysis of a single V-flame simulation with unity Lewis number Le (i.e.ratio of thermal diffusivity to mass diffusivity), which has subsequently been extended by Gao et al. [7] for the algebraic modelling of e N c for flames with a range of different values of global Lewis number Le, heat release parameter s ¼ ðT ad À T 0 Þ=T 0 (where T o and T ad are the unburned gas and adiabatic flame temperatures respectively), and turbulent Reynolds number Re t based on a priori DNS analysis.It has been demonstrated in previous a priori DNS analyses [6][7][8] that Eq. (1) can also be used for the filtered reaction rate _ w closure for most practical LES grid sizes (i.e.D > d th ), as combustion predominantly takes place at the sub-grid level.Furthermore, it has been found that the predictive capability of Eq. (1) [6][7][8].This suggests that Eq. ( 1) can be used for reaction rate closure if e N c is appropriately modelled.The SDR based reaction rate closure and the algebraic SDR model proposed by Gao et al. [7] were subsequently a posteriori assessed by Ma et al. [8] using LES of a rectangular dump combustor configuration (i.e.ORACLES burner) [9] as well as for a flame stabilised on a triangular bluff body within a rectangular channel (i.e.VOLVO rig) [10].Ma et al. [8] found that the predictions of the SDR based closure were either comparable or better than the two most promising Flame Surface Density (FSD) models [8,[11][12][13].However, the SDR based closure is yet to be analysed for LES of turbulent swirl flames involving stratified mixtures.The objective of the current analysis is to assess the quality of this newly developed SDR closure in a more complex burner configuration involving stratified mixtures, making it reasonably similar to industrial applications.The TECFLAM swirl burner of the Institute for Energy-and Power-plant Technology (EKT) at Darmstadt University offers an ideal test case for this purpose.A comprehensive database of experimental measurements already exists [14] and various computational analyses [15][16][17][18][19][20] have already been performed on this flame.
A number of experiments and simulations on the TECFLAM swirl burner configuration have been conducted at EKT [14,21]: Schneider et al. [14] recorded mean values and fluctuations of velocities as well as temporal auto-correlations, spatial cross-correlations and spectra for the reacting and non-reacting cases using Laser-Doppler Velocimetry (LDV).This velocity dataset was later on appended with experimental data on mean and fluctuating temperature profiles, obtained from Rayleigh-scattering, and major species fractions measured using Raman-scattering by Gregor et al. [21].Wegner et al. [15] performed simulations to assess unsteady Reynolds-Averaged Navier-Stokes (URANS) simulations for the isothermal flow, comparing results to the experimental measurements and LES calculations.A Direct Numerical Simulation (DNS) of the isothermal flow with a reduced Reynolds number was performed by Freitag and Klein [16] who reported excellent agreement with experimentally obtained mean velocities and two-point correlations.Freitag and Klein [16] also analysed the transport of a passive scalar using DNS and satisfactory agreement was obtained.Klein [17] subsequently focussed on the effects of turbulent inflow boundary conditions and the inclusion of the burner geometry in the computational domain.Freitag [18] was first to analyse the reactive flow in the TECFLAM burner by LES, using a level-set approach.The velocity statistics obtained by Freitag [18] were in good agreement with experimental data.More recently, Kuenne et al. [19] applied an artificially thickened flame (ATF) approach with the wrinkling factor proposed by Colin et al. [22], and the power-law formulation by Charlette et al. [23], coupled with tabulated chemistry.Kuenne et al. [19] achieved good agreement with experimental data for both velocity and temperature statistics.An Eulerian stochastic field method was used by Jones et al. [20] for the LES of TECFLAM burner, which also yielded satisfactory agreement with the velocity and temperature statistics obtained from experiments.
The SDR based reaction rate closure depends on flamelet assumption similar to the level-set modelling used by Freitag [18] and the ATF approach adopted by Kuenne et al. [19].By contrast, the Eulerian stochastic field method used by Jones et al. [20] could be used for both flamelet and non-flamelet combustion but at the expense of higher computational cost than the SDR, level-set and ATF approaches.The flame brush thickness control can be an issue for algebraic FSD and SDR closures but this problem is inherently avoided in the ATF and level-set methods.The inter-relation between the temperature field and the flame surface is not straightforward and often empirically formulated in the level-set approach, whereas the temperature and species fields are inherently interlinked by combustion thermo-chemistry in the context of FSD, SDR and ATF methodologies.The wrinkling factor models [22,23] used in the ATF methodology by Kuenne et al. [19] were proposed for unity Lewis number combustion but they do not adequately perform for non-unity Lewis number flames [24].It has been demonstrated based on a priori DNS analysis that Eq. ( 1) remains valid for a range of different Lewis numbers Le, whereas the reaction rate closure by FSD needs a correction factor, which is dependent on the global Lewis number [3].This makes the SDR based closure a promising computationally inexpensive methodology for turbulent premixed combustion modelling.
In the current analysis a SDR based reaction rate closure (i.e.Eq. ( 1)) has been used for stratified combustion in the TECFLAM burner configuration where the algebraic model for e N c proposed by Gao et al. [7] is used.In a recent a posteriori LES assessment [13] of the algebraic FSD models for the ORACLES burner [9] and the VOLVO rig [10], the models proposed by Fureby [11] and Muppala et al. [12] performed better than the alternative algebraic FSD models.The performances of the closures by Fureby [11] and Muppala et al. [12] have therefore been compared to the predictions obtained from the SDR closure too.The results from the current LES in turn have been compared to the earlier predictions from ATF [19] and Eulerian stochastic [20] closures.In this respect, the main objectives of the current analysis are: 1. To assess the SDR based reaction rate closure for a swirl flow stabilised burner involving stratified mixtures.2. To compare the predictions of the SDR based closure with the predictions of established FSD models and the results available in the literature [19,20].
The remainder of the paper is organised as follows: the model formulation and numerical implementation are presented in the next two sections.These will be followed by the presentation of results and their discussion.The main findings are then summarised and conclusions are drawn in the final section.

Model formulation
The reaction progress variable c can be defined in terms of a suitable reactant mass fraction (e.g.fuel mass fraction Y F ) in the following manner: In Eq. ( 3), Y Fu and Y Fb are the mass fractions of fuel in completely unburned and burned gases respectively; these quantities are functions of mixture fraction n, which is defined as: Here, Y F1 ¼ 1:0 denotes the mass fraction in the pure fuel stream and Y O1 ¼ 0:233 the oxidiser mass fraction in pure air, with Y O being the oxidiser mass fraction.In Eq. ( 4), s is the ratio of mass of oxidiser to the mass of fuel under stoichiometric conditions (e.g.s ¼ 4 is for methane-air mixture).The quantities Y Fu and Y Fb are given by: Interested readers are referred to Bray et al. [25] for further information in this regard.The filtered transport equations for the Favre filtered reaction progress variable c and mixture fraction ñ for globally fuel-lean turbulent stratified mixtures read: The first and second terms on the left hand side in Eqs.(5i) and (5ii) indicate the transient and resolved advection effects.The first term on the right hand side indicates the effects of molecular diffusion.The second term on the right hand side of Eq. (5i) is the filtered reaction rate _ w and the third term involves the sub-grid flux ½qu j c À qũ j c, both of which require modelling.The last term on the right hand side of Eq. (5i) arises due to cross-scalar dissipation rate of progress variable c and mixture fraction n due to mixture stratification (i.e.variation of Y Fu ¼ Y F1 n due to the change in equivalence ratio).It is worth noting the cross-scalar dissipation term (i.e. term IV in Eq. (5i)) remains negligible (at least one order of magnitude smaller in this case) in comparison to filtered reaction rate term (i.e. term II in Eq. (5i)) in predominantly premixed stratified flames [26].Furthermore, both rc and r ñ vanish as ñ approaches zero (i.e.pure air) so the term IV is exactly equal to zero for ñ ¼ 0.
It is important to note that it is necessary to solve a transport equation for ñ in order to account for the stratification of the fuel-air mixture by the co-flowing air in this configuration.
Ignoring spatial variations of mixture composition leads to spurious flame propagation into the regions characterised by equivalence ratio values below the flammability limits (e.g.combustion of the co-flowing air).This is in contrast to a previous analysis by Ma et al. [8] where the SDR closure was applied to a rectangular dump combustor (i.e.ORACLES burner) [9] and a bluff body stabilised flame in a rectangular channel (i.e.VOLVO rig) [10].In these cases no stratification occurs inside the enclosed burner configurations, and hence solving only the transport equation of c was sufficient.
Eqs. (5i) and (5ii) with closure models applied take the following form: The symbols l and r (l t and r t ) are the molecular (Eddy) viscosity and molecular (turbulent) Schmidt number respectively, and ð _ wÞ M stands for the modelled filtered reaction rate term.The dynamic viscosity l ¼ qDr in Eqs.(6i) and (6ii) varies with temperature following Sutherland's law [27], while the laminar and turbulent Schmidt numbers are kept at a constant value of 0.7.In Eq. (6i), the gradient contribution of ½qu j c À qũ j c is modelled using 6i) corresponds to the term arising from counter-gradient behaviour of ½qu j c À qũ j c.Satisfactory results have been obtained in the past [8,13] using a constant turbulent Schmidt number in the context of both FSD and SDR based closures and the same approach has been followed here.The magnitude of the unclosed turbulent transport term remains small in comparison to the magnitude of the chemical source term [13,28] but the nature of turbulent transport (i.e.gradient as opposed to counter-gradient transport) affects the flame brush thickness and flame wrinkling [13].It has been found that a change in turbulent Schmidt number from 0.7 to 0.4 did not affect the predictions of the simulations in this configuration (not shown here).A-priori assessment of the performances of the various sub-grid scalar flux models with static coefficients can be found in Gao et al. [29], which reveals that the benefit of dynamic evaluation of turbulent diffusivity is likely to be marginal for some sub-grid scalar flux models.The effects of dynamic evaluation of turbulent diffusivity on predictive capabilities of FSD and SDR based closures are beyond of the scope of this analysis, and can be a subject of a future analysis.
The cross dissipation term of Eq. ( 5i) is modelled in the following manner in [26,30] and the same approach has been adopted here: As the term IV in Eq. (5i) remains negligible in comparison to the chemical source term, the model in Eq. ( 7) does not have a major influence on the predictions of the LES simulations in the configuration considered here.
In the context of SDR based closure, ð _ wÞ M is expressed using Eq.
(1), with e N c being modelled as [7,8]: The equation includes the bridging function f b ¼ exp½À0:7ðD=d th Þ 1:7 that ensures that the LES becomes a DNS for very fine grids, as well as C Ã 3 ; C Ã 4 and b c as model parameters with the following values resulting from a priori DNS analysis [7,8]: where s 0 ¼ ðT b À T 0 Þ=T 0 is the modified heat release parameter.

Eqs. (8i) and (8ii) depend on the local sub-grid Damköhler number Da
p is the sub-grid level velocity fluctuation.It is worth noting that Eq. (8i) was obtained by extending a RANS-SDR model [4], which can be derived by assuming an equilibrium of the leading order contributors to the SDR transport for high Damköhler number (i.e.Da ) 1) combustion [4,31].A recent analysis by Gao et al. [32] demonstrated based on scaling arguments that the order of magnitudes of the leading order terms in the SDR transport equation can be used to justify the validity of Eq. (8i) for low Damköhler number (i.e.Da < 1) flamelet combustion, which was subsequently substantiated by a priori DNS analyses [7,32].This indicates that Eq. (8i) can be applied for both high and low Damköhler number combustion.Interested readers are referred to Refs.[4,[6][7][8]31,32] for the theoretical justification and a priori DNS assessment of the algebraic SDR closure given by Eq. (8i).A-priori DNS assessment of Eq. (8i) was carried out for flames for moderate values of turbulent Reynolds number [7,8,32] but this model has been a posteriori validated using LES simulations of ORACLES and VOLVO rigs in [8], and satisfactory agreement with experimental data has been obtained.Furthermore, the RANS model for SDR, which has been extended for LES to obtain Eq. (8i), was a posteriori assessed based on actual RANS simulations for several practical burners and laboratory-scale configurations in the past and satisfactory agreement with experimental observations was obtained [33][34][35].
The applicability of Eq. (8i) in the modelling of turbulent stratified flames is yet to be established.For stratified mixture combustion the effective burning velocity S 0 L and effective burned gas temperature T b due to stratification are evaluated as: In this equation, S L ðnÞ and T ad ðnÞ are the polynomial fits of the variation of laminar burning velocity and adiabatic flame temperature with mixture fraction n according to Vagelopoulos and Egolfopoulos [36].In Eq. (8iii), PðnÞ corresponds to the presumed sub-grid pdf of mixture fraction, which is approximated by the top-hat distribution, as suggested by Floyd et al. [37] for LES: The parameter P 0 is the magnitude determined by the normalisation condition R 1 0 PðnÞdn ¼ 1:0, whereas n a and n b denote the lower and upper limits of mixture fraction.The sub-grid variance of mixture fraction is estimated as: The parameter C n ¼ 0:2 is set according to Branley and Jones [38].
The thermal flame thickness d th is estimated as [39]: 7 where d z is the Zel'dovich flame thickness, which is given by: d z ¼ a T0 =S 0 L where a T0 is the thermal diffusivity in the unburned gas.In Eq. (8i) K Ã c is a thermo-chemical parameter, which for a premixed flame can be expressed as [6][7][8]31,40,41]: Here K Ã c is taken to be 0:85s 0 based on the previous findings based on a priori DNS analysis [40] for methane-air flames with an equivalence ratio close to / ¼ 0:83, which is very close to the operating condition of the TECFLAM burner.For methane-air flames with an equivalence ratio close to / ¼ 0:83, the value of c m in Eqs. ( 1) and (8ii) is found to be 0.835 based on laminar flame calculations, and this value is used for both reaction rate and SDR closure.
Ma et al. [8] used the sub-grid scalar flux model proposed by Richard et al. [42] in conjunction with the SDR based reaction rate closure, and the same approach has been adopted here, which leads to the following model expression for the unresolved flux ½qu j c À qũ j c: The parameter M i ¼ Àð@c=@x i Þ=jrcj ¼ Àð@ c=@x i Þ=jr cj is the resolved flame normal vector, and the second term on right hand side of Eq. ( 9) accounts for counter-gradient transport.The sub-grid scalar flux ½qu 1 c À qũ 1 c assumes a value equal to Àq 0 S L ð c À cÞM 1 based on mass conservation for a steady laminar 1-D flame (i.e.qu 1 ¼ q 0 S L ) moving in the negative x 1 -direction [42,43].Thus, the term Àq 0 S L ð c À cÞM 1 deterministically suggests a non-gradient transport for a planar flame propagating in the negative x 1 -direction, which is utilised by Richard et al. [42] to model the counter-gradient contribution to ½qu j c À qũ j c in Eq. (9)   as Àq 0 S 0 L ð c À cÞM j [42].The reaction source term and the molecular diffusion are commonly combined and modelled using the generalised FSD in the following manner [43,44]: The closure of Eq. ( 10) is valid in the corrugated flamelets regime of premixed combustion [45], where the flame thickness is smaller than the Kolmogorov length scale.It leads to the following form of the modelled transport equation for c for turbulent stratified combustion: The quantity N ¼ R gen =jr cj denotes the wrinkling factor and inter- ested readers can refer to previous work [13,24,46,47] for a detailed assessment of many available algebraic wrinkling factor models.
According to Fureby [11] N can be modelled as: with It is worth noting that the original version of the Fureby model [11] does not have the addition of unity in Eq. (12i), which was introduced by Ma et al. [13] to ensure that N reduces to one for regions of zero turbulence.The model given by Eq. ( 12) was used in conjunction with Eq. ( 9) in the a posteriori assessment by Ma et al. [13] and the same approach has been adopted here.
According to Weller et al. [48] the sub-grid scalar flux can alternatively be modelled in the following manner: where N M is an appropriate wrinkling factor.The first term on right hand side in Eq. (13i) accounts for gradient transport, whereas the second term takes care of counter-gradient transport.It was demonstrated by Lecocq et al. [49] and Ma et al. [13] that the turbulent transport term can be written in the following manner using Eq.(13i) under some assumptions: The combination of _ and Eq.(13ii) yields: This suggests that the combination of the source term ð _ wÞ M and counter gradient transport can also be modelled in the following manner [13,49], using an appropriate wrinkling factor and the gradient of the density weighted progress variable: In the context of Eq. ( 13) the Favre filtered progress variable transport equation takes the following form: The wrinkling factor N M according to Muppala et al. [12] is given by: The model and the sub-grid Reynolds number Re D ¼ u 0 D D=m include m, p and p o , which are the kinematic viscosity, filtered pressure and reference pressure respectively.
The unclosed terms of the transport equation of c is expressed in the form shown by Eq. ( 14) in the original formulation by Muppala et al. [12].Thus, it is only fair to consider Eq. ( 14), as a whole, for the purpose of assessment of the methodology proposed by Muppala et al. [12] because the closures of chemical source and turbulent transport terms are known to interact with each other in LES of turbulent premixed flames [13].Furthermore, the analysis by Ma et al. [13] demonstrated that the original formulation by Muppala et al. [12] predicts the experimental observation in a better manner than several other existing algebraic FSD closures for the numerical schemes used in the current analysis.A detailed analysis by Ma et al. [13] also revealed that the chemical source term closure using the wrinkling factor by Fureby [11] yields the best prediction when the sub-grid scalar flux model proposed by Richard et al. [42] (i.e.Eq. ( 9)) is used.The prediction of the LES simulations for the combination of Fureby's wrinkling factor model (i.e.Eq. ( 12)) and the sub-grid flux model by Richard et al. [42] has been found to be comparable to that of Muppala et al. [12] and better than several other alternative algebraic FSD closures in the analysis by Ma et al. [13].In this paper, the SDR based closure is compared to two best performing algebraic FSD closures [13] and thus the combinations of chemical source term and sub-grid flux, which yielded optimum performance for the models proposed by Fureby [11] and Muppala et al. [12], are retained in this analysis.A-posteriori assessment by Ma et al. [8] demonstrated that satisfactory results have been obtained when SDR based reaction rate closure (i.e.Eq. ( 1)) is used in conjunction with the sub-grid scalar flux model (Eq.( 9)) proposed by Richard et al. [42].Thus a combination of Eqs. ( 1) and ( 9) has been used here for the implementation of SDR based closure.

Flow configuration
The unconfined TECFLAM swirl burner has been studied experimentally at the Institute for Energy-and Power-plant Technology at Darmstadt University of Technology [14,21].The operating conditions of this burner ensure that a stratified mixture resulting from entrainment burns in a predominantly fuel-lean premixed mode [14][15][16][17][18][19][20][21].The swirling motion of the fluid flow is generated using a moveable block geometry upstream of the burner outlet that allows for tuning of the swirl intensity.Swirl intensity can be expressed by means of a swirl number following the original definition by Gupta et al. [50], representing the ratio of tangential and axial momentum transported by the axial flow.In the experiments the swirl number was set to 0.7, providing sufficiently high swirl intensity for the formation of a central recirculation zone due to vortex breakdown caused by a pressure gradient outweighing the overall negative streamwise pressure gradient along the jet.The region of low pressure close to the burner exit originates due to centrifugal forces, i.e. conversion of tangential momentum to radial momentum.In the reactive case, this enforces backflow of hot gases towards the nozzle exit and is responsible for flame stabilisation.Periodic distortions of this recirculation zone have been observed experimentally for the non-reactive case, leading to vortex-like structures exhibiting precessions around the centreline of the jet, also called ''precessing vortex cores'' (PVC).No precessing motion can be found for the reactive case, suggesting damping due to higher kinematic viscosities and dilatation rate effects at elevated temperatures.
The flame burns a natural gas/air mixture with an equivalence ratio of / ¼ 0:83, leaving the burner geometry through an annulus with an outer diameter of 60 mm, which surrounds a central bluff body with a diameter of 30 mm.In the reactive experiments, the bluff body was water-cooled to maintain a temperature of 353 K. Homogeneous mixing is ensured by the injection of natural gas to the air flow far upstream of the nozzle exit and has been confirmed experimentally by Schneider et al. [14].For the non-reactive experiments, the natural gas fraction was substituted by air in a way that provides the same flux of momentum at the burner outlet.
The experimental investigations at EKT addressed test cases with flames of 30, 90 and 150 kW of thermal power and the corresponding non-reactive experiments.The most important flow conditions for the 30 kW case, which is studied numerically in this work, are summarised in Table 1, where Reynolds number is calculated based on the flow in the nozzle exit plane and the diameter of the central bluff body.

Numerical implementation
A sketch of the computational domain is shown in Fig. 1.An in-house 3D low-Mach number finite volume code PsiPhi  [13,23,[51][52][53] is used, which discretises the computational domain in the form of cubic cells with uniform grid spacing.The time advancement has been carried out using a 3rd order Runge-Kutta scheme in combination with a predictor-corrector procedure.The time step width was controlled by the Courant-Friedric hs-Lewy (CFL) criterion, and a CFL number of 0.5 was used for the simulations in this analysis.The momentum equations were discretised using a second order central differencing scheme to achieve good accuracy and low numerical dissipation.A TVD scheme using the CHARM limiter function was used for the discretisation of the mass conservation equation and transport equations of c and ñ in order to avoid numerical oscillations.It was found to be necessary to apply this scheme also to the momentum equations in regions with sharp velocity gradients close to the burner exit plane, i.e. x < 60 mm.For the closure of sub-grid stresses the static Sigma model with a constant of C = 1.5 was used, as originally proposed by Nicoud et al. [54]: In this model, the magnitudes of the principal values of the resolved velocity gradient tensor @ũ i =@x j are ordered in a way that r 1 P r 2 P r 3 P 0. It is worth noting that including a part of the annular nozzle in the computational domain is necessary to reproduce flow instabilities, i.e. the precessing vortex core in the investigated case [15].In addition, Klein [17] demonstrated that underpredictions of turbulent velocity statistics and spreading of the jet were obtained when the computational inlet was located at the nozzle exit, which was also found in a trial simulation (results not shown here).However, as the PsiPhi code is based on rectangular domains, the inclusion of large parts of the burner geometry would lead to an unjustifiable expense for the computation of flow quantities in cells surrounding the configuration that are of no interest for the investigation of the flame.A short part of the nozzle with an axial extension of 10 mm, reproduced in the domain by means of immersed boundaries, proved sufficient to avoid these computational expenses but still capture important flow physics.As in this setup the computational inlet is located only slightly upstream of where experimental inflow data had been taken, experimental nozzle exit profiles were set as boundary conditions.Artificial turbulence was created by a standard inflow generator proposed by Kempf et al. [55] providing a turbulent inflow satisfying a prescribed Reynolds-stress tensor and an inflow integral length scale.Based on a detailed sensitivity analysis (not shown here) it was found that Reynolds-stresses based on experimental root-mean-square (rms) velocity fluctuations and an inflow integral length scale of 5 mm according to Schneider et al. [14] yielded the closest agreement with experimental data.A clipped Neumann condition suppressing negative axial velocities at the outlet of the domain was used as boundary conditions for momentum.Lateral boundaries were treated with a Dirichlet condition, set to the co-flow velocity.Zero gradient boundary conditions were applied to the outflow of the domain as well as the lateral boundaries for c, Favre-filtered temperature, and ñ, while a major impact isothermal reactive    of the lateral boundaries was not expected due to their distance from the area of interest.The Favre filtered mixture fraction ñ at the inlet was given a value corresponding to an equivalence ratio of 0.83, i.e. ñ ¼ 0:0454.The Favre filtered reaction progress variable c at the inlet was set to zero corresponding to the unburned state.Accordingly, the Favre-filtered temperature is taken to be the unburned gas temperature (=300 K) at the inlet.
As the flame could only be sustained after a stable recirculation zone was established, a field of fully burned products (i.e.c = 1.0) was initialised inside the annular flow field after approximately 2000 time steps.The simulations were conducted on both ''coarse'' and ''fine'' grids, with cubic cells of 2 mm and 1 mm resolution respectively, corresponding to approximately 7.5 M and 60 M cells.A combination of 450 ms for initialization and 650 ms for sampling was found to be sufficient to obtain converged statistics.The reactive simulation required 24 and 96 core hours on the coarse and fine grid respectively, using up to 192 cores in parallel.the fine grid.Vorticity has large magnitudes especially in the outer shear layer.In the isothermal case, a redistribution and transfer to smaller Eddies is observed, leading to smaller vortices that are more evenly distributed over the radius of the flame.Radial profiles of axial, azimuthal and radial velocity mean and rms values for isothermal flow conditions obtained on both grids at various axial positions are presented in Figs.3-5.Spatial coincidence of peaks in axial and azimuthal mean velocities in all measurement planes, on the fine as well as on the coarse grid, proves that the spreading of the annular jet after exiting the burner is well reproduced, while the level of agreement with experimental data deteriorates with increasing axial distance.As the internal recirculation zone is the flow feature, which plays a key role in flame stabilisation, it is crucial to reproduce the magnitude and extension of the backflow in the simulations.It can be seen from Figs. 3-5 that the strength of recirculation is mostly well captured, with only a small underestimation of the negative axial velocities further downstream in combination with a slightly short recirculation zone.

Isothermal non-reacting flow simulation
Simulations with an increased domain size have been performed but did not provide an improved prediction of the length of the recirculation zone.The agreement of radial mean velocities obtained from the simulation with the experiment is somewhat poorer than what was found for the other components, which is consistent with the results of LES simulations by Jones et al. [20] and DNS data by Freitag and Klein [16].This is perhaps not surprising given the small magnitude of the radial velocity, making small errors appear large in both simulation and experiments.The trends of turbulent fluctuations of the velocity components in all measurement planes are reasonably well captured on the fine grid except for an underestimation on the centreline of the nozzle, in particular near the burner exit plane.However, the deviations in radial and azimuthal velocity fluctuations are more noticeable than those in the axial direction.Unsurprisingly, the predictions for rms velocities obtained from the coarse grid simulation are less accurate.Compared to the fine grid results, resolved fluctuations are weaker in the coarse grid simulations, but the mean velocities are hardly affected and can be considered as converged.It is worth noting that the experimental results for the isothermal case at x = 10 mm on the centreline show significantly stronger fluctuations in radial and azimuthal direction than in the axial direction, which is in contrast to almost isotropic conditions in the reacting case.This could be due to the fact that the precessing motion, which is only present in the isothermal case, might not have been captured with sufficient accuracy to reproduce the experimentally observed velocity fluctuations in both radial and tangential directions.However, discrepancies between simulation results and  experimental data might also have been caused by the imposed turbulent inflow boundary condition and the modelling of sub-grid stresses.From the overall reasonably good agreement between the simulation and experimental data in Figs.3-5, it can be concluded that the numerical setup (e.g.inflow turbulence generation, filter width, boundary conditions and sub-grid turbulence modelling) is suitable for the simulation of the TECFLAM configuration, which is a necessary condition for a meaningful analysis of the performance of combustion models in turbulent reacting flow LES studies of this test case.

Reactive flow simulation
Figure 2 shows how the reactive flow field is qualitatively different from the isothermal one, even though the maximum values of vorticity magnitude are of the same order of magnitude.The volumetric dilatation rate r Á ũ acts to dissipate vorticity within the flame, which is aided by the augmented viscous dissipation of vorticity in the reactive case due to the higher temperature values.In fact, the dissipation of enstrophy by dilatation and viscous dissipation are found to be of same order of magnitude within the flame brush in this configuration.Instantaneous field snapshots of axial velocity, mixture fraction, temperature and calculated reaction source term are shown in Fig. 6.From the distribution of the reaction source term it is evident that combustion in the proximity to the nozzle takes place mainly on the inner shear layer of the swirl flow, i.e. in regions of almost constant value of ñ.Backflow of burnt products with approximately adiabatic flame temperature inside the jet stabilises the flame.Further downstream, however, burning occurs in more stratified mixtures and the flame might eventually move to the outer mixing layer.It is also visible from Fig. 6 that wrinkling of the flame front increases in downstream direction, while nearly no wrinkling can be observed adjacent to the nozzle.
Figures 7-9 show axial, azimuthal and radial velocity mean and rms values obtained from reactive simulations on the fine grid using the SDR model along with the measured data and the results obtained from FSD based closures proposed by Muppala et al. [12] and Fureby [11].Results produced by all three closures look very similar and it is not possible to identify a distinct trend in the differences between the velocity profiles.Overall, the agreement between simulation results and the experiment is reasonably good while, generally, recirculation is underpredicted and radial velocities are overestimated in the regions close to the nozzle exit near the centreline of the domain.A similar behaviour was observed in the LES by Jones et al. [20] and by Kuenne et al. [19] when the wrinkling factor model by Colin et al. [22] was used (Results obtained by Kuenne et al. [19] using the power-law formulation by Charlette et al. [23] do not show this behaviour, but provide less accurate estimations of azimuthal velocities).These discrepancies might have arisen from a larger extent of thermal expansion compared to the experiment.This may be due to heat losses to the water-cooled bluff body, which were not included in the simulation, so that gas temperatures near the nozzle in the centre of the jet are close to the adiabatic flame temperature of 2040 K (see Fig. 10), in contrast to roughly 1730 K in the experiment.
The wider spreading angle of the annular jet compared to the experiment is also very likely linked to the absence of heat loss.Overall, mean velocity predictions provide a comparable level of accuracy to those shown in [19,20].
Regarding velocity fluctuation levels, both qualitative and quantitative agreements were found to be reasonably good throughout the domain, with the predictions improving in the downstream direction.While it is possible that high values of viscosity due to elevated temperatures especially on the centreline close to the burner outlet are responsible for underpredictions of velocity fluctuations, it is also possible that the TVD convection scheme used in this area of high velocity gradients leads to higher levels of damping due to its dissipative nature.However, the vorticity magnitude does not show a major change in behaviour due to the change in discretisation schemes from CDS to TVD at x = 60 mm, which does not indicate a strong influence of the TVD scheme on the calculated level of turbulence.Another possible cause for the underestimation of velocity fluctuations close to the burner might also have arisen due to the neglected heat losses to the bluff body: Heat losses in the experiment probably lead to a flame lift-off with an oscillating lift-off distance inducing further fluctuations that cannot be reproduced by the simulation.
Mean and rms profiles of temperature are shown in Fig. 10.While in all other aspects the performances of FSD and SDR closures are similar, the proposed SDR model visibly captures both magnitude and position of the peak value of temperature fluctuations in measurement planes up to axial distance x = 30 mm better than the FSD models.A similar effect was found in the coarse grid simulations, as shown in Appendix A. For all simulations, up to about the axial distance x = 30 mm, the width of the high-temperature zone is significantly too large, which is consistent with the overestimation of radial velocities discussed before.This overpredicted radial flow might have originated due to the neglected heat losses again, as higher temperatures at the centreline lead to a stronger expansion of the jet, which in turn is believed to lead to higher values of ñ at large radii than those found from experiments as apparent from Fig. 11.
A radial offset between the simulated and measured temperatures can similarly be found but to a lesser extent in the results published by Kuenne et al. [19] and Jones et al. [20].However, the behaviour of the flame is as expected: combustion occurs at the inner shear layer in the regions of constant mixture fraction close to the nozzle, as can be concluded from a comparison of Figs. 10 and 11.As mentioned before, at an axial distance from the burner (x = 60 mm), high temperature gradients occur in stratified mixtures as well.
A second maximum of mean temperatures outside the outer shear layer is visible from Fig. 10 (20 mm, 30 mm) for both calculated and measured values.The temperature peak at the outer shear layer predicted by the SDR simulation shows better agreement with experimental results than previous simulations.The second temperature peak has been attributed to the formation of an external recirculation zone by Gregor et al. [21] and Jones et al. [20].From the present simulations it was found that due to the radial flapping motion of the main reactant stream, burned products are occasionally trapped on the co-flow side of the main stream where they will rest due to the low co-flow velocity.At the radial location of the main jet, however, these products are blown away by the fast main stream, leading to a smaller mean progress variable and with that a smaller mean temperature in this region than the corresponding values on the co-flow side.This behaviour is illustrated by the sequence of instantaneous temperature fields shown in Fig. 12.The extension of the internal recirculation zone and the location of approximately conical inner shear layer are visible from the instantaneous temperature fields shown in Fig. 12.The main flow separates this region from the occasionally trapped pockets of hot gas in the co-flow.
It is worth noting that results produced by the different models shown in Figs.7-11 are very similar and deviate from experimental results in the same manner.The numerical setup of the simulation including the discretisation schemes, boundary conditions, turbulence modelling etc. and the deficiencies such as ignoring the heat transfer to the bluff body leads to the same deviations from experimental data and in fact there is very little difference in the output of the models with regard to the velocity predictions on the fine grid.In terms of temperature predictions, however, both mean and especially fluctuating data from the SDR simulations are preferable to those of the two FSD models, as can be seen from Fig. 10.This is particularly true in the regions with small values of mixture fraction, i.e. in the stratified regions.This can also be 30mm seen to a much greater extent in the results from the coarse grid simulations, which are discussed in detail in Appendix A.
Figure 13 shows profiles of the resolved and the modelled part of scalar dissipation rate at various axial positions.The resolved part of the SDR (i.e. e Drc Á rc) reaches its maximum where gradi- ents of Favre-filtered values of reaction progress variable and temperature are high, i.e. in the main reaction zone at the inner shear layer.Throughout the domain the modelled part shows the same qualitative trend but is considerably larger.It is obvious that the ratio of resolved to modelled contribution will shift towards a larger resolved part with increased grid resolution, as expected.
To enable a better understanding of the flame structure, investigations into the nature of stratification have been carried out at two representative locations.Stratified flames can be back-and front-supported, where a back-supported flame travels towards leaner mixtures and a front-supported flame propagates towards richer mixtures.Samples of the relative directions of the gradients of c and ñ have been taken at a position close to the nozzle exit at the inner shear layer of the jet (x = 30 mm, r = 25 mm), i.e. inside the internal recirculation zone, results of which are shown in Fig. 14.The mean and fluctuating temperature distributions in Fig. 10 indicate that the gradients of c at this position are expected to be predominantly negative, suggesting a propagation of the flame towards larger radii.Positive gradients might appear eventually due to fluctuations in the c-profile and wrinkling of the flame front.However, fluctuations in this region are moderate and, as discussed earlier based on Fig. 6, wrinkling is weak in the inner shear layer in vicinity of the nozzle exit.Positive directions of the gradient of c (\c ¼ p=2) indicate propagation of the flame towards the centreline of the jet, which can occur locally only due to wrinkling.Therefore only few samples representing positive c-gradient directions were found here, while negative directions (\c ¼ Àp=2) are much more likely, both for the case of back-and front-support, the former being present in conjunction with negative directions of ñ-gradients and vice versa.Judging from the very flat gradient of ñ at the sampling location (see in Fig. 11), both positive and negative instantaneous directions can be expected due to weak temporal fluctuations.The existence of front-supported flames at this position, meaning that the flame propagates from a lean mixture near the centre-line towards a richer mixer further away from it, can be explained by the lean recirculated gases near the centre-line, i.e. negative mean axial velocities between x = 30 mm and x = 60 mm as can be seen in Fig. 7.The significant clustering of points in the top and bottom left quadrants of Fig. 14 implies that in this zone, both front-and back-supported flames are equally likely.
Figure 15 shows the directions of gradients of c and ñ at a position further downstream located on the outer shear layer, i.e. x = 60 mm and r = 45 mm.Wrinkling of the flame front is much more pronounced at x = 60 mm and r = 45 mm than at x = 30 mm, r = 25 mm, whereas the mean temperature gradient is still negative and the temperature fluctuations are at their maximum at this location (see Fig. 10).This leads to a larger number of samples with positive directions of c-gradients in the downstream location while flame propagation in positive radial direction remains still likely.
Gradients of c and ñ are likely to be aligned in this area so that back-supported propagation predominates in this zone.This is not surprising, as in this region being located on the outer shear layer and outside the recirculation zone, no backflow of lean gases is to be expected.The clustering in the top left quadrant observed in the previous axial location does no longer occur at this location, implying that at this point, the richer mixture is usually located nearer to the centre-line.

Conclusions
Large Eddy Simulations of the turbulent stratified TECFLAM swirl flame have been performed using a recently developed reaction rate closure based on the SDR of reaction progress variable, in addition to using two existing FSD based reaction rate models [11,12].The SDR based closure is a flamelet based methodology of turbulent premixed combustion modelling which was demonstrated to perform well based on both a priori [6,7] and a posteriori [8,[33][34][35] analyses in the past for turbulent premixed flames but its performance in turbulent stratified combustion in a complex swirl flame configuration is assessed for LES for the first time in this analysis.In this paper, the performance of SDR based closure is compared to the predictions of the algebraic FSD closures proposed by Fureby [11] and Muppala et al. [12], which were shown to perform better than several other alternative FSD models in a previous a posteriori assessment [13].The combinations of chemical source term and sub-grid flux, which yielded optimum performances in previous a posteriori assessments [8,13] have been used for both FSD and SDR closures in this analysis.A comparison of the results obtained from the simulations with experimental observations shows that the SDR closure is capable of satisfactorily reproducing the physical behaviour of the turbulent swirl flame in the TECFLAM burner configuration both qualitatively and quantitatively.The level of accuracy obtained from the algebraic SDR closure has been found to be comparable to that achieved in simulations of the same test case by other authors [19][20][21].The predictions of the FSD models by Fureby [11] and Muppala et al. [12] remain comparable but the recently developed SDR model performs equally well or better than the FSD based closures [11,12] in some respects.Especially, on coarse grids and in highly stratified conditions, the SDR model provides more accurate temperature estimations, whereas the velocity predictions with the FSD closures were still marginally superior on the coarse grid.An analysis on the nature of stratification present in the TECFLAM burner configuration shows that the character of the flame changes in streamwise direction.In the recirculation zone back-and front-support are found to be equally likely while back-support strongly predominates further downstream.The present study shows that the recently developed SDR based closure offers an inexpensive method for LES analyses of complex burner configurations involving turbulent premixed and stratified flames representing the flamelets regime of combustion.
The performance of LES simulations depends on the accuracy of both reaction rate and sub-grid scalar flux closures, which also interact with each other.Moreover, sub-grid scale velocity fluctuation u 0 D also acts as an input parameter to the SDR model and thus the modelling of u 0 D can potentially play an important role in the predictive capabilities of the SDR model.Thus, the effects of the modelling of sub-grid scalar flux and sub-grid scale velocity fluctuation u 0 D on the predictive capabilities of the SDR based closure need to be analysed in detail in the future.Furthermore, it is possible to dynamically evaluate some of the model parameters of the algebraic SDR model (i.e.Eq. (8i)), which is beyond the scope of current analysis but forms the basis of future investigations.show a similar level of accuracy for all three closures and underpredictions are more noticeable than in the corresponding fine-grid results.In terms of mean temperature profiles (see Fig. 19), results of the three closures are quite different.For all models, sensitivity towards grid resolution is stronger with regard to velocity predictions than temperature predictions, as can be deduced from a comparison between Figs. 7-10 with Figs.[16][17][18][19] respectively.However, it is evident that this sensitivity is much weaker for the SDR closure.Judging from this, it is likely that there will be very little difference between the temperature predictions of the SDR and FSD models with further grid refinement.While the FSD closures cannot reproduce the second temperature maximum at position x = 30 mm at the outer shear layer, the SDR closure follows the experimental trend better.Additionally, the SDR model provides closer agreement with the experiment at x = 60 mm at large radii.These differences are assumed to be related to the way the reaction source term is evaluated by the closure models.The mentioned positions are located in areas of high stratification, i.e. small values of equivalence ratio and laminar burning velocity.While the reaction source term for the FSD models depends proportionally on the laminar burning velocity (see Eq. ( 10)), the laminar burning velocity dependence of the SDR based closure is more complex and is expected to have a smaller impact on the calculation of the source term in this area.Furthermore, gradients of c are expected to be small here (not shown, but the qualitative trend of c is adequately represented by mean temperatures), leading to a lower estimation of the mean reaction rate by the FSD models compared to the SDR closure.Regarding temperature fluctuations, results produced by the SDR simulations are preferable compared to the FSD results.There is no major difference between the SDR and FSD models regarding the ease of implementation and the computational effort for the simulation.Thus, the results in Figs.7-10 and Figs.[16][17][18][19] indicate that the SDR model is better suited to predict the temperature distribution in highly stratified flames on coarse grids, i.e.where a large portion of the flame wrinkling appears on the sub-grid level, while being slightly less preferable for predictions of velocities.

Fig. 2 .
Fig.2.Instantaneous fields of vorticity magnitude (normalised with the maximum value outside the burner, 5.9 Â 10 3 /s) for the isothermal and reactive case.

Fig. 3 .
Fig.3.Radial profiles of mean and fluctuating axial velocities at various axial positions for the isothermal case on the coarse and fine grid compared to experimental results.

Figure 2 Fig. 4 .Fig. 5 .
Figure2shows instantaneous snapshots of normalised vorticity magnitude fields for both isothermal and reactive simulations on

Fig. 6 .
Fig. 6.Instantaneous fields of (a) Favre-filtered axial velocity ũ1 and Favre-filtered mixture fraction ñ, (b) Favre-filtered temperature T and the filtered reaction rate term according to Eq. (1) from the fine grid simulations.Note that not the entire computational domain is shown.

Fig. 7 .
Fig.7.Radial profiles of mean and fluctuating axial velocities at various axial positions for the reacting case from SDR and FSD simulations compared to experimental results.

Fig. 8 .Fig. 9 .
Fig. 8. Radial profiles of mean and fluctuating radial velocities at various axial positions for the reacting case from SDR and FSD simulations compared to experimental results.

Fig. 10 .Fig. 11 .
Fig. 10.Radial profiles of mean and fluctuating temperatures at various axial positions from SDR and FSD simulations compared to experimental results.

Fig. 12 .
Fig.12.Series of instantaneous temperature fields from the fine grid SDR simulation (10 ms between frames).

Fig. 13 .Fig. 14 .Fig. 15 .
Fig. 13.Radial profiles of mean resolved and the modelled normalised values of scalar dissipation rate at various axial positions.

Figures 16 -Fig. 16 .
show mean and fluctuating velocity profiles obtained from SDR and FSD simulations on the coarse grid.Generally, predictions are less accurate than those produced by the fine grid simulations, while in terms of mean velocities the FSD closures behave better, especially on the centerline close to the burner exit.Predictions of turbulent fluctuations of velocities

Fig. 17 .Fig. 18 .
Fig.17.Radial profiles of mean and fluctuating radial velocities at various axial positions for the reacting case from SDR and FSD simulations on the coarse grid compared to experimental results.

Fig. 19 .
Fig. 19.Radial profiles of mean and fluctuating temperatures at various axial positions from SDR and FSD simulations on the coarse grid compared to experimental results.
055 is the stoichiometric mixture fraction.Substituting Y Fu ¼ Y F1 n and Y Fb ¼ 0 for fuel-lean mixture in Eq. (4), and subsequently using the transport equations of Y F and n, yield the transport equation of c as: qðDc=DtÞ ¼ r:ðqDrcÞ þ 2ðqD=nÞrc Á rn where the cross-scalar dissipation term 2ðqD=nÞrc Á rn originates due to the variation of Y Fu ¼ Y F1 n caused by the fluctuations of mixture fraction (or equivalence ratio).

Table 1
Flow parameters in TECFLAM burner.