Numerical Investigation of Combustion Instabilities in Swirling Flames with Hydrogen Enrichment

This work presents a numerical study on technically premixed, swirl-stabilised flames in the PRECCINSTA model combustor. The employed method, BOFFIN-LES, comprises a fully compressible formulation to study unsteady combustion with thermo-acoustic instabilities. To allow for this, the iso-thermal flows are first investigated, based on which three reacting cases are established. The investigation delves into various aspects including flame topology, flow characteristics, and the related thermo-acoustic and hydrodynamic instabilities are studied and results are benchmarked against available measurement data. The dominant feedback mechanism of the observed thermo-acoustic fluctuations is identified; the evolution of the helical vortex is discussed together with the related flame stabilisation process. Furthermore, the interplay of the thermo-acoustic oscillations, helical structure, and the flame stabilisation process is summarised in the end, with the potential effect of the wall-heat transfer on them discussed. This work establishes that the Large Eddy Simulation (LES) effectively captures the iso-thermal flow dynamics and the flame topology under various operating conditions, with a good prediction of the thermo-acoustic frequencies in all the cases. The dominant driving mechanism of the observed thermo-acoustic fluctuations was identified as a combined effect of equivalence ratio and velocity fluctuations in all the cases investigated. The effect of Hydrogen enrichment on modifying the flame topology and changing the thermo-acoustic instability features are well predicted by the simulations. Moreover, different modes of the helical vortex are detected, and their periodic excitement, evolution, and effect on flame stabilisation are discussed in great detail. To conclude, this LES-based investigation offers valuable insights into the complex interplay of unsteady combustion, acoustic fluctuations, flow dynamics, and solid boundaries within swirling flames subjected to unsteady conditions.


Introduction
Fuel-lean premixed combustion, in which the fuel-air mixture contains a controlled excess of air, has the potential to fulfill both requirements of emission control and efficiency boost (Swaminathan and Bray 2011), and hence is often employed in modern combustor design.Although most combustion devices are designed to operate in stable regimes, combustors under lean premixed operations are sensitive to flow perturbations and can exhibit unexpected self-excited combustion instabilities which are typically only detected during late-stage prototype testing (Poinsot 2017).These combustion instabilities give rise to additional noise and vibration which can lead to serious problems including transient combustor issues such as local extinction and ignition, flame blow-off, flashback, increased thermal load to the burner, and even engine failure (Correa 1998).Meanwhile, the combustion of hydrogen-enriched flames has gained growing worldwide interest.Hydrogen as a fuel may offer higher energy density and lower carbon-based emission, and it can be reproduced from renewable power sources such as solar, wind, and geothermal energy (Edwards et al. 2008).Compared to CH 4 , H 2 has a higher adiabatic flame temperature (at stoichio- metric conditions in air, 2383 K and 2220 K), lower flammability limit ( = 0.14 and = 0.46) and higher maximum adiabatic laminar flame speed (320 cm/s and 40 cm/s) (Glassman et al. 2014).Hence changes in the fuel by hydrogen enrichment can drastically change the combustion characteristics, including the lean blow-off limit, auto-ignition, flashback, and combustion instabilities (Lieuwen et al. 2008).Taamallah et al. (2015a) and Shanbhogue et al. (2016) experimentally investigated the flame shape transition in a series of hydrogen-enriched swirling flames.They found that hydrogen addition brought down the critical equivalence ratio at which the flame-shape transitions occurred.Karlis et al. (2019) further examined the applicability of the extinction strain rate as a scaling parameter of the observed combustion dynamics transitions in methane hydrogen flames.
In general, combustion instabilities (often referred to as thermo-acoustic instabilities) are driven and controlled by a variety of complex, highly non-linear physical mechanisms (Lieuwen and Yang 2006).When the unsteady heat release excites the natural acoustic modes of the combustion chamber, a feedback loop between the combustion chamber acoustics and the unsteady reactive flow field is established, which generally couples the downstream flow to the upstream region.Computational Fluid Dynamics (CFD) has emerged as an increasingly sophisticated and essential technique.Complementing experimental techniques and theories, CFD is a useful tool to understand the complex processes involved in combustion dynamics as it has the ability to overcome some of the difficulties associated with measurement techniques.The balance of accuracy and computing cost required for the accurate prediction of unsteady turbulent flames has led to significant interest and development of Large-Eddy Simulations (LES) in recent years.Besides, LES has proved a useful tool to address unsteady flame dynamics due to intrinsic unsteadiness which is mandatory for the prediction of combustion dynamics.Detailed reviews of numerical prediction using LES methods for combustion instabilities can be found in the review paper by Poinsot (2017).
LES can be used for predicting combustion instabilities in two modes.The first approach is the forced open-loop LES, where incompressible LES (sometimes termed as 'low Mach number LES') is used together with low-order modeling techniques from acoustic theory (Giauque et al. 2005;Palies et al. 2011;Iudiciani and Duwig 2011;Han and Morgans 2015;Xia et al. 2019).The idea is to decompose the combustor into acoustic elements where the burner itself is seen as one part.A second method is to use a fully compressible LES formulation to compute the entire combustor as a resonator (Poinsot and Veynante 2005).The propagation of acoustic waves is therefore enabled and proper wave reflections can be included by well-defined acoustic boundary conditions.Hence, inherent coupling mechanisms between the flame, flow, and acoustics inside a combustion chamber can be directly accounted for and the LES should exhibit the self-excited modes exactly like what is found in the experiments.Despite the advantages, these self-excited LES methods face many challenges i.e. proper boundary conditions including acoustic impedance at all inlets and outlets should be prescribed and the long transition time from stable to unstable operating conditions observed in some combustors which are almost impractical to compute numerically.Further requirements for the LES study are necessary in order to describe the problem more 'accurately', e.g.An and Steinberg (2019) argued that numerical simulations hence need to be able to reproduce flame extinction due to turbulent straining and hydrodynamic instability, in order to accurately predict the swirling flame configurations.Hence, there are only a limited number of numerical works in the open literature using compressible LES to capture the self-excited modes in the combustors e.g.(Roux et al. 2005;Franzelli et al. 2012;Lourier et al. 2017;Noh et al. 2019;Tachibana et al. 2015;Chen et al. 2019;Fredrich et al. 2021a, b).
In order to predict self-excited combustion instabilities in gas turbine combustion chambers, a compressible LES code BOFFIN-LES has recently been devised (Fredrich 2019;Gong 2022).A transported probability density function (pdf) with a stochastic field method is employed for turbulence-chemistry interaction.The objective of the current work is to assess the predictive capabilities of the LES approach in the context of combustion dynamics under various operating conditions including different equivalence ratios and hydrogen addition and to investigate the flame dynamics in the swirling flames investigated.

