LES/CMC of Blow-off in a Liquid Fueled Swirl Burner

Large Eddy Simulations of two-phase flames with the Conditional Moment Closure combustion model have been performed for flow conditions corresponding to stable and blow-off regimes in a swirl n-heptane spray burner. In the case of stable flame (i.e. low air velocity), the predicted mean and r.m.s. velocities and the location and shape of the flame agree reasonably well with experiment. In particular, the presence of localised extinctions is captured in agreement with experiment. Using model constants previously calibrated against piloted jet methane flames (Sandia F) with localised extinction, we obtain that at the experimentally determined blow-off velocity of the swirling spray flame, the predicted flame also blows off, demonstrating that the LES-CMC approach can capture the global extinction point in a realistic configuration.

It requires application of sophisticated experimental equipment and advanced combustion models. In two-phase flows the situation is complicated further as the analysis must account for the multiple interactions that occur between the liquid spray and the turbulent flow field.
Experimental and numerical studies of strong unsteadiness in two-phase combustion have focused so far mainly on spark ignition and auto-ignition phenomena. The spark ignition of a heptane spray in a bluff-body configuration was examined experimentally in Ref. [1] and then it was simulated numerically in Ref. [2] using the Reynolds Averaged Navier-Stokes (RANS) approach combined with the Conditional Moment Closure (CMC) combustion model [3]. Simulations of auto-ignition of heptane sprays in a high-pressure, high-temperature closed combustion vessel has been studied applying the RANS-CMC method [4][5][6]. The application of Large Eddy Simulation (LES) for modelling of ignition together with the flame propagation in two-phase flow (kerosene-air) in an industrial combustion chamber was studied in Ref. [7] applying thickened flame combustion model and in Ref. [8] with the LES-Eulerian PDF approach. Although these studies had a rather demonstrative character, since no detailed comparison with experiment was undertaken, they showed a high potential of the LES approach to model a wide range of phenomena in two-phase combustion. Very promising results obtained using the LES-Eulerian PDF method were presented recently for a lab-scale swirl stabilised kerosene spray in a simplified combustion chamber [9,10] and good agreement with experimental data was demonstrated. Probably the first LES-CMC computations of auto-ignition in Diesel engine like conditions resulted to an acceptable agreement with experimental findings [11]. It is evident that LES, with an advanced combustion model such as CMC or the PDF method, have great potential in capturing flame transients.
Recently, an experimental analysis of spark ignition and blow-off processes of an n-heptane spray flame in a bluff-body swirl burner was done [12][13][14] and the corresponding LES-CMC of spark ignition [15] showed that the LES-CMC can capture the spatial distribution of the ignition probability. Both experiment and LES have shown a very rich and variable behaviour of the flame following spark ignition and, to a large extent, ignition probability becoming lower than unity has been attributed to localised quenching, a physical phenomenon that also lies at the heart of the global blow-off of flames.
This work focuses on capturing extinction phenomena using LES and CMC in spray flames. In the case of gaseous fuels, it was shown [16][17][18][19] that LES-CMC is able to capture correctly both the local extinction and ignition phenomena in jet flames, while the lift-off of a swirling natural gas flame (TECFLAM) has also correctly been captured [20]. However, such work has not been performed yet for flames with liquid fuels. The present paper aims to examine if LES-CMC can predict localised extinction also for spray recirculating flames and, more importantly, if the global blow-off condition can be predicted. In keeping to the spirit behind the Sandia D-F Flame research project that partly aims to inform combustion models so that they can then be used for realistic geometries, here we use the model as previously calibrated against the Sandia Flame D-F simulations [16] and apply it virtually unchanged to a flame representative of a liquid-fuelled gas turbine or industrial burner. To the authors' knowledge, this study constitutes one of the first efforts to predict global blow-off with combustion LES.

LES formulation for two-phase flows
In LES the scales of the turbulent flow are divided into the large scales, which are directly solved on a given numerical mesh, and the small scales (subgrid scales) which require modelling. This separation of scales is obtained by a spatial filtering defined as [21,22] where f stands for an arbitrary variable and G is the filter function: with a filter width , being equal to the cube root of a local mesh volume. In variable density flows, Favre filtering is applied almost without exception. It is defined as f = ρ f /ρ where ρ is the density. In this work we consider a low Mach number flow and we apply the so-called low Mach number approximation [23,24] in which the acoustic modes are eliminated from the solution. Applying the filtering procedure to the low Mach number approximation of the continuity and the Navier-Stokes equations gives: where u i are the velocity components and p is the pressure. The stress tensor of the resolved field, τ ij , and unresolved subgrid stress tensor τ sgs ij , resulting from the filtering of the non-linear advection terms, are defined as: where μ is the molecular viscosity determined from Sutherland's law. In this work the subgrid tensor is modelled by an eddy viscosity type model [22] defined as: where S ij = 1 2 ∂ũi ∂x j + ∂ũ j ∂xi and the subgrid (or turbulent) viscosity is computed according to the Smagorinsky model. The source terms appearing in Eqs. 3 and 4 are responsible for mass and momentum transfer between liquid and gas phase. The liquid phase is tracked using a Lagrangian formulation where the fuel droplets were presumed to act as point sources of mass and momentum that modify the vapour fuel distribution, depending on the gas velocity, temperature, pressure and the vapour mass fraction [25][26][27]. The source terms and NS are given as: where V is the volume of the computational cell, N d is the number of liquid droplets in that cell, the symbols m i and v i are the mass and velocity of the i-th droplet computed accordingly with formulas given in [26]. The CMC model presented in the next section belongs to the family of mixture fraction based models. Therefore, together with the fluid flow, the transport equation for the mixture fraction in two-phase flows is solved in the LES: where D = μ/ρ Sc is the molecular diffusivity and Sc = 0.7 is the Schmidt number. The term J sgs is the subgrid part modelled as: J sgs =ρ D t ∂ξ ∂xi with the subgrid diffusivity D t = μ t /ρ Sc t where the turbulent Schmidt number is assumed constant Sc t = 0.4 [28]. The CMC model also needs an estimate of the sub-grid variance of the mixture fraction and this is discussed in the next sub-section.

