2.1. Vortex breakdown

The vortex breakdown phenomenon occurs when the flow swirl velocity is sufficiently high to create an adverse pressure gradient in the axial direction while centrifuging the flow outward and to form a stagnation point at the centre of the swirling axis (Brown & Lopez Reference Brown and Lopez1990; Lopez Reference Lopez1990; Lucca-Negro & O’Doherty Reference Lucca-Negro and O’Doherty2001). This forces a sudden radial expansion of the swirling vortex, creating a depression at its centre. As the swirl velocity decays axially, so does the pressure depression at the centre of the swirling vortex, inducing a reverse axial flow and creating a CRZ. This process is described and illustrated by Syred (Reference Syred2006). The formation of a CRZ due to vortex breakdown was classified by Harvey (Reference Harvey1962) as the transition between weakly and strongly swirled flows; however, Harvey (Reference Harvey1962) reported different vortex breakdown modes depending on the level of swirl. Sarpkaya (Reference Sarpkaya1971) first identifies three types of vortex breakdown while performing experiments in laminar flows, for which the ratio of tangential to axial velocity was progressively increased at a fixed flow rate Reynolds number $({\textit{Re}}) \approx 5000$ : the first mode to appear for the lowest swirl levels is the spiral vortex breakdown, where the swirling vortex creates a helical structure. Then, if the swirl level is further increased, the spiral mode transitions to an (axisymmetric) BVB mode, creating a bubble-like CRZ. Finally, the double-helix vortex breakdown is a special mode that is observed when very slow ( $1/60^\circ /{\textrm {s}}$ ) changes in swirl velocities are applied, creating a double-helix structure. If the swirl level is further increased or if the system is disturbed, the flow transitions to the BVB mode. The last mode was observed much later when Sarpkaya (Reference Sarpkaya1995) performed new experiments to study turbulent vortex breakdowns. At ${\textit{Re}} = 200\,000$ and circulation number $\varOmega = \varGamma /U_0 D_0 = 0.50$ (where $\varGamma$ is the circulation of the vortex, $U_0$ is the cross-sectionally averaged axial velocity and $D_0$ is the tube diameter), Sarpkaya (Reference Sarpkaya1995) observed the BVB mode develop to a vortex breakdown with a fully conical shape, with linear growth from the stagnation point to the end of the wake, naming it the CVB mode. Sarpkaya (Reference Sarpkaya1995) then stated that the CVB should be considered as the fourth fundamental type of vortex breakdown, in addition to the spiral, bubble and double-helix types.

In a later study, Billant, Chomaz & Huerre (Reference Billant, Chomaz and Huerre1998) also observed the appearance of BVB and CVB in a non-confined swirling flow. They gave further insights into the flow around and inside the recirculation zones of the two vortex breakdown modes. As reported by other researchers, the BVB mode bubble exhibits a diverging flow that is followed by a contraction, resulting in the ovoid shape of the CRZ and in the creation of a toroidal vortex ring that recirculates the fluid. The flow inside the bubble is well organised and strong. In the CVB mode, which Billant et al. (Reference Billant, Chomaz and Huerre1998) were the first to characterise in detail, the flow expands from the stagnation point, but there is no contraction further downstream, as the flow recirculates inside the cone, also creating a toroidal vortex ring at the base of the cone. The conical sheet that separates the swirling jet from the CRZ is unstable and begins to distort and roll up into vortices before breaking into turbulence further downstream. The authors state that these instabilities might come from streamwise shear and/or centrifugal instabilities. From their observations of the flow in the CVB mode, Billant et al. (Reference Billant, Chomaz and Huerre1998) argue that the CVB observed by them differs from the one observed by Sarpkaya (Reference Sarpkaya1995), as the conical state reported in the latter is ‘in fact a bubble whose turbulent conical wake has come so close to the stagnation point that the whole structure looks like a cone (…) appearing as the natural evolution of the wake of a bubble towards a turbulent state as Re is increased’. Thus, since the conical state reported by Sarpkaya (Reference Sarpkaya1995) does not display any of the characteristics observed by Billant et al. (Reference Billant, Chomaz and Huerre1998) (such as the important quasi-steady, organised flow inside the cone), Billant et al. (Reference Billant, Chomaz and Huerre1998) state that the conical mode represented in their study should be regarded as a reference CVB.

2.2. The BVB-CVB vortex breakdown mode transitions

Despite the different modes being well characterised, no single mechanism fully describes all the transitions between the BVB and the CVB modes found in the literature. Several works report this transition in various configurations, discussing multiple mechanisms; however, in all cases, the swirl level and the relationship between axial and swirl velocity are key parameters for the onset of the two vortex breakdown modes, with the triggering parameter varying in each case.

Sarpkaya (Reference Sarpkaya1995), who first described the BVB-CVB mode transition, observed the change of mode increasing the Reynolds number from ${\textit{Re}} = 10\,000$ up to ${\textit{Re}} = 250\,000$ , while the circulation number was kept around $\varOmega \approx 0.6$ in a water tank. As the Reynolds number increases, the wake behind the bubble moves upstream until it merges with the bubble, resulting in a full CVB. However, as aforementioned, Billant et al. (Reference Billant, Chomaz and Huerre1998) questions the observation of a CVB by Sarpkaya (Reference Sarpkaya1995), thus the Reynolds number as a cause for the transition from the BVB to the CVB must also be regarded with caution.

Indeed, in opposition to Sarpkaya (Reference Sarpkaya1995), to trigger vortex breakdown mode transitions Billant et al. (Reference Billant, Chomaz and Huerre1998) progressively increased the swirl level (characterised by the swirl parameter, $S_{pa} = 2 V_{\theta } (R/2,x_0) / V_x (0,x_0)$ , which puts in relation the tangential velocity $V_\theta$ at the centre of the jet ( $V_{\theta } (R/2,x_0)$ ) and the axial velocity $V_x$ , both measured at the axial position of the initiation of vortex breakdown, where $R$ is the radius of the nozzle exit), but keeping the Reynolds number constant. Starting from a fully axial jet, the swirl parameter was progressively increased, causing the jet to become distorted, then helical and, eventually, an axisymmetric BVB was formed, which could be stable or later slowly evolve into an axisymmetric CVB. Thus, Billant et al. (Reference Billant, Chomaz and Huerre1998) observed that the two vortex breakdown modes, bubble and conical, coexisted for the same values of Reynolds number (in the range $300 \lt {\textit{Re}} \lt 1200$ ) and swirl parameter, creating bi-stable states. Even if the data available were not enough to propose a mechanism explaining why the change in vortex breakdown mode would sometimes be triggered, Billant et al. (Reference Billant, Chomaz and Huerre1998) found that temperature perturbations of the order of $0.1\,^\circ\rm C$ had some influence in defining which vortex breakdown mode would appear. The authors supposed that buoyancy effects could be at the origin of this influence, with higher jet temperatures (compared with the tank temperatures) favouring the appearance of CVB.

Terhaar et al. (Reference Terhaar, Reichel, Schrödinger, Rukes, Paschereit and Oberleithner2015), studying the addition of an axial air stream at the centre of the swirler as a way to mitigate flashback risks in swirling combustors, observed that increasing the axial flow rate triggered a change in vortex breakdown mode for isothermal flows. For a constant high swirl number ( $S = 0.9$ ), when the axial injection rate is increased to over $10\,\%$ of the total air injection rate, the CRZ expands and the vortex breakdown mode changes from the bubble to the conical mode. After a parametric study, Terhaar et al. (Reference Terhaar, Reichel, Schrödinger, Rukes, Paschereit and Oberleithner2015) derived a model and mapped the triggering of the mode change for the isothermal flow as a function of the level of swirl and axial flow injection rates. While the authors clearly state that the onset of vortex breakdown is controlled by the swirl number, the resulting vortex breakdown mode (BVB or CVB) is not a sole consequence of the swirl number, both appearing for various swirl numbers as for Billant et al. (Reference Billant, Chomaz and Huerre1998). However, unlike Billant et al. (Reference Billant, Chomaz and Huerre1998), who found that small temperature differences triggered the transitions, for Terhaar et al. (Reference Terhaar, Reichel, Schrödinger, Rukes, Paschereit and Oberleithner2015), the shape of the inflow velocity profile seemed to determine the selection of the vortex breakdown type. High primary swirl number at high rates of axial injection promoted CVB, whereas the same swirl number without axial injection resulted in a BVB. Even if reacting conditions were also analysed, when a methane lifted flame was present, no change in vortex breakdown mode was observed with the increase in axial flow, the vortex breakdown mode remaining the bubble one.

