Analysis of Flame Front Breaks Appearing in LES of Inhomogeneous Jet Flames Using Flamelets

Abstract

The physical mechanism leading to flame local extinction remains a key issue to be further understood. An analysis of large eddy simulation (LES) data with presumed probability density function (PDF) based closure (Chen et al., 2020, Combust. Flame, vol. 212, pp. 415) indicated the presence of localised breaks of the flame front along the stoichiometric line. These observations and their relation to local quenching of burning fluid particles, together with the possible physical mechanisms and conditions allowing their appearance in LES with a simple flamelet model, are investigated in this work using a combined Lagrangian-Eulerian analysis. The Sidney/Sandia piloted jet flames with compositionally inhomogeneous inlet and increasing bulk speeds, amounting to respectively 70 and 90% of the experimental blow-off velocity, are used for this analysis. Passive flow tracers are first seeded in the inlet streams and tracked for their lifetime. The critical scenario observed in the Lagrangian analysis, i.e., burning particles crossing extinction holes on the stoichiometric iso-surface, is then investigated using the Eulerian control-volume approach. For the 70% blow-off case the observed flame front breaks/extinction holes are due to cold and inhomogeneous reactants that are cast onto the stoichiometric iso-surface by large vortices initiated in the jet/pilot shear layer. In this case an extinction hole forms only when the strain effect is accompanied by strong subgrid mixing. This mechanism is captured by the unstrained flamelets model due to the ability of the LES to resolve large-scale strain and considers the SGS mixture fraction variance weakening effect on the reaction rate through the flamelet manifold. Only at 90% blow-off speed the expected limitation of the underlying combustion model assumption become apparent, where the amount of local extinctions predicted by the LES is underestimated compared to the experiment. In this case flame front breaks are still observed in the LES and are caused by a stronger vortex/strain interaction yet without the aid of mixture fraction variance. The reasons for these different behaviours and their implications from a physical and modelling point of view are discussed in this study.

1 Introduction

The development of efficient, ‘green’ engines requires strong understanding of the physical processes involved along the whole gas path. In modern engines, combustion is characterised by considerable levels of turbulence, which is often desirable to obtain the fast mixing needed to achieve high combustion efficiency. This mixing is imperfect in most combustion systems, thus reactants are not homogeneously burned and lean premixed combustion is only partly achieved due to short residence times and to avoid potential flashbacks (Gicquel et al. 2012; Syred and Beér 1974). These inhomogeneities cause local fluctuations that can lead to undesired local extinctions. Understanding and predicting these processes is a priority for stable and efficient operability of new generation engines.

Many experimental and numerical works have contributed to the understanding of local extinction processes. In non-premixed combustion the mixing field dominates the flame behaviour (Peters 1986) because the flame is attached to the stoichiometric mixture surface and extinctions occur due to scalar dissipation rate (SDR) peaks causing diffusion losses that exceed the local heat release (Ihme and Pitsch 2008). The mixture can then reignite due to the surrounding heat and the decrease of SDR following the local extinction (Caetano and da Silva 2015; Zhang and Mastorakos 2017; Ihme and Pitsch 2008). In premixed combustion the extinction process is complicated by the flame propagation and its coupling with mixing and turbulence. Extinctions were first studied in terms of imbalance between chemical, ignition and flow time scales (Zukowski and Marble 1986). More recent experimental (Watson et al. 1999; Stroomer et al. 1999; Lyons et al. 2005; Baudoin et al. 2013; Tuttle et al. 2013; Kariuki et al. 2015; Chowdhury and Cetegen 2018) and numerical (Zhou et al. 2015; Wang et al. 2017) works found that local extinctions occur when the flame front moves on high vorticity regions in the shear layer due to the high strain rates breaking the flame front. Whether this leads to further extinctions or reignition, however, remains less understood. The extinction process is even less understood when the combustion is partially premixed or stratified, due to the limited amount of measurements. One recent set of measurements in this direction is the Sydney/Sandia flame involving a jet with inhomogeneous inlets (Barlow et al. 2015; Cutcher et al. 2018). In these experiments the amount of local extinctions on the flame was observed to rise when the bulk velocity increases from 50% to 90% of the blow-off speed. A quasi-direct numerical simulation (DNS) was also recently performed by Zirwes et al. (2019) to provide a valuable dataset for model validations, but only the 50% blow-off speed case was simulated due to the extreme cost for higher Reynolds numbers. Thus, LES is required to analyse the higher speed cases. The significant amount of statistics measured in the extinctions region allowed an extensive process of validation of different combustion models including probability density function (PDF) approach with stochastic fields, Lagrangian particles and multi-environment approaches, diffusion flamelet based models, multiple mapping conditioning and hybrid approaches (see Barlow (2018) for a review of these numerical works and Pitsch (2006), Fiorina et al. (2015) for a review on combustion modelling). The performance of an unstrained premixed flamelets model was assessed in a recent study (Chen et al. 2020), where the predictions satisfactorily matched first and second order statistics of temperature, velocity, mixture fraction and its scalar dissipation rate from experiment in the extinction region. More interestingly, holes in the flame front along the stoichiometric line were observed and compared to the observation of spatio-temporal intermittency in the experimental OH and progress variable traces, although their extent was underestimated in the flame closer to blow-off. Whether or how this observed intermittency is related to extinguishing fluid particles in both the LES and the experiments, and to the concept of local extinction in the Lagrangian view, however, remained unclear and will become apparent in this study. For the sake of simplicity, this intermittency will be referred here regardless as a local extinction, to indicate a break on the flame front (an unambiguous definition will be given in Sect. 4). The mechanism for the so-defined extinctions to occur in first place, and how the combustion model can predict them despite its limitations, is also unclear, and is the focus of this work.

Although previous studies on local extinctions using non-premixed flamelets exist (Ihme and Pitsch 2008), the use of unstrained premixed flamelets is controversial due to the insensitivity of the so-constructed flamelet to strain effects which are of great importance in situations of imminent quench (van Oijen and de Goey 2000; van Oijen et al. 2016). Nevertheless, the classical view for which stretch is most effective at scales smaller than the Gibson scale has been contradicted in a number of recent works (Filatyev et al. 2005; Langella and Swaminathan 2016; Doan et al. 2017), where it is argued that small scales have not enough energy or lifespan to impart significant changes to the flame structure, and thus prediction of local extinctions is possible when most of the strain is resolved. The appropriateness of using an unstrained flamelets model in presence of local extinctions was discussed in Chen et al. (2020). In this work the objective is instead to shed light on the mechanism leading to the formation of the observed breaks on the flame front in inhomogeneous jet flames, and on what allows the prediction of these in the unstrained premixed flamelets formulations. The second question in particular is still unanswered within the flamelet community and requires a thorough investigation of the numerics involved. This is thus discussed in the present work along with physical insights - in the limits of the modelling approach used - of the dynamics observed. Close to blow-off flame front breaks are observed to be associated with strain and the alignment of vorticity and stoichiometric line, consistently with the above literature. At 70% blow-off speed, however, the role of strain is more moderate and aided by relatively high values of mixture fraction variance convected onto the stoichiometric line by large eddies. The prediction of extinctions is then associated with the local levels of resolved strain and mixture fraction variance combined.

This paper is organised as follows. The next section describes the experimental burner along with its experimental procedure. The LES model, SGS closures, numerical method and grid used, and boundary conditions employed are described in Sect. 3. The results are shown in Sect. 4 and discussed in Sect. 5, while the conclusions are summarised in the final section.

2 Experimental Case

The Sydney/Sandia burner with inhomogeneous inlet (Barlow et al. 2015; Cutcher et al. 2018) (Fig. 1) having an offset distance of $L_r = 75\,\hbox {mm}$ for the fuel inlet is considered in this study. The experimental campaign for the $L_r= 75\,\hbox {mm}$ cases consists of measurements of radial profiles of temperature, velocity, mixture fraction, major and minor species, and their rms values, at different axial locations. Data for three-directional scalar dissipation rate of the Bilger’s mixture fraction (Bilger et al. 1990), as well as measurements of a progress variable based on oxygen mass fraction, are also available in Cutcher et al. (2018). This experimental dataset is available for three mixture bulk velocities, $U_b = 57$, 80 and $103\,\hbox {m/s}$, corresponding to 50%, 70% and 90% of the measured blow-off velocity $U_{bo}$, which are characterised by a consistent increase in the occurrence of local extinctions from the lowest to the highest speed configuration. These three flames were simulated by Chen et al. (2020), and results were observed to satisfactorily compare to experimental data even in the regions with significant level of extinctions. The two flames with higher speed are further investigated in the present work due to the larger extent of local extinctions. Details of these configurations are summarised in Table 1.

Fig. 1

Schematic of the Sandia/Sydney burner (Barlow et al. 2015; Cutcher et al. 2018)

Table 1 Configuration parameters for the flames analysed. $L_r$ is the axial distance between the exits of the two innermost nozzles, $U_b$ is the bulk jet velocity at the burner’s mouth, $U_f$ and $U_a$ are fuel and air speeds in their respective streams, $U_{bo}$ is the blow-off velocity

3 Methodology

3.1 Combustion Modelling