CMC for two-phase flows
The CMC combustion model was formulated independently by Klimenko and Bilger in the 90's and then it was summarised in a review paper [3]. For two-phase flows the CMC equations were derived by Mortensen and Bilger [29] starting from the two-fluid approach of Kataoka [30]. In the context of LES, the CMC equations were obtained by applying the conditional filtering approach and using the Filtered probability Density Function (FDF) [28,31,32]. Following this approach the LES-CMC equations for dilute spray combustion are formulated as: ßÞ contribution from the liquid phase (10) where the contribution from the liquid phase has been explicitly pointed. The operator ( ·|η) = ( ·|ξ = η) stands for the density-weighted LES filtering conditioned on the mixture fraction ξ [28,32]. The symbol Q α = Y α |η corresponds to the conditionally filtered mass fraction of species α, and u j |η, N|η, ω α |η and |η are the conditionally filtered velocity, scalar dissipation rate, reaction rate, and mass evaporation rate respectively. The Kronecker delta, δ α, f , is equal to one for the fuel and zero otherwise. The equation for the conditionally filtered enthalpy, Q h = h|η, takes the form: ßÞ contribution from the liquid phase (11) where the radiation and pressure work terms have been neglected. The terms e α and e h in Eqs. 10 and 11 represent contributions from the subgrid scales and account for the conditional transport in physical space and also in mixture fraction space due to presence of evaporation source terms (discussed later). The conditionally filtered variables are related to the LES filtered variables as, where P is the Filtered probability Density Function [31,33]. In this work we assume that P(η) has a β-function shape parametrized by the resolved mixture fraction ξ and its subgrid variance ξ 2 . In simulations of evaporating sprays the mixture fraction varies between zero and its local saturation value ξ s that, except for high temperature regions, is usually low. As pointed out in Refs. [34][35][36], a particular treatment may be needed in evaluating the beta-function as it must be properly defined in the range 0 < ξ < ξ s (although, for boiling droplets, ξ s = 1 and so the standard usage of the beta-function may be sufficient). Ge and Gutheil [35] proposed a scaled β-function PDF defined as: where (x) is the Gamma function and α, β are defined as: α = ξ( ξ(1 − ξ )/ ξ 2 − 1) and β = α(1 − ξ )/ ξ and the parameters η max,min have to be adjusted according to flow conditions. In this work, we use Eq. 13. For the lower limit we assume η min = 0, whereas the upper limit was set η max = 0.7 which was estimated based on a solution with a fully developed flame (see Section 4) and examining the maximum droplet temperature (and hence saturation mixture fraction reached). Note that for η max = 1 and η min = 0, Eq. 13 reduces to the classical β-function PDF.
The terms e α , e h , ω α |η, u j |η, N|η as well as the source terms h |η and |η are unclosed and require modelling. In the presence of spray the term e α is defined as: where the second term on the right hand side involves the subgrid conditional joint fluctuations of the evaporation rate and species. This term is analogous to the conditional fluctuations in the RANS-CMC approach defined as u j Y α |η [6,29]. In [6] this term was neglected, hence, in the LES where the subgrid fluctuations are expected to be small, we proceed in the same way, i.e. the second term on the right hand side of Eq. 14 is neglected. The joint conditional fluctuations of the velocity and species are modelled with a gradient model [28,32]: and a similar expression is assumed for the joint conditional enthalpy-velocity fluctuations u j h|η − u j |ηQ h appearing in the definition of e h in which the subgrid conditional joint fluctuations of the evaporation rate and enthalpy are also neglected. The conditionally filtered turbulent diffusivity is modelled as D t |η ≈ D t [28]. The chemical reaction rate is modelled by first order closure [3]: where n is the number of reacting scalars. The chemical scheme used is discussed later. Finally, the conditional velocity may be modelled by a linear model [3] or by explicit conditional filering [32] or assumed constant [28]. Here we follow the simplest approach and assume the conditional velocity to be uniform in η-space and equal to the local filtered velocity, i.e. u j |η ≈ u j . The most debatable sub-model at this stage of maturity of the LES-CMC approach for sprays seems to be the modelling of the conditional source terms h |η, |η, which are responsible for the mass and heat exchange between the gas and the spray. Due to lack of consensus how to model these terms and arguing that they are small, a common practice has been to neglect them together with the whole contribution from spray, see Eqs. 10 and 11. This took place for instance in RANS simulations [4,5,37], but also in recent LES-CMC [11]. In Ref. [11], studying the auto-ignition of heptane spray in a bomb, it was observed that the auto-ignition times were over-predicted for high temperature cases. Borghesi et al. [6] found that inclusion of the spray terms affected the solution in the region where the evaporation was the strongest and that these terms delayed the auto-ignition slightly, which suggests that the spray terms in the CMC equation could probably improve the results obtained in Ref. [11]. The terms h |η, |η have to be related to their unconditional counterparts. The latter is connected with defined by Eq. 7, and the former is related to: where the symbols θ i and Cp L are the temperature of the i-th droplet and the heat capacity of the liquid fuel. Modelling of |η in the context of RANS approach was proposed in Ref. [6] and it was further extended to LES for modelling spark ignition of liquid-fuelled burners [15]. A very similar approach was later used in [38] for modelling of LES-CMC of an acetone flame. The difference between the model proposed in [15] and the formulas used in [38] and [6] is the method of conditioning of |η. In Ref. [6] it is done based on the fuel mass fraction on the droplet surface with |η computed individually for every droplet. In Ref. [38] the conditioning is based on the filtered mixture fraction as obtained from the LES solver which is usually lower than at the surfaces of evaporating droplets. For their flame, the authors claimed that this approach gives more accurate results than the ones obtained with the model proposed in [6]. In the present work and in [15] the evaporation rate is based on a cell-mean value of the fuel mass fractionξ s at saturation conditions. The calculation ofξ s is performed by averaging over the droplets residing in particular cells. Thus, the model for |η is defined as: We assume that for each droplet ξ s is equal to the surface fuel mass fraction computed based on the saturation pressure and the droplet temperature [6,39]. Equation 17 ensures that |η satisfies the constraint that the convolution integral of conditionally filtered variable with P(η) leads to the unconditionally filtered value as in Eq. 12.
Modelling for h |η is done analogously to Eq. 17. The conditional scalar dissipation rate is computed with the Amplitude Mapping Closure (AMC) model used widely in RANS-CMC [40][41][42] defined as: where erf is the error function. The filtered scalar dissipation rate N is computed as the sum of the resolved and subgrid part [16,17,32]: (19) where C N is the model constant which is discussed later. The model for the mixture fraction variance is formulated following Refs. [15,43,44]: where we substituted the model (17) for |η. Modelling of ξ asξ s is similar to the model proposed in Ref. [45] where ξ = 1 ρV È Nd i=1 ξ s, iṁi was suggested. Actually, when ξ s, i in a given cell is the same for all droplets, the model for ξ used here and the one from Ref. [45] are equivalent. The advantage of the present approach relies on the fact that it is consistent with the definition of the conditionally filtered evaporation rate |η.
The constant C V in Eq. 20 is equal to 0.1 as suggested in [46] for gaseous flows. More questionable is the value of the constant C N in Eq. 19 as there is no clear guidance for its value, especially for sprays. The literature suggests that depending on the flow problem, the value of this constant varies considerably. For instance, in the LES-CMC simulations of hydrogen auto-ignition [17,47] and of methane bluff body stabilized flames [48] C N was assumed of the order of O(1), whereas in [16] the value of C N = 42 was used as the result of calibration against experimental data of scalar dissipation rate from Sandia Flame D. The same value was used in [19] for prediction of the local extinction for the Delft III flame with satisfactory results. In addition, it was found that using the constant as calibrated against the Sandia Flame D reproduces flame lift off in the natural gas swirling TECFLAM flame [20]. In the present paper, it is demonstrated that C N may have also a crucial influence on the flame behaviour in two-phase flow flames, however it is also shown that the value of C N = 42 selected by tuning against the Sandia Flame D does indeed give reasonable results also for the present swirling spray spray flame and, in particular, it predicts global extinction at the experimentally-determined blow-off condition.
The CMC model is very expensive computationally and therefore a typical approach uses two separate meshes, i.e. one for the flow field (called CFD mesh) and the second (coarser) for the CMC equations (called CMC mesh). Hence, the conditional variables ( f |η) that appear as coefficients of the CMC equations and that have been computed on the CFD mesh must be transferred to the CMC mesh where the CMC equations are solved. Details of implementation issues of the CMC model together with various strategies for linking CFD and CMC meshes are presented in Ref. [28]. In the present work, the simplest approach was used in which the conditional variables on the CMC mesh (denoted as f |η * ) are determined by the mass weighted volume integral within the CMC cells (V CMC ) defined as: Thus, the conditionally filtered variable f |η * corresponding to the CMC cell is common for a group of the CFD nodes embedded in that CMC cell. The resolved variables f are computed from Eq. 13 in which f |η is replaced by f |η * .

