Large eddy simulation and low-order modelling of azimuthal thermoacoustic

instability in a gas turbine model annular combustor

Zhi X. Chen,1, a) Nedunchezhian Swaminathan,2 Marek Mazur,3 Nicholas A. Worth,4 Guangyu Zhang,5 and Lei
Li5
1)
State Key Laboratory of Turbulence and Complex Systems, Aeronautics and Astronautics, College of Engineering,
Peking University, Beijing 100871, P.R. China
2)
Department of Engineering, University of Cambridge, Trumpington Street, Cambridge CB2 1PZ,
UK
3)
CORIA-UMR 6614 Normandie Université, CNRS-Université et INSA de Rouen, Campus Universitaire du Madrillet,
Saint Etienne du Rouvray, France
4)
Department of Energy and Process Engineering, Norwegian University of Science and Technology,
N-7491 Trondheim, Norway
arXiv:2111.00731v1 [physics.flu-dyn] 1 Nov 2021

5)
Research Institute of Aero-Engine, Beihang University, Beijing 100191, P.R. China
Self-excited spinning mode azimuthal instability in an annular combustor with non-swirling flow is investigated
using large eddy simulation (LES). Compressible Navier-Stoke equations are solved with a flamelet combustion
model to describe the subgrid chemistry−turbulence interactions. Two flamelet models, with and without heat
loss effects, are compared to elucidate the non-adiabatic wall effects on the thermoacoustic instability. The
azimuthal modes are captured well by both models with only marginal differences in the computed frequencies
and amplitudes. By comparing with the experimental measurements, the frequencies given by the LES are
approximately 10% higher and the amplitudes are well predicted. Further analysis of the experimental and
LES data shows a similar dominant anti-clockwise spinning mode, under which a good agreement is observed
for the phase-averaged heat release rate fluctuations. Dynamic mode decomposition (DMD) is applied to
shed more light on this spinning mode. The LES and experimental DMD modes reconstructed for their
azimuthal mode frequencies agree very well for the heat release fluctuations. The DMD mode structure for
the acoustic pressure from the LES shows a considerable non-zero profile at the combustor outlet, which
could be essential for azimuthal modes to establish in this annular combustor. Finally, a low-order modelling
study was conducted using an acoustic network combined with the flame transfer function extracted from
LES. The results show that the dominant mode is associated with the plenum showing a first longitudinal and
azimuthal mixed mode structure. By tuning the plenum length to match the effective volume, the predicted
frequency becomes very close to the measured value.

I. INTRODUCTION on supercomputers, accompanied with a vast increase of

computational cost2 . Due to the complex burner ge-
Azimuthal thermoacoustic instabilities remain an un- ometry and boundary conditions, analytical and low-
resolved issue for gas turbine (GT) industries1,2 . This order network models still play essential roles in study-
instability arises when the unsteady combustion of mul- ing azimuthal modes, see Refs.9,10 for extensive reviews
tiple circumferentially aligned flames couples with the on these models. To reduce the computational bur-
acoustic resonant modes of the annular chamber, which den, hybrid approaches combining incompressible LES
often dominates over other modes exhibiting higher os- and acoustic solvers have also been developed11,12 , where
cillation frequencies and amplitudes3 . As pressure waves the connection between solvers, i.e., the flame response
propagate not only in the longitudinal direction but to acoustic perturbations is considered through a flame
also in the transverse direction along the circumfer- transfer or describing function (FTF or FDF). The ad-
ence, azimuthal modes involve more complex dynam- vantage of this approach is that the flame calculation, e.g.
ics such as flame/flame interaction4 , modal dynamics5 , reactive flow LES, is decoupled from the system acoustics
azimuthal/longitudinal mode coupling (e.g., the slanted calculations so that it can still be performed using the
mode)6 , etc. These complexities make azimuthal modes computationally much cheaper low-Mach formulations.
very difficult to understand, predict and control. The major difficulty, however, is how to apply appropri-
To simulate and predict thermoacoustic instabilities, a ate forcing to the flame when the nozzle inlet involves
range of approaches have been developed over the past complex geometry with multiple swirlers, or the pertur-
several decades, from hand-calculated original analyti- bations need to be applied in transverse directions for
cal7,8 models, then desktop running low-order network azimuthal instabilities. With the rapid advances in high-
models, to large eddy simulation (LES) now performed performance parallel computing over recent years, com-
pressible LES for full annular lab-scale burners has now
become computationally feasible13,14 , offering a promis-
ing alternative to investigate azimuthal instabilities and
a) Authorto whom correspondence should be addressed: chen- the complex underlying physical mechanisms for modal
[email protected] (Z.X. Chen). dynamics and structures.
 2