The interaction between combustion and turbulence at the subgrid scales is taken into account using a presumed PDF approach, where the thermochemistry is modelled using a database of unstrained premixed flamelets (van Oijen et al. 2016). The underlying flamelet assumption states that the smallest turbulent eddies are not small enough, or do not last long enough, to impart significant changes to the flame reaction structure and thus the turbulent flame can be seen as an ensemble of thin, one-dimensional flames (flamelets). This approach has been validated in a number of previous studies, including the present configuration (Chen et al. 2020; Langella and Swaminathan 2016; Langella et al. 2017, 2016a) (see also Pitsch (2006) for a review). Transport equations for a filtered progress variable $\widetilde{c}$ and mixture fraction $\widetilde{\xi }$ are needed to track the reaction progress and mixing level respectively, which can be expressed as:

$$\begin{aligned}&\overline{\rho }\frac{\mathrm {D}\widetilde{c}}{\mathrm {D}t} = \nabla \cdot \left[ \left( \overline{\rho \mathcal {D}} + \frac{\mu _t}{\mathrm {Sc}_t} \right) \nabla \widetilde{c} \right] + \overline{\dot{\omega }_c} \end{aligned}$$

(1)

$$\begin{aligned}&\overline{\rho }\frac{\mathrm {D}\widetilde{\xi }}{\mathrm {D}t} = \nabla \cdot \left[ \left( \overline{\rho \mathcal {D}} + \frac{\mu _t}{\mathrm {Sc}_t} \right) \nabla \widetilde{\xi } \right] \end{aligned}$$

(2)

In the above expressions $\mathrm {D}/\mathrm {D}t$ indicates the total derivative, $\overline{\rho }$ is the filtered density, $\mu _t$ is the sub-grid viscosity and $\overline{\rho \mathcal {D}}\approx \overline{\mu }/\mathrm {Sc}$ is the filtered molecular diffusivity, assumed to be the same for all scalars. The Schmidt number is $\mathrm {Sc} = 0.7$, and unity Lewis number is assumed. The filtered dynamic viscosity $\overline{\mu }$ is computed from the Favre-filtered temperature using Sutherland’s law, while the sub-grid Schmidt numbers is $\mathrm {Sc}_t = 0.7$. The Bray number (Veynante et al. 1997) of both flames has been evaluated and is less than one in the region of interest, indicating that counter-gradient transport is negligible for this configuration. The progress variable c is defined from the mass fractions Y of CO and $\hbox {CO}_2$ (Langella et al. 2017; Fiorina et al. 2003) as $c=\psi /\psi _b$, where $\psi = Y_{{\mathrm{CO}}}+Y_{{\mathrm{CO}_{2}}}$ and the subscript b refers to burnt conditions. This scaled definition is preferred to an unscaled version of the progress variable, as it was observed to better capture the local extinction process (Chen et al. 2020).

The filtered reaction rate $\overline{\dot{\omega }}_c$ in Eq. 1 is modelled as:

$$\begin{aligned} \overline{\dot{\omega }}_c \approx \overline{\dot{\omega }}_{\mathrm {fp}} + \overline{\dot{\omega }}_{\mathrm {np}} \end{aligned}$$

(3)

where $\overline{\dot{\omega }}_{\mathrm {fp}}$ represents the contribution due to stratified premixed combustion mode and $\overline{\dot{\omega }}_{\mathrm {np}}$ represents a correction due to the non-premixed combustion mode contribution, and can be positive or negative. The additional terms describing the effects of SGS cross-correlation of mixture fraction and progress variable on reaction rate are not included as they were shown to be negligible for well-resolved LES (Bray et al. 2005; Ruan et al. 2014). The stratified premixed term is modelled as:

$$\begin{aligned} \overline{\dot{\omega }}_{\mathrm {fp}} = \overline{\rho } \int _0^1 \int _0^1 \frac{\dot{\omega }_{\mathrm {fp}} \left( \zeta , \eta \right) }{\rho \left( \zeta , \eta \right) } P\left( \zeta , \eta \right) \mathrm {d}\zeta \mathrm {d} \eta \end{aligned}$$

(4)

where $\dot{\omega }_{\mathrm {fp}}$ is the precomputed reaction rate of a laminar, one dimensional freely-propagating planar unstrained flame at a fixed mixture fraction, and P is the joint-PDF of progress variable ($\zeta$) and mixture fraction ($\eta$). The joint-PDF is modelled using the Bayesian decomposition $P(\zeta ,\eta )= P(\zeta |\eta )\, P(\eta )$, where for simplicity $P(\zeta |\eta )=P(\zeta )$ under the assumption of statistical independence of $\zeta$ and $\eta$, which is a reasonable assumption in well-resolved LES following past studies (Chen et al. 2018; Bray et al. 2005; Ruan et al. 2014). $P(\zeta )$ and $P(\eta )$ are presumed using beta-functions, which require an expression for the SGS variances of progress variable, $\sigma _{c\text {,sgs}}^2=\widetilde{c^2}-\widetilde{c}^2$, and mixture fraction, $\sigma _{\xi \text {,sgs}}^2=\widetilde{\xi ^2}-\widetilde{\xi }^2$. It should be stressed that the choice of beta-functions is known to yield an overprediction of the reaction rate if the mesh size is larger than the laminar flame thickness and if the flame is close to laminar conditions. In these cases, alternative choices for the PDF are possible (Mahdi Salehi et al. 2013). For the present flame configuration beta-functions are an adequate choice, since grid resolution is sufficiently smaller than the laminar flame thickness, as will be highlighted later in Sect. 3.2. $\sigma _{c\text {,sgs}}^2$ and $\sigma _{\xi \text {,sgs}}^2$ are computed in this work using their transport equations:

$$\begin{aligned}&\overline{\rho }\frac{\mathrm {D}\sigma _{c\text {,sgs}}^2}{\mathrm {D}t} = \nabla \cdot \left[ \left( \overline{\rho \mathcal {D}} + \frac{\mu _t}{\mathrm {Sc}_t} \right) \nabla \sigma _{c\text {,sgs}}^2 \right] \nonumber \\&\quad +2 \frac{\mu _t}{\mathrm {Sc}_t} \left( \nabla \widetilde{c} \cdot \nabla \widetilde{c} \right) - 2\chi _{c\text {,sgs}} + 2 \left( \overline{c \dot{\omega }_c} - \widetilde{c}\overline{\dot{\omega }_c} \right) \end{aligned}$$

(5)

$$\begin{aligned}&\overline{\rho }\frac{\mathrm {D}\sigma _{\xi \text {,sgs}}^2}{\mathrm {D}t} = \nabla \cdot \left[ \left( \overline{\rho \mathcal {D}} + \frac{\mu _t}{\mathrm {Sc}_t} \right) \nabla \sigma _{\xi \text {,sgs}}^2 \right] + 2 \frac{\mu _t}{\mathrm {Sc}_t} \left( \nabla \widetilde{\xi } \cdot \nabla \widetilde{\xi } \right) - 2\chi _{\xi \text {,sgs}} \end{aligned}$$

(6)

In the above equations, $\chi _{\xi \text {,sgs}}$ and $\chi _{c\text {,sgs}}$ are the subgrid scalar dissipation rates (SDR) of mixture fraction and progress variable respectively, and require modelling. The former is closed using a linear relaxation model, $\overline{\rho } \chi _{\xi \text {,sgs}} = C_\zeta \left( \mu _t/\varDelta ^2 \right) \sigma _{\xi \text {,sgs}}^2$, where $\varDelta$ is the LES filter size. The commonly used value of the model constant $C_\zeta =2$ (Pitsch 2006) is also chosen for this work following previous LES studies showing satisfactory results for a range of different configurations (e.g. see Langella et al. (2017); Chen et al. (2020); Langella et al. (2020)), although a recent DNS investigation (Sitte et al. 2020) has pointed out that this value might be an order of magnitude different, at least for the conditions examined in that study. The SDR of the progress variable must instead balance both reaction and turbulence production in Eq. 5, and thus the model proposed by Dunstan et al. (2013) and validated a posteriori in a number of studies (Langella et al. 2017, 2016a, b), is used instead:

$$\begin{aligned} \chi _{c\text {,sgs}} = \mathcal {F} \left[ 2K_c \frac{s_L}{\delta _\mathrm {th}} + \left( C_3-\tau C_4\mathrm {Da}_\varDelta \right) \left( \frac{2{u^\prime _\varDelta }}{3\varDelta } \right) \right] \frac{\sigma _{c\text {,sgs}}^2}{\beta _c} \end{aligned}$$

(7)

where ${u_\varDelta ^\prime }=(2k_\mathrm {sgs}/3)^{1/2}$ is the SGS velocity scale and the SGS kinetic energy $k_\mathrm {sgs}$ is computed with a transport equation. The scale-dependent parameter $\beta _c$ is obtained with a dynamic procedure following previous studies (Langella and Swaminathan 2016; Langella et al. 2015; Gao et al. 2014, 2015). The reader can refer to the LES works above for details about the other model parameters.

The reaction term in Eq. 5 is modelled in a manner consistent to Eq. (4). The non-premixed combustion correction in Eq. 3 is modelled as

$$\begin{aligned} \overline{\dot{\omega }}_{\mathrm {np}} \approx \overline{\rho } \widetilde{c} \widetilde{\chi _{\xi }} \int _0^1 \frac{1}{\psi _b \left( \eta \right) } \frac{\mathrm {d}^2 \psi _b \left( \eta \right) }{\mathrm {d}\eta ^2} P\left( \eta \right) \mathrm {d}\eta \end{aligned}$$

(8)

where $\widetilde{\chi _{\xi }}$ is the sum of resolved and SGS SDR and the term within the integral is also precomputed and stored in the look-up table (Ruan et al. 2014). This modelling strategy for reaction rate is essential to correctly capture the multiple combustion modes present in this flame configuration. Further details about the model can be found in Chen et al. (2020).