More recently, Santhosh, Miglani & Basu (Reference Santhosh, Miglani and Basu2014) and Rajamanickam & Basu (Reference Rajamanickam and Basu2018) observed the formation of a CVB in an isothermal flow with two coaxial jets (an internal central axial jet and an external tangential one) while varying the swirl number for various Reynolds numbers. Vortex breakdown forms when the swirl number attains a certain value, but only in cases where the tangential jet dominates the flow. In these two works, this regime is characterised by a modified Rossby number ${\textit{Ro}}_m = |\Delta u_y|/u_{\theta } \lt 1$ , where $\Delta u_y$ is the difference of the bulk velocity based on the mass flow rate and tube surface section between the two jets and $u_{\theta }$ is the tangential velocity of the tangential jet. As in Billant et al. (Reference Billant, Chomaz and Huerre1998), first a BVB is observed, but this time, when the swirl number is further increased (keeping ${\textit{Ro}}_m \lt 1$ ), the vortex breakdown mode changes into a CVB, the change being clearly triggered by the increase in swirl number. However, increasing the swirl further makes the vortex breakdown mode re-transition into a bubble mode. Thus, Rajamanickam & Basu (Reference Rajamanickam and Basu2018) stresses that CVB exists only for a very narrow band of swirl number $S_G$ . Additionally, they note that CVB was observed for vane swirlers, but only with coaxial jets. A similar statement is made by Terhaar et al. (Reference Terhaar, Reichel, Schrödinger, Rukes, Paschereit and Oberleithner2015), where the authors claim the CVB mode could be linked with a simultaneous high level of swirl and axial flows in the jet, creating a specific and propitious inlet velocity profile that allows the stabilisation of the CVB mode. Comparing the results, it is possible to infer that the CVB observed by Rajamanickam & Basu (Reference Rajamanickam and Basu2018) is a similar state as the one observed by Billant et al. (Reference Billant, Chomaz and Huerre1998), characterised by the quick opening of the swirling jet and by a slow, organised and quasi-steady flow inside the CRZ, as illustrated by the streamlines for the $S_G = 2.53$ CVB case.

In a subsequent group of studies, Moise & Mathew (Reference Moise and Mathew2019) used direct numerical simulations (DNS) to study the effect of increasing swirl number on the topology of laminar swirling flows, confirming past experimental observations and the transition from spiral vortex breakdown to BVB, and finally to CVB. Afterwards, Moise (Reference Moise2020) discussed the co-existence of both BVB and CVB for the same swirl number, highlighting the importance of hysteresis effects on swirling flow topologies. Following these studies, Moise & Mathew (Reference Moise and Mathew2021) used LES to investigate the impact of turbulence on vortex breakdown modes. They observed a similar progression in the turbulent regime – BVB to CVB with increasing swirl – but additionally reported the emergence of a wide-open CVB structure, likely associated with a Coandă-type attachment of the flow to the wall. Moise & Mathew (Reference Moise and Mathew2021) thoroughly characterised each of the modes, showing the flow fields and discussing the exclusive characteristics of each of the modes.

In parallel, Douglas et al. (Reference Douglas, Emerson and Lieuwen2022) used space-state global bifurcation analysis (in a different approach to linear stability analysis of Terhaar et al. Reference Terhaar, Reichel, Schrödinger, Rukes, Paschereit and Oberleithner2015; Dagan et al. Reference Dagan, Arad and Tambour2016) to provide a stability map of different swirling flow modes of a laminar annular jet as a function of swirl number, Reynolds number and centre-body diameter. They identified the three modes: (i) wake, a pre-breakdown flow; (ii) vortex breakdown, a BVB-like isolated central recirculation bubble; and (iii) wall jet, which can be divided into two regimes: a CVB-like regime where the recirculation bubble expands radially due to the influence of wall-jet entrainment, but with the jet still separated from the wall, and a wide-open CVB-like configuration (as per the Moise & Mathew Reference Moise and Mathew2021 definition) where the jet attaches along the outer wall due to a full Coandă effect. These transitions were characterised through a sequence of bifurcations, including symmetry-breaking Hopf and saddle-node bifurcations, along with pronounced hysteresis loops. They retrieve the classic onset of BVB as swirl increases enough to generate a stagnation point and reverse flow on the axis; however, Douglas et al. (Reference Douglas, Emerson and Lieuwen2022) show that the vortex breakdown can occur at lower swirl levels caused by flow perturbations induced by a PVC, for example. The transition to CVB involves the merging of this bubble with the wake behind the bluff body and the radial expansion of the bubble, but more importantly, the study concluded that the transition to CVB is governed by a nonlinear exchange of dominance between the central low-pressure region associated with the vortex breakdown recirculation zone and the outer low-pressure region generated by wall-jet entrainment. At even higher swirl, the Coandă effect drives the formation of a wall-attached jet analogous to the wide-open CVB of Moise & Mathew (Reference Moise and Mathew2021). Another important result from Douglas et al. (Reference Douglas, Emerson and Lieuwen2022) for the context of this work is the conclusion that a transition to or from the wall-jet (or CVB) regime is unlikely to occur spontaneously due to flow and swirl level fluctuations, normally requiring a significant change in bulk flow swirl level. The authors argue that a BVB-CVB transition due to flow fluctuations would require strong perturbations, such as external forcing or special initial conditions, that could cause the flow to ‘jump’ spontaneously between the breakdown (BVB) and wall-jet (CVB) solution manifolds. This conclusion is similar to that of Mesquita et al. (Reference Mesquita, Vié, Zimmer and Ducruix2021a ), which demonstrated how flow oscillations induced by a thermoacoustic instability led to a transition from the CVB mode to the BVB mode. Douglas et al. (Reference Douglas, Emerson and Lieuwen2022) also found that larger centre-body diameters and higher Reynolds numbers lower the critical swirl needed for these transitions, due to increased azimuthal momentum and reduced entrainment, respectively. Their analysis reveals that even in laminar flows, these mechanisms produce a richly structured and multi-stable parameter space.

Finally, Keeton et al. (Reference Keeton, Carpio, Nomura, Sánchez and Williams2022) employed theoretical and numerical approaches to study vortex breakdown in variable-density laminar swirling jets. They assessed the effect of increasing the swirl number for cases with different density ratios between the swirling jet and the atmosphere. This work provides a first insight into how swirling flow topologies could behave under reacting conditions. The authors demonstrate that altering the jet-to-ambient density ratio impacts the critical swirl number that triggers the BVB and CVB modes, revealing a significant decrease with increasing density ratio for thin initial shear-layer thicknesses. However, for thicker initial shear layers, the critical swirl numbers showed little dependence on the density ratio.

These few works addressing the transition of BVB to CVB mode highlight the complexity of the problem. In particular, Rajamanickam & Basu (Reference Rajamanickam and Basu2018) emphasise the scarcity of available studies on CVB and the limited understanding of the problem. Despite the investigation of various parameters triggering the BVB-CVB transition (Reynolds number, temperature, inlet velocity profile, geometry and swirl number) and the fact that all the works that observed a CVB (Sarpkaya Reference Sarpkaya1995; Billant et al. Reference Billant, Chomaz and Huerre1998; Santhosh et al. Reference Santhosh, Miglani and Basu2014; Terhaar et al. Reference Terhaar, Reichel, Schrödinger, Rukes, Paschereit and Oberleithner2015; Rajamanickam & Basu Reference Rajamanickam and Basu2018; Moise & Mathew Reference Moise and Mathew2019; Moise Reference Moise2020; Moise & Mathew Reference Moise and Mathew2021; Douglas et al. Reference Douglas, Emerson and Lieuwen2022; Keeton et al. Reference Keeton, Carpio, Nomura, Sánchez and Williams2022) showed that it evolves from an already stabilised BVB, it is still an open question if a general mechanism controlling the change of a BVB into a CVB exists. Moreover, Moise (Reference Moise2020), Moise & Mathew (Reference Moise and Mathew2021) and Keeton et al. (Reference Keeton, Carpio, Nomura, Sánchez and Williams2022) stress the importance of further studying the effect of temperature on BVB to CVB transitions.

Furthermore, when considering turbulent reacting flows, the study of BVB to CVB transitions, and in particular their impact on the change in flame shape, is even less documented in the literature. Even if vortex breakdown mode transitions are likely to occur simultaneously with the flame shape changes (and vortex breakdown mode changes are certain to impact the flame shape), the existing studies of transitions and bifurcations in reacting flows focus only on flame shape changes as a function of operating parameters, rather than the flow topology. Nevertheless, some works where these phenomena might be happening simultaneously are discussed in the following paragraphs, as well as other works that provide relevant context for this study.