Due to the high cost of running annular test-rigs, few ing the number of burners (i.e., flame-to-flame spacing)
laboratory annular combustor experiments are available and operating flow conditions, various azimuthal modal
(e.g.,5,15 ). It has been shown in these studies that the dynamics were observed involving complex coupling and
swirl-generated azimuthal mean flow can introduce com- interplay among different physical processes, most of
plex time- and operating condition-dependent switching which still remain to be understood. To simplify the sys-
among spinning, standing and mixed modes. To bet- tem and thus better understand the individual processes,
ter understand the modal behaviour of spinning mode this combustor was recently modified by Mazur et. al 16
without azimuthal bulk flow, Mazur et. al 16 removed at NTNU by removing the swirlers in the 12 flames con-
the swirlers from the 12-burner annular rig studied in4,5 figuration (largest flame spacing). As such, azimuthal
and showed well-characterised spinning modes occurring modes could be established in the absence of azimuthal
for high amplitude instabilities over a wide range of op- mean flow and flame−flame interactions, in a geometri-
erating conditions. Therefore, these experiments pro- cally simpler flow.
vide an ideal test case for LES of azimuthal instabilities. The experimental geometry has been described in de-
The main objective of this study is to simulate the non- tail in the previous work5 , and is only summarised briefly
swirling annular premixed combustor of Mazur et. al 16 here together with important changes. A schematic of the
and explore the LES capabilities of capturing spinning modified burner is shown in Fig. 1 with all the dimen-
modes observed in the experiments. Moreover, previous sions in mm. A perfectly premixed ethylene-air mixture
LES studies have been performed over relatively short is passed through a cylindrical plenum chamber, before
run times (tens of ms, much shorter than experiments). being divided around a central bluff body into 12 injector
The present study also investigates the effect of running tubes. Each injector tube holds a central rod, to which
for extended times in order to evaluate the time varying a conical bluff is mounted with a diameter of 13 mm.
mode development. In the present work, a bulk flow velocity of 24.5 m/s
Amongst the available subgrid combustion models, and a equivalence ratio of 0.95 are employed. Further
flamelet approaches are preferable for complex geome- to Ref.16 , the geometry was slightly modified in order
tries because of their robustness and low computational to more closely replicate the simulation boundary condi-
cost. However, compared to models that compute chem- tions. In particular, all flow conditioning elements were
istry explicitly (not necessarily refer to computing chem- removed from the plenum and a high level of uniformity
istry using resolved quantities known as no model or was demonstrated suggesting that the conditioning pro-
quasi-laminar chemistry, but also in thickened flame and vided by the plenum is sufficient. An additional plate was
partially stirred reactor models, where local temperature mounted inside the inner annular wall cylinder (marked
and pressure effects are inherently considered in reac- as ‘blocking plate’ in Fig. 1) in order to prevent the en-
tion rate closure), tabulation models require additional trainment of ambient air through the central opening
non-trivial modelling treatment to account for heat loss during the experiment. Measurements16 were performed
effects, which were shown to have a strong impact on the using dynamic pressure sensors (fs = 51.2 kHz, sam-
FTF/FDF12 . More recent studies by Kraus et al.17 and pled over ts = 10 s) at 3 equi-spaced angular locations,
Agostinelli et al. 18 employed advanced conjugated heat as shown in Fig. 1. Pressure ports were flush-mounted
transfer method to account for the two-way coupling be- in the injector tubes, at a distance 45 mm upstream of
tween the flame LES and burner walls. They showed that the dump plane. Heat release rate imaging is performed
the heat loss effects are important for the swirl burners using a high speed intensified camera (fc = 10 kHz, sam-
considered and need to be taken into account in the LES pled over tc = 2 s, spatial resolution 3.2 pixels/mm).
in order to correctly capture the thermoacoustic instabil-
The computational domain includes the full geome-
ity. Thus, a non-adiabatic flamelet approach is adopted
try of the combustor including the upstream plenum, 12
from Massey et al.19 for this work and the heat loss effects
premixers (feeding tubes), combustion chamber and the
on the instability simulation are investigated.
downstream atmospheric buffering domain. Note that
This paper is organised as follows. The target experi-
the conical volume at the bottom of the plenum was
mental configuration is described in Section II, followed
omitted and only the cylindrical part (marked by green
by the LES modelling framework and numerical setup in
in Fig. 1) is considered to simplify the inflow conditions.
Section III. The LES results are discussed in Section IV
However, this simplification has an considerable impact
and the low-order modelling is presented in Section V.
on the predicted mode frequencies, which will be dis-
The conclusions are summarised in Section VI.
cussed in detailed in Section V. A typical snapshot of
LES undergoing an anti-clockwise (ACW) spinning mode
is also shown in Fig. 1. Chamber pressure contours are
II. EXPERIMENTAL TEST CASE plotted with the arrows indicating the direction of the
travelling azimuthal wave. The oscillating flame surfaces
The target annular combustor was first experimentally are marked using volume rendering for the reaction rate
studied at Cambridge by Worth and Dawson4,5 to char- of progress variable. The axial velocity contours in the
acterise and understand azimuthal modes in laboratory- bluff-body mid-plane are shown for two opposite burn-
scale swirled burners at atmospheric pressure. By vary- ers. From one flame to another along the annulus, the
 3

Fig. 1. Mid-plane cut view of the experimental configuration16 (left) and typical LES snapshot undergoing ACW spinning
mode (indicated by white arrows). Zoomed view of the mesh near bluff-body is also shown.

flame length varies in response to the acoustic pertur- and see details in21 and Refs. therein.
bations resulting in azimuthal heat release fluctuations, Compressible 3D Navier-Stokes equations are solved
which give rise to further acoustic perturbations in the for filtered mass and momentum in a pressure-based form
feedback loop. To capture this phenomenon, an appro- for the low-Mach number flows considered in this study:
priate LES framework considering compressibility effects
is required along with correct numerical boundary condi- ∂ρ  
tions for the acoustics. In addition, a robust and efficient + ∇ · ρU
e = 0, (1)
∂t
SGS combustion model is also essential and these are de-
scribed next.
∂ ρUe  
+ ∇ · ρU
eUe = −∇ p + ∇ · τeff , (2)
III. LES FRAMEWORK AND NUMERICAL SETUP ∂t

where the modified filtered pressure p is computed by

The LES modelling approach used for this study in- solving a Poisson equation. The filtered mixture den-
volves the flamelet chemistry of 1D freely propagat- sity, ρ, is obtained using the ideal-gas equation of
ing laminar premixed flames and a presumed probabil- state. The effective shear stressi tensor is given by
ity density function (PDF) to account for the subgrid
h
τeff = 2ρ ( νe + νt ) Se − (∇ · U
e ) I/3 , where Se is the fil-
scale (SGS) turbulence−chemistry interactions. This ap-
proach has been previously tested and validated exten- tered strain rate, νe and νt are the filtered molecular
sively for a broad range of flames and configurations. In and subgrid eddy viscosities respectively. The latter is
particular, it was shown in20,21 that the unstable flame computed using the Smagorinsky model22 . To describe
dynamics and a longitudinal (Helmholtz) mode thermoa- SGS turbulent premixed combustion with heat losses,
coustic instability in a dual swirl GT burner were well the non-adiabatic flamelet thermochemical states are pre-
captured. Thus, it is of interest here to further assess tabulated using three control variables: filtered reaction
progress variable, e c, its SGS variance, cf002 , and total en-
this approach for azimuthal modes in annular combus-
tors. The modelling framework is briefly described below thalpy, h (sensible+formation), which are solved using
e
 4