3.2 Numerical Setup

The Favre-filtered LES equations for mass, momentum and total enthalpy, as well as for the four additional scalars required by the combustion model described in Sect. 3.1, are solved using the OpenFOAM libraries (Weller et al. 1998) with a pressure-based approach for low-Mach flows. The PISO loop (Bressloff 2011) is used for pressure-velocity coupling, and an additional, external loop with 5 to 10 iterations is used to further enforce convergence of the flow and combustion equations at each time step. A second-order central scheme is used for the velocity convective fluxes,

which is not uncommon in reacting LES (see for example (Butz et al. 2015; Zhang and Mastorakos 2017; Popp et al. 2021; Yu et al. 2020)) and was considered sufficient on the basis of the validation and mesh sensitivity analysis carried out in Barlow (2018) and Chen et al. (2020). Although the use of higher order schemes is possible and could yield improved results (Zhou et al. 2015; Rochette et al. 2018), this choice is known for causing instabilities in the flame region due to the strong gradients, unless grid resolution across the flame front approaches DNS requirements. An implicit Eulerian scheme is used for time advancement following previous works (Chen et al. 2020; Langella et al. 2017) and a constant time step is chosen to ensure a CFL number below 0.3 everywhere in the flow. The sub-grid viscosity is modelled using a transport equation for the sub-grid kinetic energy $k_\mathrm {sgs}$ as in previous studies (Langella et al. 2017; Chai and Mahesh 2012), $\mu _t = C_k \overline{\rho } \varDelta k_\mathrm {sgs}^{1/2}$, where $C_k=0.1$ is a model constant. The boundary layer is treated with a two-layer wall model (Piomelli and Balaras 2002). The temperature field is computed from the filtered enthalpy as $\widetilde{T} = T_0 + (\widetilde{h} - \widetilde{\varDelta h}_f^0) / \widetilde{C}_{p,\mathrm {eff}}$, where $T_0=298.5$ K is a reference temperature, $\widetilde{\varDelta h}_f^0$ is the enthalpy of formation of the mixture, and the effective heat capacity of the mixture is calculated as $\widetilde{C}_{p,\mathrm {eff}} = (\int _{T_0}^T C_p dT')/(T-T_0)$. The filtered density is computed from the operative pressure and filtered temperature via the state equation. The filtered thermochemical quantities are pre-computed for each flamelet and pre-integrated using an equation consistent with Eq. (4). Outside the flammability limits, these are linearly interpolated between the lean (rich) limit value and air (fuel) conditions.

Zero and unity mixture fraction values are set at the inlets for the air and fuel streams respectively, while a stoichiometric value $\xi _\mathrm {st} = 0.055$ is set at the pilot inlet. The progress variable is set to unity at the pilot inlet to represent combusting gases, and to zero at the other inlets. Sensitivity to these values is further discussed in Chen et al. (2020). Air and fuel inlets in the numerical domain are located $100\,\hbox {mm}$ ahead of the burner mouth, which allows the two streams to develop before entering the mixing section. Inlet turbulence is imposed with a synthetic eddy method (Kornev and Hassel 2007) requiring non-uniform velocity root-mean-square (rms) velocity values and turbulence length scales, through an auxiliary RANS (Chen et al. 2020). This approach was shown to be effective in Chen et al. (2020) to handle the inhomogeneous fuel/air mixing upstream of the main jet exit, which has a pronounced impact on the downstream local extinction behaviours. Inlet temperature is set equal to ambient temperature ($300\,\hbox {K}$) for fuel and air streams following the experiments, and to the adiabatic flame temperature at stoichiometric conditions ($2196\,\hbox {K}$) for the pilot stream. Far field and pressure-outlet boundary conditions are used at the sides and outlet of the domain respectively.

The presence of the pilot stream poses the issue of potential spurious flames at the pilot/coflow interface when using a scaled formulation for the progress variable. This is addressed with the introduction of a fuel-tracking variable ($Y_\mathrm {CH_4}$), whose transport equation is similar to Eq. (1), for which the fuel reaction rate is also tabulated. The progress variable reaction rate is computed following Eq. (4) for positive values of fuel mass fraction and is zero otherwise. Correspondingly, the thermodynamic quantities are taken from the table or interpolated with air and fuel conditions depending on the local value of fuel mass fraction. The alternative unscaled definition of progress variable, despite providing similar predictions for the profiles of temperature, major and minor species, was observed to underpredict the levels of local extinction (Chen et al. 2020). The progress variable normalisation required before accessing the look-up table in this case can introduce considerable numerical errors close to the flammability limits, where $\psi _b$ is small, yielding incorrect values of important parameters, such as reaction rate, temperature and effective heat capacity.

A structured grid of 3.2M elements is used, with local refinement at the walls and in the shear layer. The radial element size spans from 0.1 to $0.6 \delta _\mathrm {th,st}$ for $x/D<30$, where $\delta _\mathrm {th,st}$ is the laminar flame thickness at stoichiometric conditions. This mesh was shown to be adequate (Chen et al. 2020) to predict the expected level of mixing and satisfies the criterion of Pope for turbulent kinetic energy (Pope 2000), with $k_\mathrm {sgs}$ observed to be smaller than 5% over most of the domain, including the region of interest, at any time. Hence, turbulent fluctuations mostly occur at the resolved scales and therefore the LES solution will be considered as an accurate representation of the instantaneous unfiltered flow. Mesh sensitivity analysis carried out on three grids having 2.2, 3.2 and 8M elements showed no significant difference in the results obtained with the current 3.2M setup and the finest configuration. The laminar computations were performed using PREMIX (Kee et al. 1985) with the GRI3.0 mechanism for methane/air combustion. The pre-computed table is discretised using 101 and 51 equidistant points respectively for the progress variable and its SGS variance spaces. The mixture fraction space is discretised using 51 points within the flammability limits, while its SGS variance is discretised using 15 exponentially spaced points between 0 and $\xi (1-\xi )$. Each simulation took about 24 hours on 240 cores to compute 12 flow-through times, defined as 30$D/U_b$. Statistics were collected over the last 6 flow-through times.

4 Results

The results from the numerical simulations of the Sydney/Sandia flame with inhomogeneous inlets have been extensively validated by Chen et al. (2020), showing good agreement with the experimental measurements of temperature, Bilger’s mixture fraction and major species, for both mean and rms values. Mean velocity and its rms values were also presented in Barlow (2018) for both reacting and non-reacting cases showing also good agreement with the experimental data for different inlet flow speeds. Some overprediction of reaction rate was reported in Chen et al. (2020) for $x/D>30$, causing an excess consumption of reactants and slightly higher values of flame temperature, which was argued to be due to the dominant non-premixed combustion mode at downstream locations not accurately interpreted by the premixed flamelets approach. Detailed comparison with the experiments, including that for minor species, can be found in the supplementary material to Chen et al. (2020). The experimental data further showed presence of local extinctions as discontinuities in the traces of progress variable and OH mass fraction (Barlow et al. 2015; Cutcher et al. 2018), which were also captured in the numerical results (Chen et al. 2020). This spatio-temporal intermittency was also highlighted by the significant scatter observed near the stoichiometric value in the temperature scatter plots. Most significantly, the LES was shown to predict the correct levels of mixture fraction SDR in the flammable region, which was experimentally observed to play a crucial role in the formation of extinctions (Barlow et al. 2015; Cutcher et al. 2018). Despite the computed probability distributions showed slightly narrower ranges of SDR compared to the experiment, thus hinting at a slightly smaller amount of local extinctions predicted by the LES, the increasing trend of extinction phenomena with bulk speed was also captured, at least qualitatively. Nonetheless, neither the experimental (Barlow et al. 2015; Cutcher et al. 2018) nor the numerical (Chen et al. 2020) study discuss the mechanism of formation and temporal evolution of these extinction holes on the flame front. Moreover, whether the observed intermittency in the progress variable is related to the concept of local extinction in the classical view (burning fluid particles that experience quenching) or by other meanings, is not clear. Keeping this in mind, for the sake of simplicity, this intermittency of the flame front will be referred to as a local extinction for the remainder of the paper. The study of the numerical mechanism that enables the unstrained, steady flamelet model to predict these extinctions, as well as the interpretation of the possible physics behind it, are discussed in this section using Lagrangian and Eulerian analyses.

4.1 Lagrangian Analysis

Fig. 2

Location of the injection points for the passive flow tracers originating from the pilot (red) and main jet (blue) with respect to the mean contours of progress variable for the two flames of Fig. 4

Passive flow tracers are injected every 200 time steps, corresponding to $\varDelta t = 0.04\,\hbox {ms}$ and are tracked during the post-processing stage for $2.8\,\hbox {ms}$ for each simulation. This allows a consistent comparison between the two cases of Table 1 since the thermochemistry tables, thus the chemical time scales, are the same. The choice of passive tracers is necessary for our purposes, as it allows carrying out a Lagrangian inspection of flame front breaks, without any additional Lagrangian modelling of the physics and chemistry that would make the comparison with Chen et al. (2020) inconsistent. It is worth noting that tracking passive tracers during post-processing is conceptually the same as using a Lagrangian particle tracking algorithm during the LES with no Eulerian-Lagrangian coupling. For this reason “flow tracers”, “Lagrangian particles” or simply “particles” will be used interchangeably in the current work. The injection period of $\varDelta t = 0.04\,\hbox {ms}$ is smaller than the eddy-turnover time associated to the smallest integral length scales in the LES, so ensuring appropriate time resolution. A sensitivity analysis (not shown), conducted using a sampling rate of 400 LES time steps, indeed shows that the statistics presented in this section are hardly affected by this change in time resolution. A uniform polar grid with 21 radial and 44 tangential injection locations is defined, resulting in about 900 particles generated at each injection step, of which about 400 originate from the main jet stream and 500 from the pilot stream. This grid (Fig. 2) is aligned with the burner lip and has the outermost injection location equal to the pilot radius. The behaviour of each particle is analysed in time and categorised according to: origin (pilot/main jet); presence of fuel, to distinguish mixing of particles from the pilot with fresh reactants, to that with the co-flowing air ($\widetilde{Y}_\mathrm{CH4}=0$); mixture fraction (particles within/outside the flammability limits); reaction rate, to distinguish between mixing and burning particles; and progress variable, to detect particles that stop burning due to an extinction and due to mixing with cold gases.

Quenching events are characterised by a sudden drop of reaction rate which could be used in principle to identify an extinction. However, the turbulent nature of the flow would make it difficult to distinguish between particles that stop burning due to a local extinction and those that are in a fully burnt or unburnt status. Instead, the quenching event is defined here as a sudden drop of the progress variable of a particle near the stoichiometric line from a burning state to a value below $c^*=0.4$. This arbitrary threshold is supported by the fact that the reaction rates are observed to drop below 5% of the peak value for $c<0.4$, as shown in Fig. 3. A similar approach was used in Mitarai et al. (2003) to detect extinction and re-ignition in DNS data of a diffusion flame.

Using the aforementioned criteria at any time t after injection, the passive flow tracers have been classified in the following separate groups:

• (a) particles originating in the pilot and mixing with the coflow, not burning;

• (b) particles originating in the pilot and mixing with the fresh air-fuel mixture, burning;

• (c) particles originating in the pilot and mixing with the fresh air-fuel mixture that start and then stop burning before complete combustion;

• (d) particles originating in the air-fuel jet that cannot burn because outside flammability limits;

• (e) particles originating in the air-fuel jet, burning;

• (f) particles originating in the air-fuel jet, with mixed behaviour (c)-(d);

• (g) particles originating in the air-fuel jet and within flammability limits that stop burning before complete combustion.

Fig. 3

Contour plot of the flamelet reaction rate normalised by its maximum value in the mixture fraction - progress variable space

Fig. 4

Classification of the flow tracers by origin and behaviour, for the two configurations of Table 1, based on the analysis over a time window of $2.8\,\hbox {ms}$ from the time of injection

The proportion of particles lying in each group is shown, for both flames of Table 1, in Fig. 4. In both cases a sample of 2500 particles for the main and 2000 particles for the pilot streams is used. More than 90% of the particles (groups b and c) originating from the pilot enter the mixing region with the fresh reactants and start burning because of the increase of fuel due to diffusion processes (both resolved and SGS). About 25% and 58% of these (group c) experience quenching for the Lr75-80 and Lr75-103 cases respectively. The increasing trend is due to the higher momentum received by the particles in the mixing region in the Lr75-103 case, which promotes the entrainment of co-flowing air at downstream locations and forces burning particles to quickly cross the lean flammability limit. The behaviour of the particles originating from the reactants stream is less straightforward. About 7% of these particles experience quenching, and this proportion does not increase with the blow-off rate. In fact, in this case the amount of particles crossing the lean flammability limit is marginal, due to their more central location that limits the interaction with co-flowing air. Interestingly, part of the particles that would either burn or are too rich to burn (groups e and d) in case Lr75-80, does not end up in group g (i.e., start and then stop burning) when the bulk speed increases. The decreased number of particles in these two groups is balanced by an increase of particles with mixed behaviour (group f). These particles cluster near the rich limit (similarly to group d) but can occasionally cross the flammability line. Nevertheless, relatively high values of mixture fraction SGS variance, low values of progress variable and/or small residence times within the flammability range prevent sustained combustion (similarly to group c). It is worth noting that particles of group e start burning due to diffusion of progress variable in the region of mixing with the pilot.

Groups c and g include any particle with value of mixture fraction within the flammability limits and not just stoichiometric values. Local extinctions occurring near the stoichiometric line are therefore a subset of these two groups and account to 1.5% of the total number of particles (from both pilot and main jet). Stoichiometric quenching does not show a sensible increase with blow-off rate, consistently with the observations made in Chen et al. (2020). This aspect will be further investigated in the next sections.

4.1.1 Lagrangian Behaviour of Fluid Particles at 70% Blow-off Rate

The behaviour of the Lagrangian particles is analysed for the Lr75-80 case for different instants of time in Figs. 5, 6, showing scatter plots of normalised filtered temperature, $T^* = (\widetilde{T}-T_\infty )/(T_a - T_\infty )$, reaction rate, $\dot{\omega }_c^* = (\overline{\dot{\omega }}_c/\overline{\rho })/(s_L/\delta _\mathrm {th})_\mathrm {st}$ (the subscript st indicating stoichiometric conditions), progress variable diffusion term ($\mathcal {D}_c = \nabla \cdot \left( \overline{\rho \mathcal {D}} \nabla \widetilde{c} \right)$, see Eq. (1)), $\mathcal {D}_c^* = (\mathcal {D}_c/\overline{\rho })/(s_L/\delta _\mathrm {th})_\mathrm {st}$, and particles location, $\left( x/D,r/D\right)$. From the Lagrangian perspective, burning particles can only extinguish due to diffusive processes, which in the studied configuration are caused either by regions of fast mixing with cold reactants (progress variable dropping at fixed mixture fraction) or with air and fuel (mixture fraction crossing the flammability limits). Particles of groups a and d of Fig. 4 never start burning either because diffusion processes are too slow (group d) or because they never mix with the fuel within the time window explored (group a). Their behaviour is thus irrelevant for the analysis in this work and not explored further. Although the behaviour of group f is intriguing, it also does not add significant insight for the Lagrangian analysis presented here, and thus it is not shown. The behaviour of particles from the pilot is qualitatively similar to that from the main jet and thus scatter plots are also not shown. Scatter plots for groups a to d and f can be found in the Supplementary Material.

Fig. 5

Scatter plots of normalised filtered temperature, reaction rate, progress variable diffusion term and particles location, for group e of Fig. 4 and for the Lr75-80 case at three instants of time from the particles injection. Colours represent the radial position particles had at time of the injection

Fig. 6

Scatter plots of normalised filtered temperature, reaction rate, progress variable diffusion term and particles location, for group g of Fig. 4 and for the Lr75-80 case at three instants of time from the particles injection. Colours represent the radial position particles had at time of the injection

The behaviour of particles originating from the main jet that start burning is shown in Fig. 5 (group e). These particles are prevalently above the rich limit at the time of injection and move within the flammability limits by diffusion with the pilot stream as observed for $t=0.4\,\hbox {ms}$ in Fig. 5. At following times these particles heat up due to diffusion processes (represented in this context by $\mathcal {D}^*_c$) with the surrounding flame front until they start burning ($t=1.4\,\hbox {ms}$). The relatively large scatter in radial positions observed at following times suggests strong turbulent mixing with entrained air. Part of the particles are further shifted towards the lean limit, suggesting that they may have crossed the flame front. However, the decrease of temperature (and progress variable) due to laminar and SGS diffusion processes is not fast enough to quench them ($\widetilde{c}<0.4$) and thus the reaction rate is always non-zero within the flammability limits. Note that particles from the pilot (group c – scatter plots available in the Supplementary Material), injected at stoichiometric conditions, have a similar behaviour to those of group e except for they exhibit an initial drop of progress variable and temperature from burnt condition, and an increase of the amount of fuel, due to diffusion processes in the region of mixing with the fresh reactants.

The behaviour of quenching particles is discussed next (groups c and g). Since the qualitative behaviour between groups c and g is similar, scatter plots are only shown for group g in Fig. 6, while those for group c can be found in the Supplementary Material. The initial behaviour of the particles in group g is similar to those in group e, at $t=0.8\,\hbox {ms}$, and consists of particles crossing the rich flammability limit and then burning. However, at following times these particles experience a strong decrease of temperature ($t=2.2\,\hbox {ms}$), which is mainly due to about 37% of the particles in this group quickly crossing the lean flammability limit and 46% crossing the rich flammability limit. This suggests the presence of large vortices that enhance mixing with both air and fresh reactants by rapidly casting the Lagrangian particles outwards in the coflow region or inwards towards the burner axis before combustion is complete. This behaviour is further supported by the large scatter in radial positions of the particles and the reduction in diffusion of progress variable observed for $t>2.2\,\hbox {ms}$. Only about 17% of the particles in group g, corresponding to the 1.5% of the total amount of particles from the main jet, experience quenching at near-stoichiometric conditions. Since local extinctions were described in Barlow et al. (2015), Cutcher et al. (2018) and Chen et al. (2020) as discontinuities in the OH concentration or progress variable field occurring at near-stoichiometric conditions, Lagrangian particles crossing these extinction holes are further investigated in order to shed light on the above behaviour. It is worth mentioning that, even if only a few of these particles are shown next in Figs 7, 8 for illustrative purposes, these are representative of typical reoccurring behaviours found for all particles quenching along the stoichiometric line.

Figure 7 (top) shows contours of progress variable on top of isolines of stoichiometric mixture fraction for different instants of time when extinction holes are observed. These are clearly marked in the figure by the $\widetilde{c}=0.4$ isolines. The trajectory of four characteristic Lagrangian particles crossing the stoichiometric line within the extinction hole at the time represented by the contours is also shown in the same figure, and coloured by normalised reaction rate. Particles crossing an extinction hole can be distinguished in two groups: those reacting before impinging on the stoichiometric line (Fig. 7a, b) and those not reacting (Fig. 7c, d). Particles of the first group stop burning after impinging on the extinction hole and thus also correspond to subsets of groups c and g in Table 1 where quenching occurs close to stoichiometric conditions. Furthermore, particles of both groups can either remain on the flame front (Fig. 7b, d) or move away towards radially outward locations (Fig. 7a, c). Note that in the particular case shown in the figure, the particles of the second group (Fig. 7c, d) are injected from different locations but impinge on the same extinction hole at about the same time. It is interesting to observe that a portion of stoichiometric line within an extinction hole is always found in a region of relatively intense strain. This is observed in Fig. 7 (bottom) plotting the Q-criterion at the same times of the progress variable contours. The Q-criterion is defined for the purposes of this analysis from the second invariant of the velocity gradient tensor as

$$\begin{aligned} Q = \frac{1}{2}\left( \left\Vert \mathbf {\varOmega }^2 \right\Vert - \left\Vert \mathbf {S}^2 \right\Vert - \nabla \cdot \mathbf {u} \right) \end{aligned}$$

(9)

where $\mathbf {\varOmega }$ and $\mathbf {S}$ are the vorticity and strain rate tensors, respectively (Hunt et al. 1988). The contribution of the dilatation term, $\nabla \cdot \mathbf {u}$, was observed to be well below 1% in the region of interest and thus was considered to be negligible. The Q field suggests that the passive flow tracers move from a region of intense vorticity right before crossing the extinction hole, to a region of high strain right after. This characteristics may play a role in the formation of the hole itself due to flame strain caused by the flame-vortex interaction, and will be further discussed in Sect. 5.

Fig. 7

Contours of filtered progress variable (top) and Q-criterion (bottom), and pathlines for four passive flow tracers crossing an extinction hole. The pathlines are coloured by normalised reaction rate. The location of the particles at the time of the contour (annotated in the box) is indicated with a red marker. The white isoline indicates the stoichiometric line, the black isoline indicates the $\widetilde{c}=0.4$ threshold. Level saturation is used to highlight the regions with positive and negative Q-criterion

Fig. 8

From top to bottom: time history of normalised filtered progress variable; its normalised diffusion term $\mathcal {D}_c^*$ and reaction rate $\dot{\omega }_c^*$; mixture fraction and fuel mass fraction; and mixture fraction SGS variance for the four particles of Fig. 7. The horizontal lines mark the values of $\widetilde{c}=0.4$ and stoichiometric mixture fraction. The red lines mark the time the particle first impinges on the stoichiometric line within the extinction hole. They grey area in the third row indicates the region outside the flammability limits

Figure 8 shows the time evolution of the four particles examined in Fig. 7. In order to eliminate the high-frequency turbulent content, the time history of each particle is smoothed with a Savitzky-Golay filter fitting a first-order polynomial over a window length of $\varDelta t = 0.12\,\hbox {ms}$. This window is still sufficiently smaller than the characteristic flame time scale ($\sim 1\,\hbox {ms}$) so as to reveal the temporal behaviour of the reacting particles. The behaviour of the particles that were burning before crossing the extinction hole is discussed first. For the case of Fig. 8a the particle originates from the edge of the main jet (low $\widetilde{\xi }$ value) and first experiences a fast diffusion of mixture fraction from the jet centre due to the strong inhomogeneity of the mixture at upstream locations, so that this quantity at $t\approx 0.8\,\hbox {ms}$ is rich, consistent with what observed for Fig. 6. Consequently, the fuel mass fraction also increases. However, the slower increase of fuel mass fraction as compared to mixture fraction suggests that part of the fuel is burning, which is due to the process at the start of combustion described in Fig. 5. The simultaneous increase observed in the progress variable is thus partly due to the burning fuel, and partly to diffusion with fluid from the pilot. For $1< t < 1.7\,\hbox {ms}$ the particle continues its trajectory moving radially outward, experiencing a decrease of mixture fraction as it approaches the stoichiometric line, and an increase in reaction intensity with peak $\dot{\omega }_c^* = 1.5$. At time $t\approx 1.9\,\hbox {ms}$ the particle impinges on the stoichiometric line within the extinction hole as observed in Fig. 7 and its value of progress variable quickly drops below the threshold $c^*=0.4$. Correspondingly, the reaction rate drops to zero, supporting the earlier argument that the extinction of the Lagrangian particles (subset of group g of Table 1) is linked with the passage through the extinction hole. The sudden drop of progress variable and reaction rate is associated to an increase of mixture fraction variance as shown in the bottom plot of Fig. 8, suggesting that SGS mixing processes between fuel and air are strong in this region. This is surprising, as values of $\sigma ^2_{\xi ,sgs}$ above $10^{-4}$ are mainly observed up to $x/D<10$ in the central unmixed stream of reactants, while the particle is quite off-centred and at $x/D\approx 13$ at this time. A local increase of mixture fraction and fuel is also observed in correspondence of the peak $\sigma ^2_{\xi ,sgs}$ value, suggesting that these peaks, and thus the prediction of the extinction, may be due to diffusion with cold and richer reactants present nearby. At following times the particle moves away from the flame front, as can be inferred by the value of mixture fraction becoming significantly leaner. It is worth noticing that diffusion processes due to mixing between burnt and unburnt fuel (term $\mathcal {D}_c^*$ in the figure) do not seem to play an important role near the extinction hole and thus the Lagrangian particle extinction is only associated to presence of high gradients of mixture fraction, consistently with past literature (Watson et al. 1999; Lyons et al. 2005), and its variance.

The behaviour in Fig. 8b is similar to that discussed for Fig. 8a. In this case the flow tracer originates from the pilot stream and therefore belongs to group c of Fig. 4. The progress variable carried by this particle first gradually decreases due to diffusion with cold reactants from the main jet, and its reaction rate remains zero up to $t \approx 1.5\,\hbox {ms}$, despite its trajectory following the stoichiometric line. This is because the particle from the pilot does not have enough fuel to burn yet. This particle then begins to burn, which is shown by the increase of reaction rate, and the increase in progress variable. Local extinction occurs at $t \approx 2.2\,\hbox {ms}$ with a sudden drop of progress variable and reaction rate on the stoichiometric line, which is again accompanied by a sudden increase of mixture fraction variance and correspondingly of fuel. Differently from the behaviour in Fig. 8a, at further times the particles remains on the stoichiometric line, thus on the flame front, and restarts burning. This extinction-reignition process occurs in less than half millisecond and during the same time the extinction hole has propagated downstream of a distance of about 2D. Note that, despite the Lagrangian extinction occurs at a slightly later time in respect to that observed in Fig. 8a, the extinction hole is located significantly more upstream in this case ($x/D<10$), which is due to the weaker momentum of particles from the pilot.

The characteristics behaviour of flow tracers passing through an extinction hole without previously burning is analysed in Fig. 8c, d. These particles are observed to originate near the jet axis and start heating up by diffusion as they enter the flame front from the rich side. However, differently from the behaviour observed for for Fig. 8a, these particles carry quite high values of mixture fraction SGS variance, which slows down reaction on the flame front. Moreover, they move radially outward relatively fast and reach the stoichiometric line at $t\approx 1.5\,\hbox {ms}$, earlier than what was observed for the particle in Fig. 8a. Since the values of $\sigma ^2_{\xi ,\text {sgs}}$ carried by these particles at the time they impinge on the stoichiometric line are similar to those in the peaks observed in Fig. 8a, b, it is possible that these particles are contributing to the formation of the extinction hole itself, while the flow tracers examined in Fig. 8a, b are quenched by it. This behaviour will be further analysed in the next section. Similarly to that discussed earlier, the Lagrangian particle at following times can remain near the stoichiometric line and start burning (Fig. 8c), or move outward and interact with co-flowing air (Fig. 8d).

The analysis in this section has revealed different scenarios for the quenching of the Lagrangian particles (groups c and g of Table 1). Most of these stop burning as they quickly move radially outward or inward towards the jet axis, thus crossing the lean or rich flammability limit respectively. A smaller part quenches at near stoichiometric condition, which is caused by the particles crossing a relatively cold region with high SGS mixing due to the formation of an extinction hole on the flame front. This in turn might be generated and sustained by other particles stemming from the centre of the jet and impinging on the stoichiometric line at relatively cold conditions and with relatively high values of mixture fraction variance. However, these characteristic behaviours only partly explain the dynamics of the formation of an extinction hole, and do not fully address the role of the large flow scales in the interaction between flame and oncoming turbulent flow of fresh reactants. Thus, an Eulerian analysis is carried out next to shed light on these extinction dynamics.

4.2 Eulerian Analysis of the Flow at 70% Blow-off Rate

The Lagrangian analysis highlights some of the complex dynamics in local extinctions, and shows that regions of relatively strong strain and high SGS variance of mixture fraction can play a role in the formation of an extinction hole. An analysis is then performed from an Eulerian point of view, to shed light on the above mechanisms using the control volume approach.

Fig. 9

Contours of progress variable, Q-criterion, mixture fraction SGS variance and normalised reaction rate for two consecutive steps before and after the formation of an extinction hole in the Lr75-80 flame. The black isolines mark the threshold $c^*=0.4$, the white isoline is the stoichiometric line. The annotations indicate the locations of vortex interaction with the flame (a), induced strain (b) and mixture fraction SGS variance dragged from the main jet onto the stoichiometric line (c)

Figure 9 shows the state of the flow right before and after the formation of a local extinction on the flame front. Vortical structures are first formed in the shear layer between main jet and pilot streams, and are then convected downstream. These vortices have thus a preferential orientation in the flame tangential direction, orthogonal to the flow direction, and tend to align along the flame front near the stoichiometric line before the extinction takes place (see circular arrows in Fig. 9), as also observed in past studies (Tuttle et al. 2013; Kariuki et al. 2015; Chowdhury and Cetegen 2018; Zhou et al. 2015; Wang et al. 2017). The initial interaction between the vortex and the flame front is characterised by relatively modest levels of strain and a local reduction of the turbulent flame thickness (Fig. 9a, points A and B). During this initial process the value of progress variable on the stoichiometric line is observed to decrease from $\widetilde{c}\approx 0.8$ to $\widetilde{c}\approx 0.6$ in the region indicated in the figure, and thus it remains above the threshold of 0.4. Moderate levels of reaction rate are still observed along the stoichiometric line at this time, but the front is continuous. A local peak of mixture fraction SGS variance is observed to cross the flame front simultaneously (point C). About $0.4\,\hbox {ms}$ later the same amount of mixture fraction SGS variance is observed along the stoichiometric line (point C in Fig. 9b), and an extinction hole is formed. This variance is observed to be convected from the inhomogeneous centre of the jet by large eddies, and it exacerbates the aforementioned effect of resolved strain by inducing a further decrease of the reaction rate by a factor 10 to 20. Consequently, the progress variable drops and the flame front breaks. This physical mechanism is consistent with the Lagrangian behaviour observed for Fig. 8c and d. However, it does not fully explain how the unstrained flamelets model used for this study is able to predict it. Since the reaction rate in the flamelets model can only vary due to a variation of its controlling parameters, $\widetilde{c}$, $\sigma ^2_{c\mathrm {,sgs}}$, $\widetilde{\xi }$ and $\sigma ^2_{\xi \mathrm {,sgs}}$, these are investigated next.

The SGS variance of progress variable, $\sigma ^2_{c\mathrm {,sgs}}$, is observed to be of order of $10^{-5}$ for the mesh used in this study, and does not change significantly near a local extinction. A dominant role of $\widetilde{\xi }$ in the local extinction mechanism is also to be excluded, as the extinction phenomenon is observed in this analysis along the stoichiometric line. The ability of the model to capture local extinctions therefore must depend on progress variable and mixture fraction variance, whose time variation is represented respectively by Eqs. 1 and 6. For convenience of discussion, these equations are rewritten as:

$$\begin{aligned}&\frac{\partial \overline{\rho } \widetilde{c}}{\partial t} = -\mathcal {C}_c + \mathcal {D}_c + \overline{\dot{\omega }}_c \end{aligned}$$

(10)

$$\begin{aligned}&\frac{\partial \overline{\rho }\sigma ^2_{\xi \mathrm {,sgs}}}{\partial t} = -\mathcal {C}_{\sigma \xi } + \mathcal {D}_{\sigma \xi } + \mathcal {P}_{\sigma \xi } - \mathcal {X}_{\sigma \xi } \end{aligned}$$

(11)

where $\mathcal {C}$, $\mathcal {D}$, $\mathcal {P}$ and $\mathcal {X}$ indicate convection, diffusion, turbulent production and dissipation, respectively. The behaviour of these terms along the stoichiometric line before and after the extinction is shown in Fig. 10 for the same characteristic extinction and the two time steps of Fig. 9, but the same qualitative behaviour is observed for any extinction hole on the flame front. The balance of the terms in Eq. 10 reveals that the initial decrease of progress variable from $\widetilde{c} \approx 0.8$ to $\widetilde{c} \approx 0.6$ is strongly correlated to an increase in magnitude of the convective term $\mathcal {C}_c$. Two peaks of reaction rate are observed aside the local minimum. Simultaneously, no significant variation of SGS mixture fraction variance is observed at the same location, although a peak is observed slightly downstream, also corresponding to a increase for the convection term $\mathcal {C}_{\sigma \xi }$. It is worth noting that the diffusion terms of $\widetilde{c}$ and $\sigma ^2_{c\mathrm {,sgs}}$ do not seem to play an important role during this process (not shown for the SGS variance).

Fig. 10

Spatial behaviour along the portion of stoichiometric line within an extinction hole of $\widetilde{c}$, $\sigma ^2_{c\mathrm {,sgs}}$, and the terms on the right-hand side in Eqs. 10 and 11 (normalised by $\delta _\mathrm {th})/(\overline{\rho } s_L$) for the two time steps of Fig. 9. The curvilinear distance $\ell$ is centred on the extinction

At the subsequent time the progress variable further drops to about $\widetilde{c}=0.3$, the reaction rate at the corresponding location decreases of about an order of magnitude, and the extinction hole is formed, as discussed for Fig. 9. This event is accompanied by a further increase in the magnitude of $\mathcal {C}_c$, and a simultaneous increase of SGS mixture fraction variance. Moreover, a sharp increase in the SDR of mixture fraction right ahead of the location of the progress variable minimum, and a strong increase of the convective term $\mathcal {C}_{\sigma \xi }$ immediately downstream, are also observed. Two specific points are worth noting here. First, presence of SDR peaks right upstream the extinction locations was also observed experimentally (Cutcher et al. 2018). Given the local peak of variance at the same position, it is expected that this SDR peak is a consequence of the quench (due to lower temperature) and not the cause, as also argued in Chen et al. (2020). Second, the peak values of SDR and turbulent production in Eq. (11) are significantly smaller than the peak values of the convective terms. One can recognise that a negative peak in the convective term implies a positive (or incoming) flux of variance by Gauss theorem. By looking at contours of $\mathcal {C}_c$ and $\mathcal {C}_{\sigma \xi }$ at the time of the extinction, one can notice an alignment of positive and negative convective flux respectively of progress variable and mixture fraction SGS variance along the stoichiometric line (indicated with A in Fig. 11). This indicates an outgoing flux of $\widetilde{c}$ and an incoming flux of $\sigma ^2_{\xi \mathrm {,sgs}}$ at this location, suggesting that the formation of the extinction hole is associated to mixture of cold reactants with relatively high mixture fraction variance convected onto the stoichiometric line. This mixture must be travelling from regions near the jet centre and not along the stoichiometric line, consistently with the observations made for the Lagrangian analysis. On the other hand, this alignment of convective terms with the stoichiometric line is also observed at locations where extinctions do not occur. An extinction hole is only observed when this alignment is complemented by relatively high values of resolved strain, as can be seen in Fig. 11 using the Q-criterion, and is related to the time scales involved during this process. This aspect will be further discussed in Sect. 5.

Fig. 11

Contours of Q-criterion and convection of c and $\sigma ^2_{\xi \mathrm {,sgs}}$ at $t=2.04\,\hbox {ms}$, right before an extinction hole is formed. White markers indicate the eddies location, black arrows their orientation. Dashed circles indicate regions of alignment between the peaks of the two convective terms. White and black isolines mark respectively stoichiometric line and $\widetilde{c}=0.4$. Level saturation is used to highlight positive and negative regions

4.3 Analysis of the Flow at 90% Blow-off Rate

The unstrained flamelets model with the mechanism described in the previous section was observed to underestimate the amount of local extinctions as compared to experiments for the 90 % blow-off speed configuration in Table 1 (Chen et al. 2020). The Lagrangian and Eulerian analyses performed for the 70% blow-off case are then repeated here to shed light on possible phenomenological differences in the process of local extinction between the two cases.

From the Lagrangian point of view, the temporal evolution of groups a to g of Fig. 4 is qualitatively analogous to the one observed for the Lr75-80 flame (Figs. 5 and 6). The representative scatter plots for group g only at three instants of time are shown in Fig. 12 for comparison with Fig. 6. The evolution of the flow tracers is qualitatively identical to that described for the Lr75-80 case, except for the faster time scales (note the different times chosen for the plots) due to the increased bulk jet speeds: particles first enter the flammability limits from the rich side due to strong diffusion with air and start to burn ($t = 0.63\,\hbox {ms}$). About $1\,\hbox {ms}$ later most of these particles reach the lean flammability limit ($t=1.71\,\hbox {ms}$) and stop burning after crossing it, indicating that they are moving radially outward. The increased bulk speed leads to stronger mixing in the region of the jet/pilot interface and thus stronger diffusion of temperature (represented in this context by the diffusion of progress variable) and earlier combustion of the particles as compared to Fig. 6. Nevertheless, the amount of Lagrangian particles exhibiting quenching is similar to that observed for case Lr75-80 (see also Fig. 4), although a significant increase of mixed behaviour (group f) is observed, as also shown in Fig. 4. A similar qualitative behaviour is also observed for particles in group c (not shown), except for the increase in the amount of particles experiencing quenching due to crossing the lean flammability limit as discussed for Fig. 4.

Fig. 12

Scatter plots of normalised filtered temperature, reaction rate, progress variable diffusion term and particles location, for group g of Fig. 4 and for the Lr75-103 case at three instants of time following the particles injection. Colours represent the radial position particles had at time of the injection

Fig. 13

Contours of progress variable, Q-criterion and mixture fraction variance for an extinction in the Lr75-103 flame. Lines and annotations as in Fig. 9

Fig. 14

Balance of the terms in Eqs. 10 and  11. Plotted along the stoichiometric line for the areas in the dashed squares and at the time steps of Fig. 13. Legend as in Fig. 10

From the Eulerian point of view, the dynamics leading to the formation of an extinction hole on the flame front is also similar to that described for the Lr75-80 case. The characteristic alignment of a vortical structure with the stoichiometric line is shown for a time immediately after the rupture of the flame front in Fig. 13. Compared to the Lr75-80 case, a significant increase of vortex intensity is observed near the stoichiometric line, amounting to an average 30% increase in mean vorticity magnitude and 50% increase in peak vorticity when a hole is formed. This increase is a direct consequence of the increased Reynolds number and higher radial velocity gradients.

According to the mechanism explained for the Lr75-80 case, vortical structures bring fresh reactants onto the flame front and thus the increased vorticity is expected to correspond to higher convective fluxes. This is observed in Fig. 14, showing the spatial variation, along the portion of the stoichiometric line within an extinction hole, of the terms in Eqs. (10) and (11) (normalised) at the time the extinction hole is formed. The progress variable within the extinction hole falls below the $\widetilde{c}=0.4$ threshold as expected. Negative peaks of convective flux of progress variable $\mathcal {C}_c$ (fresh fuel coming in) are on average 30% higher than in the 70% blow-off case. On the other hand, the peak levels of variance which are cast onto the stoichiometric line are now decreased by about one third, which can also be observed in the contour of Fig. 13. The behaviour of these terms highlights two main differences with the 70% blow-off case. First, a well-defined peak in $\sigma ^2_{\xi ,\mathrm {sgs}}$ is not observed anymore at the time the extinction hole is formed, and seems to be compensated by the stronger influx of fresh reactants. Second, turbulence production and dissipation terms are now observed to be predominant as compared to the convection of variance (see Fig. 14). This increase in turbulent production and SDR, which is observed over the entire flame, is a direct consequence of the higher radial gradients of mixture fraction and its variance (not shown). As for the Lr75-80 case, local peaks of SDR are observed within the extinction hole, supporting the argument that they are a consequence of the vanishing reaction rate following the formation of the extinction.

5 Discussion

5.1 Extinction Mechanism at 70% Blow-off Rate

The Lagrangian and Eulerian analyses have shed some light on the mechanism of formation of flame front breaks/local extinctions for the Sydney/Sandia flame with inhomogeneous inlets. The analysis of the scatter plots for the flame at 70% blow-off rate indicates that combustion is predominantly rich in the stratified premixed zone ($x/D<5$), but reaction takes place across the whole flammable range further downstream ($x/D>8$), as can be observed from Figs 5 and 6. The inspection in particular of the flow tracers in groups c and g of Fig. 4 reveals a correspondence between a subset of these particles and the appearance of extinctions holes along the flame front in the contours of $\widetilde{c}$ (Figs. 8 and 7). These extinctions form due to the impingement of particles carrying cold reactants and high values of mixture fraction variance from the main jet onto the stoichiometric surface. The Eulerian analysis indicates that this motion of fluid to and from the stoichiometric line is driven by the dynamics of the large eddies within the shear layer, as shown in Fig. 11 and consistently with past studies (Watson et al. 1999; Stroomer et al. 1999; Lyons et al. 2005; Baudoin et al. 2013). These vortices are convected downstream by the bulk flow and occasionally align with the stoichiometric line as shown schematically in Fig. 15. This alignment produces convective fluxes in and out, as indicated in the sketch. In particular, the flux near the stoichiometric line destabilises the flame as it brings fresh reactants ($\mathcal {C}_c>0$) and high mixture fraction variance ($\mathcal {C}_{\sigma \xi }<0$), which both weaken the reaction rate. However, as the same vortex moves downstream and away from the stoichiometric line, the inflow of cold reactants terminates and the flame is able to restart burning due to high reactivity of the surrounding fluid. From the Lagrangian perspective, when the cold particles from the main jet move towards outward radial locations, they are preheated due to diffusion with the surrounding hot fluid, which eventually leads to combustion. Reaction onset happens with a delay, which is comparable to the chemical time scale of the flamelet $t_c = \delta _\mathrm {th}/s_L=1.2\,\hbox {ms}$ at stoichiometric conditions, and to the delay observed in Fig. 5, where the Lagrangian particles start reacting at $t \approx 1.4\,\hbox {ms}$ after a gradual increase of progress variable. This time scale is not significantly smaller than the mean flow-through timescale based on 30 diameters and the bulk flow velocity, which is $t_{tf}=30D/U_b=2.8\,\hbox {ms}$. Although most particles in Fig. 5 start burning before reaching the stoichiometric line, some of them are accelerated by the large eddies so that they impinge on the stoichiometric surface before combustion can be initiated. This corresponds to the characteristic behaviour observed for Fig. 8c, d, showing that the preheating time (roughly estimated as the time needed by the particle to either cross the stoichiometric line or achieve $\widetilde{c}=0.4$ from the beginning of the of sharp increase in $\widetilde{c}$) for these particles is only about $0.5\,\hbox {ms}$, which is significantly smaller than the typical flame time scale $t_c$.

Fig. 15

Schematic of the effect of turbulent eddies on convective terms of progress variable and mixture fraction variance in the generation of extinction holes. The legend is as in Fig. 9

The analyses conducted in the previous sections further show that the extinction of Lagrangian particles does not contribute to the formation of the extinction holes on the flame front, but is rather a consequence of diffusion processes when these particles cross the extinction hole itself (see Fig. 8a, b). This occurs for only 1.5% of the particles, while most of the Lagrangian particles (remaining 23.5% in group c and 4.5% in group g of Fig. 4) stop burning as they rapidly cross the flammability limits. The number of Lagrangian extinctions is evidently too small to justify the significant amount of extinctions observed on the flame front both experimentally (Barlow et al. 2015; Cutcher et al. 2018) and numerically (Chen et al. 2020). This is indeed predictable for a model that does not take into account the effects of subgrid strain and non-premixed combustion mode in the look-up table. Nonetheless, good mesh resolution ensures that the LES is able to resolve most of the strain and capture its large-scale effects on the flame, as will be discussed next.

5.2 Role of Resolved Strain

Resolved strain plays an important role in terms of time scales required for the formation of an extinction hole on the flame front. In Fig. 7 it has been shown that the flow within an extinction hole is dominated by strain, highlighted by the negative values of Q-criterion. The Eulerian analysis further indicates that this resolved strain can curve the flame front and influence the way it interacts with the cold reactants. In non-extinguishing circumstances, the fraction of time spent by a vortex approaching the stoichiometric line, destabilising the flame by the effect of the convective fluxes ($\mathcal {C}_c>0$ and $\mathcal {C}_{\sigma \xi }<0$), is not sufficient and the flame is immediately re-energised by the change of sign of the fluxes ($\mathcal {C}_c>0$ and $\mathcal {C}_{\sigma \xi }<0$) due to diffusion and flame propagation discussed in Sect. 5.1, although the flame front is observed to become thinner. In this condition, the amount of fresh reactants entering the flame region results in a cooling of the flame front due to diffusion of heat, but this process is slower than the time required by the cold particles to start burning. On the other hand, the large eddies can occasionally stretch the flame in such a way that a larger portion of the stoichiometric line is aligned to a region with prevailing destabilising fluxes, as the event marked with A in Fig. 11. Under this circumstance, a higher volume of cold reactants can impinge on the stoichiometric line, thus offsetting the preheating effect and leading to the formation of an extinction hole. This mechanism, also observed experimentally in diffusion jet flames (Lyons et al. 2005), is of primary importance for the prediction of extinctions with LES, where large scale interactions between strain and flame are directly computed. This is also backed by recent works arguing that most influential flame strain effects are at resolved level (Filatyev et al. 2005; Langella and Swaminathan 2016; Doan et al. 2017). Nonetheless, possible SGS strain effects in the specific configuration investigated here cannot be excluded with certainty, and this can be quantified only using DNS, which is beyond the reach of the present analysis.

The mechanism of alignment between the vortical structures and the stoichiometric line is driven by the strain due to the flame-turbulence interaction as described next. Strain is related to vorticity via its transport equation:

$$\begin{aligned} \frac{\mathrm {D}\omega _i}{\mathrm {D}t} = \omega _j S_{ij} - \omega _i S_{kk} + \frac{\epsilon _{ijk}}{\rho ^2}\frac{\partial \rho }{\partial x_i} \frac{\partial p}{\partial x_k} - \frac{\epsilon _{ijk}}{\rho ^2}\frac{\partial \rho }{\partial x_j}\frac{\partial \tau _{kl}}{\partial x_l} + \frac{\epsilon _{ijk}}{\rho }\frac{\partial ^2 \tau _{kl}}{\partial x_j \partial x_l} \end{aligned}$$

(12)

where $\omega _i$ indicates the vorticity component, $S_{ij}$ is the strain component and $\epsilon _{ijk}$ is the Levi-Civita permutation parameter (Langella et al. 2018). The effects of density and viscous stress gradients (fourth term on the RHS), as well as the combined molecular diffusion and dissipation term (fifth term on the RHS), were observed to be of second-order importance in the current flame and thus can be omitted for simplicity following previous works (Langella et al. 2018; Hamlington et al. 2011). The latter in particular was observed to become more relevant in regions of less intense flame dilatation; however, this does not lead the generalities of the discussion presented below. Equation (12) identifies the mechanism in which vorticity not aligned to the flame normal (second term on the RHS) is suppressed by the flame due to the dilation effect proportional to $S_{kk}$. Since most vorticity in the cases under investigation is tangential to the flame (see for example the vortex sketch in Fig. 15), this is suppressed, as can also be observed in the contours of Q-criterion of Fig. 9 before and after the extinction hole is formed. The highly anisotropic eddies observed at $5<x/D<30$ thus tend to disappear and exhibit no preferential alignment at downstream locations. This vorticity suppression mechanism competes with the process of formation of a local extinction. It is observed that on average about $0.4\,\hbox {ms}$ for the Lr75-80 case is needed by a large eddy aligned to the stoichiometric line to create a hole within the flame front ($0.3\,\hbox {ms}$ for the Lr75-103 case due to the higher vorticity and stronger convective fluxes), estimated as the time the progress variable decreases from $\widetilde{c}=0.8$ to the threshold value of 0.4 (see Fig. 8). This time is evidently not sufficient for vortex suppression to occur before it can trigger a flame front break. One can make an order of magnitude estimation of the suppression time by solving Eq. (12) for a stationary vortex and only retaining the second term on the RHS, as it is of leading importance. Assuming $S_{kk}$ is constant during the process, the solution is exponential as $\omega _i(t) = \omega _0 e^{-S_{kk}t}$, where $\omega _0$ is the planar vorticity before the suppression action. Assuming $S_{kk}\approx (\varDelta u_n/\delta _{th})_{st}$ ($\varDelta u_n$ being the velocity variation in the direction normal to the flame) and using the flamelet’s continuity and state equations one can easily write $S_{kk} \approx (\tau s_L/\delta _{th})_{st}$. Thus, 90% of planar vorticity would get suppressed in about $0.37\,\hbox {ms}$, and the 99% of it in about $0.74\,\hbox {ms}$. Despite the errors introduced by the approximations, these timescales are similar to those observed directly from the LES results in Sect. 4. Most importantly they are slower than – or at most comparable to – those observed for the formation of the extinction hole (0.3 to 0.4 ms), suggesting that the balance between these times scales indeed plays in favour of a vortex causing an extinction before it gets suppressed.

5.3 Role of Mixture Fraction SGS Variance

A crude approximation of the isolated effect of $\sigma ^2_{\xi \mathrm {,sgs}}$ on the generation of a local extinction is also provided to quantify the importance of this term in the process of formation of an extinction hole in the flame front. At a given time step n, the transport equation for $\widetilde{c}$ can be discretised as:

$$\begin{aligned} \widetilde{c}^{n+1} \approx \widetilde{c}^{n} + \varDelta t \left( -\mathcal {C}_c^n+\mathcal {D}_c^n + \dot{w}_c^n\right) \end{aligned}$$

(13)

where $\dot{w}_c^n=\overline{\dot{\omega }}^n/\overline{\rho }^n$ and convection and diffusion terms are already divided by density. If we assume no extinction has taken place and that there is an instantaneous balance between convection, diffusion and reaction rate, $\widetilde{c}$ does not change after a small $\varDelta t$. This gives $c^{n+1}\approx c^n$ and $-\mathcal {C}_c^n+\mathcal {D}_c^n \approx - \overline{\dot{w}}_c^n$. One can then estimate how much the sudden change in mixture fraction variance observed during the extinction process contributes to the decrease of reaction rate. Considering the characteristic behaviour of Fig. 10, $\sigma ^2_{\xi \mathrm {,sgs}}$ increases on average from $2\times 10^{-4}$ at $\widetilde{c}=0.6$ to $3\times 10^{-4}$ at $\widetilde{c}=0.4$. By accessing the stoichiometric flamelet in the precomputed table with a value of $\widetilde{c}=0.6$ and $\sigma ^2_{c\mathrm {,sgs}}=0$, it is observed that the above variation causes a 30% drop of reaction rate. Assuming that this is a sudden change, and that convection and diffusion do not have time to adjust to the new level of reaction rate, one can write $\dot{w}_c^{n+1} \approx 0.7\dot{w}_c^{n}$ and $-\mathcal {C}_c^{n+1}+\mathcal {D}_c^{n+1} \approx -\dot{w}_c^{n}$. Equation (13) becomes:

$$\begin{aligned} \widetilde{c}^{n+2} \approx \widetilde{c}^{n} + \varDelta t \left( -\dot{\omega }_c^{n} + 0.7\dot{w}_c^{n} \right) = \widetilde{c}^{n} - 0.3\, \varDelta t \;\dot{\omega }_c^{n} \end{aligned}$$

(14)

If the variance was the sole cause for a change in c of about $\varDelta c = 0.33$ (similar to that observed for Fig. 10), the action time of $\sigma ^2_{\xi \mathrm {,sgs}}$ in these “quasi-frozen” conditions would be:

$$\begin{aligned} \varDelta t \approx \left| - \frac{\varDelta \widetilde{c} }{0.3\, \dot{w}_c^{n}}\right| \approx 1\,\hbox {ms} \end{aligned}$$

(15)

The analysis of the extinction over different snapshots performed in Sect. 4.2 and Fig. 10 shows that the estimated value is about 2.5 times the actual time needed for c to drop from 0.63 to 0.3, and thus this effect would not be able to cause local extinction alone. Although this time scale is considered as independent of strain effect in the above approximation, it suggests that role of $\sigma ^2_{\xi \mathrm {,sgs}}$ in the local extinction is of the same order of the vortex-induced incoming flux of cold reactants acting as a strongly coupled effect. This result is indeed interesting and shows that the formation of flame holes predicted by LES with unstrained flamelets is not solely a cold mixing effect due to the resolved strain, but relies on a simultaneous and equally important contribution due to SGS unmixedness, modelled through mixture fraction variance.

5.4 Extinction Mechanism at 90% Blow-off Rate

Unlike in the experimental observations (Barlow et al. 2015; Cutcher et al. 2018), the amount of flame extinctions does not increase in the numerical simulation for the 90% blow-off case, as also observed in Chen et al. (2020), despite the stronger destabilising convective fluxes observed in Fig. 14, due to the stronger flame-tangential vorticity. This is due to the fact that the stronger destabilising flux is balanced out by significantly lower peaks of mixture fraction SGS variance observed to be brought onto the flame front, which is due to higher SDR values observed on the flame front at the higher speed condition. Thus, the extinction at 90% blow-off rate must be mostly controlled by the strain mechanism discussed in Sect. 5.2. As also argued in Chen et al. (2020), it is possible that near blow-off, the precomputed thermochemistry (not considering subgrid strain and cross premixed/diffusion mode) may be insufficient, which is difficult to verify given the modelling uncertainty and limitations for strain effect and multi-mode combustion at subgrid level. The analysis conducted here further indicates that in this case the limitation of the unstrained flamelets model would lie in the over-prediction of the flame-tangential vorticity suppression process ($S_{kk}$ in Eq. (12) does not decrease with strain at flamelet level), which either affects the time scale needed by the vortices to be effective, or inhibits their ability to bring high variance from the centre of the jet (e.g. by rapidly changing the vortices orientation). This gives indication for modelling of flame extinctions in future studies.

6 Conclusions

In this work a mechanism that leads to formation of extinction holes on the flame front in LES with unstrained flamelets and presumed PDF closure is investigated, and its physical implications are discussed in the context of piloted jet flames with inhomogeneous inlets. In particular two flames are examined, respectively having 70% and 90% of the measured blow-off speed. To the best knowledge of the authors, this is the first study aiming to explain how the coupled effect of resolved large flow structures and modelled subgrid mixing can predict local extinction-like flame breaks in LES with flamelet combustion modelling.

Eulerian and Lagrangian analyses have identified the predominant role of the local convective fluxes of progress variable and SGS mixture fraction variance, corresponding to cold particles from the inhomogeneous jet centre that impinge on the flame front carried by large eddies. The shear-induced eddies have been observed to align to the stoichiometric mixture fraction for both flames and thus induce strain that weakens the flame front, as also observed in many previous studies. For the 70% blow-off case the eddies are observed to carry relatively high values of mixture fraction variance from the jet core onto the stoichiometric line. This subgrid-scale mixing effect, coupled with the resolved-scale strain, yields an additional one-third drop in reaction rate leading to the formation of local extinction holes. As the strain effect is moderate and mainly at resolved level, the unstrained flamelets model can predict the extinctions with good accuracy. In the 90% blow-off case, the strain effects are stronger but the effect of the variance in breaking the flame front is not observed. The combustion model in this case can still predict the extinction due to the impingement of cold reactants on the flame front. However, the amount of local extinctions is not observed to increase as compared to the lower speed case. It is argued that the flame-tangential vorticity is too-quickly suppressed by the heat release, and SGS strain modelling is in this case relevant to reduce this effect.

From the Lagrangian perspective, about 4.5% of the particles from the main jet experience quench due to strong mixture fraction gradients as they move radially inwards or outwards (crossing the flammability limits), while only about 1.5% of the particles quench at near-stoichiometric conditions. Similarly, 24 to 55% of the particles from the pilot experience quenching due to interaction with either co-flowing air or fuel-rich gases, while still only 1.5% of the quenching events occur along the stoichiometric line. The analysis conducted in this work suggests that the Lagrangian extinction (particle quenching along fluid path) does not lead to the breakage of the flame front, but rather these particles quench when occasionally crossing the local extinction sites.