2.3. Flame shape transitions

2.3.4. Fuel stratification effects

In staged fuel injectors, some studies discuss how flame shape transitions are caused by a change in the ratio between fuel injected by each stage. Hermeth et al. (Reference Hermeth, Staffelbach, Gicquel, Anisimov, Cirigliano and Poinsot2014) studied the occurrence of bistable flames in a two-staged gaseous industrial burner and the triggering of the transition caused by the change in fuel staging. The authors showed that two different flame shapes were possible for the same operating conditions: an attached M flame (similar to figure 2 f, seemingly a CVB/A-ISL/A-OSL flame, even if the vortex breakdown mode is not identified) and a lifted M flame (similar to figure 2 g, seemingly a BVB/L-BLE/A-OSL flame, even if the vortex breakdown mode is also not identified). To achieve this, a different path must be taken while changing the operating conditions to allow the stabilisation of each of the flames. To obtain the attached M flame at the operating point of interest, the operating parameters (namely, pressure, compressor temperature, air and fuel mass flow rates at each stage) must be increased slowly in a quasi-steady manner, allowing the flame to maintain its original attached shape. When operating conditions are changed abruptly, the flame detaches, stabilising in the lifted M shape. The authors observed significantly different HRRs, velocity, and temperature fields for each of the flames, with a PVC only present in the lifted M flame. In particular, the flow topology differs, seemingly indicating a change in vortex breakdown mode as the flame changes, although this is not discussed by the authors. Additionally, the acoustic activity of each flame was investigated through the construction of flame transfer functions for each flame, revealing that the detached flame exhibited a higher response amplitude for all studied frequencies, indicating that it was more acoustically active than the attached one. Finally, Hermeth et al. (Reference Hermeth, Staffelbach, Gicquel, Anisimov, Cirigliano and Poinsot2014) also studied other ways of triggering the flame bifurcations after they were stabilised on the operating conditions of interest. The first transition studied was from the attached M flame with a large CRZ to the lifted M flame, where a bubble-like CRZ is present. The first way to trigger this transition was to modify the pilot fuel ratio (defined as the ratio between the current pilot fuel mass flow rate and the nominal pilot fuel mass flow rate). This ratio controls the local equivalence ratio in different parts of the flame. The authors observe that the flame detaches when this ratio is reduced below $0.5$ or increased beyond $4.0$ , highlighting that there is an ideal range of local equivalence ratio where an attached flame is possible. The second transition, from the lifted M flame to the attached M one, is triggered by introducing a forced modulation of the flow. The flame leading point of the lifted M flame flashes back when the forcing amplitude is higher than $30\,\%$ of the mean diagonal inlet velocity. Each of these mechanisms only affects the flames in one direction, meaning that forcing the attached flame at the studied amplitudes or changing the pilot fuel mass flow rate inside the same range of values for the attached flame did not trigger the flashback of the lifted flame. This study highlights the role of the pilot injector and of flow modulations on the motions of the central portion of the flame and its stabilisation mechanism.

In a separate family of studies on gaseous staged combustion, Han et al. (Reference Han, Laera, Morgans, Lin and Sung2018) and Han et al. (Reference Han, Laera, Morgans, Sung, Hui and Lin2019) demonstrate that different flame shapes can be achieved by varying the stratification ratio or the global equivalence ratio. However, the flow field is not significantly impacted by the changes in flame shape. Moreover, vortex breakdown mode changes are neither discussed nor seem to be linked with the flame shape changes, except in very extreme conditions where some of the swirler passages are closed, but these flow fields are not representative of real operation.

Finally, the BIMER combustor family of works brings the flame transition studies closer to aeronautical applications by introducing liquid fuel. Providakis (Reference Providakis2013) characterised the impact of the staging factor $\alpha$ (ratio of fuel injected by the pilot injector, over the total fuel injected, as defined in figure 2) on the flame shape, while also studying the two-phase flow combustion physics – spray atomisation, droplet evaporation and fuel mixing. Providakis (Reference Providakis2013) observed a stable V flame from a pilot-only condition ( $\alpha =100\,\%$ ) down to $\alpha = 60\,\%$ (similar to figure 2 d, so a CVB/S-ISL/L-OSL flame). Further decreasing the staging factor, Providakis (Reference Providakis2013) observed the development of a strong thermoacoustic instability at around $\alpha = 40\,\%$ with high levels of sound pressure level (SPL). This led to an unstable behaviour of the flame, which no longer presented a V shape, sometimes seeming to switch to a tulip shape (similar to figure 2 c, so a BVB/S-ISL/A-OSL flame). Even further reducing the staging factor to $\alpha =20\,\%$ , the flame bifurcates and stabilises in an M shape (similar to figure 2 e, so a BVB/L-BLE/A-OSL flame), increasing the SPL in the chamber. When the staging factor was re-increased from $\alpha =20\,\%$ to $\alpha =60\,\%$ , the flame retained its M shape instead of reverting to the original V shape. Thus, one of the main findings of Providakis (Reference Providakis2013) was to identify this bi-stable behaviour of the burner and also the conditions under which bi-stability was observed. However, no discussion was made concerning vortex breakdown modes or correlations between vortex breakdown mode changes and flame shape changes. Additionally, Providakis (Reference Providakis2013) created a database for validating numerical models of the BIMER combustor, with flow velocity and droplet diameter distributions available for various conditions.

Following the experimental works of Providakis (Reference Providakis2013), LES using an Euler–Euler approach to model the disperse liquid phase have been performed by Cheneau et al. (Reference Cheneau, Vié and Ducruix2019a ). Cheneau et al. (Reference Cheneau2019 a) started his numerical cycle with the ignition of the chamber at $\alpha =100\,\%$ (air mass flow rate $\dot {m}_{\textit{air}} = 53\,\rm g\,s^-{^1}$ and global equivalence ratio of $\phi = 0.6$ ), resulting in a tulip flame (similar to figure 2 a, so a BVB/S-ISL flame). When changing to $\alpha = 60\,\%$ to compare the numerical flame with the experimental one, the changing in the spatial distribution of fuel triggers a slow modification of the flame and flow: the flame moves upstream and widens passing from a tulip to a V flame after $\approx 80\,\textrm{ms}$ (similar to figure 2 d, so a CVB/S-ISL/L-OSL flame). Two paths are then followed after this bifurcation: either the staging factor can be re-increased to $\alpha =100\,\%$ , where the flame keeps its V shape, or the staging factor can be further decreased, following the procedure of Providakis (Reference Providakis2013). When the staging factor is further decreased, Cheneau et al. (Reference Cheneau, Vié and Ducruix2019a ) showed that the numerical simulations were not able to capture the strong thermoacoustic instability observed by Providakis (Reference Providakis2013) around $\alpha =40\,\%$ and, thus, the flame kept its V shape throughout the whole path until $\alpha =0\,\%$ . However, after approximately $100\,\textrm{ms}$ , a perturbation suddenly led to the V flame transitioning to an M flame (similar to figure 2 e, so a BVB/L-BLE/A-OSL flame). From the M flame obtained at $\alpha =0\,\%$ , the staging factor is then re-increased up to $\alpha =100\,\%$ , showing that once the flame stabilises in an M-shape, it is no longer able to switch shapes. With this thorough study, Cheneau et al. (Reference Cheneau, Vié and Ducruix2019a ) was able to map all possible flame shapes as a function of the staging factor and also demonstrated that the $\alpha =100\,\%$ case is a multi-stable point, where the three flames are possible and stable. Cheneau et al. (Reference Cheneau, Vié and Ducruix2019a ) also performed a thorough validation of the LES against the experimental data, showing a good match of velocity profiles and of flame shapes. However, no discussion was held concerning vortex breakdown modes, their relationship with flame shapes, or the mechanisms driving the transitions.