their transport equations: speed drops smoothly. Figure 3 compares the measured
and computed mean flame shapes using time-averaged
∂ ρec  
OH∗ chemiluminescence from the experiments and heat
+ ∇ · ρU c = ∇ · ( ρDeff ∇ e
e e c ) + ω̇c , (3)
∂t release rate from LES with non-adiabatic (NAD) and adi-
abatic (AD) flamelet models. Note that since the OH∗
002
image is a line-of-sight visualisation, the LES results are
∂ ρ cf 
002
 
002

integrated in the transverse direction. It is seen that the
+ ∇ · ρU
e cf = ∇ · ρDeff ∇ cf − 2ρχ esgs
∂t flame shape and length are well captured in the LES.
νt
c |2 + 2 c ω̇c − e


+ 2ρ |∇ e c ω̇c , (4) In comparison with the AD model, the NAD case shows
Sct weaker burning outer branches with extinguished flame
roots near the dump plane, which is similar to the ex-
    perimental image. Bauerheim et. al 12 argued that the
∂ ρe
h   νe νt Dp
+ ∇ · ρU h = ∇ · ρ
e e + ∇h +
e , lifted flame base is more prone to acoustics-induced os-
∂t Pr Prt Dt cillations, which may have a impact on the flame response
(5) for azimuthal modes. This is further discussed later in
Section IV B.
where the laminar and turbulent Prandtl numbers, Pr
and Prt , are both taken to be 0.7. The effective diffusiv-
ity Deff is modelled using ( D
e + νt /Sct ) with a turbulent
Schmidt number Sct = 0.423,24 , and D e is the filtered AD-LES NAD-LES Expriment
molecular diffusivity. The progress variable is defined
using the sum of CO and CO2 mass fractions and it is
normalised to vary from 0 to 1 between the fresh and
burnt gases. The unclosed reaction rate and dissipation
rate terms are modelled through a revised flamelet for-
mulation as detailed in Refs.21,25
Laminar 1D premixed flame are calculated using Can-
tera for a range of heat loss levels. The full USC mech- Fig. 3. Mean shape of single flame: line-of-sight integrated
anism (C1−C3)26 is used for ethylene/air combustion LES heat release rate using adiabatic (left) and non-adiabatic
at 300 K and 1 atm. The simulation results are plot- (middle) flamelets, compared with the experimental16 OH∗
ted against the experimental measurements27 in Fig. 2. image (right).
Heat loss effects are implemented in a similar fashion to

The LES modelling framework described above is im-

plemented using the OpenFOAM 2.3.0 libraries. The
computational domain includes the upstream tubes and
plenum, and a fixed mass flow rate condition (u0 = 0)
is employed for the inlet boundary. To mimic the open
downstream boundary, a large hemispherical domain is
added and the far-field boundary is specified as a pressure
outlet (p0 = 0). A fully unstructured grid with about 18
million (∼1.5M per burner) tetrahedral cells is used and
a zoom-in view of the mesh near the bluff-body is shown
in Fig. 1. Mesh sensitivity was investigated for a single
sector using a finer mesh and similar results were found
for the mean fields and also the FTF (see Appendix A).
Isothermal condition of 1000 K12,14 is given for the com-
bustion chamber and bluff-body top walls and adiabatic
condition is specified for all other walls. The Time-step
Fig. 2. Laminar flame speed versus equivalence ratio for ethy- is determined such that CFL < 0.5 for the entire domain
lene/air mixture at 300 K and 1 atm. Symbols: measure- and the typical size is between 2 − 5 × 10−7 s. A fully
ments27 ; Thick line: adiabatic simulation; Thin lines: simu- developed cold flow is first simulated for about 200 ms
lations with heat losses. before ignition and then there is a transient period for az-
imuthal modes to settle and subsequently establish limit-
that used in Ref.19 by scaling down the heat release rate cycle oscillations. The length of this transient period
in the energy equation. As can be seen in the figure, ranges between tens and hundreds of ms, highly depend-
the computed adiabatic flame speeds agree well with the ing on the ignition method used and this is discussed
measurements and as the heat loss increases, the flame next.
 5

1000
longitudinal transition azimuthal
500

-500
FI case P1 P2 P3

-1000
0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34
1000 500

0
500
-500
0.33 0.334 0.338
0
500
-500
SI case 0
-1000
0.2 0.22 0.24 0.26 -500
0.258 0.26 0.262

Fig. 4. Pressure time series for two different ignition cases (for adiabatic simulations). Probes P1, P2 and P3 are marked in
Fig. 1 (upper tapping position).

IV. RESULTS AND DISCUSSION nition to fully excited azimuthal mode is observed in the
LES. Shortly after numerical ignition at t = 0.195 s, sig-
nificant longitudinal oscillations (>500 Pa) appear and
A. Onset of computed azimuthal modes: effect of the
ignition method the three pressure signals collapse onto one another. This
period is quite long for the FI case and the large oscilla-
tions occur several times spanning over almost 0.1 s. For
For LES of multi-sector combustors, ignition is com- the SI case, the longitudinal oscillations reach a much
monly conducted by duplicating the reactive fields of a higher amplitude, about twice as that for the FI case,
separate single-sector simulation to reduce the cost14 . which consequently initiate the transition period at a
This differs from the experiments16 where an ignition much earlier time about 0.015 s after ignition. During
source was used to light the burner from downstream the transition period, mixed longitudinal and azimuthal
and subsequently the flame propagated upstream until oscillations co-exist and the length of this period is quite
stabilised around the bluff-bodies. Since the ignition pro- similar between the two ignition cases, although the am-
cess (known as light-around ) involves an azimuthal mo- plitudes are generally higher for the SI case. Finally,
tion of the flame, which gives rise to azimuthal pressure dominant azimuthal modes are established in both cases
waves, the ignition transients may have an impact on the showing very similar frequency and amplitude as can be
start-up of azimuthal instabilities by acting as an initial seen in the zoomed signals in Fig. 4. This ignition sen-
transverse forcing. sitivity test suggests that although azimuthal instability
To shed light on this, two separate LES cases, referred can be reproduced using either ignition methods, the on-
to as full ignition (FI) and spark ignition (SI) respectively set of the azimuthal modes is strongly influenced by the
hereafter, are performed based on the same initial fields. transient ignition processes and the spark ignition ap-
Figure 4 presents the computed time series of the pressure proach followed by a natural light-around process leads
fluctuation at probes P1, P2 and P3, which are located in to a significantly faster establishment of the instability.
the injector tubes, 120◦ apart in the azimuthal direction This has a significant implication for LES of azimuthal
and 45 mm upstream of the dump plane. For the FI case, instabilities as the computational cost can be substan-
the entire annular chamber is ignited by duplicating the tially reduced.
single-sector solutions. The ignition location for SI is on
the centreline of one of the bluff-bodies and the height is
the same as the inner wall (see marked in Fig. 1). For
both the FI and SI cases, a transient evolution from ig-
 6