LES solver and solution strategy
The CMC model was implemented in a 2nd-order finite-volume code named PRECISE (unstructured version) [49] with an implicit 2nd-order integration method in time. The CFD and CMC meshes consisted of 3 × 10 6 and 54,000 hexahedral cells, respectively. Inside the combustion chamber the CFD mesh was stretched radially and axially towards the bluff body and refined near the air inlets. The CMC mesh in the region near the bluff body (radially: The CMC equations were solved applying the operator splitting approach where the transport in physical space and transport in mixture fraction space are solved separately and sequentially. Additionally, in mixture fraction space the Strang splitting method was used with respect to the source terms and diffusive terms. The second-order terms in physical and mixture fraction spaces were discretized with 2nd-order central differences. The TVD (Total Variation Diminishing) scheme with van Leer limiters [50] was applied for the convective terms in physical space. The numerical time step was assumed constant and equal to t = 5.0 × 10 −6 s, which on the CFD mesh resulted in CFL number around 0.2. The same time step was used for droplets tracking and for the CMC model. The CMC mesh is coarser than the CFD one and this yields even smaller CFL number on the CMC mesh. Therefore, the CMC equations in physical space were integrated in time applying a first-order explicit Euler method. However, expecting high velocity fluctuations during possible local extinctions or blow-off, the explicit integration procedure was supplemented with an automatic time step reduction procedure. In the case of sudden jump of the velocity resulting in CFL number larger than 0.5 the time step for the CMC equations was properly decreased and inner loops were performed. It is necessary to note that such situations took place only during flame initiation where very large density and temperature gradients occurred. The numerical mesh in mixture fraction space consisted of 51 nodes stretched in the vicinity of the stoichiometric mixture fraction, which is equal to ξ st = 0.062 for heptane. In mixture fraction space the time integration was fully implicit. The 2nd-order discretization applied for the difusive terms resulted in the tri-diagonal system that was solved with the TDMA (Tri-Diagonal Matrix Algorithm) algorithm. The last step of the splitting procedure was the solution of the chemical terms of the CMC equations, and at this stage the VODPK (Variable-coefficient Ordinary Differential equation solver with the Preconditioned Krylov method [51]) solver, which is suitable for stiff systems, was applied.

Modelling the chemical reaction
The chemical reaction was modelled by a modified one-step chemistry for heptane, following the method proposed in [52], which is based on tuning the heat release rate as a function of the local equivalence ratio to give approximately the correct flame speed and temperature across the whole flammable range. The scheme for heptane used here has been developed by Richardson [53].
In this work, we are focusing on extinction and so it is important to consider how this scheme may reproduce the extinction conditions of a non-premixed flamelet. Simplified analysis may be performed using the present CMC code by removing the physical transport and spray terms in Eqs. 10 and 11. Computations with use of simplified CMC equations (usually denoted as "0D-CMC") can be thought of as laminar transient flamelet calculations with unity Lewis number and a prescribed scalar dissipation.
Solutions of the 0D-CMC equations allow to determine the level of scalar dissipation rate at which the flame extinguishes. The so-called critical value N 0,crit is estimated by solving the 0D-CMC equations starting from a steady burning solution for small N 0 . Then by increasing N 0 progressively, successive steady solutions are obtained up to the moment when N 0,crit is reached. When this occurs the temperature and species jump to their inert profiles-this process is shown in Fig. 1. The present mechanism gives that N 0,crit = 273 s −1 , which compares well with a detailed chemistry calculation [54] that gave N 0,crit = 285 s −1 . We note that the solutions obtained in this way cannot be regarded as a complete proof of the accuracy of the modified one-step scheme. Nevertheless, these results certainly have a demonstrative character illustrating the performance of the present scheme in an idealised steady situation. Good agreement of one-step and detailed chemistry in such conditions is a basic condition needed to correctly capture extinction. More elaborate study in this context would require analysis of a dynamical behaviour of the schemes, for instance how they respond to a time-varying scalar dissipation rate. A parametric study taking into account the frequency and amplitude of scalar dissipation rate variations would allow for more detailed comparisons and conclusions.

Flow configuration and initial and boundary conditions
The computational domain is shown in Fig. 2. The flow has been developed in an experimental investigation to study spray flame blow-off [14]. It consisted of a circular duct of inner diameter D = 37 mm fitted with a conical bluff body of diameter D b = 25 mm. A 60 • swirler was located 40 mm upstream of the bluff body. Analysed cases included stable combustion regimes (in the following denoted as SWH1) with air bulk velocity at the inlet U bulk = 10.73 m/s, and blow-off conditions (denoted as SWH3) with U bulk = 13.95 m/s that corresponded to the experimental blow-off point [14]. Due to the conical shape of the bluff-body the air flowing to the combustion chamber accelerates and at the bluff-body edge reaches a velocity U b = 14.3 m/s for flame SWH1 and U b = 18.5 m/s for SWH3. In the computations the boundary conditions for the gaseous phase were: (i) inlet velocity with constant temperature; (ii) no-slip adiabatic walls; (iii) convective outflow. A hollow-cone atomiser was placed inside the bluff body and the liquid fuel (n-heptane) was injected through a slot of 0.15 mm diameter. The fuel flow rate was equal to 0.12 g/s and was kept constant for all air flow velocities. In the simulations the fuel injection was modelled by 64 discrete injection points uniformly distributed along a circle with 0.15 mm diameter. The temperature of the droplets  [13] was 300 K, their velocity were estimated based on the mass flow measurement and was 9.9 m/s, and the droplets were injected in a deterministic manner without any velocity fluctuations. We assumed ten classes for the droplet diameters, which were computed from a modified Rossin-Ramler distribution with the Sauter mean diameter equal to 30 μm and with the exponential factor q = 3 in the Rossin-Ramler formula [39]. The boundary conditions for the liquid phase were: (i) deterministic inlet velocity, injection angle and location (ii) the droplets hitting the walls were reflected. Additionally, the minimum droplet diameter was limited to 1 μm, and smaller droplets were regarded as evaporated. The maximum droplet temperature was limited to the boiling temperature for heptane, T bn = 371 K.
The conditional variables in mixture fraction space were initiated as inert. The temperature at the boundary η = 0 was assumed 300 K and the species mass fraction corresponded to pure air. The boundary condition at η = 1 is more problematic as it cannot be a priori said what value of the heptane mass fraction will result from the droplets and what will be the temperature of the evaporated fuel. For safety and simplicity we assumed that at η = 1 the heptane mass fraction Y C 7 H 16 | η=1 = 1 and additionally we assumed that this may happen only when the fuel orginates from the droplets having the boiling temperature, i.e. in mixture fraction space we have T η=1 = T bn . In the physical domain the boundary surfaces for the conditional variables at η = 0 and η = 1 were kept constant, whereas in the range 0 < η < 1 they were: (i) assumed as inert at the inlet plane and kept constant for the entire simulation time; (ii) computed from the Neumann condition of zero gradient at the walls and outflow.
The computations were performed for conditions SWH1 and SWH3 with two values of the constant C N in Eq. 19, i.e. with C N = 2 and C N = 42 as these values are used in the cited literature. To verify the accuracy of the PRECISE code and to check numerous computational settings, simulations for the cold flow conditions without the fuel flow were also performed. In the experiment [14], the blow-off condition SWH3 was obtained by gradually increasing the air flow from the steady flame condition SWH1. Due to the unsteady nature of the blow-off phenomenon and particularly from the point of view of the correct prediction of the blow-off time (i.e. the duration of the blow-off event), it is very important to consider the initial conditions before the velocity is altered. In the present work, seven simulations were run, summarised in Table 1, to which we will refer in the next sections. A comparison between the stable flames Case-A1 and Case-A2 will inform us about the sensitivity of the predictions to the constant C N , while a comparison between Case-A1 with Case-B1 and Case-A2 with Case-B2 will inform about the effect of the flow velocity. The approach to blow-off was done experimentally by starting from the stable flame and increasing the air velocity, a procedure that is reproduced in simulations B1 (for a low value of the constant C N ) and B3 (for the value of C N deemed appropriate for the Sandia Flame D-F simulations [16]).

Velocity Measurements
The mean and fluctuating axial and swirl velocities have been measured for cold flow (i.e. with no fuel flow) and for the stable flame SWH1 with Laser-Doppler Anemometry. The details of the LDA set-up are given in Ref. [55]. The measurements used 2000 data points, which resulted in a statistical uncertainty of 1-3 % depending on the location. For the cold flow data, which are used in an attempt to validate the basic LES, the measurements are relatively straightforward. For the spray flame, the uncertainties in the region where droplets exist are larger, as it is not clear if the velocity recorded is due to the tracking particles seeded with the air or due to the spray droplets, since no amplitude discrimination was made. A comparison of these measurements with droplet-only velocity measurements [55] and use of the Mie scattering images [14] suggests that: (i) at axial locations farther downstream than about 40 mm from the bluff body, there are no more droplets so the LDA signal comes mostly from the seed particles, hence truly representing the air flow; (ii) at radial locations corresponding to the air annular jet, the droplets are captured by the air flow and the droplets and air have almost the same velocities. Therefore, the cold flow data are unambiguous to interpret everywhere, while the spray flame data are unambiguous to interpret only for z > 40 mm, with earlier profiles used for indicative purposes only.

LES verification of the cold flow
Prior to the analysis of the flame behaviour in the reacting cases, the accuracy of the LES code, the boundary conditions, and the mesh density were verified for the   Table 1) and SWH3 (Case-B in Table 1)     .75 the experimental data are almost twice as large as the LES results. Nevertheless, the trend is well predicted, and hence, one may conclude that in the cold flow conditions the aerodynamics in the burner are captured satisfactorily.