The most recent experimental and numerical coupled studies on the BIMER combustor were conducted by Renaud (Reference Renaud2015) and Mesquita (Reference Mesquita2021), respectively. Renaud (Reference Renaud2015) concentrated on further characterising the liquid-phase behaviour of each injection stage, the combustor stability and the flame shapes by varying the staging factor from $\alpha =100\,\%$ to $\alpha =0\,\%$ and modulating the flow with a siren. Moreover, the air mass flow rate and total power were reduced (respectively to $\dot {m}_{\textit{air}}=[32, 43] \,\rm g\,s^-{^1}$ and $\textit{Power} = [45, 73]\,\rm \,\rm kW$ ) to allow the closure of the hysteresis cycle first observed by Providakis (Reference Providakis2013). The first part of Renaud (Reference Renaud2015) study was dedicated to the complete characterisation of the liquid phase in terms of diameter and velocity distributions (through PIV and Phase Doppler Anemometry measurements). He first analysed the pilot-only case ( $\alpha =100\,\%$ ), to characterise the pilot spray, and then a Lean Premixed Prevaporised (LPP) representative point ( $\alpha =15\,\%$ ), to characterise the spatial distribution of fuel in a case dominated by multipoint injection. These were then studied numerically in Mesquita, Vié & Ducruix (Reference Mesquita, Vié and Ducruix2018), where a numerical model was developed to match the spray injection. Renaud (Reference Renaud2015) moved, then, to reacting cases, igniting the flame with a pilot-only injection ( $\alpha =100\,\%$ ), obtaining a V flame (figure 2 b, so CVB/S-ISL flame). The influence of the ignition procedure and operating conditions was then studied numerically in Mesquita et al. (Reference Mesquita, Vié and Ducruix2021b ), where it was shown that the flame shape after ignition can have a V, an M (figure 2 g, so a BVB/L-BLE/A-OSL flame) or a tulip (figure 2 a, so a BVB/S-ISL flame) shape depending on how the flame is ignited. Mesquita et al. (Reference Mesquita, Vié and Ducruix2021b ) also showed that the M and tulip flames retain the non-reacting flow field and BVB mode topology, while the V flame changes the vortex breakdown mode to a CVB. Renaud (Reference Renaud2015) then reduced the staging factor to $\alpha =15\,\%$ , where the stable V flame (figure 2 d, so a CVB/S-ISL/L-OSL flame) was studied in these mostly premixed conditions, and then to $\alpha =0\,\%$ (figure 2 f CVB/L-ISL/A-OSL flame), where, after the development of a thermoacoustic instability, the flame transitioned from a V shape to an M one (similar to figure 2 e, so a BVB/L-BLE/A-OSL flame). This transition was studied numerically by Mesquita et al. (Reference Mesquita, Vié, Zimmer and Ducruix2021a ), who demonstrated that the thermoacoustic instability triggers an intricate vortex collective interaction phenomenon that destabilises the flow field and causes the vortex breakdown mode to transition from CVB to BVB, resulting in a change in the flame shape as a consequence. Renaud et al. (Reference Renaud, Ducruix and Zimmer2017) showed that this V-to-M transition could also be triggered by forcing the $\alpha =15\,\%$ V flame with a siren. From the M flame at $\alpha =0\,\%$ , the staging factor was progressively increased and around $\alpha =20\,\%$ , Renaud et al. (Reference Renaud, Ducruix, Scouflaire and Zimmer2015b ) observed a flashback of the M flame, where its central portion moved upstream and reattached at the spray in a tulip shape (similar to figure 2 c). In a numerical study, Mesquita et al. (Reference Mesquita, Vié and Ducruix2023a ) detailed the mechanism by which the combustion chamber temperature, the thermoacoustic cycle and the change in staging factor all contribute to causing the M flame to flashback to a tulip flame. Renaud (Reference Renaud2015) showed that, in the conditions just after the transition, this tulip flame proved to be metastable, and after some milliseconds, it transitioned again to a V shape, closing the hysteresis cycle. With this comprehensive screening of the staging factor, Renaud (Reference Renaud2015) characterised the V, M and tulip flames under different operating conditions. These were carefully studied numerically in Mesquita et al. (Reference Mesquita, Vié and Ducruix2023b ), showing the effect of the staging factor on each flame shape, the CVB stabilisation mechanisms for the V flame and the BVB one for both M and tulip flames.

As shown in the literature review above, works from different groups suggest that the vortex breakdown mode changes simultaneously with flame shape transitions, even if this is not explicitly discussed by the authors. Furthermore, if we infer from the presented literature that BVB to CVB mode changes are indeed happening simultaneously with the flame shape transitions, results of Chterev et al. (Reference Chterev2014), Malbois et al. (Reference Malbois, Salaun, Frindt, Cabot, Renou, Grisch, Bouheraoua, Verdier and Richard2017), Taamallah et al. (Reference Taamallah, Dagan, Chakroun, Shanbhogue, Vogiatzaki and Ghoniem2019) and Cheneau et al. (Reference Cheneau, Vié and Ducruix2019a ) suggest that an increase in the flame reactivity drives BVB to CVB transitions, results of Providakis (Reference Providakis2013), Terhaar et al. (Reference Terhaar, Oberleithner and Paschereit2014), Hermeth et al. (Reference Hermeth, Staffelbach, Gicquel, Anisimov, Cirigliano and Poinsot2014) and Cheneau et al. (Reference Cheneau, Vié and Ducruix2019a ) suggest that a decrease in the flame reactivity drives CVB to BVB transitions, and results of Hermeth et al. (Reference Hermeth, Staffelbach, Gicquel, Anisimov, Cirigliano and Poinsot2014), Renaud et al. (Reference Renaud, Ducruix and Zimmer2017) and those of Mesquita et al. (Reference Mesquita, Vié, Zimmer and Ducruix2021a ) demonstrate flow modulations triggering vortex breakdown mode transitions. All these triggering mechanism and transitions are summarised in table 2. These results, in particular those suggesting a BVB to CVB transition with increasing equivalence ratio, provide quite a relevant context for this work, which demonstrates that the BVB transitions to a CVB simultaneously with the tulip flame's transition to a V-shaped flame as the equivalence ratio is increased. Therefore, not only do we show that the equivalence ratio is a clear trigger for this transition, but also that it is more appropriately identified as a BVB/S-ISL to CVB/S-ISL flow--flame transition rather than a simple change in flame shape.

Table 2.

Summary of the literature review on flame shape transitions likely occurring simultaneously with transitions in vortex breakdown mode.

4.1. Sub-grid closure models for LES

The sub-grid terms represent the motion of the fluid for the turbulent scales that are filtered and, therefore, not resolved on the LES grid. Normally, the resolved LES scales comprise the majority of the turbulence spectrum, and the modelled sub-grid scales can be assumed to be statistically homogeneous, isotropic and responsible for the viscous dissipation of energy, which simplifies their modelling.

For this work, we model the sub-grid scales with the functional wall-adapting local eddy viscosity (Nicoud & Ducros Reference Nicoud and Ducros1999). It was developed to reproduce the asymptotic behaviour of the turbulent viscosity in function of $(y^+)^3$ near the wall. It does not need any modelling of the boundary layer, nor any turbulent wall law. Thus, the turbulent viscosity is defined as

(4.11) \begin{equation} \nu _t = C^2_W \bar {\Delta }^2 \dfrac {(s^d_{\textit{ij}} s^d_{\textit{ij}})^{3/2}}{(\tilde {S}_{\textit{ij}} \tilde {S}_{\textit{ij}})^{5/2} + (s^d_{\textit{ij}} s^d_{\textit{ij}})^{5/4}} ,\end{equation}

with $s^d_{\textit{ij}}$ the traceless symmetrical part of the square of the velocity gradient tensor

(4.12) \begin{equation} s^d_{\textit{ij}} = \dfrac {1}{2}(\tilde {g}^2_{\textit{ij}} + \tilde {g}^2_{ji}) - \dfrac {1}{3}\tilde {g}^2_{kk}\delta _{\textit{ij}} ,\end{equation}

where $\bar {\Delta }$ is the characteristic size of the LES filter and $\tilde {g}_{\textit{ij}} = {\partial \tilde {u}_{\!j}}/{\partial x_i}$ the gradient of the resolved velocity. Thanks to numerical simulations of the isotropic homogeneous turbulence, the constant $C_W$ has been estimated as $0.5$ by (Nicoud & Ducros Reference Nicoud and Ducros1999).

The sub-grid terms of (4.2), (4.3) and (4.4) are approximated here as follows (Boileau Reference Boileau2007):

(4.13) \begin{align} &\tau ^{\textit{sgs}}_{\textit{ij}} = \bar {\rho } (\tilde {u_i u_{\!j}} - \tilde {u}_i \tilde {u}_{\!j}) \approx 2 \bar {\rho } \nu _t \left ( \dfrac {\partial \tilde {u}_i}{\partial x_{\!j}} + \dfrac {\partial \tilde {u}_{\!j}}{\partial x_i} \right ) - \dfrac {2}{3} \bar {\rho } \nu _t \dfrac {\partial \tilde {u}_k}{\partial x_k} \delta _{\textit{ij}} ,\end{align}