B. Comparison with the experiments Figure 5a compares the measured pressure power spec-
tra and those obtained from LES with the AD and NAD
approaches described earlier in Section III. It can be seen
The statistics obtained from experiments and simula- that both models are able to capture the azimuthal insta-
tions for the spinning azimuthal modes are compared in bility although the predicted frequencies are 100−150 Hz
this subsection. It should be noted that the start-up pro- higher than the measured value at around 1800 Hz. The
cess, i.e., from ignition until limit-cycle oscillation, was reason for this difference is that the acoustic mode struc-
observed to be quite long (up to almost 10 s), although ture dictating the excited mode lies in the upstream
the azimuthal instability already appears in the first sev- plenum, having one wavelength in both the longitudinal
eral hundreds of milliseconds. During this process, the and azimuthal directions. Therefore, omitting the conical
oscillation frequency remains nearly the same from the volume (mentioned earlier in Fig. 1) effectively resulted
beginning, whereas the amplitude ramps up gradually. in a shorter longitudinal length and hence a higher funda-
As simulating the entire process for 10 s is computation- mental frequency. More details on this will be discussed
ally impractical, the same sampling period, 400 to 500 ms later in Section V.
after ignition, was chosen for both the experimental and
LES pressure data to make an appropriate comparison. The amplitude in Fig. 5 is also under-predicted by both
LES models by 5−8%. By including the heat loss ef-
fects, the NAD model gives a slightly lower frequency and
(a) 300 smaller amplitude, but the overall difference is marginal
compared to the adiabatic approach. This suggests that
NAD-LES unlike the compact swirling flames12 , the bluff-body sta-
AD-LES bilised flames are much longer and their response to
|pb′ | [Pa/Hz]

200 Exp. acoustic perturbations are less affected flame root dy-
namics induced by the wall heat losses. This is also con-
sistent with previous numerical studies on the FTF/FDF
of a single bluff-body flame28 , see e.g. Refs.29,30 , where
100 the non-linear flame response was well captured by adia-
batic models. As the two LES models give similar results,
further analyses hereafter are performed only for the adi-
abatic case for simplicity.
0 To compare the azimuthal mode characteristics, the
1000 1500 2000 measured and computed pressure data at probes P1, P2
and P3 (see Fig. 1) are processed to obtain the super-
f [Hz] position of azimuthal wave amplitudes following the
least-square fitting procedure used in Refs.5,16 . The am-
1 0.8 0.6 0.4 0.2 0
PDF [1/Pa2] x10-3
plitudes of anti-clockwise and clock-wise spinning waves
are denoted as A+ and A− respectively. Based on these
ACW-spinning amplitudes, the spinning and standing modes can be well
600 defined and the joint PDFs of A+ and A− obtained from
the measurement and LES are compared in Fig. 5b. A
Exp. Standing good agreement can be observed with both exhibiting
an evident ACW-spinning mode. While the experimen-
tal PDF shows a considerable scatter, the LES results
400 fall very close to the vertical axis with high PDF values,
A+ [Pa]

0 0.04 0.08 suggesting a stronger preference in the ACW spinning

1 direction. This may be indicative of less noise in the
LES31,32 . It is also to be noted that in this symmetric
0
SR

200 Exp. burner configuration (without bulk swirl), there should

LES LES be equal probabilities for CW or ACW modes. Some
-1 level of asymmetry due to manufacturing or other un-
certainties resulted in a preferences for ACW observed
CW-spinning over many experimental runs. However, the geometry
0 is perfectly symmetric in the LES, which should prefer
0 200 400 600
A- [Pa]
neither directions, and the ACW mode observed is likely
to be a coincidence for a single LES run. Using the def-
inition given in15 , the spin ratio can be calculated as:
Fig. 5. Comparison of measured and computed (a) pressure SR = (A+ − A− )/(A+ + A− ). The measured and com-
power spectra and (b) joint PDF of the azimuthal wave am- puted SR are plotted in the inset of Fig. 5b. Consistently
plitudes with the spin ratio (SR) plotted in the inset. with the joint PDFs, both the LES and experiment give
 7

-0.05 0 0.05

Exp.

LES

Fig. 6. Comparison of phase-averaged fluctuating heat release rate (FHRR) obtained from overhead OH∗ chemiluminescence
measurements (top) and line-of-sight integration of LES results (bottom).