Low velocity flame, sensitivity to C N
The simulation for Case-A1 started from the cold flow conditions (Case-A), then the liquid fuel was injected and after a while the spray spread in the whole domain. The combustion was initiated by a numerical spark located at (z, r) = (13, 10) mm. The spark was modelled by prescribing a "burning flamelet solution" (Fig. 1 for N 0 = 5 1/s) extending over 3 × 3 × 3 CMC cells, from which the flame developed in the chamber. Figure 8 shows the heat release value (60 MJ/m 3 s) in selected time instants during the flame expansion. We can see that starting from the region where the spark was located the flame grows successively and finally, in about 20 ms, the flame stabilizes and fluctuates remaining attached to the bluff-body.
The flame exhibits almost no signs of extinction. However, this is not fully consistent with the experiment because, even in the stable SWH1 conditions, observations with OH-PLIF show that the flame is locally torn that could be interpreted as local instantaneous extinction [14]. Further analysis showed that the flame behaviour changes depending on the way of computing the conditional scalar dissipation rate (SDR). Alternative tests included volume averaging weighted with the PDF function (see [28]) and also application of the AMC model directly on the CMC mesh. Here one should note that applying the AMC model to compute conditional SDR always leads to a bell-shape distribution scaled by some maximum value. Indeed, it turned out that this maximum was the most influencing  [16] in the case of the LES-CMC simulations of methane diffusion flames. There, it was shown that the local extinction phenomenon is sensitive to the scalar dissipation level. Garmory and Mastorakos [16] managed to capture the extinctions by adjusting the model constant for the subgrid SDR defined by Eq. 19. By setting C N = 42, instead of the generally assumed C N = 2, they obtained a higher level of total SDR which agreed reasonably well with experimental data from Sandia D, and this then resulted in relatively good prediction of the amount of localised extinction for Sandia F. According to the AMC model the increased level of the subgrid SDR influences on N|η level. In the present configuration the experimental data do not include measurements of N|η nor N, nevertheless because of lack of any visible sign of extinction in the simulations it seemed reasonable to modify the value of C N . Hence, without having experimental reference point concerning the level of SDR the next computations were performed with C N = 42, although the flow problem considered here and in [16] are drastically different. In the computations, the constant C N had to be increased progressively, i.e. (c) OH-PLIF experiment Fig. 9 a, b Contours of the instantaneous heat release for simulations SWH1. The black line corresponds to the stoichiometric mixture fraction. c Instantaneous snapshots of OH from Planar Induced Fluorescence from the experiment [14] for flame SWH1 a few hundreds of time steps had to be performed with C N = 20, then with C N = 30, and finally we could set C N = 42 (this procedure was necessary for stability).
The results obtained with C N = 2 (Case-A1) and C N = 42 (Case-A2) are compared in Fig. 9 showing a few instantaneous contours of the heat release in the central cross-sectional plane of the combustion chamber. The black line represents the stoichiometric mixture fraction where nominally the heat release should be high, if localised extinction had not occurred. Additionally, in Fig. 9c the experimental data showing instantaneous OH snapshots are presented for comparison. For Case-A1 the selected time instants correspond to the results presented in Fig. 8, whereas in the simulations Case-A2 the time distances are counted from the time when C N was set to 42. Clearly, in the simulations of Case-A2 with the high value of the constant C N the flame is much shorter and localised extinctions are visible, however, the flame is overall stable and certainly far from global extinction. The flame shape and the presence of localised breaks in the flame sheet is fully consistent with the experiment (Fig. 9c). Figure 10a shows instantaneous contours of the mixture fraction together with the liquid droplets (in the figure, the droplet sizes are proportional to the diameter in the range 1-100 μm). As one may see in Fig. 10a, many droplets survive in the region of the flame and, having a large inertia, they leave the recirculation zone-some of them move directly towards the outflow while the others hit the walls. The contours of the mixture fraction show that the flame attaches to the bluff body and its maximum value is around ξ max = 0.57 and does not change with time considerably. The temperature contours in Fig. 10b are shown for the range > 301 K, and it is evident that the cold flow (seen as the holes in the contours close to the bluff body edge) penetrate the combustion chamber for about 0.5D b .
The results showing the time averaged heat release contours and experimental data from OH-PLIF measurements overlapped with Mie-scattering (white lines) data are shown in Fig. 11. The symbols P1, P2, P3, P4 in Fig. 11 denote locations of the nodes of the CMC mesh where the solution will be analysed in mixture fraction space. The time averaging was performed over 0.15 s and this took about 400 h of continuous simulations using 32 processors. Despite this nominally long averaging time, it is seen that the averaged heat release values are not fully converged to the (a) mixture fraction and droplets (b) temperature statistically steady state as the contours are not fully symmetric as one could expect. Nevertheless, they certainly allow an assessment of the flame and droplet position with respect to the experimental data. The white isolines in the computational results (Fig. 11, left) represent the contours of the liquid volume concentration, which indicate the location of the droplets. Numerical results show that the inertia of the droplets is very high and their spreading rate is very small. This results that the droplets follow trajectories along the injection angle for a long distance. Although a very similar behaviour is seen in the experimental data, the computational results seem to be more narrow. This may result from the fact that in the computations the droplets are injected without any randomness, neither in location nor in velocity. In general, analysing Fig. 11 one may conclude that the shape of the flame agrees reasonably well with the experiment. The regions of high heat release are seen downstream of the bluff body (around 1 D b -1.5 D b ) and also in the shear layer of the recirculation zone where the evaporated fuel mixes with the fresh air. The double-conical reaction sheet is therefore captured. The accuracy of the LES results is further assessed by comparison of the axial and tangential air velocity components with the experimental data. The contours of the velocity obtained from LES are shown in Fig. 12. The region where the droplets influence the flow field is clearly seen in the figure with the axial velocity contours. Close to the bluff-body and along the axis (x = 0) the velocity increases due to drag forces exerted by the injected droplets. Apart from that, it is noticeable (a) LES -axial velocity (b) LES -swirl velocity Fig. 12 Case-A2-contours of the mean axial and swirl velocities in a cross-sectional plane that the velocity field exhibits much more symmetrical behaviour than the contours of the heat release presented in Fig. 11. This means that the quantities related to the chemistry are more unsteady than the flow field variables. Most likely this is a consequence of the dependence on the solution in mixture fraction space. The radial profiles of the mean and RMS values of the velocity at various locations from the bluff-body are shown in Figs. 13 and 14. It is seen that they agree with the measured   We conclude that the simulation of Case-A2 seems closer to the experiment than Case-A1 from the point of view of flame shape and presence of localised extinctions, and that there is good quantitative agreement of the predicted mean and RMS air velocities with the experiment. Therefore, the stable flame SWH1 has been captured well by the LES and now the behaviour of the simulation at the higher velocities corresponding to blow-off can be explored.