(4.14) \begin{align} &J^{\textit{sgs}}_{i,k} = - \bar {\rho } ( {u_{\!j} \tilde{Y}_k} - \tilde {u}_{\!j} \tilde {Y}_k) \approx -\bar {\rho } \left ( D_t \dfrac {W_k}{W} \dfrac {\partial \tilde {X}_k}{\partial x_i} - \tilde {Y}_k \tilde {V}_{\!j}^c \right ) ,\end{align}

(4.15) \begin{align} &q_i^{\textit{sgs}} = - \bar {\rho } \left ( \tilde {u_{\!j} E} - \tilde {u}_{\!j} \tilde {E} \right ) \approx - \lambda _t \dfrac {\partial \tilde {T}}{\partial x_i} + \sum _{k=1}^N J^{\textit{sgs}}_{i,k} \tilde {h}_{s,k} , \\[9pt] \nonumber\end{align}

with $\nu _t$ the sub-grid scale turbulent viscosity, $D_t = {\nu _t}/{Sc_t}$ the turbulent molecular diffusion coefficient and $\lambda _t = \bar {\rho } \nu _t \bar {C}_{\!p} Pr_t$ the turbulent thermal diffusion coefficient. Also, $Sc_t$ and $Pr_t$ are respectively the turbulent Schmidt and Prandtl numbers.

4.2. Two-phase flow modelling: Lagrangian discrete particle simulation approach

One of the key points to correctly capture the governing phenomena in two-phase combustion is the numerical modelling of the liquid phase. As liquid fuel is injected to feed the combustion process, accurately representing the liquid phase and simulating the two-phase flow dynamics are essential for accurate predictions. The liquid fuel is injected at high pressure, and the liquid sheet is destabilised, disintegrating into filaments until it is broken down to form a mist of droplets with diameters of the order of $1$ to $100\, \unicode{x03BC} {\rm m}$ , which evaporate to feed the flame in gaseous fuel. Excepting the region near the liquid injection, where primary atomisation, coalescence and secondary atomisation take place, the majority of the flow in aeronautical combustion chambers is located in the diluted suspension region, where the mist of droplets is located.

Since resolving numerically primary atomisation is currently still out of reach for various reasons in the case of spray combustion simulations (Fiorina et al. Reference Fiorina2016), a spray of droplets is commonly used to model the injection. Considering the context and application of this work, we rely on the Lagrangian point-particle strategy to model the liquid phase. In particular, a deterministic discrete particle simulation (DPS) is employed, where each numerical particle corresponds to a single physical particle. Thus, the fuel droplets are considered as point particles that are also local source terms for the gas under the following hypotheses:

1. (i) The spray can be considered as diluted, meaning that the volume occupied by the liquid and the collisions between the droplets are negligible. Moreover, the interactions between neighbouring droplets without contact are sufficiently small to also be neglected.

2. (ii) The droplets are spherical as a consequence of their small size. Indeed, for the small droplets, the surface tension forces are preponderant compared with other exterior forces that may modify the droplet shape. Thus, the equilibrium shape for the droplet is a sphere.

Therefore, with these hypotheses, a single droplet can be considered as a discrete point in space that is affected by the surrounding conditions, which will be evaluated at the particle’s position (often indexed by $f@p$ , meaning fluid at the particle’s position). The evolution of a particle’s state can be individually traced using the equations of motion that will be further discussed. The coupling with the gas phase is done through source terms of mass, momentum and energy (already presented in (4.1), (4.2), (4.3) and (4.4)). For this work, a two-way coupling is considered, where the liquid-phase influence over the gas phase ( $l\rightarrow g$ ) and the gas-phase influence over the liquid phase ( $g\rightarrow l$ ) are both modelled. However, this choice of modelling implies the need for another hypothesis: the considered particle perturbations to the flow field at $f@p$ are negligible compared with the perturbations caused by the other surrounding droplets to the flow field. This hypothesis is needed to consider that the flow at $f@p$ , modelled as the far-field environment for the macroscopic models, is not altered by the analysed particle. The error introduced by this hypothesis is studied by Sirignano (Reference Sirignano2010).

In the Lagrangian DPS approach, the trajectory, the equations of mass, momentum and energy conservation are solved for each droplet individually with the following set of equations:

(4.16) \begin{align} \dfrac {{\rm d} X_{p,i}(t)}{{\rm d}t} &= u_{p,i} ,\end{align}

(4.17) \begin{align} \dfrac {{\rm d} m_{\!p}(t)}{{\rm d} t} &= \dot {m}_{\!p} ,\end{align}

(4.18) \begin{align} \dfrac {d m_{\!p} u_{p,i}(t)}{{\rm d}t} &= F_{i,p}^{\textit{ext}} ,\end{align}

(4.19) \begin{align} \dfrac {{\rm d} m_{\!p} h_{s,p}(t)}{Dt} &= \dot {\varPhi }_{\!p}. \end{align}

With the index $p$ representing the particle, $i$ the corresponding spatial coordinate, $X$ the position, $m$ the mass, $F_{i,p}^{\textit{ext}}$ represent the external forces acting on the droplet, $\dot {m}_{\!p}$ is the mass variation and $\dot {\varPhi }_{\!p}$ is the droplet’s internal energy variation.

The coupling with the gas phase is done through the conservative transfer of the mass flux, momentum and energy flux of each droplet $n$ to the gaseous LES equations by means of the source terms $s^{l\rightarrow g}$

(4.20) \begin{align} s_m^{l\rightarrow g} &= \dfrac {1}{\Delta\! V} \sum _{n=1}^{N} \varPsi _n (x_{p,n})\dot {m}_{p,n} ,\end{align}

(4.21) \begin{align} s_{{\textit{mo}},i}^{l\rightarrow g} &= \dfrac {1}{\Delta\! V} \sum _{n=1}^{N} \varPsi _n (x_{p,n})(-m_{m,p} F_{i,p,n}^{\textit{ext}} + \dot {m}_{p,n} u_{p,n,i}) ,\end{align}

(4.22) \begin{align} s_{\textit{en},i}^{l\rightarrow g} &= \dfrac {1}{\Delta\! V} \sum _{n=1}^{N} \varPsi _n (x_{p,n}) \left ( -m_{p,n} \boldsymbol{F}_{p,n} \boldsymbol{\cdot }\boldsymbol{u}_{p,n} + \dfrac {1}{2} \dot {m}_{p,n} \left \| \boldsymbol{u}_{p,n} \right \|^2 - \dot {\varPhi }_{p,m} \right ) ,\end{align}

(4.23) \begin{align} s_{\textit{fuel}}^{l\rightarrow g} &= \dfrac {1}{\Delta\! V} \sum _{n=1}^{N} \varPsi _n (x_{p,n}) \dot {m}_{p,n}. \end{align}

Where the index $\textit{fuel}$ on the source term $s_{\textit{fuel}}^{l\rightarrow g}$ for the conservation of species equation indicates that only the fuel species will have a source term coming from the evaporation of the droplets. Here, $\Delta \!V$ is the control volume of the flux transfer (generally, the nodal volume, i.e. the fictive volume around a node), $N$ the number of particles inside that volume and $\varPsi _n(x_{\!p})$ an interpolation function depending on the mesh. For unstructured meshes, the interpolation is done by weighting the inverse of the particle’s distance to the considered node ( $d_{\!j}$ ) with the sum of the inverse of the distances of the particle to each node of the cell containing the particle ( $d_l$ )

(4.24) \begin{equation} \varPsi _{n,j}(x_{\!p}) = \dfrac {\dfrac {1}{d_{\!j}}}{\sum _{l=1}^{N_v} \dfrac {1}{d_l}} ,\end{equation}

where $j$ is the index for the nodes of the containing cell and $N_v$ is the number of nodes in the cell.

The Lagrangian DPS approach enables modelling the spray as if it were composed of a group of isolated droplets, allowing for the subsequent analysis of the spray as the collective behaviour of these isolated droplets. The necessary closure models are then presented in the following sections.

4.3. External forces acting on an isolated droplet

The study of the dynamics of an evaporating droplet often neglects the effects of evaporation on the drag forces, considering only the motion of a non-evaporating rigid sphere in a flow. This reasonable hypothesis simplifies the modelling, and the forces acting on the surface of the rigid, spherical droplet can be divided into steady and unsteady forces.