SR values close to 1 throughout the sampling period, fea- is applied here to further investigate the spinning mode
turing a reasonably stable ACW spinning mode. Mode azimuthal instability.
switching is not observed for this case due to the rel- The DMD algorithm implementation of33 is followed
atively high equivalence ratio and bulk velocity, which to process 516 snapshots collected from continuous mea-
was shown to prefer spinning modes16 . surements and simulations spanning over about 30 ms.
The strength and distribution of the fluctuating heat This corresponds to approximately 50 cycles of the az-
release rate (FHRR) and its response to the spinning az- imuthal mode. Note that the experimental snapshots are
imuthal waves are of prime importance for understanding 2D line-of-sight OH∗ images taken from the top of the
azimuthal instabilities. Thus, phase-averaged normalised burner, and the LES snapshots are 3D heat release rate
FHRRs for six angles (0◦ , 60◦ , ..., 300◦ ) obtained from field for the entire annular combustion chamber. Prior to
the experimental OH∗ images and LES results are com- the DMD computation, the LES results are integrated in
pared in Fig. 6. The maximum heat release fluctuation in the streamwise direction and then projected to the same
the LES is about 7%, which agrees well with the exper- 2D plane as the OH∗ images. The DMD mode growth
imental value of 6%. The previously reported16 crescent rates for the peak frequency (about 1800 and 1900 Hz
shaped structures around both sides of the flame is also for the experiment and LES respectively) are near-zero
observed here in both the measurements and simulations. negative values, suggesting that limit-cycle oscillations
Despite the considerable noise in the experimental data, have been established before the sampling period. It is
a quite symmetric FHRR distribution is seen with half of worth remarking here that in the DMD spectrum for the
the burners forming a positive semicircle and the other LES, two eigen-frequencies at about 900 and 1100 Hz
half forming a negative one. This structure spins as one also showed considerable amplitudes, which are consis-
moves along the cycle from left to right in the figure. tent with the two small peaks seen in Fig. 5 and the low
The above comparisons for the pressure and heat re- frequency content seen on the time series in Fig. 4. These
lease rate show an overall good agreement between the are associated with the longitudinal modes of the com-
experiment and simulation, and this suggests that the bustor geometry. However, both of these DMD modes
spinning mode and its major characteristics are well cap- have large negative growth rates suggesting that they will
tured in the LES. eventually be damped out in a longer simulation. Never-
theless, the azimuthal mode of interest is well established
already in the present simulation duration and hence is
C. DMD analysis of the azimuthal mode discussed next.
The DMD reconstructions of the heat release rate at
Dynamic mode decomposition (DMD)33 is a powerful the azimuthal oscillation frequencies are shown in Fig. 7
technique to study coherent spatial-temporal structures for three phase angles. The reconstructed pressure field
in dynamic systems. Compared to proper orthogonal de- (not available from the experiments) is also shown for
composition (POD), DMD generates a set of modes as- the LES. Note that these phases here do not necessar-
sociated with their own oscillation frequency, amplitude ily correspond to those shown in Fig. 6, because there
and growth/decay rates, which are desired for a better is an arbitrary phase shift during the DMD process. It
interpretation of combustion dynamics34 . Hence, DMD is seen that the reconstructed modes show relatively less
 8

Min Max

Exp.

-1 0 1
LES
Normalised

Fig. 8. DMD mode structure of the acoustic pressure recon-

structed at the azimuthal mode frequency from the LES.
Min Max

Fig. 7. Phase-shifted DMD reconstruction of the FHRR at for details) and the sketch here only shows a half cross-
the azimuthal mode frequency for three phase angles. The
sectional view of the 3D model. We consider the pressure
corresponding combustion chamber pressure mode is shown
for the LES.
disturbance in the form of
p0 = Re[p̂ eiωt ], (6)
noise compared to the simple phase averages, particu- where the operator Re[·] takes real part of the complex
larly for the LES. Again, good agreement is obtained value. The premixers are so thin that the waves inside
here for the DMD modes for FHRR, which further con- can be assumed to be planar waves, written as
firms that the spinning mode behaviour is well captured
−iωx iω(x−L1 )
in the LES. It is also seen that the DMD modes for pres- + c1 +u1N
p̂1N (x) = P1N e −
+ P1N e c1 −u1N , (7)
sure and FHRR seem to be synchronised in the annu-
lar chamber. Finally, the DMD mode structure for the
LES domain is shown in Fig. 8. The mode structure 1 −iωx iω(x−L1 )
+ c1 +u1N −
features a first longitudinal/first azimuthal mixed mode û1N (x) = [P1N e − P1N e c1 −u1N ], (8)
ρ1 c1
in the upstream plenum, a first longitudinal mode with
the pressure node at the middle of the premixer, and where P1N+ −
and P1N denote the complex amplitudes of
a quarter-wave longitudinal/first azimuthal mixed mode the downstream and upstream propagating waves respec-
in the combustion chamber. These are found to be quite tively, and N denotes the sequence number of the pre-
similar to that computed in12 using a Helmholtz solver. mixers. The axial coordinate x is relative to the left
However, evident differences are observed here for the end of each network element (0 for plenum, 1 for pre-
downstream combustion chamber outlet (top of the outer mixer and 2 for combustion chamber). The local speed
wall), where the LES acoustic pressure retains a consid- of sound and mean velocity are denoted by c and u re-
erable value compared to the p0 = 0 condition imposed in spectively, with the subscript corresponding to the given
the Helmholtz solver. This downstream mode structure network element. The waves in the plenum and combus-
may be essential for self-excited azimuthal modes to oc- tion chamber are three-dimensional, and can be obtained
cur. It also implies that downstream acoustic boundary from the wave equation using common notations:
condition could have a major role on the overall stabil-
m=+∞ +∞
ity and azimuthal modes, and hence should be handled X X
appropriately in the analytical and low-order modelling p̂j (r) = ψjm (zjmn r)× (9)
approaches. m=−∞ n=1
+ −
−
+
[Pjmn eikjmn x + Pjmn eikjmn (x−Lj ) ],