High velocity flame, Case-B1
The first simulation at the blow-off condition was Case-B1, i.e. starting the simulation from the results obtained for the fully developed stable flame computed with C N = 2. The flow conditions corresponding to the experimental blow-off regimes were achieved by increasing smoothly the inlet air velocity from U bulk = 10.73 m/s (SWH1) up to U bulk = 13.95 m/s (SWH3). The change of the boundary conditions was done in a short time, i.e. 0.002 s, which is quicker than in the experiment that used a 2 % rise in velocity over 20 s [14]; such a long transient would be prohibitively expensive in the context of LES.
The time variations of the conditional temperature T|η and conditional SDR (N|η) at the stoichiometric mixture fraction in the CMC nodes in points P1, P2, P3, P4 are presented in Fig. 15 Sample instantaneous results are shown in Fig. 16 presenting the isosurfaces of the stoichiometric mixture fraction colored by the resolved temperature. The "holes" in these isosurfaces correspond to the temperature below 400 K, and for better visibility the holes were surrounded by the black line. Clearly, there are very few places with the temperature falling below 400 K, the flame is stable and not close to global extinction.

High velocity flame, Case-B2
The differences in the results obtained for the SWH1 condition for Case-A1 and Case-A2 indicated the importance of the constant C N in the sub-grid scalar dissipation model. Here, we extend the investigation of this sensitivity for conditions in which the flame should extinguish according to experiment. Intuitively one should  Fig. 17 present T|η and N|η in the same points as in Fig. 15. Except for location P4, the effect of alteration of C N is evident, with T|η and N|η now varying substantially with time. One may notice that the "inertia" of the CMC solution is very low, i.e. a change of N|η causes almost immediate change of the conditional temperature. This is so not only for the values of T|η at the stoichiometric mixture fraction, but the entire mixture fraction space.   [16,19].
Such behaviour in mixture fraction space has direct consequences for the resolved temperature and also the density, which almost immediately influences the flow. Figure 19 shows the time variations of the total (volume integrated) heat release (Q [kJ/s]) and the maximum temperature (T max [K]) in the combustion chamber. The values of T max obtained for Case-B1 are virtually constant for the entire simulation time, and the total heat release also stays at a high level, with slow oscillations due to the unsteady flow. The results for Case-B2 are significantly different. The arrow points to the moment of time when the constant C N has been increased. One may see that alteration of C N has no sudden impact on T max that remains high and almost constant for 20 ms. The lowering of the temperature starts to be seen only after that time, but even in this case in the analysed period of time T max dropped by 20-30 K only. Persistence of T max at a high level for a long time was reasonable and expected. This is because somewhere in the flow domain there was a place with a very slow scalar dissipation and hence high temperature even if the flame has begun to extinguish.
The flame shrinking is confirmed by analyzing the variations of the total heat release (Fig. 19 on the right hand side). In this case one may see that Q changes almost immediately with modification of C N and at the end of the simulation time it decreased to about one third of the value compared to the case with C N = 2. This certainly had to be related to the flame shape, and indeed, the instantaneous 3D results showed that the flame exhibits very unstable behavior and as the time progress it seems to vanish. This is clearly seen in Fig. 20 showing the isosurfaces of Fig. 20 Case-B2: isosurface of the stoichiometric mixture fraction coloured by the resolved temperature the stoichiometric mixture fraction-it is evident that the flame extinguishes, and in this example the extinction starts from the left side of the bluff body with a very large extinction region.
The simulation was stopped at the last instant shown in Fig. 20 because this simulation was meant only as a preliminary indication on the effect of the scalar dissipation model. Properly simulating the blow-off transient, where the velocity increases from low to high values, and using the more appropriate model for the scalar dissipation constant, is described in the next sub-section.

