Azimuthally-driven subharmonic thermoacoustic instabilities in a swirl-stabilised combustor

A joint experimental and computational study of thermoacoustic instabilities in a model swirl-stabilised combustor is presented. This paper aims to deliver a better characterisation of such instabilities through the examination of measurements in conjunction with the numerical results obtained via Large Eddy Simulation (LES). The experimental configuration features a cylindrical combustion chamber where the lean premixed methane/air flame experiences self-sustained thermoacoustic oscillations. The nonlinear behaviour of the acoustically excited flame is experimentally investigated by broadband chemiluminescence and dynamic pressure measurements. In LES, the Eulerian stochastic field method is employed for the unknown turbulence/chemistry interaction of the gas-phase. Comparisons of the predicted frequency spectrum and phase-resolved flame structure with measurements are found to be in good agreement, confirming the predictive capabilities of the LES methodology. The presence of Period 2 thermoacoustic oscillations, as a consequence of period-doubling bifurcation, is also confirmed in LES. Through the application of nonlinear analysis both to the experimental and numerical acoustic fluctuations, it is highlighted that the nonlinear behaviour of combustion instabilities in the burner under investigation follows a pattern typical of Period 2 oscillations. Furthermore, the current work demonstrates a useful approach, through the use of dynamic mode decomposition (DMD), for the investigation of unstable flame modes at a specific frequency of interest. The experimental and numerical DMD reconstructions suggest that hot combustion products are convected azimuthally at a rate dictated by the subharmonic frequency. © 2018 The Authors. Published by Elsevier Inc. on behalf of The Combustion Institute. This is an open access article under the CC BY license. ( http://creativecommons.org/licenses/by/4.0/ )


Introduction
Lean premixed combustion is a promising way to reduce NO x emissions in gas turbine engines. Unfortunately, the development of swirl-stabilised combustors operating at this regime is hindered by their susceptibility to the onset of high amplitude thermoacoustic instabilities when the thermal and acoustic fields oscillate in phase. Such instabilities have been studied under various scopes to optimise combustor design parameters. Effects of flow-flame interactions [1] , flame shape transitions [2] and alteration of combustor geometrical parameters [3] on unsteady dynamics have been investigated to increase our understanding of the underlying instability mechanisms. Linear analysis [4] provides insight into the frequency, the growth rate and the unstable oscillation modes. In this linear framework, the feedback coupling between acoustics and the thermal field is described through flame describing functions between unsteady heat release and velocity. However, oscillations may exhibit nonlinear behaviour when more than one frequency is dominant. Nonlinear dynamic techniques can therefore help describe the complex unstable flame behaviour. In [5] the rich nonlinear flame behaviour is investigated by altering the heat release source location along a duct enclosing a laminar flame. The first bifurcation in the system, i.e., a subcritical Hopf bifurcation, is found to cause departure from stable operation into limit cycle oscillations. It is also observed that self-excited oscillations undergo transition into quasi-periodic oscillations via a secondary Hopf bifurcation. The complex behaviour of the system is studied by reconstructing the pressure time-series in phase space. In [6] the authors report a period-doubling bifurcation through the gradual emergence of a hydrodynamic frequency in the acoustic spectra of an unsteady bluff body stabilised combustor. The dynamics depart from a single period state (Period 1 oscillations) while the equivalence ratio is decreased, and gradually lock in Period 2 oscillations prior to lean blowout. In the experimental work [7] where a swirl-stabilised flame is forced, an abrupt period-doubling bifurcation is shown to take place when an excitation threshold is exceeded in association with saturation of the CH * − u transfer function.
Further to experimental approaches, thermoacoustic instabilities can be characterised through high-fidelity simulations such as Large Eddy Simulation (LES). The approach has shown its predictive capabilities for a wide range of turbulent flames and can be directly extended to investigate the flame/acoustics interaction. In the work of Selle et al. [8] , Thickened Flame Model integrated within a compressible LES solver is applied to capture the presence of natural instability modes, such as precessing vortex core (PVC), for non-reactive flows and its suppression due to combustion. The importance of resolving the fuel/air mixing prior to the combustion chamber is highlighted in [9] as a self-excited unstable mode cannot be captured under the assumption of perfect premixing. This argument supports the experimental hypothesis of Meier et al. [10] under the same operating condi-tions. To the best of our knowledge, there have been only few joint experimental and numerical works that attempt to address thermoacoustic instability mechanisms, see for example [11] . On the other hand, various stand-alone computational investigations [12,13] have been performed in the past to study thermoacoustic phenomena experimentally observed. In spite of ongoing efforts, the mechanisms of self-sustained oscillations are not fully understood [14] .
The present work combines experimental and LES investigations into a unique approach, to study the nonlinear behaviour of self-excited instabilities in a model swirl-stabilised configuration. The experimental combustor consists of an axial swirler with a transparent cylindrical combustion chamber. The lean premixed methane/air flame is studied under conditions with φ = 0 . 65 and Re = 17,000, i.e., a condition that demonstrates self-excited thermoacoustic oscillations. The computational investigation is performed using LES wherein the unknown interaction between turbulent motions and chemical reaction is closed by the sub-grid probability density function ( pdf) approach [15] . The operating condition is chosen to jointly investigate the fundamental instability as well as the influence of the subharmonic on the flame behaviour. First, the LES results are compared with measurements in terms of the dynamic spectra from acoustic fluctuations and the flame shape evolution during a period of oscillations. Nonlinear time-series analysis is then used to describe the self-excited flame behaviour which exhibits Period 2 limit cycles with the presence of period-doubling bifurcation. In addition, the current work demonstrates a systematic approach through the application of dynamic mode decomposition (DMD) [16] which allows for the extraction of information on the spatial modes by the frequency of oscillations and the growth/decay rate. To the authors' knowledge, this is the first work that attempts the DMD analysis based on the experimental and numerical snapshots in order to achieve the description of the flame dynamics at the frequency of interest.

Experimental configuration
An overview of the swirler cross section that is used in this study is shown in Fig. 1 . We provide compressed air at 4 bars and monitor the flow rate through a thermal mass flow meter (M+W Instruments, Mass Stream D-6280) and control it through a thermal mass flow controller (M+W Instruments, Mass Stream D-6383). We regulate the flow through the following equipment. We install a critical Venturi nozzle (Cussons Technology Ltd, BS/ISO 9300:2005) to choke the flow and provide a sonic boundary. Upstream, a pressure reducer is installed (Backpressure Regulator, Mackenberg UV 5.1). The flow profile is conditioned before the nozzle through a tube bundle containing a wire mesh gauge and through a flow straightener. We measure the mass flow rate of fuel gas through a thermal mass flow meter (M+W Instruments, Mass Stream D-6250, measurement uncertainty ± 0.625 g/s, 3% factory accuracy at full scale) and control it through a thermal mass flow controller (M+W Instruments, Mass Stream, D-6320, measurement uncertainty ± 0.01 g/s, 1% factory accuracy at full scale). Methane is provided using a fuel injector at the midway of the centrebody, as its cross section is highlighted in Fig. 1 , and is delivered out of ten equally spaced holes into the swirling annular air flow. The major amount of air is provided to the swirler and a smaller percentage to the centreline of the centrebody to ensure adequate cooling. The resulting swirl number is approximately 0.7. The fresh air/fuel mixture is expanded through a diffuser. We ignite the mixture through a spark plug igniter located at the end of the quartz section. The quartz section length is 450 mm with a diameter of 70 mm. Dynamic pressure measurements are accommodated through the installation of a semi-infinite tube (SIT) configuration 60 mm downstream of the far end of the quartz section. The microphone (Kulite MIC-190M, nominal sensitivity: 0.0035 mV/Pa at 160 dBa excitation) is installed 130 mm over the bore of the exhaust and acoustic waves from the combustor are measured through a 1 mm hole drilled on the SIT configuration. Chemiluminescence signals are acquired using a photomultiplier. Data in each experimental case are collected for 10 s at a sampling rate of 32,768 Hz. 5000 high speed broadband luminosity images (bandpass interference filter: 350 ± 40 nm) of the flame are captured at 1000 frames per second. The camera (Photron ICCD) is triggered synchronous to the dynamic pressure signal. Images are acquired at a resolution of 880 × 512 pixels each pixel resolving 0.1367 mm of the imaging area.

LES solver
The numerical investigation is performed using a compressible flow version of the in-house LES solver BOFFIN [17] . The pressure-based LES code solves the governing equations for compressible flows based on the use of structured, boundary conforming, multi-block grids. To include compressibility effects, an equation for the total enthalpy is solved. In the LES equations, there are several unknown terms arising due to the application of a density-weighted filtering for the purpose of scale separations between large-and smallscale turbulent motions. In the filtered momentum equation, the unknown residual stress tensor is represented by a dynamically calibrated version of the Smagorinsky model [18] . Following Gao and O'Brien [15] , the sub-grid pdf approach is employed to represent the unknown interaction between turbulence and chemistry once again caused by the spatial filtering in LES. The solution to the temporal evolution of a joint one-point, one-time pdf of relevant scalars, namely the total specific enthalpy and species mass fractions, is obtained through the Eulerian stochastic field method [19] . The LESpdf formulation has proven its predictive capabilities for turbulent premixed flames, e.g., an industrial gas turbine combustor [20] and a premixed propane flame stabilised on a bluff body [21] . The methane chemistry is characterised by the reduced reaction mechanism of Lu and Law [22] involving 19 chemical species and 15 reaction steps.

Inlet/boundary conditions
The computational domain ( Fig. 2 ) with around 7 million cells and 380 sub-domains has been created using ICEM-CFD to represent all the possible sources for the onset of thermoacoustic oscillations in the combustor under investigation. The minimum and maximum cell sizes are 0.075 and 46 mm, respectively. The grid refinement is made in regions where the premixing of air and fuel, the flame stabilisation and the pressure reflection at the end of the open pipe take place. As suggested in [9] , where the resolution of incomplete mixing between fuel and air was found necessary to numerically capture self-excited unstable flames, the air passage including the axial swirler as well as all the injection holes are discretised in order to explicitly compute the extent of premixing prior to the entry to the combustion chamber with its length of around 1 m. The domain also includes a large cylindrical box such that the reflection of pressure waves at the end of the exhaust pipe open to the atmosphere is computationally captured. A non-reflecting outlet boundary condition at atmospheric pressure can therefore be applied along the outer boundaries of the cylindrical box. This solution approach, i.e., the explicit computation of the pressure reflection at the end of the chamber open to the atmosphere, is an essential feature that triggers the acoustic/flame interaction and is carried out in a manner consistent to the limited examples in literature [9,11] . The inlets (air supply and fuel injection holes) as well as the outlet are described by the Navier-Stokes Characteristic Boundary Conditions (NSCBC) [23,24] . Furthermore, no-slip wall conditions are imposed along the surface of the entire combustor. The velocity components of air and fuel are determined by conserving the measured mass flow rates while fresh gases are injected at 293.15 K. The computational time step is found to vary between 1 and 2 . 5 × 10 −7 s using a constrained CFL number well below 1 based on the speed of the travelling pressure waves. The numerical results presented in the following are those obtained using 30 flame cycles. It should be noted that the quality of all the results is found to be almost insensitive to the number of flame cycles after around 20.

Validation of the numerical results
This section compares the measured and predicted results in terms of the dynamic spectra from pressure fluctuations and the phase-resolved flame structure.

Analysis of natural acoustic modes
The dynamic pressure signals acquired from measurements and LES along with the respective frequency spectra are shown in Fig. 3 . It is worth noting that the experimental time-series represents only a small fraction of the entire measurement duration such that it becomes clear for the purpose of comparison to that predicted by the LES simulation. The frequency resolution of the experimentally and numerically calculated dynamic spectra is 2 . 8 × 10 −6 Hz and 4.85 Hz, respectively. The fundamental frequency, 170 Hz, corresponds to the first quarter wave of a closed to open tube with its length equal to that of the present combustor at the adiabatic flame temperature for the methane/air mixture ( φ = 0 . 65 ). The LES prediction reveals a value (around 170 Hz) close to the measured fundamental frequency. The amplitude of the simulated pressure fluctuations is also in good agreement with that measured experimentally at roughly 155 dBa. In a previous work on the same geometry [25] , it is reported that, while increasing the equivalence ratio above 0.55 for a given Reynolds number, a subharmonic frequency equal to half the fundamental acoustic frequency gradually emerges through a period-doubling bifurcation. The amplitude of the measured subharmonic frequency is found to be around 135 dBa. In the LES computation, the appearance of the subharmonic is clearly evident and its amplitude (140 dBa) is reasonably well captured. However, the predicted subharmonic frequency is shifted towards a lower value of around 65 Hz in comparison to 85 Hz measured in the experiment. Rather than a single, sharp peak as in the experimental spectrum, the computed spectrum features a broader peak due to its poorer frequency resolution; however, the experimental subharmonic frequency of 85 Hz falls within this range. In addition, the amount of the energy contained in the range of frequencies around the subharmonic peaks of LES and measurements ( ± 30 Hz) is investigated. It is interesting that both spectra are found to contain a nearly identical fraction of around 1.6 % of the total energy. Overall, the present dynamic spectra feature an acoustically related fundamental frequency as well as harmonics of the fundamental. Secondary peaks at a subharmonic frequency and at linear combinations of the subharmonic with the fundamental harmonics can also be observed. These observations indicate the establishment of Period 2 oscillations [26] .

Phase-averaged global flame dynamics
To achieve a better representation of the interaction between the acoustic and thermal fields, the measured broadband flame chemiluminescence as well as the simulated heat release rate (line-of-sight integration) are conditionally averaged with respect to the phase angle ψ of the pressure oscillations obtained from measurements and LES, respectively, and averages are then performed for 6 phases (60 °± 30 °, 120 °± 30 °, ..., 360 °± 30 °). The experimental and computational results are presented in Fig. 4 where the phase ψ = 360 • corresponds to the maximum peak of the pressure fluctuations. The experimental flame shape is an ensemble of broadband CO * 2 since in lean methane flames the most significant contribution in chemiluminescent emissions is provided by broadband CO * 2 and by OH * . The relation between broadband flame luminosity and heat release is investigated in [27] . The authors suggest that CO * 2 attains a monotonic relationship with heat release and the choice of imaging radical to be indicative of spatial distribution of heat release is therefore supported. In general, the numerical prediction of the phase-resolved flame structure in response to the acoustic perturbations reveals most of the experimentally observed phenomena. The local extinction of the flame between ψ = 60 • and ψ = 120 • is both experimentally and numerically observed. This extinction occurs as a result of the expected pressure/heat release coupling. From ψ = 180 • until ψ = 240 • , the flame propagates along the combustor walls as the swirling jet entrains the chamber. At ψ = 300 • , the flame lifts off the inlet as captured by LES and the experimental imaging and is further stabilised along the shear layers of the internal recirculation zone (IRZ). At ψ = 360 • , the experimental phaseaveraging suggests that the flame extinguishes from the inlet area and the maximum heat release is locally situated downstream the IRZ. Although the location of the maximum heat release is numerically represented in agreement with measurements, the flame transition to "V" shape is not evident in the computation.

Characterisation of thermoacoustic instabilities
Following the good agreement between LES and experimental data in Section 3 , the nonlinearity in combustion oscillations, namely the emergence of the subharmonic, is characterised through nonlinear analysis. The oscillatory behaviour is described by the reconstruction of attractors in phase space while recurrence plots are employed to quantify the recurrence of phase-space trajectories. The DMD analysis is then performed to reveal the flame dynamics at the subharmonic frequency.

Nonlinearity of thermoacoustic oscillations
The focus here is given on the use of nonlinear analysis techniques to effectively visualise the thermoacoustic oscillations present in the combustor; please refer to [28] for a detailed mathematical description. Following the embedding theorem suggested in [29] , the embedding dimension (hereafter m ) and the optimal time delay ( τ ) are determined for the reconstruction of the dynamics in phase space out of the experimental and numerical dynamic pressure signals. Figures 5  and 6 represent global phase-space portraits highlighting the evolution of the experimental and numerical attractors, respectively, as well as their incremental evolution. From the global attractors in Figs. 5 a and 6 a, it can be observed that the dynamics exhibit a toroidal transition between two manifolds, i.e., from a low amplitude to a higher amplitude limit cycle or vice versa. The presence of Period 2 oscillations is further exemplified by the incremental attractor evolution. For instance, both the experimental and numerical trajectories in Figs. 5 d and 6 c are shown to form a double-looped attractor in phase space instead of a single closed loop in the case of limit cycle oscillations. The double-looped attractor in phase space appears due to the emergence of a subharmonic [30] .
Recurrence plots [31] provide a two-dimensional graphical representation that indicates whether two points on the trajectory are close to each other in phase space. The technique has been used in [32] to effectively define the number of time scales involved in thermoacoustic oscillations. Period 2 and higher oscillations are characterised by the same number    of distinct diagonal line patterns. From the knowledge of the phase-space attractor using the embedding dimension and the optimal time-delay determined earlier for phase-space reconstruction, the recurrence of the phase-space trajectory can be visualised on a two-dimensional map. The threshold for recurrence is defined as a fraction of the maximum Euclidean distance between two points. The fraction is subject to optimisation as a high threshold will obscure the recurrence. The threshold used in this work is 10% of the attractor diameter. Phasespace portraits with m = 3 may downproject the dynamics, whereas recurrence plots reveal further information because they are constructed by the true dimension of the dynamics in phase space. Figure 7 shows the experimental and numerical recurrence plots. The observed white bands generally indicate the excursion of the trajectory away from the regions with a high density of recurrence points in phase space and they correspond to the transitional period between the two manifolds. Typically, the appearance of only equally spaced diagonal lines parallel to the main diagonal in recurrence plots gives an indication for periodic behaviour with a single dominant frequency, i.e., Period 1 oscillation. In comparison, the measured recurrence plot features a pattern with intermittently interrupted diagonal lines (refer to the experimental time scale between 0.15 and 0.2 s in Fig. 7 ) that is indicative of the dynamical motion on a Period 2 attractor [32] . The numerical recurrence plot appears to demonstrate similar patterns to the experimental one including the repeated appearance of white areas and the presence of Period 2 oscillations in the numerical time scale between 0.105 and 0.135 s.

Flame dynamics at the subharmonic frequency
The flame structure that gives rise to the subharmonic frequency is investigated using the DMD analysis on the experimental and computational heat release images as shown in Fig. 8 . The experimentally defined mode at 82 Hz has a near zero negative growth rate equal to −0 . 12 s −1 as one would expect for a periodic process. The reconstruction of the mode depicts the evolution of azimuthally convected hot combustion products around the centreline of the combustor. The burning intensity increases as the hot gas propagates towards the centreline and disappears once it is convected across the opposite side of the combustor. It is worth noting that this phenomenon is not observed in Period 1 oscillations since it is a characteristic dynamic mode that is accompanied with the emergence of Period 2 oscillations. A similar observation has been made in [33] by conducting the DMD analysis on the PIV snapshots of a swirl-stabilised flame that undergoes transitions intermittently between a "V" and "M" shape. A stationary vortex convected azimuthally at a subharmonic frequency is detected therein and its appearance is related to the intermittent transition of the flame towards the outer recirculation zone from the inner shear layer stabilisation zone. The numerical DMD reconstruction exhibits features similar to those from measurements. The emergence of hot combustion products along the right side of the combustor, their azimuthal propagation around the centreline and their eventual disappearance on the opposite side of the chamber are qualitatively captured in the numerically predicted spatial mode at 64 Hz (this frequency corresponds to the numerical subharmonic frequency as shown in Fig. 3 ). To further appreciate the azimuthal propagation of burnt products, another numerical DMD analysis is performed on a transverse plane at an axial location of 100 mm as represented in Fig. 9 . The hot gas initially residing on the right side of the burner propagates around the centre and reaches the other side. It eventually vanishes as the burning intensity decreases. Fig. 9. The DMD reconstruction based on the simulated heat release rate which reveals the azimuthal propagation of hot combustion products around the centre of the combustor. Please note that this plane is normal to the window in Fig. 8 and the quantity is not post-processed using line-of-sight integration.

Conclusions
A lean premixed methane/air flame which exhibits self-sustained combustion instabilities is studied through a joint numerical/experimental investigation in the present work. The characteristic of thermoacoustic instabilities experimentally identified as Period 2 oscillations is compared with that of the numerical pressure fluctuations from the compressible LES computation. The LESpdf formulation is found to capture both the frequency content as well as the amplitude of the acoustic oscillations. Further validation is achieved through the comparison of the phase-averaged flame shape in response to acoustic perturbations. The underlying characteristics of the spatial heat release distribution are found to agree qualitatively well between the experimental and numerical results. It is therefore suggested that the fuel/air premixing chamber as well as the presence of the atmosphere should be taken into account in order to numerically capture the flame/acoustics interaction. This suggestion is in line with other few numerical works [9,11] .
The phase-space reconstruction and recurrence plots are used to characterise the nonlinear behaviour of the thermoacoustically excited flame based on the experimental and numerical results. Both of these representation techniques reveal similar patterns indicative of Period 2 oscillations. Specifically, the numerical phase-space portraits are able to demonstrate the intermittent injection of the dynamics between two attractive manifolds as well as the toroidal transition between the two. The experimental and numerical DMD modes that represent the flame structure at the subharmonic frequency indicate the azimuthal transport of hot combustion products around the centreline of the combustion chamber. The decomposition technique is found capable of extracting the spatial mode at the subharmonic frequency. Future work may include the extension of the DMD analysis to relate the flame structure with the flow field at the subharmonic frequency.