Experimental Set-ups
The target swirling flame configuration is the well-known PRECCINSTA burner at DLR (Meier et al. 2007), which can be operated under a 'technically' premixed (TP) mode and a perfectly premixed (PP) mode.In reality, however, the premixing of fuel and oxidiser must occur inside the combustor for safety reasons, creating only partially premixed reactants because of the limited space and time available for mixing, which is referred to as the TP mode.In the TP mode, the fuel is injected into the swirler and does not premix perfectly with the air before reaching the flame.Compared to PP flames, thermo-acoustic oscillations in TP flames exhibit an additional complexity as they are driven not only by velocity fluctuations but also by equivalence ratio fluctuations in the unburned mixtures entering the flame zone (Lieuwen and Zinn 1998;Schuermans et al. 2010;Hermeth et al. 2013).Meier et al. (2007) identified a fluctuating fuel supply, subject to a convective time delay, as the main mechanism of the thermo-acoustic fluctuation in pulsating flame F2.
Recently Chterev and Boxx (2019) and Kushwaha et al. (2021) extended the range of operating conditions of the TP flames: P th from 10 kW to 30 kW, from 0.65 to 1.05 and hydrogen concentration from 0% to 40%.The flame and flow behaviour were found to vary significantly across different operating conditions, and the impact of H 2 enrichment on the flame dynamics also strongly depended on the operating conditions.The effect of hydrogen enrichment was also investigated by (Chterev and Boxx 2021) at elevated pressure up to 5 bar with a slightly increased burner length design.
Figure 1 shows a schematic of the PRECCINSTA burner, with four sections from left to right: the air plenum, injection unit, combustion chamber, and exhaust pipe.The pressure signals are recorded to provide information on combustion dynamics using 2 microphone probes: one in the plenum (Probe PL) and one in the chamber (Probe CH) as shown in Fig. 1.The flames are visualised by line-of-sight OH chemiluminescence (OH-CL) measurement which is a common and reliable marker for the heat release rate (Chterev and Boxx 2021).Velocity fields are measured using 10 kHz stereo partial image velocimetry (PIV) with the optical view marked in Fig. 1.The swirl number is fixed and is estimated at a value of 1 in the experiments.It is worth noting that the swirl number here is referring to the 'geometry swirl' number which is a predetermined value using the geometry dimensions of the swirler.The effective swirl number S describing the degree of swirl in a swirling flow originally proposed by Chigier and Beér (1964) is defined as: where Gy is the axial flux of the azimuthal momentum, Gz is the axial flux of axial momentum and R n is the burner nozzle equivalent radius.This effective S estimated using axial and azimuthal velocities was found to differ significantly from the theoretical estimate, which is smaller than the corresponding geometrical one in most cases (Palies et al. 2010).
Among all operating conditions, the present work focus on five cases as summarised in Table 1, with the burner operated under TP mode.Constant and uniform velocities, consistent with the mass flow rates, were applied at both the air and fuel inlet boundaries.All cases have a similar thermal power ∼ 20 kW.C1 features a lower fuel/air equivalence ratio = 0.65, while in C2 = 0.85 .In C1 and C2 the fuel consists of pure methane, while in C3 the fuel is enriched with 38% hydrogen by volume.Two iso-thermal cases N1 and N2 feature the same operating conditions as C1 and C2, respectively.All selected cases are performed under atmospheric pressure, and the fuel and air mixture are slightly preheated to reach temperatures between T ∼ 310 K and T ∼ 380 K.
The reason for such a selection is the rich physics observed in these cases.In C1, an oscillating flame with strong self-excited thermo-acoustic oscillations has been detected; the flame presents an M shape, which is closely related to the precessing vortex core (PVC) (Oberleithner et al. 2015) in the flow.C2 represents an unstable case with weak (1) Fig. 1 The experimental configuration of the PRECCINSTA burner (Chterev and Boxx 2019) thermo-acoustic oscillations: the dominant thermo-acoustic fluctuation amplitude is about 20% of the one in C1 in Pascal.Besides, the flame remains attached to the burner's centre body, which features a V shape.C3 aims at investigating the effect of H 2 enrichment compared to C2, where the thermo-acoustic oscillations found in C2 are suppressed and the flame retains a V shape.When compared with previous work on the PRECCINSTA flames, flame C1 has great similarity to the pulsating flame F2 in Meier et al. (2007), which numerically has been investigated by Franzelli et al. (2012); Lourier et al. (2017);Fredrich et al. (2021bFredrich et al. ( , 2021a)).The present work is one of the first numerical works to investigate the C2 flame and C3 flame (or any other stable flame) using a fully compressible LES method.