The resulting total force ( $F_{i,p}^{\textit{ext}}$ ) from (4.18) for this work finally considers the drag (the particle’s weight is neglected, as in this application is much inferior to the other forces)

(4.25) \begin{equation} F_{i,p} = \dfrac {1}{\tau _{\!p}}(u_{\! f@p} - u_{p,i}) ,\end{equation}

and the characteristic time $\tau _{\!p}$ with the Schiller & Naumann (Reference Schiller and Naumann1935) correction is defined as

(4.26) \begin{equation} \tau _{\!p} = \dfrac {\rho _l d_{\!p}^2}{18 \mu _g f({\textit{Re}}_{\!p})} \quad \text{with } f({\textit{Re}}_{\!p}) = \left(1+0.15{\textit{Re}}_{\!p}^{0.687}\right)\! .\end{equation}

4.4. Evaporation source terms

The evaporation model used in this work is the one developed by Abramzon & Sirignano (Reference Abramzon and Sirignano1989). The authors proposed to modify the Sherwood and Nusselt numbers to consider the effect of the convective boundary layer on the evaporation model, by introducing two correction factors $F_M$ and $F_T$ that are defined as

(4.27) \begin{align} &F_M = (1+B_{\!M})^{0.7} \dfrac {\text{ln}(1+B_{\!M})}{B_{\!M}} ,\end{align}

(4.28) \begin{align} &F_T = (1+B_T)^{0.7} \dfrac {\text{ln}(1+B_T)}{B_T} ,\end{align}

where $B_{\!M}$ and $B_T$ are Spalding mass and heat transfer numbers, respectively. The corrected Sherwood and Nusselt numbers are then

(4.29) \begin{align} &\textit{Sh}^* = 2 + \dfrac {(\textit{Sh} - 2)}{F_M} ,\end{align}

(4.30) \begin{align} & \textit{Nu}^* = 2 + \dfrac {(\textit{Nu} - 2)}{F_T} ,\end{align}

where $Sh$ is the Sherwood number and $Nu$ the Nusselt number. This model allows for a significant improvement in the numerical representation of droplet evaporation. The final equations for the droplet mass flux $\dot {m}_{\!p}$ , gaseous conductive flux $\varPhi _g^{\textit{cd}}$ and droplet enthalpy flux $\varPhi _{\!p}$ are then

(4.31) \begin{align} &\dot {m}_{\!p} = -\textit{Sh}^* \pi d_{\!p} [\rho D_F] \text{ln}(B_{\!M} +1) ,\end{align}

(4.32) \begin{align} &\varPhi _g^{\textit{cd}} = \textit{Nu}^* \pi d_{\!p} \dfrac {\lambda }{C_{\!P}}(C_{\!P}(T_{\zeta })T_{\zeta } - C_{\!P}(T_{\infty })T_{\infty })\dfrac {\text{ln}(B_T + 1)}{B_T} ,\end{align}

(4.33) \begin{align} &\varPhi _{\!p} = \dfrac {d}{{\rm d}t}(m_{\!p} h_{\!s,l} (T_{\!p})) = -\varPhi _g^{\textit{cv}} - \varPhi _g^{\textit{cd}} ,\end{align}

where $d_{\!p}$ is the droplet diameter, $D_F$ is the fuel’s diffusivity, $\lambda$ is the thermal diffusivity, $C_{\!P}$ is the constant pressure heat capacity of the mixture, $T_{\zeta }$ is the temperature at the surface of the droplet, $T_{\infty }$ is the far-field temperature, $\varPhi _g^{\textit{cd}}$ and $\varPhi _g^{\textit{cv}}$ are respectively the gaseous conductive and the convective fluxes and $h_{s,l}$ is the sensible enthalpy of the droplet.

4.5. Modelling the pilot spray: the fuel injection model by upstream reconstruction (FIMUR) for pressure-swirl atomisers

The internal geometry of a typical pressure-swirl simplex atomiser, representative of the BIMER pilot injector. The liquid enters the atomiser through a tube perpendicular to the profile cross-section, generating a swirling motion of the liquid inside the atomiser cavity. The fuel is centrifuged to the walls of the atomiser, while the swirling movement causes the entrainment in a cylindrical region of radius $R_a$ at the axis of rotation inside the injector. The discharge channel leads to an orifice of radius $R_0$ , from where fuel is then injected inside the chamber. The angle formed by the axis of rotation and the droplets is called the spray half-angle and defined as $\theta _S$ . This angle is mostly a function of the injector geometry and of the size of the injection hole, while the droplet size depends mostly on the thickness of the liquid sheet inside the atomiser (affected by the fuel characteristics), the fluid velocity (controlled by the injection pressure) and interaction with the surrounding air (Rizk & Lefebvre Reference Rizk and Lefebvre1983).

Using the assumption of a frictionless flow inside the atomiser and knowing the spray angle $\theta _S$ , it is possible to use empirical relationships to determine the contraction factor at the discharge orifice X (Rizk & Lefebvre Reference Rizk and Lefebvre1985)

(4.34) \begin{equation} X = \dfrac {R_{a}^2}{R_0} = \dfrac {\sin^2\theta _S}{1+\cos^2\theta _S} .\end{equation}

Knowing the ratio X allows us to define the discharge coefficient $C_D$ (Rizk & Lefebvre Reference Rizk and Lefebvre1985)

(4.35) \begin{equation} C_D = 1.17 \sqrt {\dfrac {(1-X)^3}{1+X}} .\end{equation}

The FIMUR (Sanjosé Reference Sanjosé2009) utilises these geometric parameters ( $ R_a$ , $ R_0$ and $C_D$ ), the fuel mass flow rate $\dot {m}_l$ and the spray characteristics ( $\theta _S$ ) to determine the velocity profiles of the injected droplets. The model defines $R^0_S = (R_0 + R_a)/2$ and supposes that the axial velocity is constant at the plane $x = x_0$ and that, as a consequence of the swirling motion, the tangential velocity evolves linearly with the radius, having its maximum near the injector walls. A solid-body rotation is assumed for the liquid inside the atomiser cavity, which leads to a zero radial velocity. Finally, to determine the zone of particle injection, it is supposed that the liquid volume fraction is one between $[R_a, R_0]$ and zero at the central air entrainment zone. These assumptions lead to the following set of equations imposing the velocity profiles at the injector exit ( $x=x_0$ ) and the region of injection:

(4.36) \begin{align} u^0_{l,x}(\theta ,r_0) &= \dfrac {\dot {m}}{\rho _l \pi R_0^2 (1-X)} ,\end{align}

(4.37) \begin{align} u^0_{l,r}(\theta ,r_0) &= 0 ,\end{align}

(4.38) \begin{align} u^0_{l,\theta }(\theta ,r_0) &= \dfrac {\dot {m}}{\rho _l A_{\!p}} \dfrac {r_0}{R^0_S} ,\end{align}