High velocity, Case-B3
The analysis of variations of the integrated heat release and behaviour of the flame shape for Case-B1 and Case-B2 as well as the previous results for stable conditions SWH1 (Case-A1 and Case-A2) showed that C N , and hence N|η, have a crucial impact on the occurrence of the local extinction. The increase of C N leads to the local extinctions that could possibly lead to global blow-off.
The last simulation, Case-B3, was performed as follows: (i) the solution started from the results obtained previously for the stable flame condition SWH1 (Case-A2; the results of which were discussed in Section 4.2); here, for comparison we additionally analysed how the total heat release and the maximum temperature behave in time; (ii) at some moment in time the inlet velocity was increased smoothly up to the SWH3 condition, which corresponds to the experimental blow-off point. Figures 21-23 show the time variation of the maximum temperature and the maximum velocity modulus |u| = √ u 2 + v 2 + w 2 in the chamber, and the total heat release. The lines with symbols correspond to the SWH1 condition (Case-A2) and the lines without symbols to SHW3 (Case-B3). The arrow points to the time instant (about 0.08 s) when the flow rate switched to SHW3. The maximum temperature for the low-velocity case shown in Fig. 21 remains high for the entire simulation time and its changes are small even if the total heat release oscillates substantially (see Figs. 9 and 23; recall that the flame for Case-A2 remained stable for the total of the simulation time.) The results in the Case-B3 show a considerably different scenario. First, when the inlet air flow increases the maximum velocity in the combustion chamber also increases, which is seen in Fig. 22. We remind that switching from SWH1 to SWH3 conditions was done smoothly by altering the inlet velocity values. This was performed in a short time (2 ms) but the flow needed more time to adjust to the higher air flow rate, and the stabilisation of the flow conditions (assessed using the maximum velocity modulus) required at least 0.01 s. The evident change in the flow conditions did not influence the maximum temperature for a quite long time, with T max staying virtually constant for about 0.025 s and only after that time it dropped down suddenly. This has to be related to a decrease of the integrated heat release Q, which, contrary to the time behaviour of T max starts to change immediately after the high velocity condition was reached, as shown in Fig. 23. After 0.025 s, Q reaches very low levels that correspond to the fast lowering of T max (see Fig. 21). Soon after this, the flame totally vanishes and T max slowly decreases. This is illustrated further in Fig. 24 showing the isosurfaces of the stoichiometric mixture fraction coloured by the temperature. The presented time sequence of the instantaneous solutions cover the period preceding the blow-off and also the moment of time when the flame fully extinguishes. It is seen that initially the isosurface shrinks, then it separates from the bluff-body, and finally completely vanishes as the temperature decreases and the droplets do not evaporate. From the point of view of the blow-off event duration, the experiments have provided an average value by performing many realisations of blow-off [14]. It was found that the blow-off event is a relatively slow process and that the average blow-off event duration was about 12 ms, defined as the time taken for the integrated heat release rate to drop from 90 % of its pre-extinction (stable) value to 10 %, on its way to zero. In the present computations and using the same definition on the data shown in Fig. 23, the blow-off event took approximately 20 ms, which agrees reasonably well with the measured value. We conclude that, by starting from a good LES solution of a stable flame, and using values of the scalar dissipation constant previously calibrated against jet flames for which scalar dissipation data exist, when the velocity increases to the experimentally-determined blow-off velocity, the simulated flame also fully extinguishes. Therefore the simulation captures the global blow-off point of this swirling spray flame, which brings confidence to the use of the LES-CMC model for flame stability problems. The key to this success seems to be the accurate prediction of the presence of localised extinction even in the stable flame at conditions approaching blow-off.

Conclusions
The results obtained here show that the LES-CMC approach offers a way to capture the localised extinction and the global blow-off phenomenon in swirling spray flames. Using model constants for the scalar dissipation rate previously calibrated against piloted jet methane flames with localised extinction, we obtain that the predicted heptane spray flame is extinguished at the experimentally determined blow-off velocity. We also obtain that the flame extinction process during a blow-off event lasts a time duration consistent with measurements. Further evaluation of the accuracy of the LES-CMC for capturing quantitatively realistic flame extinction needs to consider simulating more blow-off points. The present simulations demonstrate a promising degree of maturity of the LES-CMC approach for capturing the blow-off limits of realistic combustors.