Numerical Set-ups
In order to investigate the thermo-acoustic response of the system, compressible LES, BOFFIN-LES is performed.In the present compressible LES scheme, the filtered Navier-Stokes equations together with the modeling strategy for the sgs stresses have the same implementation as those in previous works and details can be found in Gong (2022).The details of the combustion model are presented here for the completeness of the work.
Due to the difficulties encountered in evaluating the filtered values of the chemical source terms, a joint sgs (or filtered fine-grained) pdf evolution equation formulation is adopted.An exact equation describing the evolution of the pdf can be derived by standard methods, e.g.Gao and O'Brien (1998).Following the approach of Brauner et al. (2016) the resolved part of the convection and 'molecular' mixing are added to both sides of the equation so that the modeled form of the equation becomes: (2) where is the sample space for the scalar .ω () is the net formation rate through a chemical reaction in the case of chemical species, the pressure derivative Dp Dt for total enthalpy.The first term on the right-hand side of (2) represents the closed sgs-pdf transport using a Smagorinsky type gradient model with sgs and sgs representing the sub-grid scale diffusion coefficient and Schmidt number respectively; and sgs is assigned the value 0.7.The last term on the right-hand side represents sgs micro-mixing which is closed by the Linear Mean Square Estimation (LMSE) closure (Dopazo 1975;Dopazo and O'Brien 2003), and the mixing time-scale is given by sgs = Δ 2 ∕ sgs .The sgs mixing constant C d is assigned the value of 2.0 (Jones and Prasad 2010).The main advantage of the formulation is that the chemical reaction terms and the pressure derivative appear in a closed form and that no combustion regime-dependent modeling assumptions are invoked.
The equation describing the evolution of the pdf, (2) is solved using the Eulerian stochastic field method.Psgs () is represented by an ensemble of N stochastic fields for each of the N s + 1 scalars namely n (x, t) for 1 ≤ n ≤ N, 1 ≤ ≤ N s + 1 .In the present work, the Itô formulation of the stochastic integral is adopted.The stochastic fields evolve according to: where dW n i represents increments of a (vector) Wiener process, different for each field but independent of the spatial location x.
The reduced chemical kinetics mechanism employed for the oxidation of methane and hydrogen is the 15-step reduced mechanism by Lu and Law (2008).GRI-Mech 3.0 has been validated in many experimental works for CH 4 ∕H 2 ∕air flames and very good agree- ment was obtained for key combustion attributes including laminar flame speeds, ignition delay time, profiles of major species, extinction strain rates, etc Ren et al. (2001Ren et al. ( , 2002)); Halter et al. (2005); Chaumeix et al. (2007); Hu et al. (2009).Reduced mechanisms developed from GRI − 3.0 (Hawkes and Chen 2004;Safta and Madnia 2006) and from GRI − 2.11 (Jackson et al. 2003) were also examined for CH 4 ∕H 2 flames, where good agreement over a wide range of mixture conditions were achieved.Hence, it is convincing that the reduced mechanism compares well with the full GRI mechanism to capture the qualitative effect of H 2 addition.
The computational domain, as shown in Fig. 2 includes the combustion chamber, swirler, and the air plenum to resolve the fuel-air mixing since the perfect premixing assumption is not appropriate when predicting self-excited instabilities Franzelli et al. (2012).The computational mesh consists of about 2.7 million grid cells, which are clustered in the reaction zone and the fuel-air mixing region as shown in Fig. 2. It can be seen that in the swirler and nozzle region, the resolution of the mesh ranges from 0.1 to 0.8 mm in all three spatial directions; the smallest dimension of the swirler is about 6 mm, which makes it well resolved by the mesh; in the downstream region of the nozzle and the bottom half of the combustion chamber where the reaction takes place, the size of the grids range from 0.1 to 0.4 mm, which is comparable to a Kolmogorov length scale on the order of 0.1 mm estimated in the experiment by Meier et al. (2007) and (3) provides a sufficient resolution for the flame front.Besides, the mesh is refined towards the domain walls resulting in a dimensionless wall-normal distance y + between 20 and 30 in the swirler and nozzle regions, and below 10 in the combustion chamber.A meshindependent study has been done with the isothermal flow and the current mesh has been shown to reproduce the main flow motion and fluctuations with satisfactory accuracy as shown in Sect.4.2.1.Non-reflective outlet boundary condition (Yoo et al. 2005;Lodato et al. 2008) is employed, with a relaxation coefficient of 0.28 as recommended by Yoo et al. (2005), to minimise flow distortion and numerical reflection at the exit.The reflective inlet boundary condition is imposed at the inlet of the air plenum which should have a minimal effect in the combustion chamber since the swirler and nozzle form almost a fully reflective boundary for the chamber inlet.
Non-adiabatic combustion chamber walls with fixed temperature (Dirichlet boundary conditions) are implemented on all the chamber walls and the centre body to include the wall heat transfer effect.Ideally, wall temperatures under each operating condition should be measured and used as boundary conditions for the simulations, which, however, is not available for many experiments.In the PRECCINSTA burner, the surface temperatures were measured only by Yin et al. (2017) with the operating condition P th = 20 kW and = 0.7 , which is different from the operating condition selected for the current study but close to C1 with a slightly lower equivalence ratio ( = 0.65 ).The chamber side walls were measured to be around 1400 K and the base plate wall around 700 K, and these values are used as fixed wall temperatures in the current simulations.The resulting preheated air and fuel temperature are estimated as 320 K in the simulations based on the measurements of the F1 flame.The rationality of these estimations is confirmed by the experimentalists (Chterev 2020).The simulations were carried out with a constant time step of Δt = 1 × 10 −7 s .The largest CFL number is monitored instantaneously and limited below 0.25.All simulations were set to develop for about 6 combustor flow-through times based on the bulk flow velocity through the combustion chamber.Statistics were subsequently collected over a further 6 combustor flow-through times.The total computational cost to perform an LES/PDF simulation for each case amounts to about 1,500 CU (1 ARCHER2 CU = 1 Node hour = 128 CPUhs) on the ARCHER2 UK National Supercomputing Service (https://www.archer2.ac.uk).

Flow Field Statistics
For the characterisation of the flow field, the non-reacting case, N1 and N2 are studied in the first place.The time-averaged flow field of N1 is shown in Figs. 3 and 4 as an example.LES predictions are compared with measurements in terms of mean axial (W), radial (U) and azimuthal (V) velocities, and rms axial ( rms w ), radial ( rms u ) and azimuthal ( rms v ) velocities, respectively.For all the images in the present work, the chamber dump plane is defined as z = 0 mm.Contours on a mid-plane of the combustor are shown with radial profiles at four downstream positions ( z1 = 7 mm, z2=15 mm, z3 = 25 mm, z4 = 35 mm) presented for a quantitative comparison.
As can be seen in Fig. 3, there is an artifact in the velocity field measurement out to about 7 mm from the dump plane in the PIV figures from the experiments.This is caused by a strong reflection of the laser at the burner surface and has little influence on the regions of interest.
Although the centerbody is tapered, at the current high swirl number the flow separates along the centerbody resulting in three regions (Meier et al. 2007): the swirling conical jet; the outer recirculation zone (ORZ), which is a toroidal recirculating zone Fig. 3 Comparison of mean axial (W), radial (U) and azimuthal (V) velocities.Top: Velocity contours; Bottom: velocity profiles caused by the sudden area expansion from the nozzle to the chamber and the confinement of the chamber walls; the inner recirculation zone (IRZ), which is created by the swirl-induced vortex breakdown.The comparison of all the contours and profiles shows that the overall flow field is well reproduced.The jet axial velocity reaches a peak ∼ 25 m/s at the nozzle exit, while the highest backflow velocity in the IRZ is ∼ −10 m/s at about z2.The mean velocities are well predicted for all velocity compo- nents including the azimuthal velocity V, which indicates the swirl level in the LES is in good agreement with the measured one.According to the definition of the swirl number, as shown in Eq. 1, the swirl number at the dump plane is found to be S ∼ 1.06 in the LES and remains nearly constant.The width of the IRZ from z = 35 mm to further downstream is slightly under-predicted, which indicates a smaller spreading angle of the swirling jet.The flow spreading angle is strongly affected by the flow separation at the nozzle lip, which requires a very fine near-wall mesh to resolve.The slightly narrower ISL suggests that the flow separation is early in the simulation, resulting in a slight under-prediction of the jet spreading angle.This is a well-known limitation of most LES in combustion chambers but has a limited influence on present work which focuses on combustion dynamics and thermo-acoustic oscillations (Roux et al. 2005).
The rms velocities predicted by LES are also in excellent agreement with the experimental results, which suggests the fluctuations in the flow are well predicted.There are high levels of axial fluctuations ( rms w ) in the IRZ close to the nozzle outlet around the centreline.This is contributed by the turbulent fluctuation as well as the coherent structures in the flows, which is discussed in the analysis of the transient behaviour of the flow field in the next section.The time-averaged flow field (mean and rms) of N2 shows similar results as N1, except for a lower velocity magnitude for all components due to the lower total mass flow rate as listed in Table 1.Overall, the mean flow field and the velocity fluctuations in the non-reacting cases are well-reproduced by the LES.

Dynamic Analysis
In swirling flows as the swirl number increases, the flow undergoes a strong vortex break down, which can induce a helical vortex structure in the vicinity of the ISL as a result of a super-critical Hopf bifurcation (Liang and Maxworthy 2005).These are known as helical instabilities in swirling flows (Gallaire et al. 2006;Oberleithner et al. 2011).An azimuthal wave number m is associated with the helical instabilities, which reflects its angular periodicity.The most common mode is the m =1 mode, (referred to as a PVC), has been found in many swirling flames (Syred 2006;Wang et al. 2005;Stöhr et al. 2011;Moeck et al. 2012;Terhaar et al. 2015).Very rarely, in some instances, a double helix structure associated with m = 2 is observed in experimental studies of swirling flames (Syred 2006;Vigueras-Zuñiga et al. 2012;Vignat et al. 2021); this type of breakdown is highly sensitive to disturbances.In the present work, the terminology 'PVC' is only used to describe the m = 1 helical mode while the m = 2 mode is referred as 'a double helical instability' for clarity.The role of the PVC in isothermal swirling flows is discussed in the reviews of (Lucca-Negro and O'Doherty 2001), and in reacting flows in the reviews of Syred (2006); Candel et al. (2014); Manoharan et al. (2020).The emergence of the PVC in iso-thermal PREC-CINSTA flows was first reported by Roux et al. (2005), where low-pressure iso-surfaces were used to visualise the PVC structure.Large scale vortical structures in the flow are often visualised using low pressure iso-surfaces (Poinsot and Veynante 2005) since coherent structures are found to be associated with local minima of the pressure field (Jeong and Hussain 1995) which are strongly coupled with high levels of velocity fluctuations.Hence in the present work, coherent structures in the simulations are also identified using the same method with low-pressure iso-surfaces.
As shown in Fig. 5a, a helical structure can be found in the vicinity of the ISL, originating from the region in the nozzle below the dump plane.The vortex core grows anticlockwise from the top view (same as the swirl direction), while it precesses clockwise, which is typical for PVC to precess in the same direction as the flow but is winding in the opposite direction (Oberleithner et al. 2011).Figure 5b shows the flow field on a 2D slice perpendicular to the y-axis as shown by the red dashed line in Fig. 5a, where a distinct zig-zag vortex pattern can be found.The PVC features a precessing frequency related to the vortex core precession f pvc which is a global frequency and is typically the same throughout the combustor (Syred 2006).In the present burner f pvc can be deter- mined using velocity probes in the ISL or the chamber acoustic signals (Boxx et al. 2012).Therefore, velocities are monitored at three probes within the IRZ ( w IRZ ), ISL ( w ISL ) and Jet ( w Jet ), with the locations shown in the frame in Fig. 5.
The existing modes in these signals are retrieved via frequency spectra obtained using fast Fourier transformation (FFT).Figure 6 compares the axial velocity flucutaion spectra of w ′ IRZ , w ′ ISL and w ′ JET achieved from the LES and the experiment results in N1.It should be noted that the experimental FFTs are built for a duration of 1 s, while the numerical spectra are calculated over a time of 0.1 s.Despite the shorter physical time in LES, the results have been found to be almost insensitive to the length of duration after about 0.05s.The sampling frequency of the velocity and pressure is 10 kHz and 100 kHz, respectively in the experiments, and 1000 kHz in the LES.Frequencies of interest are below 1000 Hz in the PRECCINSTA configuration (Meier et al. 2007).Hence, the waves of interest are well resolved in both the experiments and in the LES.
The w ′ ISL spectra in Fig. 6 show a distinctive single peak at 479 Hz in the LES and 443 Hz in the experiment, which is also found in the spectra of w ′ JET in the jet region with a relatively lower magnitude but not seen in the w ′ IRZ .This indicates the existence of strong periodic hydrodynamic fluctuations in the ISL and the jet with a fixed frequency, which corresponds to the observed PVC precession.The predicted dominant f pvc is hence 479 Hz in the LES which is close to the measured 443 Hz in the experi- ments (about 8% difference).In isothermal flows, the frequency of the PVC can be read- ily characterised by the swirl number S and a Strouhal number defined as St = fl∕u , where f is the frequency of vortex shedding, l is the characteristic length and u is the flow velocity.Hence, given the accurate simulation of the velocity field, the slight discrepancy in the frequency of the PVC is likely due to a mismatch of the predicted vortex shedding details in the vicinity of ISL, the resolution of which may require a much finer mesh.
Due to the fact that no measurements are available for N2, the LES results of N2 are compared with N1 as shown in Fig. 7.The frequency of the PVC in N2 is found to be 431 Hz compared to the 479 Hz found in N1, which is consistent with the claims by Syred (2006), Stöhr et al. (2012) and Massey et al. (2022) that the PVC frequency increases linearly with the flow rate.In the previous experimental study of the PRECCINSTA burner, a PVC is always found active in the non-reacting flows, while in reacting flows the PVC was found to be suppressed under certain conditions (Roux et al. 2005;Oberleithner et al. 2015).It is also found in the present work that the isothermal flow field is dominated by a consistent PVC in each case investigated, which causes the expansion and contraction of the IRZ at the same frequency.
To summarise, the flow field and hydrodynamic features of the iso-thermal cases are well predicted by the simulations, which suggests the LES framework and sub-models are capable of capturing the complex flow dynamics in this combustor.

Flame Topology and Flow Field
In order to visualise the flame topology and its location, the time-averaged heat release rate (HRR), OH mass fraction ( Y OH ), and axial velocity (W) on the mid-plane are shown in Fig. 8 from left to right.The W fields are overlaid on top by the contours of a binarised Y OH using a threshold value.The threshold value of Y OH is set as 80% of the maximum value in each case, with which the key features of the flames are reserved.From top to bottom Fig. 7 Comparison of axial velocity spectra in N1 and N2 are the results of C1, C2 and C3 respectively.Flame regions are characterised using HRR images and the Y OH field showing the presence of hot burned products.The results present two main types of flames in terms of the flame topology: in C1 the flame features an M shape which is lifted above the centre cone; while in C2 and C3 a V shape flame stabilised along the ISL can be found.ISL can start upstream of the dump plane, and so is ISL stabilised flame in C2 and C3 where the flame is attached to the centre cone and origins below the dump plane.The different flame shapes are mainly caused by the higher in C2 and C3 and hydrogen addition in C3 which increase the extinction stretch rate of the flame base.Hence, C2 and C3 flame can withstand higher hydrodynamic strain in the ISL and are likely to remain attached, compared to C1 flame.The flame lifting process is also closely related to the coherent structure: when the flame cannot withstand the local strain in the ISL, the flame detaches partially from the ISL which is then aided by helical vortex structure until the formation of a detached M-shaped flame (Oberleithner et al. 2015).The flame lift-off process in the present work will be discussed in detail together with the coherent structure in Sect.4.3.3.
Besides, C2 and C3 present a blurry reaction region in the upstream region in the nozzle in C2 and C3, which is absent in C1.This indicates the flame stabilisation location move in time during the observed strong thermo-acoustic oscillations.The flame flashbacks and other transient motions will be discussed further in Sect.4.2.4.Another observation is the increasing HRR level from C1 to C3.For C1 and C2, the equivalence ratio increases from 0.65 to 0.85 leading to higher HRR in C2.For C2 and C3, despite the same , hydrogen enrichment results in a higher flame temperature and HRR due to the higher heat value of the fuel.Besides, with hydrogen addition, the flame becomes more compact and shorter in length and stabilises closer to the combustor nozzle, which suggests a higher thermal loading to the combustor hardware.This is likely due to the increased flame speed, which enables the flame to stabilise in the regions with higher flow speed, i.e. the region close to the nozzle exit.
Figure 9 compares the simulated and measured mean flame topology.The measured averaged OH-CL images are compared with the line-of-sight integrated mean HRR in the LES.It can be observed that the flame form in each condition is well predicted in LES; an M shape lifted flame in C1 and a V shape flame in C2 and C3.The effect of hydrogen enrichment on shortening the flame and making it more compact is found both in the LES and experiment.However, it can be found in Fig. 9, the most active region with the highest HRR signal exists in the region above z = 30 mm in C2 and above z = 20 mm in C3 in the experiment, while in the LES the flame is more active in the region close to the centre bluff body.Moreover in C2 and C3, the flame is found to be longer in the longitudinal direction and spread less in the transverse direction, compared to the OH-CL from the measurements.This is due to the coherent structures found in the LES, which interact with the ISL and significantly modify the flame stabilisation.
The flow fields are also significantly affected by the observed strong thermo-acoustic oscillations and hydrodynamic instabilities.Figures 10 and 11 show the radial profiles of the mean axial velocity W and the rms axial velocity rms w for C1, C2 and C3 at the same downstream locations as presented in the non-reacting cases.For C1, the mean axial velocity magnitude and the jet spread angle are well-reproduced in the LES despite a slight overprediction of the jet velocity.However, in C2 and C3, the difference between the measurement and simulation is obvious: LES predicts a narrower jet and higher jet velocity, as well as a higher reverse flow velocity in the IRZ.The predicted rms velocities for the pulsating flames with the PRECCINSTA burner are checked numerically for the first time, which are shown in Fig. 11.In the experiments, the rms velocity in C1 is found to be the highest compared to those in C2 and C3, which is consistent with the thermo-acoustic features detected as discussed in Sect. 2. In LES, strong axial velocity fluctuations are found in all three cases, especially in C2 and C3 where rms w can reach as high as 40 m/s.In each LES case, the IRZ exhibits the most violent fluctuations, followed by the jet region.The flows are not only subjected to turbulent fluctuations but to the coherent structures and strong thermo-acoustic oscillations, which all contribute to the high rms velocities in all reacting cases investigated.The discrepancy between the predicted and measured flow fields is consistent with the flame topology comparison as shown in Fig. 9, which are related to the flame dynamics and will be discussed in Sect.4.3.3.

Limit-Cycle Oscillations
Strong flame oscillations are observed in the LES, which has a significant impact on the flame topology and velocity moments as shown above.Previous studies on the PRECCI-NSTA burner identified the main driving mechanism of the thermo-acoustic instabilities for this geometry is the convective coupling between the Helmholtz modes of the burner and the heat release via mass flow and equivalence ratio fluctuations (Meier et al. 2007;Franzelli et al. 2012).The air supply line has high acoustic impedance but the pressure in the large plenum volume responds to the acoustics resulting in air mass flow fluctuations.In the TP mode when air and fuel are injected separately, the high impedance of the fuel line leads to fuel accumulation in the swirler, creating equivalence ratio disturbances subject to the convective delay.A similar driving mechanism of thermo-acoustic instabilities is detected in the currently investigated cases, with detailed discussions provided in this section.
The dynamic behaviour of the three cases is monitored using the signals of the pressure in the plenum p ′ pl , chamber p ′ ch and the axial velocity at the ISL probe w ′ ISL as done already in the discussion of non-reacting cases.The air mass flow rate ṁtot and fuel mass flow rate ṁfuel at the dump plane are also monitored.In addition, the volume-integrated (or global) HRR in each case is monitored.In the system stability analysis, a widely used approach to determine whether the coupling between the acoustic field and the unsteady reaction rate leads to instability is the global Rayleigh criterion (Rayleigh 1878) defined as: where p ′ and ̇q′ are the local pressure and HRR fluctuations, and T and Ω are the time dura- tion and volume of interest (flow domain in the current case).Equation 4states that the relative phasing of the acoustic pressure and the HRR fluctuation would determine whether positive feedback exists in the thermo-acoustic cycle.
Figure 12 shows the simulated signals of the pressure fluctuations of p ′ pl and p ′ ch , the HRR, the mass flow rates ṁtot and ṁfuel , and the axial velocity at P2 of C1.A period of 4 ms contains about 10 oscillation cycles.It can be found that a period of 4 ms contains about 10 oscillation cycles during which all the signals oscillate with nearly constant amplitudes.The amplitudes of all the signals are nearly constant, indicating the state of zero growth rate of the acoustic energy: the acoustic energy source from the unsteady heat release and (4) the acoustic energy losses are balanced (Poinsot and Veynante 2005).This is commonly referred to as limit-cycle oscillation and is controlled by a variety of complex, highly nonlinear physical mechanisms introduced in Sect. 1. Limit-cycle oscillations are also found in the C2 and C3 cases, suggesting all three reacting cases in the LES are subject to thermoacoustic instabilities.In contrast to the measurements where C3 is found to be stable, it

Self-Excited Instabilities
Fig. 13 compares the measured and computed PSD of the chamber pressure and HRR in all three cases.The duration and resolution of the signals are kept the same as those for the iso-thermal case.It is worth noting that in the experiments, the acoustic mode peaks are significantly damped and broadened by various effects including acoustic damping due to loosely mounted and not completely air-tight optical windows (Chterev and Boxx 2021).
Hence, the flame is labeled as thermo-acoustically unstable if there is a distinct peak in the heat release spectrum.During self-excited thermo-acoustic instabilities, one or more acoustic modes can be excited and the heat release may couple to one or more of these modes.
In C1 for the experiment p ′ ch and HRR spectra, the highest PSD is found at a frequency of around 224 Hz, which represents the dominant thermo-acoustic oscillations.The second peak at about the first harmonic frequency of 449 Hz is detected, which is also found in the HRR PSD, indicating its thermo-acoustic nature.In LES, these two thermo-acoustic modes are also found at around 240 Hz and 480 Hz respectively, with additional harmonics at multiples of the dominant mode frequency.The peak amplitudes measured to be 132 dB and 116 dB, respectively, are denoted by the blue star marks.
In a recent study, (Lourier et al. 2017), the loosely mounted quartz windows used in the experiment as the chamber side walls have been found to cause acoustic damping on the chamber pressure PSD by about 10 dB with the frequencies unaltered.With this estimation, the two dominant mode peaks are marked by red stars, and good agreements are achieved in the predicted and measured peak amplitudes (within 3 dB difference).One possible reason is acoustic damping, and as discussed by (Lourier et al. 2017), the main contributor to damping is the loss of acoustic energy at the boundaries of the system.The previous estimation of the damping effect of the quartz windows is obtained from the oscillating flame of Meier et al. (2007) which features an M shape and resembles the C1 flame.However, the cases C2 and C3 exhibit different flame configurations and dynamic features, and the effect of the acoustic damping caused by the quartz windows remains unclear.
In the PSD of the experiment p ′ ch in C2 there are also two dominant modes, with frequencies of 303 Hz and 607 Hz respectively; in the PSD of the experiment HRR, only the first mode is detected, indicating the thermo-acoustic nature of this dominant mode.In the LES, the frequencies of these two modes are predicted to be 305 Hz and 610 Hz, which are in excellent agreement with the measured ones.As for the PSD of HRR, a dominant peak at the first mode frequency is found, with harmonics existing at its multiples (610 Hz and 915 Hz, etc).Over predictions are found for the amplitudes by about 10 dB for the first and second peaks in p ′ ch PSD with the estimated damping effect taken into account.The discrepancies become more pronounced in the C3 case, where the dominant peak is found to be an acoustic mode in the experiment but a thermo-acoustic mode in the LES.That is, in C3 the HRR fluctuations are coupled with the acoustic fluctuations in the LES, which is not found in the experiment.Besides, an underestimation of about 22 Hz of the acoustic frequency is found in the LES.
These thermo-acoustic instabilities are highly impacted by thermal boundaries and initial conditions.In general, the thermo-acoustic frequency is determined by the natural frequency of the active acoustic mode of the system, which is related to the geometry and speed of sound (Chterev and Boxx 2021).The geometry is kept unchanged in the current study leaving the speed of sound the key parameter, which is closely related to the temperature field.As discussed in Sect.3, the wall temperature and the inlet air temperature are fixed according to the available measurements, which are the same for all three cases.However, these temperatures should definitely change across cases with different heat release distributions in the burner with different flame topologies.Figure 14a compares the time-averaged temperature field of C1 and C3, and (b) the corresponding radial profiles along the bottom of the burner at z = 0.5 mm and two upstream positions in the nozzle at z = −4 mm and z = −15 mm.Due to the axial symmetry, only half of the profiles are pre- sented.It can be found that the temperature field presents a significant difference in the two cases: in C1 when the flame is lifted, the temperature in the region near the bluff body is around 800 K, which reduces gradually upstream and reaches about 320 K at z = −15 mm; while for the attached V-flame in C3, a temperature ∼ 2000 K is reached near the centre cone tip on the internal side.This is not unexpected, as a high value of OH mass fraction is found in the region just downstream of the center body tip in C2 and C3 as shown in Fig. 8.The region in the upstream nozzle is also greatly heated in C3 which is related to the flame flashbacks.Hence it can be expected, the different flame structures lead to differences in the solid wall temperatures especially for the centre body, which make the same prescribed wall temperatures in all cases not accurate.
The slight mismatch of the measured and predicted acoustic frequencies is closely related to these imposed thermal boundaries.The various near-wall temperature fields shown above would definitely lead to different wall temperatures, which in turn, control the heat transfer between the flow and the solid walls and as a result, affect the preheated temperature of the inlet gas (instead of being fixed at 320K in all cases).These changes in the temperature field modify the local speed of sound (Poinsot 2017), which affects the acoustic frequency directly.The preheated gas temperature also affects the laminar flame speed, thereby influencing the flame topology and anchoring position, which is also directly related to the wall-heat transfer, adding to the complexity of the problem.
The overestimation of the fluctuation amplitudes and the thermo-acoustic instabilities found in LES for C3 involves more complex reasons.One possible reason is acoustic damping, and as discussed by Lourier et al. (2017), the main contributor to damping is the loss of acoustic energy at the boundaries of the system.The previous estimation of the damping effect of the quartz windows is achieved for the oscillating flame by Meier et al. (2007) which features an M shape and resembles the C1 flame.However, C2 and C3 exhibit different flame configurations and dynamic features, and the effect of the acoustic damping caused by the quartz windows remains unclear.In LES, the outlet boundary conditions are treated with the NSCBC formulation as discussed in Sect. 3 with a prescribed fixed pressure relaxation factor based on the optimal value found in previous work.This coefficient controls the 'stiffness' of the boundary which may explain LES/EXP discrepancies.Hence to better represent the acoustic impedance at the exit of the combustion chamber, in a separate study the atmosphere is included in the computational domain with a large extended domain.However, despite a slight reduction in the pressure amplitude, this does not suppress the observed thermo-acoustic mode.Hence an underestimation of the acoustic damping may not be the reason for the thermo-acoustic oscillations observed in the LES.The second and more important reason is again related to the fixed wall temperature as the thermal boundary conditions.
In a recent LES study on the impact of wall heat transfer on the flame dynamic in the PRECCINSTA burner by Agostinelli et al. (2021), the simulation using isothermal boundary conditions with fixed wall temperatures were found to predict thermo-acoustic oscillations which are not observed experimentally.The numerically predicted instabilities were found to be supported by equivalence ratio fluctuations caused by fuel pockets periodically pushed into the chamber, and could not be suppressed by changing the acoustic impedance of the inlet boundary conditions.Instead, they found, with the conjugate heat transfer (CHT) approach where the wall temperature was solved and coupled to the flow computation, the flame response to the disturbances was slightly but sufficiently modified and the flame was found to be stable as in the experiments.Improvement in the prediction of combustion instabilities with the CHT method has also been claimed by Kraus et al. (2018) recently in an LES study of a swirling flame, where the frequency and amplitude of the oscillations were better predicted due to the better computed preheating temperature of the gases in the plenum.
Figure 15 further compares the predicted thermo-acoustic modes in C1, C2 and C3, in order to show the effect of increasing equivalence ratio and hydrogen addition.When comparing C1 and C2, where increases from 0.65 to 0.85 with a slight shift of the thermal power, the dominant frequency increases in the LES.C2 predicts higher peak amplitude than C1, which contrasts the experiment results as shown in Fig. 13.This is due to different flame shapes (hence HRR distribution) and HRR levels, resulting in significant differences in heat transfer distributions to the combustion hardware (centre body, combustor liner, dump plane) (Chterev et al. 2014).In V-shape flames with the flame root attached to the center body tip, the regions where the flame tips reach the wall are very hot as shown in Fig. 14; while in reality and in the CHT study by Agostinelli et al. (2021), the temperature field and hence the sound speed in that region are more homogeneous due to the heat conduction in the solid, and a globally lower sound speed is found.This non-homogeneous heat transfer with fixed wall temperature was found to have been closely related to the numerically predicted thermo-acoustic instabilities (Agostinelli et al. 2021), which is less significant in the M-shape C1 flame where the flame is lifted and the temperature field in the region near the centre body is already more homogeneous.
The effect of H 2 addition under the condition = 0.85 and p th = 22.8 kW is shown when comparing C2 and C3.As shown before by the equation c = √ RT , H 2 enrichment increases the flow temperature, leading to a higher speed of sound and a higher thermoacoustic frequency in C3 compared to C2.Moreover, with more H 2 in the fuel, there will be more H 2 O in the hot product in which acoustic wave propagates faster compared in CO 2 due to a smaller molar mass.Besides, as shown in Fig. 15, H 2 addition decreases the acous- tic and the HRR peak amplitude by about 30%.This trend is in good agreement with the experimental findings that the H 2 addition damps the acoustic fluctuations under these con- ditions.This is closely related to the driving mechanism of the thermo-acoustic instabilities in the present work (details shown in the next section).As introduced at the beginning of this section, the convective delay between this equivalence ratio disturbance being created and reaching the flame determines the phase difference between the acoustic oscillations and HRR oscillations, which according to the Rayleigh criterion, determines the coupling between the acoustic and HRR oscillation, e.g.phase difference must be within 90 • for a thermo-acoustically unstable system, and a smaller value leads to a stronger coupling.This time delay is a function of parameters such as flame stand-off location, flame length, and speed of sound (Chterev et al. 2014).As discussed in the recent experimental work by (Chterev and Boxx 2021), H 2 enrichment was found to increase the phase delay between the pressure and heat release at the flame location.This is caused by decreasing the pressure wave travel times in the plenum and swirler sections of the burner and by decreasing the flame length while simultaneously increasing the pressure wave speed within the burner, and the balance between the thermo-acoustic period and the total travel time to the flame leads to a longer phase delay.Hence, according to the Rayleigh criterion, the thermoacoustic amplitude is decreased in C3 when the initial phase delay is positive in C2.

Thermo-Acoustic Feedback Loop
The driving mechanism and feedback loop of the thermo-acoustic instabilities are discussed in this section, together with the temporal evolution of the flow and reaction field.
Figure 16a compares the phase-resolved p ′ pl , p ′ ch and the integrated HRR in a full period of oscillation in C1, C2 and C3.The pressure drop between p ′ pl and p ′ ch is computed as p diff = p � pl − p � ch and plotted in Fig. 16b together with ṁtot and ṁfuel simultaneously.These mass flow rates are evaluated in the dump plane at z = 0 mm.The equivalence ratio of the unburned gas ch is computed according to ṁtot and ṁfuel and present for the interval when ṁtot is positive in Fig. 16c.It is worth noting when the flame flashbacks with the thermo- acoustic oscillations as shown in Fig. 17, the equivalence ratio of the reactant mixture can not be evaluated at the z = 0 mm plane.Hence instead, ch is computed as the plane z = −15 mm to monitor the fuel accumulation effect in all cases.
= 0 • is defined at the maximum of the HRR, which is used throughout the rest of the analysis for the reacting cases.The corresponding flow evolutions are shown in Fig. 18 with the instantaneous axial velocity, mixture fraction, and HRR at different phases angles from = 0 • to = 360 • .Several observations can be made when comparing these signals.First, as shown in Fig. 16a, the phase difference between the p ′ ch and p ′ pl peaks in the three cases are roughly 56 • , 59 • and 63 • respectively, corresponding to a propagation time delay of 0.65 s, 0.55 s, and 0.51 s.With the fixed positions of the acoustic probes, this decline in  16b.A phase lag between the peak of ṁtot and ṁfuel can be found, which reduces from 31 • in C1 to 14 • in C3.This phase lag is mainly caused by the different acoustic impedances of the fuel and air supply lines, which lead to fuel accumulation in the swirler unit.The resulting ch fluctuates periodically in all three cases, as shown in Fig. 16c.It can be found, in C1 ṁfuel as well as ch exhibits rapid growth and decline, while in C2 and C3 ṁfuel and ch drop mildly in a similar way as the total mass flow rate ṁtot .ch in C1 fluctuates between a wide range from 0.38 to 1.3 (over Δ ch = 0.92 ), while this range difference reduces to Δ ch = 0.4 in C2 and C3.This indicates that in C1 fuel accumulation has a distinct and more pronounced effect.Besides, the mass flow rates in C2 and C3 exhibit a nearly constant level, before a sudden increase at around 180 • .This is caused by the violent flashback of the flame in C2 and C3 as can be seen in Fig. 18 at 90 • − 135 • for C2 and C3.As a con- sequence, the reactant mixture is consumed before reaching the z = 0 mm plane.
Figure 18 shows the instantaneous axial velocity and mixture fraction at different phases from 0 • to 360 • .Iso-contours of a constant-value HRR are overlapped on top to show the flame region in the chamber.In the swirler inlet region, the slices are taken along the passage in order to show the flow locally as illustrated in Fig. 18 on the top.In C1 at 0 • when the HRR reaches its maximum, the local high HRR generates pressure disturbances (acting like an acoustic source).Due to the propagation delay, p ′ ch increases ahead of p ′ pl , and hence p diff has a negative value, leading to the drop of w in the nozzle and thereby less sup- ply of fresh reactants into the flame region.Hence the HRR drops afterwards ( 0 • to 180 • ).With the increase of p ′ pl , p diff grows and turns positive at around = 60 • , after which the axial velocity in the nozzle begins to increase ( 135 • to 180 • ).The velocity within the fuel channels where fuel is injected into the swirler, by contrast, remains relatively constant due to the high impedance at the fuel inlet.A significant amount of fuel accumulates within the swirler ( 90 • to 135 • ) when the nozzle velocity is declining.When the axial velocity grows, the fuel-rich pockets of gas are convected downstream and reach the chamber inlet at around 180 • .With these large quantities of reactants transferred into the reaction region, combustion is greatly enhanced as shown by the climbing HRR at 270 • , which continues to increase until the next maximum is reached when the feedback loop is closed.
A similar oscillation behaviour is found in C2 and C3, with some notable differences.When the flame flashes back ( 90 • ), in C1 the flame reaches the tip of the nozzle, while  in C2 it extends into the centre body halfway, and in C3 is even further upstream.This is mainly caused by the increasing flame speed from C1 to C3. Besides, the flame in C2 and C3 remains attached to the centre body throughout the cycle, while the flame is lifted off in a fairly long range in the C1 cycle.The flame lift-off is closely related to the hydrodynamic structure in the flow and will be discussed in detail in the following section.
In summary, a combination of mass flow rate fluctuations and equivalence ratio fluctuations is identified as the dominant feedback mechanism for thermo-acoustic fluctuations in all three cases.The flame experiences violent, periodic oscillations in the axial direction, correlated with thermo-acoustic fluctuations.A similar oscillating behaviour of the flame was reported in the experimental work by Temme et al. (2014) where a swirl-stabilised, liquid fuel gas turbine combustor featured strong self-excited combustion instabilities driven by a comparable equivalence ratio oscillation mechanism.Tachibana et al. (2015) also identified periodical flame flashback into the pre-mixing annulus passage in a liquidfuel aero-engine combustor both experimentally and numerically, as a result of different responses of the fuel and airflow rates to the pressure fluctuations.In the recent LES study by Fredrich et al. (2021a) of the F2 flame periodic flame lift-off and flashback were reported together with fuel supply oscillations in the system.

Hydrodynamic Instabilities
In this section, the hydrodynamic instabilities resulting from flame-vortex interaction are investigated.The formation and evolution of the vortical structures and their potential impact on the flame are also discussed.

Swirl Number Fluctuation
Periodic excitation and suppression of large-scale vortical structures have been detected in all three cases.This is linked to the transient interaction between the swirler and the strong acoustic fluctuations.
The mechanism describing the interaction between the swirler and incident acoustic perturbations was characterised analytically by Palies et al. (2010).As illustrated in Fig. 19 when a swirler unit is impinged by acoustic waves under low Mach number conditions, axial ( w ′ ) and azimuthal ( v ′ ) velocity perturbations are induced.v ′ is found to have the same order of magnitude as the incoming acoustic disturbance (Palies et al. 2010), and is caused by vorticity waves ( ) generated at the trailing edge of the swirler blades which is subject to a mode conversion process (Cumpsty et al. 1977;Komarek and Polifke 2010).The vorticity wave has a convective nature and propagates downstream at the bulk flow velocity.These velocities waves impact the system mainly in two ways: the axial velocity fluctuations solely can generate toroidal vortices at the dump plane, rolling up the flame surface and thereby modulating the overall HRR in the burner; the combined effect of the axial and azimuthal velocities fluctuations gives rise to the fluctuations of the swirl number.The swirl number controls the motion and strength of the IRZ, modifies the flame spreading angle and the distance between the nozzle and the flame front, and controls the formation, geometry, and frequency of the PVC, hence having a significant impact on the flame dynamics.Periodic PVC deformation was reported experimentally by Caux-Brisebois et al. ( 2014) in a PP PRECCINSTA flame and was identified to be induced by swirl number fluctuations at the thermo-acoustic frequency.Fredrich et al. (2021a) first identified a periodic excitation and suppression of a PVC and flame angle oscillations induced by Fig. 19 The interaction between the swirler and incident acoustic perturbations.The blues arrows represent the acoustic waves propagating at the speed of sound c, and the black arrows denote the jet flow rate at which the vorticity waves travel downstream the swirl number fluctuations in the F2 flame, which were found to be directly linked to the dominant thermo-acoustic mode.

Temporal Evolution of Vortical Structures
Helical vortex structures are observed with both M-and V-shape flames in C1 to C3.The helical vortex existence was found to be strongly dependent on the density gradient in the wave-maker region (the separating ISL from the center cone) (Oberleithner et al. 2015), which is hence closely related to the heat transfer at the centre body surface and the preheated air temperature.In previous studies with the PRECCINSTA burner, a PVC was found only in the M-flame (Steinberg et al. 2011;Oberleithner et al. 2015).In the recent experiments by Chterev and Boxx (2021) with elevated pressure and preheated air from 613K to 663K, PVC existed in all reacting cases with M-and V-flames.
Figure 20a shows the axial ( w ′ ) and azimuthal ( v ′ ) velocity fluctuation in one acoustic cycle in C1.The amplitude of the v ′ wave is about half of the w ′ wave, with a phase lag of ∼ 33 • corresponding to 0.4 ms.This phase lag together with the different wave amplitude leads to fluctuating momentum flow rates G z and G y , and hence a fluctuating swirl number S as computed by Eq. 1, which are shown in Fig. 20b.The resulting flow field evolution is shown in Fig. 21.Longitudinal and transverse planes of HRR field are visualised which are overlapped with pressure contours to show the interaction of the flame and the vortical structures.The red dashed lines in the longitudinal planes mark the position of the transverse planes.On the bottom, the 3D vortical structures are shown by different levels of low-value pressure iso-surfaces: blue contours are for the vortex ring (VR); green contours describe the vortex cores (VC) and the edge of the vortex structure; red contours represent the kernel (core) of the vortex with a lower value pressure iso-surface.
From 0 • to 108 • , G z and G y are both low due to a low mass flow rate ṁtot , and no distinct hydrodynamic structures exist in the field.At 126 • , a vortex ring forms at the outer rim of the nozzle at the chamber inlet as a result of the increased G z , which is convected down- stream along the OSL before breaking down from 162 • and vanishing afterward.As the swirl number S grows, a central vortex core (CVC) is excited at around 162 • as a result of the rapidly increasing azimuthal and axial velocities.PVC with a single helical vortex is then excited along the ISL with a high instantaneous S at 180 • , leading to a characteristic deformation of the IRZ and the jet.The helical vortex rotates anti-clockwise from the top view while precessing clockwise and continuously growing in size from 180 • to 270 • .The resulting flow field and flame topology show great asymmetry, and the flame is lifted off when the PVC is excited ( 198 • to 360 • ), which will be discussed further in the following section.It is difficult to identify the specific event that triggers the formation of the PVC since the flow is at the border of instability and thus even minor turbulent fluctuations of velocity or density may initiate the formation.The single helix starts to break down when S drops and finally disappears at the end of the cycle.The existence of the PVC is closely related to the instantaneous swirl number and is found to exist when S is higher than a geometrically dependent critical value, which is found to be S = 0.6 in the work by Chigier and Chervinsky (1967) and 0.55 in the DRL dual swirl.This critical value of S is not examined in the current work, however, it can be readily found from Fig. 20b and the temporal field evolution that the CVC is excited at around 162 • where S is very close to the value 0.6 as marked by a red dashed line in Fig. 20b This periodic evolution of the vortices is observed in every thermo-acoustic cycle, strongly coupled with the dominant thermo-acoustic mode.
The PVC with a single helical vortex structure discussed above is observed to be the most common structure in C1 and is most often present in each cycle.However, a secondary helical structure with double helices is also detected in some of the cycles, which is shown in Fig. 22.The corresponding swirl number S fluctuations do not show a distinct difference from the one during a single helix evolution, and the excitation of this double helical structure occurs in a random manner.Unlike the single helical structure where a zig-zag vortex pattern is presented as shown in Fig. 21, this double helical structure leads Fig. 22 A double helical vortex cycle in C1 to a nearly axisymmetric distribution of the local flow field and the flame topology on a 2D vertical plane.Another notable difference is instead of being lifted off, the flame is attached when the double helical vortices are excited ( 198 • and 252 • ).A central vortex core forms at 180 • with two distinct branches, which both co-rotate around the centreline clockwise and grow into two helices ( 198 • and 252 • ).This double helical structure grows along the ISL and gains in size, before collapsing into a merged vortex core ( 360 • ), when the flame is slightly lifted with this single helix structure, similar to the flame detachment shown in Fig. 21.Finally, the vortices break into small eddies and vanish ( 18 • ) as a result of the declining momentum fluxes and S.
In C2 and C3, a periodic evolution of the vortical structures linked to velocity fluctuations is also detected with a representative cycle as shown in Fig. 23.The phase difference between the v ′ wave and the w ′ wave increases to ∼ 44 • , and the peak amplitudes of G z and G y both reduce due to a lower mass flow rate.The resulting swirl number S fluctuates between 0.1 to 0.9 and remains above 0.6 during about 2/3 of the cycle, suggesting a less oscillating hydrodynamic field and longer existing time of the helical structure compared to C1.In C2 and C3, instead of being bi-stable, a double helical structure is found to be excited and suppressed in every thermo-acoustic cycle, with the structure similar to the one displayed in Fig. 22 for C1.
The evolution of vortical structures in the cycle shown in Fig. 23 is presented in Fig. 24.A central vortex core forms at 180 • , which grows into two branches both co-rotating around the centreline of the burner, with the flame remaining attached throughout the cycle.The double helical structure rotates along the ISL before breaking down as shown at 360 • .Compared to the double helices detected in C1, this double helical vortices in C2 expand less in the radial direction, which is likely due to the lower swirl number magnitude when the helical structure is excited as shown in Fig. 23.
As discussed in Sect.4.1, helical instabilities are characterised by the wave number m, and m = 1 mode is the most common mode which is known as PVC.A recent experimental study by Vanierschot et al. (2020) identified the m = 2 mode as an independent global mode instead of a second harmonic of the m = 1 mode in a swirling iso-thermal flow.Vignat et al. 2021 investigated the transient PVC dynamics in a swirled spray flame under stable and unstable conditions.In the stable flame, the helical instability took the form of Fig. 23 Spatially integrated variables at the dump plane in C2 a single helix with brief intermittent switching to a double, followed by a triple helical (m = 3) geometry, and these modes corresponded to flashback events that occur in a random way.In the unstable flame, the double helix was found to be the dominant PVC mode, which was modulated at the thermo-acoustic frequency due to the oscillating mass flow rate at the injector.Only very few numerical works have reported the double helix structure.Zhang and Vanierschot (2021) used incompressible LES to predict single and double PVC in the iso-thermal flow investigated by Vanierschot et al. (2020).Mansouri et al. (2016) studied a double helical vortex structure and its interaction with a premixed propane flame in a swirl-stabilised burner using Detached Eddy Simulation.The vortex core was found to be characterised by solid-body rotation, which contains only fresh gases with no flame detected in the core region.
In a recent experimental work of a perfectly premixed PP PRECCINSTA flame, Yin and Stöhr (2020) used multi-resolution POD (MRPOD) to study the PVC and thermo-acoustic instabilities in a bi-stable flame, where the flame alternates irregularly between an M-shape and a V-shape.The PVC-related dynamics were found to have spatial and temporal bimodal behaviours: the dominant PVC was observed to only prevail during M-shape flame, while its harmonics, 2PVC and 3PVC, exist in a bi-modal trend: they could possess either a single or double helical structure and they could be active during either V or M flame periods with the same characteristic frequencies.Although under different operating conditions, these experimentally observed bi-modal PVCs can be compared to the one predicted here in C1 with the flame shape and position periodically oscillating.In addition, in the present work, the mode of the helical instabilities is found to be related to the flame shape and fuel composition.In C2 and C3, a V-shape flame exists with a higher equivalence ratio = 0.85 and hydrogen addition in C3, the single PVC mode is never observed and the flame is always attached to the centre body.This can be attributed to the fact that the higher equivalence ratio and hydrogen addition increases the extinction stretch rate of the lean premixed flame (Taamallah et al. 2015b), which determines the localised flame extinction in the turbulent combustion (Sarli and Benedetto 2007).A recent experimental study (Datta et al. 2021) on hydrogen addition and PVC in TP PRECCINSTA flames shows H 2 enrichment results in the suppression of flame lift-off and the PVC, which is consistent with this analysis.Stochastic extinction near the flame base reduces the radial density gradient, which was found to be a favorable condition for the formation of PVC (Oberleithner et al. 2015).The C3 flame is thereby less likely to blow off at the flame base, leading to an unfavorable condition for the single PVC.However, in the simulation, instead of being suppressed, the helical structure presents a double-helix mode.Since the helical vortex formation was found to be strongly dependent on the density gradient region near the center cone, it is, again, closely related to the temperature distribution there and hence is subject to the complex interaction of the thermo-acoustic and wall heat transfer interaction.The determinant for the formation of the double helical structure remains unclear and a rigorous analysis of it is beyond the scope of the present work.

Helical Vortex Structure and Flame Stabilisation
As shown above in Figs.21, 22 and 24, flame attachment and detachment are closely correlated to the helical structure in the flow and will be discussed in detail in this section.
The dynamics of the single PVC formation and its interaction with the flame is illustrated in Fig. 25.The regions where the vortices interact with the flame are highlighted with dashed boxes and the enlarged views are shown on the right.First, as described in Sect.4.3.2, the azimuthal momentum flux increases followed by the formation of a single helix structure as shown at 198 • with a typical zig-zag vortex pattern.The flame moves upstream as a result of the high-speed reverse flow, deflected by the vortices, before colliding with the incoming flow on the rhs.A stagnation point (SP) and a stagnation line occur in the flow, as marked by the black dots and black dashed lines, where the stagnation line represents the zone where the incoming flow and the recirculating flow from the IRZ collide.The SP blocks the reverse flow of hot burned gas toward the nozzle and thus stabilises the flame base near this point.When the helical vortex grows and rotates around the central axis, the vortices move downstream in the 2D plane, and another vortex appears on the right-hand side as shown at 216 • .The SP also moves downstream due to its correlation with the vortex motion, resulting in the reaction zone being lifted.In a recent experimental study on a bi-stable flame by Stöhr et al. (2012); Oberleithner et al. (2015); Stöhr et al. (2018), a single PVC was also found to cause the flame detachment by creating an unsteady lower stagnation point whose dynamics coupled to the vortex motions.As described in Stöhr et al. (2012), the unsteady SP is an important feature of the PVC as it determines the position of the flame root.At 234 • when the helical vortex further grows, the flow field becomes more complex due to the vortex motions, and an additional stagnation point is found at the top rhs of the enlarged view.The local collision of the opposed flows of burned and unburned gas with the thin flame brush between them forms a local turbulent counter-flow flame (TCF) near the two SPs.Along these stagnation lines, the burned gas supplies heat and radicals to the fresh gas and which are favourable for local re-ignitions.
Figure 26 shows the temporal evolution of a double helical vortices and the resulting flame and flow field.At 180 • , a pair of symmetric vortices form above the dump plane.Stöhr et al. (2018) detected a symmetric vortex pair above the inlet in a transient flame in 2D PIV data and identified that they were caused by the pumping motion of the TA  and 24, these symmetric vortices correspond to a double helix structure in the present work.The double helices grow along with the formation of a continuous reverse flow in the IRZ, causing the flame to move upstream ( 198 • ).Instead of being deflected towards one side by the single PVC, the flame propagates axisymmetrically along the centreline due to the symmetry of the flow field.As a consequence, a local stagnation point is absent above the dump plane and the flame is further convected upstream and touches the centre cone at 234 • .As one of the helices breaks ( 360 • ), the flame is lifted as a result of the remained sin- gle helix as a result of the SP begins to form at the lower right corner.In C2 and C3 where a double helical structure is repeatedly excited in every thermo-acoustic cycle, the flame attachment is enhanced by the back-flow in the IRZ and the vortex motions coupled to the double-helix structure.

Combined Effects
In addition to the helical structure effect on the flame stabilisation, other interactions of the helical vortex and the flame can also be detected.V-shape flames stabilise along the ISL, and helical vortex interacts with and distorts the ISL significantly, hence PVC contributes to the observed narrower flame as shown in Figs. 8 and 9.In addition, by changing the flame topology and position as described above with single or double helical vortex, the flame length and anchoring position changes accordingly, which, as discussed in Sect.4.2.3, impacts the phase difference between the acoustic oscillations and HRR oscillations that determines the thermo-acoustic intensity.
The helical structure also changes the unsteady heat release locally.As shown in Figs. 25 and 26, the flame front is locally rolled up and elongated by the large-scale vortices, thereby the helical instabilities directly modulate the overall flame 'surface area' and generates unsteady heat release (Fredrich et al. 2021a).The helical instabilities affect the combustion by enhancing the local mixing of the fuel and oxidiser and by enhanced mixing of burned and unburned gas as shown in Figs. 25 and 26, as reported in numerous studies (Stöhr et al. 2011;Froud et al. 1995;Stöhr et al. 2015).Stöhr et al. (2012), PVC could affect the location of the regions of high compressive strain rates which highly coincided with the reaction zone.It has also been found that the combustion is enhanced along the helix edge in some instances as shown in Figs. 25 and 26, and explains why the integrated HRR shown in Fig. 9 has a high value in the region along ISL close to the centre body, where PVC affects most.
It is now clear that both single and double helical structures are coupled with the local HRR via the mechanism discussed above.Prior studies on consistently active PVC have suggested that PVC precession around the burner did not contribute to fluctuations of the global HRR: the associated fluctuations in local heat release have previously been shown to cancel out in space over the oscillation period Moeck et al. (2012).However, when the PVC is periodically modulated due to the swirl number fluctuation, a phase-specific influence on the global HRR was observed (Caux-Brisebois et al. 2014;Fredrich et al. 2021a).The flow-flame interactions associated with deformation or excitation of the PVC at the thermo-acoustic frequency couple the HRR fluctuations with the acoustic oscillations.The change in total HRR over the thermo-acoustic cycle is not due to the PVC's precession around the burner, which occurred at a different frequency, but rather at the thermo-acoustic frequency.The hydrodynamic evolutions amplify the global instability by increasing HRR in phase with the dominant thermo-acoustic fluctuations, thereby feeding acoustic energy into the system.Hence, the feedback loop of the self-sustained strong thermoacoustic fluctuations in all the cases investigated in the present work is a combined effect of the equivalence ratio fluctuations and the hydrodynamic instabilities.
Figure 27 briefly summarises the interplay of the thermo-acoustic oscillations, helical vortex structure, and the flame stabilisation process, with the pre-mentioned wall transfer effect on them in the LES study.Thermo-acoustic instabilities interact with the swirl unit, periodically exciting or suppressing the helical vortex structure, which, in turn, increases HRR in phase with the thermo-acoustic oscillations, hence amplifying the thermo-acoustic fluctuations.Single or double helical vortices determine the flame stabilisation, leading to different flame topology and thereby changing the phase difference between the acoustic oscillations and HRR oscillations that impacts the intensity of thermo-acoustic instabilities.Wall heat transfer directly impacts the temperature field in the near wall region, which impacts the flame stabilisation; it affects the density field close to the center body, which plays an important role in the helical instabilities.Finally, fixing the wall temperature was found closely related to the simulated thermo-acoustic instabilities with high non-homogeneous heat transfer near the centre body, which results in the strong helical vortex modulation which closes the feedback loop.

Conclusion
The self-excited combustion dynamics in swirling flames under different operating conditions are investigated using LES.The targeted technically premixed PRECCINSTA flames have previously been shown to exhibit combustion-driven instabilities under various operating conditions.In the iso-thermal flows, the simulated and measured velocities are overall in good agreement.A continuously active single PVC is observed in all the iso-thermal flows, and the predicted PVC frequencies are in good agreement with the value observed experimentally.In the reacting cases, a lifted M shape flame is detected in C1 while an Fig. 27 The interplay of the thermo-acoustic oscillations, helical structure, and flame stabilisation, with the wall-heat transfer effect on them potentially attached V shape flame is found in C2 and C3, which are consistent with the experimental observations.Self-sustained, limit cycle oscillations are observed in all three cases, which are identified as thermo-acoustic fluctuations.The predicted dominant acoustic frequencies are in overall good agreement with the measured ones, while the fluctuation amplitudes are over-estimated and strong coupling between the heat release rate fluctuations and the acoustic fluctuations are found in cases that are experimentally observed stable.These numerically predicted thermo-acoustic oscillations are closely related to the predicted wall heat transfer using fixed wall temperatures.The hydrogen addition is found to affect the thermo-acoustic instabilities by affecting the flame shape and the heat release level, leading to a higher dominant frequency and lower peak amplitude.
The driving mechanism of the observed thermo-acoustic fluctuations is identified as a combination of equivalence ratio fluctuations and mass flow rate fluctuations in all reacting cases (C1-C3).The acoustic waves interact with the burner unit, generating swirl number fluctuations which leads to a periodic excitation and suppression of the helical structure.The helical instability in C1 exhibits a bi-modal behaviour: these instabilities attain either a single or double helical structure in a random manner, while in C2 and C3 a double helical structure is always excited.The single helical vortex (m=1 mode) is found to be related to the flame detachment by creating a dynamic stagnation point due to the resulting asymmetric flow field, which is absent when the double helical vortex structure (m=2 mode) is excited and the flame remains attached as a consequence.The helical vortex structure modulates the flame surface area and local heat release, amplifying the thermo-acoustic fluctuations, which completes the feedback mechanism of the predicted self-sustained instabilities.

Fig. 2
Fig.2A view of the computational domain (cut by half) and the mesh details

Fig. 5
Fig. 5 Coherent structure visualisation: a 3D low pressure iso-surfaces; b In-plane velocity vectors and the pressure contours overlapped on the velocity magnitude

Fig. 8
Fig. 8 Time-averaged mid-plane images of HRR (left), Y OH (middle), and W fields (right) for C1 (top), C2 (middle) and C3 (bottom).The W fields have been overlaid on top by the contours of the mean binarised Y OH images in red ▸

Fig. 9 Fig. 10
Fig. 9 Experimental mean OH-CL images (top) compared to line-of-sight integration of the LES mean HRR for C1, C2 and C3

Fig. 14 a
Fig. 14 a Mid-plane contour of the predicted time-averaged temperature field in C1 and C3.b Comparison of the temperature profiles

Fig. 15
Fig. 15 Power spectra of the predicted chamber pressure and heat release rate in C1, C2 and C3.Amplitudes are shown in normalised values

Fig. 16
Fig. 16 Temporal signal evolution in a full acoustic cycle for C1, C2 and C3: (a) Air plenum pressure p⅊϶∨̂‶∨̄϶ pl , chamber pressure p ′ ch and integrated HRR.(b) Pressure difference p diff between p ′ pl and p ′ ch and mass flow rates ṁtot and ṁfuel in the same cycle.(c) Equivalence ratio of the unburned gas entering the chamber

Fig. 17
Fig. 17Flame flashback and planes on which mass flow rate and equivalence ratio are monitored Fig. 17Flame flashback and planes on which mass flow rate and equivalence ratio are monitored

Fig. 18
Fig. 18 Periodic evolution of the flow in C1, C2 and C3 at six phases from = 0 • to = 360 • : Left: Mixture fraction; Right: Instantaneous axial velocity.Iso-contours of instantaneous HRR are overlapped on top to show the flame topology

Fig. 20
Fig. 20 Spatially integrated variables at the dump plane in C1. a axial ( w ′ ) and azimuthal velocity fluctuation ( v ′ ).b axial and azimuthal momentum flow rate ( G z and G y ) and the corresponding swirl number (S) oscillations

Fig. 24
Fig. 24 Temporal evolution of vortices structures and their impact on the flame in C2

Fig. 25
Fig. 25Instantaneous interaction of a single PVC with the flame.Left: CH 4 mass fraction overlapped with HRR iso-contours in red with a value of 80MW∕m 3 which is about 5% of the maximum value.Right: inplane velocity vectors coloured by velocity magnitude overlapped on the HRR colour plot.The green isocontours represent the vortical structures