(4.39) \begin{align} \alpha ^0_l(\theta ,r_0) &= \left \{\begin{matrix} 0 \, \text{if} \, r_0 \in [0, R_a],\\[3pt] 1 \, \text{if} \, r_0 \in [R_a, R_0]. \end{matrix}\right . \end{align}

Particles are injected from a ring with radii varying in the range $[R_a, R_0]$ and with an accordingly assigned velocity. The profiles resulting from these equations are illustrated in Sanjosé (Reference Sanjosé2009).

Regarding the droplet diameter, the FIMUR does not utilise any correlation; instead, the user is responsible for inputting the distribution from which the diameter will be randomly sampled. This means in particular that the velocity and diameter attribution to a droplet are independent. The FIMUR has been validated and applied with success for various aeronautical applications (Sanjosé Reference Sanjosé2009; Sanjosé et al. Reference Sanjosé, Senoner, Jaegle, Cuenot, Moreau and Poinsot2011; Hannebique et al. Reference Hannebique, Sierra, Riber and Cuenot2013) and on the BIMER combustor (Cheneau Reference Cheneau2019; Cheneau et al. Reference Cheneau2019a , Reference Cheneau, Vié and Ducruixb ; Mesquita et al. Reference Mesquita, Vié and Ducruix2017, Reference Mesquita, Vié and Ducruix2018, Reference Mesquita, Vié, Zimmer and Ducruix2021a , Reference Mesquita, Vié and Ducruixb , Reference Mesquita, Vié and Ducruix2023a , Reference Mesquita, Vié and Ducruixb ). For the BIMER pilot spray, the injected droplet properties (i.e. velocity, diameter distribution) have been extracted from the experimental studies of Renaud et al. (Reference Renaud, Ducruix, Scouflaire and Zimmer2015b ) and the injected Probability Density Functions (PDFs) correspond to the measurements. This injection set-up has been compared with the experiments and validated in the dedicated study of Mesquita et al. (Reference Mesquita, Vié and Ducruix2018). Due to the absence of measurements for different staging values, the diameter distribution is assumed to be constant, and the velocity is scaled according to the fuel mass flow rate.

4.6. Modelling the multipoint injection

The injection of fuel through the multipoint system is much simpler compared with the pilot pressure-swirl atomiser. In this system, fuel is injected through 10 equally spaced holes with only axial velocity, creating a jet-in-cross-flow that promotes atomisation. The numerical strategy is similar to this procedure: droplets are directly initialised over the surfaces of the injection holes, and the axial velocity is imposed directly from the fuel mass flow rate. This simple numerical modelling, despite ignoring some of the physics of jet-in-cross-flow atomisation, such as the column or surface breakups, has shown a satisfactory performance on the numerical simulations of the BIMER combustor by Cheneau et al. (Reference Cheneau, Vié and Ducruix2019a ,Reference Cheneau, Vié and Ducruix b ); Mesquita et al. (Reference Mesquita, Vié and Ducruix2017, Reference Mesquita, Vié and Ducruix2018, Reference Mesquita, Vié, Zimmer and Ducruix2021a , Reference Mesquita, Vié and Ducruixb , Reference Mesquita, Vié and Ducruix2023a , Reference Mesquita, Vié and Ducruixb ), but also in other similar applications (e.g. Jaegle Reference Jaegle2009). Similarly to the pilot spray, the injected droplet properties (i.e. velocity, diameter distribution) for the multipoint injection have also been modelled from the experimental studies of Renaud et al. (Reference Renaud, Ducruix, Scouflaire and Zimmer2015b ), after a detailed study and validation in Mesquita (Reference Mesquita2021). Thus, injected PDFs correspond to the measurements, with the diameter distribution considered constant for all staging values and the velocity is scaled according to the fuel mass flow rate.

4.7. Chemistry modelling: the 2S_KERO_BFER reduced scheme

The fuel used in this work is kerosene. The primary factor in selecting a chemical scheme for this work was the balance between computational cost and accuracy, as many operating points needed to be simulated, and the use of Lagrangian modelling for the liquid phase already implies a high computational cost. Cheneau et al. (Reference Cheneau, Vié and Ducruix2019b ) and Mesquita et al. (Reference Mesquita, Vié, Zimmer and Ducruix2021a , Reference Mesquita, Vié and Ducruixb , Reference Mesquita, Vié and Ducruix2023a ,Reference Mesquita, Vié and Ducruix b ) compared the simulations with the experiments and showed that modelling kerosene with the 2S_KERO_BFER reduced scheme developed by Franzelli et al. (Reference Franzelli, Riber, Sanjosé and Poinsot2010) (assuming a unity Lewis number) produces good results for the LES of the BIMER combustor, validating its choice. Not only is the power of the burner well-captured numerically in all conditions and with all flame shapes, but dynamic phenomena, as ignition and flame shape transitions, are also well captured, as shown by Mesquita et al. (Reference Mesquita, Vié, Zimmer and Ducruix2021a ,Reference Mesquita, Vié and Ducruix b , Reference Mesquita, Vié and Ducruix2023a ,Reference Mesquita, Vié and Ducruix b ).

The kinetic scheme 2S_KERO_BFER is based on 6 species and 2 reactions, the oxidation of kerosene and the $\textit{CO}-\textit{CO}_2$ equilibrium

(4.40) \begin{align} \textit{KERO} + 10 O_2 &\rightarrow 10 \textit{CO} + 10 H_2 O ,\end{align}

(4.41) \begin{align} {\textit{CO}} + 0.5O_2 &\leftrightarrow {\textit{CO}}_2. \end{align}

4.8. Combustion modelling: thickened flame model for LES

The thickened flame model for LES (TFLES) coupled with the 2S_KERO_BFER reduced chemistry is a combustion model that has already proven to produce accurate simulations in the context of turbulent two-phase multi-regime combustion (Shum-Kivan et al. Reference Shum-Kivan, Marrero Santiago, Verdier, Riber, Renou, Cabot and Cuenot2017; Cheneau et al. Reference Cheneau2019b ; Mesquita et al. Reference Mesquita2021a , Reference Mesquita, Vié and Ducruixb , Reference Mesquita, Vié and Ducruix2023a ,Reference Mesquita, Vié and Ducruix b ).

The principle behind the thickened flame model, firstly proposed by Butler & O’Rourke (Reference Butler and O’Rourke1977), is to divide the reaction rate and multiply the species and thermal diffusion by a factor $\mathcal{F}$ , aiming to increase the flame thickness by the same factor $\mathcal{F}$ , allowing it to be resolved on the LES mesh, and keeping unaltered the flame speed. Recalling that for a laminar premixed flame the flame thickness $\delta _L$ and the flame speed $s_{\!L}$ are defined as (Kuo & Acharya Reference Kuo and Acharya2012)

(4.42) \begin{align} s_{\!L} \propto \sqrt {\dot {\omega }_T \mathcal{D}_{\textit{th}}} \qquad \text{and} \qquad \delta _L \propto \dfrac {\mathcal{D}_{\textit{th}}}{s_{\!L}} = \sqrt {\dfrac {\mathcal{D}_{\textit{th}}}{\dot {\omega }_T}} ,\end{align}

where $\mathcal{D}_{\textit{th}} = \lambda _{\textit{th}}/\rho C_{\!p}$ is the thermal diffusivity and $\dot {\omega }_T$ is the HRR. The thickened quantities become then:

(4.43) \begin{align} \bar {\dot {\omega }}_k = \dfrac {\dot {\omega _k}}{\mathcal{F}} \qquad &\text{and} \qquad \bar {\dot {\omega }}_T = \dfrac {\dot {\omega }_T}{\mathcal{F}} ,\end{align}

(4.44) \begin{align} \mathcal{D}_k \rightarrow \mathcal{F}\mathcal{D}_k \qquad &\text{and} \qquad \mathcal{D}_{\textit{th}} \rightarrow \mathcal{F}\mathcal{D}_{\textit{th}}. \end{align}

Leading to

(4.45) \begin{align} \delta _L^{\textit{TF}} = \mathcal{F} \sqrt {\dfrac {\mathcal{D}_{\textit{th}}}{\dot {\omega }_T}}. \end{align}

However, the artificial thickening of the flame modifies the flame–turbulence interactions, a modification that is visible through the consequent reduction by the same factor $\mathcal{F}$ of the Damkhöler number, which compares the characteristic time of turbulence $\tau _t$ with the characteristic time for chemistry $\tau _c$

(4.46) \begin{align} \textit{Da} = \dfrac {\tau _t}{\tau _c} = \dfrac {l_t s_{\!L}}{u' \delta _L} \longrightarrow {\textit{Da}}^{\textit{TF}} = \dfrac {l_t s_{\!L}}{u' \delta _L \mathcal{F}} .\end{align}

The thickened flame is then less sensitive to turbulence and the wrinkling of the flame front is underestimated, thus reducing the flame surface and, consequently, the reactants consumption rate (Angelberger et al. Reference Angelberger, Veynante, Egolgopoulos and Poinsot1998). For this reason, the effect of turbulence on the flame wrinkling was studied with DNS, and Colin et al. (Reference Colin, Ducros, Veynante and Poinsot2000) proposed an efficiency function $\mathcal{E}$ to account for the loss of wrinkling due to thickening and for the sub-grid flame wrinkling, aiming to improve the modelling of flame–turbulence interactions

(4.47) \begin{align} \bar {\dot {\omega }}_k = \dfrac {\mathcal{E} \dot {\omega }_k}{\mathcal{F}} \qquad &\text{and} \qquad \bar {\dot {\omega }}_T = \dfrac {\mathcal{E} \dot {\omega }_T}{\mathcal{F}} ,\end{align}

(4.48) \begin{align} \mathcal{D}_k \rightarrow \mathcal{E} \mathcal{F} \mathcal{D}_k \qquad &\text{and} \qquad \mathcal{D}_{\textit{th}} \rightarrow \mathcal{E} \mathcal{F} \mathcal{D}_{\textit{th}}. \end{align}

The TFLES model is then summarised in table 3.

Table 3.

Summary of the TFLES model impact on combustion.

Many efficiency functions are available in the literature. The chosen efficiency function for this work is the one proposed by Charlette, Meneveau & Veynante (Reference Charlette, Meneveau and Veynante2002)

(4.49) \begin{equation} \mathcal{E} = \varXi \left ( \dfrac {\Delta _e}{\delta ^0_L} \right ) = \left ( 1 + \textit{min} \left [ \dfrac {\Delta _e}{\delta ^0_L} , \varGamma \dfrac {u^{\prime}_{\Delta _e}}{s^0_L} \right ] \right ) ^\beta \!,\end{equation}

with $\beta = 0.5$ a constant of the model, $u^{\prime}_{\Delta _e}$ the sub-grid scale turbulent velocity and $\varGamma$ the spectral efficiency function estimated in Charlette et al. (Reference Charlette, Meneveau and Veynante2002).

The thickened flame model was initially designed for premixed flames and must be modified to be applicable in partially premixed cases, as the thickening in the original model is applied uniformly, modifying the species and thermal diffusion even for the species mixing in the non-reacting zones and overestimating diffusion fluxes by the factor $\mathcal{F}$ . Aiming to extend the flame thickening approach to non-premixed combustion, Legier, Poinsot & Veynante (Reference Legier, Poinsot and Veynante2000) proposed a dynamic version of the model, where the thickening factor is only applied on the reacting zones. Thus, a flame sensor $\mathcal{S}$ was developed to identify the reaction zones and the thickening is only applied where the sensor is activated

(4.50) \begin{align} &\mathcal{S} = \text{tanh} \left ( \beta ' \dfrac {\varOmega }{\varOmega _0} \right ) \text{ with } \beta '=50 ,\end{align}

(4.51) \begin{align} &\varOmega = Y_{F}^{\nu ^{\prime}_{F}} Y_{O}^{\nu ^{\prime}_{O}} \text{exp} \left (\! - \varGamma _T \dfrac {E_a}{\mathcal{R} T} \right )\! ,\end{align}

(4.52) \begin{align} &\mathcal{F} = 1 + (\mathcal{F}_{\textit{max}} - 1)\mathcal{S} ,\end{align}

where $\varOmega _0$ is the maximum value of $\varOmega$ in a one-dimensional laminar non-thickened flame, indices $F$ and $O$ correspond, respectively, to $Fuel$ and $\textit{Oxidant}$ and $\varGamma _T$ is the parameter used to activate the thickening function before the reaction, with a typical value of $\varGamma _T = 0.5$ .

It is common for simulations of complex geometries to encounter a wide range of mesh sizes. As a consequence, a single value for the filter size $\varDelta _e$ will not be adapted, as the number of points to resolve the flame front will vary. To counter this practical issue, the maximum value of the thickening factor is weighted by the local mesh size $\varDelta _x$ and by $n$ , the user-defined number of points used to resolve the thickened flame front

(4.53) \begin{equation} \mathcal{F}_{{\textit{max}}} = n \dfrac {\delta ^0_L}{\varDelta _x} .\end{equation}

Generally, $n=5$ points are sufficient to correctly solve the flame front, a value adopted for the simulations performed in this work.

For numerical simulations of two-phase combustion, modelling the interactions between the liquid phase and the flame front is central, as the liquid spray reaches the reaction zone in most industrial applications. Thus, intense evaporation of fuel occurs, imposing a strong coupling between the evaporating liquid phase and the gaseous chemical reaction, with the evaporation potentially controlling the local reaction rate. Even if the TFLES model was developed for gaseous premixed flames only, modifications were introduced to account for liquid fuel combustion and to consider flame–liquid-phase interactions, thereby correctly representing the behaviour of the non-thickened flame front. In the BIMER combustor, despite being an LPP burner (and thus with an intended operation in lean-premixed regime), sometimes part of the flame burns in a non-premixed regime, as the liquid phase does not always pre-vaporise and premixes entirely before reacting. Nevertheless, these regions remain marginal when compared with the total flame, which burns in a premixed regime in the majority of cases considered in this work. Thus, even if using TFLES introduces an error when these non-premixed regions are present, this error has been shown to be marginal for the LES of the BIMER combustor and not to affect the combustion dynamics in several comparisons with experiments. Indeed, Cheneau et al. (Reference Cheneau, Vié and Ducruix2019b ) and Mesquita et al. (Reference Mesquita, Vié, Zimmer and Ducruix2021a , Reference Mesquita, Vié and Ducruixb , Reference Mesquita, Vié and Ducruix2023a ,Reference Mesquita, Vié and Ducruix b ) have shown that the use of TFLES coupled with 2S_KERO_BFER produced very good results for the BIMER combustor when compared with the experiments for pilot-only to multipoint-only operation, even in cases where the combustion was not completely premixed. Thus, this choice of modelling is kept in this work.

For a Lagrangian modelling of the liquid-phase, the approach used for the two-phase TFLES (TP-TFLES) is the one derived by Paulhiac (Reference Paulhiac2015), who proposed an analysis of the thickened flame model as a modification of time (t) and space (x) through the thickening factor

(4.54) \begin{align} &x^{\textit{TFLES}} \approx \dfrac {x}{\mathcal{F}} ,\end{align}

(4.55) \begin{align} &t^{\textit{TFLES}} \approx \dfrac {t}{\mathcal{F}}. \end{align}

Thus, to keep consistency, the same scale transformation must be applied to the liquid phase. This is achieved by applying the same scale transform as in (4.54) and (4.54) to the Lagrangian derivative ( $({D}/{Dt}) = ({\partial }/{\partial t}) + u_{p,i}({\partial }/{\partial x_i})$ ) and to the Lagrangian motion equations of particle $p$ ((4.16), (4.18), (4.17), (4.19)), leading to the modification of the drag force, evaporation rate and heat transfer as the particles’ characteristic time $\tau _{\!p}$ (4.26), mass variation $\dot {m}_{\!p}$ (4.31) and internal energy variation $\varPhi _{\!p}$ (4.33) are modified by the thickening factor

(4.56) \begin{align} \tau ^{\textit{TFLES}}_{\!p} &= \mathcal{F} \tau _{\!p} ,\end{align}

(4.57) \begin{align} \dot {m}^{\textit{TFLES}}_{\!p} &= \dfrac {\dot {m}_{\!p}}{\mathcal{F}} ,\end{align}

(4.58) \begin{align} \varPhi ^{\textit{TFLES}}_{\!p} &= \dfrac {\varPhi _{\!p}}{\mathcal{F}}. \end{align}

This implies that in the thickened space, the droplets evaporate $\mathcal{F}$ times slower and experience a drag force that is $\mathcal{F}$ times lower than the actual drag force. The calculations made by Paulhiac (Reference Paulhiac2015) demonstrate that these modifications ensure the TP-TFLES modelling for Lagrangian particles accurately captures the non-thickened spray flame behaviour.

4.9. Numerical set-up

The LES is performed using the AVBP code, co-developed by CERFACS and IFPEN (Wolf et al. Reference Wolf, Balakrishnan, Staffelbach, Gicquel and Poinsot2012), using a two-step Taylor–Galerkin (TTGC) numerical scheme (third order in space and time), along with the artificial diffusion implemented by Colin & Rudgyard (Reference Colin and Rudgyard2000) to treat numerical perturbations and stiff gradients. The Navier–Stokes characteristic boundary conditions developed by Poinsot & Lele (Reference Poinsot and Lele1992) are used to impose the gaseous inlet and outlet boundary conditions. Walls are considered non-slipping and adiabatic, as no experimental measurements of the wall temperatures were available, except for the two water-cooled walls, which were treated as isothermal at $373\,\textrm {K}$ , the measured exit temperature of the cooling water.

An unstructured tetrahedra mesh of 32 million elements is used (figure 3). This mesh is a result of several mesh convergence studies and improvements presented in Cheneau, Vié & Ducruix (Reference Cheneau, Vié and Ducruix2015), Cheneau (Reference Cheneau2019) and Mesquita (Reference Mesquita2021), which resulted in several successful validations against experimental data (e.g. gaseous velocity, liquid velocity and diameter, flame shape, flame dynamics, SPL), as shown in Cheneau et al. (Reference Cheneau, Vié and Ducruix2019b ) and Mesquita et al. (Reference Mesquita, Vié and Ducruix2018, Reference Mesquita, Vié and Ducruix2021b , Reference Mesquita, Vié and Ducruix2023a , Reference Mesquita, Vié and Ducruixb ). In the reaction region, the mesh keeps the flame thickening factor between 3 and 7.