V. LOW-ORDER MODELLING +
+∞
1 X
+
kjmn +
Pjmn eikjmn x
ûj (r) = − ψjm (zjmn r)[ + (10)
In order to further understand the azimuthal mode pre- ρj n=1 ω + kjmn uj
dicted in the LES and its dependency on the upstream −
− −
plenum geometry, a low-order modelling analysis is per- kjmn Pjmn eikjmn (x−Lj )
formed. The model combustor geometry is simplified to + − ],
ω + kjmn uj
establish an analytical network acoustic model, as shown
in Fig. 9. Note that the network model considers the sys- where j = 0, 2 and ψ denotes the eigenfunction of the
tem acoustics in the three-dimensional space (see Ref.35 circular duct (j = 0) or annular duct (j = 2), m the
 9

Closed inlet Open outlet

u = 0 p = 0
+
P 2mn
+
P0mn +
P1N
Dp P1N−
R2o −
P2mn
−
R0 P0mn Flame
R2i

0 1 2
L0 L1 Lc

Fig. 9. Sketch of the half cross-section view of the simplified combustor. The sketch is not to scale.

azimuthal mode number and n the radial mode number.

± TABLE I. Model geometrical parameters.
The wavenumbers kjmn can be obtained by dispersion
equations35 . Parameters Symbols Values
Since the combustion zone is nearly isobaric, the flame Length of plenum L0 0.2 m
can be seen as compact and is located at the dump Radius of plenum R0 0.106 m
plane10 . The convective terms involving vorticity and Length of premixer L1 0.155 m
entropy are neglected. Hence, the energy conservation Cross-section area of premixer Ap 0.00015 m2
across the flame is reduced to Length of combustion chamber Lc 0.3 m
γ ∂ Outer radius of combustion chamber R2o 0.106 m
(pûj ) = q̂. (11) Inner radius of combustion chamber R2i 0.0635 m
γ − 1 ∂xj

The unsteady heat release rate term q̂ using the linear

n − τ flame transfer function (FTF) is obtained from the TABLE II. Frequency variation against the plenum length.
LES (full-annular adiabatic case), with the form L0 /m 0.12 0.14 0.16 0.18 0.20 0.22 0.24
f /Hz 3009 2637 2386 2115 1900 1817 1712
q̂N û1N −iωτ
= kf e , x = L0 + L1 , (12) Growthrate/(rad/s) × × × × 409.2 × ×
qN u1N

where the gain kf = 1.1 and the phase difference ωτ =

1.77. The reference point for the velocity fluctuations The list of model parameters used in the above low-
was taken to be 60 mm upstream of the dump plane order calculation are given in Table I. The length of the
to keep clear of the inlet turbulence disturbances. It is plenum is varied while keeping all other parameter un-
worth remarking that the choice of this location has a changed to investigate its effect on the instability mode
quite considerable effect on the above FTF values due characteristics. Note that the length used for the LES
the standing mode inside the premixers. Also note that mesh is 0.2 m (see Fig. 1). The modelling result for the
the FTF is only reliable at the self-excited frequency of frequency and growth rate of the corresponding mode
the LES results. Equation (11) is integrated over the versus the plenum length is presented in Table II. The
(r, θ) plane36 . Then we have the closed equations that frequency of 1900 Hz for the plenum length equals 0.2
characterize the disturbances in the model combustor m agrees very well with LES result of 1910 Hz. Since
the FTF is only reliable at the excited mode frequency,
X(ω)P = 0. (13) only the growth rate at 1900 Hz is presented. The pre-
dicted growth rate is positive implying that the system is
By solving det[X] = 0, we have the complex eigenfre- unstable, which is again consistent with the LES. It can
quency ω = ωr + iωi of the system, where ωr represents be also seen that the predicted frequency monotonically
the natural frequency, and −ωi represents the growth decreases with the increasing plenum length, showing a
rate. The positive growth rate means the system is un- strong longitudinal dependency. This provides an expla-
stable, while the negative growth rate means stable. Re- nation for the over-predicted instability frequency in the
fer to Ref.35 for elaborate modelling details. LES because an effectively shorter length is considered
 10

without including the bottom conical inlet volume. The θ p'

height of this cone is about 60 mm, which would give an 2π 1
effective length of 20 mm if converted to a cylindrical vol-
ume with the same circular area. Consequently, the cor-
responding plenum length L0 = 0.22 m gives a frequency
of 1817 Hz, which is close to the measured value around
1800 Hz. This observation suggests that one should be
0
extremely careful while simplifying the burner geometry
for either LES or low-order modelling as small details
may have significant impact on the system instability be-
haviour.
The predicted mode shape by the low-order model is
shown in Fig. 10. The overall distribution also matches -1
well with that obtained by the LES (see Fig. 8). Figure 11 0 L0 L1 L2 x
depicts the axial pressure mode shape at the mean radius
of the model combustor when L0 = 0.2 m and θ = 0◦ .
It is seen that one wavelength exists in both the plenum Fig. 10. Pressure mode shape at the mean radius of the model
and premixer, while the pressure amplitude in the pre- combustor when L0 = 0.2 m.
mixer is larger than that in the plenum. Although the
combustion chamber occupies the largest length scale, 5
the axial pressure variation there is the smallest due to
the higher temperature and hence larger speed of sound
and wavelength. A notable difference between the mode
shapes predicted by LES and low-order model is that the

'
LES gives a higher amplitude in the plenum compared 0
to that in the combustion chamber, whereas it is the op-
posite in the low-order modelling. It is speculated that
this is caused by the simplified outlet boundary geom-
etry condition, i.e. a flat open outlet shown in Fig. 9,
whereas in the present LES and previous Helmholtz sim- -5 Plenum Premixer Combustion chamber
ulation12 the outer wall is much higher than the inner 0 1
wall (see Fig. 8). The former is fully reflective while the
latter is not. As a result, the amplitude in the combustion
chamber given by the low-order model is relatively higher
Fig. 11. Axial pressure mode shape at the mean radius of the
than that obtained from LES and Helmholtz solver. Un-
model combustor when L0 = 0.2 m and θ = 0◦ .
fortunately, experimental measurement of the pressure
amplitudes in the plenum and combustion chamber is
unavailable for validation. Nevertheless, the LES and
proach to capture self-excited azimuthal modes is demon-
low-order modelling results presented here provide use-
strated for the first time in this work. Although the
ful information for the pressure mode structures inside
amplitudes are predicted well, the computed frequencies
the combustor, which complements the experimental ob-
are higher than the measured value (about 1800 Hz) by
servations with useful physical insights. Future works
8−12%, which is due to the simplified upstream plenum
will be directed towards understanding and predicting
geometry. The subsequent low-order modelling results
the various modal dynamics and explore the effects of
show that when the effectively plenum length is matched
the downstream boundary conditions.
the frequency agrees very well with the measurement.
Analysis of the pressure data from experiments and
LES indicates an anti-clockwise spinning mode with the
VI. CONCLUSIONS spin ratio close to positive unity. Comparison of phase-
averaged heat release rate fluctuations shows a good
A comprehensive numerical study including LES and agreement for this spinning mode. DMD is applied to fur-
low-order modelling of the azimuthal mode thermoacous- ther investigate the modal behaviours at the azimuthal
tic instability in an annular combustor with non-swirling mode frequency. A non-zero combustor outlet mode
flows is presented. Two LES flamelet models, with and structure is observed for the pressure in the LES, which
without considering the heat losses, are used to inves- suggests that the downstream geometry may play an im-
tigate non-adiabatic effects on the azimuthal instability. portant role in establishing the azimuthal modes in this
The azimuthal modes are captured by both models and combustor. The LES observation is further supported
only marginal differences are observed for the frequency by the low-order modelling results showing very similar
and amplitude. The ability of flamelet-based LES ap- mode shapes in the plenum and premixers, but differ-
 11

ent in the combustion chamber due to a simplified out- imuthal mode based on the perimeter of annular cham-
let boundary. The higher amplitude in the combustion ber is considered. Also, as the amplitude of pressure
chamber indicates that a fully reflective open boundary oscillation is below 1% of the ambient pressure12 , the
condition may be inadequate and possible improvements two relatively low forcing amplitudes, 5% and 10% of the
will be explored in future works. mean flow velocity respectively, are used for the forcing.
Figure A.14 presents the time series (after transition)
of the normalised fluctuating velocity at the bluff-body
APPENDIX A. MESH SENSITIVITY exit and heat release response to the excitation for the
two meshes (see Fig. A.12 for the formulae used). The
To avoid the high computational cost of repeatedly overall fluctuating behaviour of the two signals are quite
running full annular combustor simulations, the mesh similar despite some slight differences in the heat release
sensitivity is investigated using a single sector of the 12 amplitude and the time delay (within 10%). To show
burner configuration16 . A schematic of the simulation this in clearer manner, the gain and phase of the transfer
setup is shown in Fig. A.12. The bulk velocity at the functions are computed as shown in Fig. A.15. Again,
bluff-body exit is Ub = 27 m/s and the equivalence ra- similar values are obtained from the two meshes for both
tio is 0.9 for the premixed ethylene/air mixture at 300 K the gain and phase. From 5% to 10% forcing amplitude,
and 1 atm. Two meshes having 1.2 and 2.1 million tetra- the gain slightly decreases and the phase increases, and
hedral cells are tested here and the typical grid spacing this trend is also captured in both the 1.2M and 2.1M
near the shear layer is about 0.6 and 0.3 mm respectively. cases.
The above grid sensitivity study on the single sector
suggests that the smaller 1.2M mesh is sufficient to cap-
ture the reactive flow fields and also the flame response to
external forcing excitations. Thus, this grid distribution
is applied to all 12 sectors for the annular combustion
LES shown in the main text of this paper.

ACKNOWLEDGEMENTS

ZXC and NS acknowledge the support of Mitsubishi

27 Heavy Industries, Takasago, Japan. This work used
the ARCHER UK National Supercomputing Service
(http://www.archer.ac.uk) with the computational time
provided by the UKCTRF.
Fig. A.12. Schematic of the single sector simulation with up-
stream velocity forcing. 1 T. Lieuwen and V. Yang, eds., Combustion instabilities in gas
turbine engines: Operational experience, fundamental mecha-
nisms, and modeling (AIAA Inc., Reston, VA, 2006).
First, the unforced reacting flow is simulated for the 2 T. Poinsot, “Prediction and control of combustion instabilities in
1.2M and 2.1M meshes. After the flame is fully estab- real engines,” Proc. Combust. Inst. 36, 1–28 (2017).
3 J. O’Connor, V. Acharya, and T. Lieuwen, “Transverse com-
lished, statistics are collected for about 4 flow-through-
time, which is calculated based on the combustion cham- bustion instabilities: Acoustic , fluid mechanic , and flame pro-
cesses,” Prog. Energy Combust. Sci. 49, 1–39 (2015).
ber length of 80 mm and the bulk velocity. Figure A.13 4 J. Dawson and N. Worth, “Flame dynamics and unsteady heat
compares the computed mean axial velocity and temper- release rate of self-excited azimuthal modes in an annular com-
ature radial profiles at axial distances of x/D = 0.5, 1 bustor,” Combust. Flame 161, 2565–2578 (2014).
5 N. Worth and J. Dawson, “Modal dynamics of self-excited az-
and 2. The bluff-body diameter is D = 13 mm. As
imuthal instabilities in an annular combustion chamber,” Com-
can be seen, the two simulations give very similar results
bust. Flame 160, 2476–2489 (2013).
for the mean axial velocity field at all three streamwise 6 J. Bourgouin, D. Durox, J. Moeck, T. Schuller, and S. Candel,
positions except for some considerable differences in the “A new pattern of instability observed in an annular combustor:
bottom of the recirculation zone. This is probably due to The slanted mode,” Proc. Combust. Inst. 35, 3237–3244 (2015).
7 L. Crocco, “Aspects of combustion instability in liquid propellant
the relatively short statistics sampling time for the 2.1M
rocket motors. Part I: Fundamentals. Low frequency instability
case as indicated by the asymmetry. The difference in with monopropellants,” J. Am. Rocket Soc. 21, 163–178 (1951).
the temperature fields is almost negligible suggesting a 8 L. Crocco, “Aspects of combustion instability in liquid propellant
good grid convergence. rocket motors. Part II: Low frequency instability with bipropel-
Further grid sensitivity is assessed by performing LES lants. High frequency instability,” J. Am. Rocket Soc. 22, 7–16
with forced inlet velocity (See Fig. A.12) to compare the (1952).
9 M. Bauerheim, F. Nicoud, and T. Poinsot, “Progress in analyt-
flame transfer functions given by the two meshes. As ical methods to predict and control azimuthal combustion insta-
the azimuthal modes are of interests here, a fixed forc- bility modes in annular chambers,” Phys. Fluids 021303 (2016),
ing frequency of 1800 Hz corresponding to the first az- 10.1063/1.4940039.
 12

Fig. A.13. Mean axial velocity and temperature comparisons for the unforced case.

Fig. A.15. Time series of the fluctuating velocity at the bluff-

body exit (black) and the heat release response (red) to the
excitation for the 1M (solid) and 2M (dashed) meshes.

3413 (2012).
14 N. Zettervall, N. Worth, M. Mazur, J. Dawson, and C. Fureby,
Fig. A.14. Time series of the fluctuating velocity at the bluff- “Large eddy simulation of CH4-air and C2H4-air combustion in
body exit (black) and the heat release response (red) to the a model annular gas turbine combustor,” Proc. Combust. Inst.
37, 5223–5231 (2019).
excitation for the 1M (solid) and 2M (dashed) meshes. 15 J.-F. Bourgouin, D. Durox, J. Moeck, T. Schuller, and
S. Candel, “Self-sustained instabilities in an annular combus-
tor coupled by azimuthal and longitudinal acoustic modes,”
in Proc. ASME Turbo Expo, GT2013-95010 (Texas, USA,
10 A. Dowling and S. Stow, “Acoustic analysis of gas turbine com- 2013) https://asmedigitalcollection.asme.org/GT/proceedings-
bustors,” J. Propul. Power 19, 751–764 (2003). pdf/GT2013/55119/V01BT04A007/4223970/v01bt04a007-
11 X. Han, J. Li, and A. Morgans, “Prediction of combustion insta- gt2013-95010.pdf.
bility limit cycle oscillations by combining flame describing func- 16 M. Mazur, H. Nygard, J. Dawson, and N. Worth, “Characteris-

tion simulations with a thermoacoustic network model,” Com- tics of self-excited spinning azimuthal modes in an annular com-
bust. Flame 162, 3632–3647 (2015). bustor with turbulent premixed bluff-body flames,” Proc. Com-
12 M. Bauerheim, G. Staffelbach, N. Worth, J. Dawson, L. Gic- bust. Inst. 37, 5129–5136 (2019).
quel, and T. Poinsot, “Sensitivity of LES-based harmonic flame 17 C. Kraus, L. Selle, and T. Poinsot, “Coupling heat transfer and

response model for turbulent swirled flames and impact on the large eddy simulation for combustion instability prediction in a
stability of azimuthal modes,” Proc. Combust. Inst. 35, 3355– swirl burner,” Combust. Flame 191, 239–251 (2018).
3363 (2015). 18 P. Agostinelli, D. Laera, I. Boxx, L. Gicquel, and T. Poinsot,
13 P. Wolf, G. Staffelbach, L. Gicquel, J.-D. Müller, and T. Poinsot, “Impact of wall heat transfer in Large Eddy Simulation of flame
“Acoustic and Large Eddy Simulation studies of azimuthal modes dynamics in a swirled combustion chamber,” Combust. Flame
in annular combustion chambers,” Combust. Flame 159, 3398– 234, 111728 (2021).
 13

19 J. Massey, Z. Chen, and N. Swaminathan, “Modelling heat loss 28 R. Balachandran, B. Ayoola, C. Kaminski, A. Dowling, and
effects in the Large Eddy Simulation of a lean swirl-stabilised E. Mastorakos, “Experimental investigation of the nonlinear re-
flame,” Flow, Turbul. Combust. 106, 1355–1378 (2021). sponse of turbulent premixed flames to imposed inlet velocity
20 Z. Chen, N. Swaminathan, M. Stöhr, and W. Meier, “Interaction oscillations,” Combust. Flame 143, 37–55 (2005).
between self-excited oscillations and fuel-air mixing in a dual 29 X. Han and A. Morgans, “Simulation of the flame describing

swirl combustor,” Proc. Combust. Inst. 37, 2325–2333 (2019). function of a turbulent premixed flame using an open-source LES
21 Z. Chen, I. Langella, N. Swaminathan, M. Stöhr, W. Meier, and solver,” Combust. Flame 162, 1778–1792 (2015).
H. Kolla, “Large Eddy Simulation of a dual swirl gas turbine 30 S. Ruan, T. Dunstan, N. Swaminathan, and R. Balachandran,

combustor: Flame/flow structures and stabilisation under ther- “Computation of forced premixed flames dynamics,” Combust.
moacoustically stable and unstable conditions,” Combust. Flame Sci. Technol. 188, 1115–1135 (2016).
203, 279–300 (2019). 31 G. Ghirardo and M. Juniper, “Azimuthal instabilities in annular
22 J. Smagorinsky, “General circulation experiments with the prim- combustors: standing and spinning modes,” Proc. R. Soc. A 469,
itive equations,” Mon. Weather Rev. 91, 99–164 (1963). 20130232 (2013).
23 Z. Chen, S. Ruan, and N. Swaminathan, “Large Eddy Simulation 32 G. Ghirardo and F. Gant, “Background noise pushes az-

of flame edge evolution in a spark-ignited methane-air jet,” Proc. imuthal instabilities away from spinning states,” arXiv preprint
Combust. Inst. 36, 1645–1652 (2017). arXiv:1904.00213 (2019).
24 H. Pitsch and H. Steiner, “Large-eddy simulation of a turbulent 33 P. Schmid, “Dynamic mode decomposition of numerical and ex-

piloted methane/air diffusion flame (sandia flame D),” Phys. Flu- perimental data,” J. Fluid Mech. 656, 5–28 (2010).
ids 12, 2541 (2000). 34 C. Huang, W. Anderson, M. Harvazinski, and V. Sankaran,
25 I. Langella, Z. Chen, N. Swaminathan, and S. Sadasivuni, “Analysis of self-excited combustion instabilities using decom-
“Large-Eddy Simulation of reacting flows in industrial gas tur- position techniques,” AIAA J. 54 (2016).
bine combustor,” J. Propul. Power 34, 1269–1284 (2018). 35 G. Zhang, X. Wang, L. Li, and X. Sun, “Effects of perforated
26 H. Wang and A. Laskin, “A comprehensive kinetic model of ethy- liners on controlling combustion instabilities in annular combus-
lene and acetylene oxidation at high temperatures,” (1998). tors,” AIAA J. 58, 3100–3114 (2020).
27 C. Xu and A. Konnov, “Validation and analysis of detailed ki- 36 S. Akamatsu and A. P. Dowling, “Three dimensional thermoa-

netic models for ethylene combustion,” Energy 43, 19–29 (2012). coustic oscillation in a premix combustor,” in ASME Turbo Expo
2001: Power for Land, Sea, and Air (American Society of Me-
chanical Engineers, 2001) pp. V002T02A001–V002T02A001.