<figure><img src="image_1.png" /></figure>

HAL
open science

# Locking of azimuthal modes by breaking the symmetry in annular combustors

José G Aguilar, J.R. Dawson, Thierry Schuller, D. Durox, K. Prieur, S. Candel

► To cite this version:

José G Aguilar, J.R. Dawson, Thierry Schuller, D. Durox, K. Prieur, et al.. Locking of azimuthal modes by breaking the symmetry in annular combustors. Combustion and Flame, 2021, 234, pp.111639. 10.1016/j.combustflame.2021.111639. hal-04537486

HAL Id: hal-04537486

https://ut3-toulouseinp.hal.science/hal-04537486

Submitted on 8 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

<figure><img src="image_2.png" /></figure>

# Locking of azimuthal modes by breaking the symmetry in annular combustors

<figure><img src="image_3.png" /></figure>

José G. Aguilar<sup>a,*</sup>, J.R. Dawson<sup>a</sup>, T. Schuller<sup>b</sup>, D. Durox<sup>c</sup>, K. Prieur<sup>c</sup>, S. Candel<sup>c</sup>

## ARTICLE INFO

**Article history:**
Received 11 March 2021
Revised 21 July 2021
Accepted 22 July 2021
Available online 16 August 2021

**Keywords:**
Combustion dynamics
Limit cycle oscillations
Symmetry breaking
Flame describing function

## ABSTRACT

In this paper we investigate the effects that breaking the symmetry in an annular configuration has on the features of standing, spinning and slanted limit cycle oscillations of thermoacoustic instabilities coupled by azimuthal modes. Experiments are carried out in the MICCA combustor, an annular system comprising sixteen identical matrix injectors. The symmetry of the system is broken by partially blocking the injectors or by changing their geometry. The resulting pressure fields are analysed using a newly introduced pressure ansatz which serves to characterize the nature of the instability (standing or spinning) and the pressure antinodes. Two distinct scenarios are identified. First, when the distribution of the unsteady rate of heat release is slightly changed, for example by partially blocking a single injector, the nature of the limit cycle remains the same. A spinning/standing mode in the symmetric combustor remains spinning/standing in the combustor in which the symmetry is weakly broken. However, it is also found that the orientation of the pressure antinodal line does change and that it is controlled by the pattern of asymmetries introduced in the system. Second, when the distribution of unsteady rate of heat release is changed to a larger extent, for example by partially blocking several burners and by changing their geometries, the nature of the limit cycle can change. A dynamical system based on the pressure ansatz is used to gain insight into the experimental observations. Using a model flame describing function it is shown that one can retrieve many of the features observed in the experimental data, including the nature of the instabilities and the pressure antinodal position as a function of the injection asymmetries. The dynamical system is also used to assess the strength of inhomogeneity introduced in the system. As in the experiments, it is found that in the model starting from a symmetric combustor and then partially blocking a single burner (by modifying its FDF) one is able to alter the pressure antinode. However if in the model one starts from a system with a strong variation of the unsteady heat release rate in the azimuthal direction, partially blocking one injector does not change perceptively any of the features of the predetermined limit cycle.

© 2021 The Author(s). 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/](http://creativecommons.org/licenses/by/4.0/))

## 1. Introduction

Annular combustors are generally used in high power density systems that are prone to thermoacoustic instabilities [1]. These instabilities are classified as longitudinal, standing, spinning, mixed or slanted modes [2]. In recent years research efforts have been focused on the field of modal dynamics, to analyse the onset and nature of these instabilities. In practice these instabilities need to

be avoided because they induce vibrations, structural fatigue and in extreme cases mechanical damage and failure. An understanding of these dynamical processes is required to devise mitigation and control strategies [3] but this is not an easy endeavour. The main challenges arise from the presence of a turbulent combustion process [4] that is coupled by transverse acoustic waves [5]. Using large eddy simulation (LES) for reacting flows Wolf et al. [6] confirmed that at a given operating condition standing and spinning modes could switch randomly. This process was also observed in the experimental data recorded in laboratory scale combustors [7,8] and has been attributed to a stochastic forcing associated with the high level of background noise in the system [9],

<footnote><sup>a</sup>Institute for Energy and Process Engineering, Norwegian University of Science and Technology (NTNU), Strømningsteknisk Kolbjørn Hejesvei 2 7491 Trondheim, Norway</footnote>

<footnote><sup>b</sup>Institut de Mécanique des Fluides de Toulouse, IMFT, Université de Toulouse, CNRS, Toulouse, France</footnote>

<footnote><sup>c</sup>Laboratoire EM2C, CNRS, CentraleSupélec, Université Paris Saclay, 3 Rue Joliot Curie, Gif-sur-Yvette France</footnote>

<footnote>* Corresponding author.
E-mail address: jose.aguilar@ntnu.no (J.G. Aguilar).</footnote>
<figure><img src="image_3.png" /></figure>

# Locking of azimuthal modes by breaking the symmetry in annular combustors

<figure><img src="image_4.png" /></figure>

José G. Aguilar<sup>a,*</sup>, J.R. Dawson<sup>a</sup>, T. Schuller<sup>b</sup>, D. Durox<sup>c</sup>, K. Prieur<sup>c</sup>, S. Candel<sup>c</sup>

## ARTICLE INFO

**Article history:**
Received 11 March 2021
Revised 21 July 2021
Accepted 22 July 2021
Available online 16 August 2021

**Keywords:**
Combustion dynamics
Limit cycle oscillations
Symmetry breaking
Flame describing function

## ABSTRACT

In this paper we investigate the effects that breaking the symmetry in an annular configuration has on the features of standing, spinning and slanted limit cycle oscillations of thermoacoustic instabilities coupled by azimuthal modes. Experiments are carried out in the MICCA combustor, an annular system comprising sixteen identical matrix injectors. The symmetry of the system is broken by partially blocking the injectors or by changing their geometry. The resulting pressure fields are analysed using a newly introduced pressure ansatz which serves to characterize the nature of the instability (standing or spinning) and the pressure antinodes. Two distinct scenarios are identified. First, when the distribution of the unsteady rate of heat release is slightly changed, for example by partially blocking a single injector, the nature of the limit cycle remains the same. A spinning/standing mode in the symmetric combustor remains spinning/standing in the combustor in which the symmetry is weakly broken. However, it is also found that the orientation of the pressure antinodal line does change and that it is controlled by the pattern of asymmetries introduced in the system. Second, when the distribution of unsteady rate of heat release is changed to a larger extent, for example by partially blocking several burners and by changing their geometries, the nature of the limit cycle can change. A dynamical system based on the pressure ansatz is used to gain insight into the experimental observations. Using a model flame describing function it is shown that one can retrieve many of the features observed in the experimental data, including the nature of the instabilities and the pressure antinodal position as a function of the injection asymmetries. The dynamical system is also used to assess the strength of inhomogeneity introduced in the system. As in the experiments, it is found that in the model starting from a symmetric combustor and then partially blocking a single burner (by modifying its FDF) one is able to alter the pressure antinode. However if in the model one starts from a system with a strong variation of the unsteady heat release rate in the azimuthal direction, partially blocking one injector does not change perceptively any of the features of the predetermined limit cycle.

© 2021 The Author(s). 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/](http://creativecommons.org/licenses/by/4.0/))

## 1. Introduction

Annular combustors are generally used in high power density systems that are prone to thermoacoustic instabilities [1]. These instabilities are classified as longitudinal, standing, spinning, mixed or slanted modes [2]. In recent years research efforts have been focused on the field of modal dynamics, to analyse the onset and nature of these instabilities. In practice these instabilities need to

be avoided because they induce vibrations, structural fatigue and in extreme cases mechanical damage and failure. An understanding of these dynamical processes is required to devise mitigation and control strategies [3] but this is not an easy endeavour. The main challenges arise from the presence of a turbulent combustion process [4] that is coupled by transverse acoustic waves [5]. Using large eddy simulation (LES) for reacting flows Wolf et al. [6] confirmed that at a given operating condition standing and spinning modes could switch randomly. This process was also observed in the experimental data recorded in laboratory scale combustors [7,8] and has been attributed to a stochastic forcing associated with the high level of background noise in the system [9],

<footnote><sup>a</sup>Institute for Energy and Process Engineering, Norwegian University of Science and Technology (NTNU), Strømningsteknisk Kolbjørn Hejesvei 2 7491 Trondheim, Norway</footnote>

<footnote><sup>b</sup>Institut de Mécanique des Fluides de Toulouse, IMFT, Université de Toulouse, CNRS, Toulouse, France</footnote>

<footnote><sup>c</sup>Laboratoire EM2C, CNRS, CentraleSupélec, Université Paris Saclay, 3 Rue Joliot Curie, Gif-sur-Yvette France</footnote>

<footnote>* Corresponding author.
E-mail address: jose.aguilar@ntnu.no (J.G. Aguilar).</footnote>
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 1.** (a) Photograph of the MICCA combustor equipped with matrix injectors. (b) Schematic of the experimental set up displaying the microphone locations around the plenum. (c) Excerpt of the pressure signals recorded for the experimental case P1A-01 with $\phi = 0.88$, $\bar{u} = 1.16$ m/s (Spinning mode).</figcaption></figure>

as experimentally confirmed by Faure-Beaulieu et al. [10]. One aspect that deserves special attention is that of defining the nature of limit cycles coupled by azimuthal modes and the possible effects of symmetry breaking.

To learn about the onset of a limit cycle it is common to look into the linear regime of the instabilities. This type of analysis can predict whether a specific mode is stable or unstable to acoustically driven perturbations together with the structure of the mode [11,12]. In symmetric combustors, a mode is said to be degenerate [13] if at the same frequency there exist two linearly independent mode structures. For degenerate modes, a linear framework cannot provide information about the amplitude or the nature of the limit cycle, i.e., whether a given azimuthal mode will saturate into a standing mode or a spinning mode. However, if the symmetry is broken, and the degeneracy is resolved, a linear framework can predict the nature of the oscillation, but the amplitude remains undetermined. Nonlinear analysis using the describing function framework [14–16] is better suited to that purpose. One of the key factors in this method is the introduction of a saturation model for the unsteady rate of heat release that allows a finite amplitude response of the perturbations. Such models have been proposed based on phenomenological arguments [14,15], physics informed models such as those based on flame kinematics [17] or experimental observations [16,18,19].

It is shown by Schuermans et al. [15] that the saturation model applied to an annular system leads to a selection between spinning and standing modes. When the amplitude saturates giving rise to a limit cycle, they found that only spinning modes are stable. Using a cubic saturation model for the heat release rate as a function of the acoustic pressure amplitude, Noiray et al. [18] describe the annular combustor dynamics in terms of Van der Pol oscillators and conclude, again, that in rotationally symmetric configurations, the system saturates into a spinning mode. It is also indicated that if the unsteady rate of heat release of azimuthal mode of order $n$ has a nonuniform azimuthal distribution, the magnitude of the $2n$ Fourier component, determines whether this mode saturates into a mixed mode or, if the magnitude is above a critical value, into a standing mode at limit cycle. Addressing the same issues, Ghirardo and Juniper [20] extended the previous formulation and considered only a symmetric distribution of heat release rate, but introduced transverse velocity fluctuations and used the strength of their coupling with the heat release rate as another bifurcation parameter to control whether the system saturates to a standing or a spinning mode. The previous predictions are limited to systems featuring a cubic saturation of the heat release rate as a function of pressure. For a generic time invariant flame describing function (FDF) Ghirardo et al. [21] derive a dynamical system and provide general rules to test the stability of standing and spinning modes. These conditions are used by Laera et al. [22] to test two different operating conditions in the MICCA burner that yield a spinning and a standing mode respectively.

Some of the aforementioned literature concerns cases where the rotational symmetry is broken. As briefly summarized by Poinot [3] breaking the symmetry of a combustor constitutes an appealing way to mitigate azimuthal instabilities. In practice this can be done by staging the injection units in a combustor [18]. However, depending on the way the symmetry is broken, this may either mitigate the instabilities (e.g., [23]) or if a pair of degenerate modes splits, the system may become more unstable [13,24]. In the latter case the system will saturate into a limit cycle [18,25].

The present study concerns features of limit cycles under broken symmetry conditions. Experiments are carried out in the MICCA combustor from EM2C [26,27] equipped with matrix injectors which have been demonstrated to excite well defined spinning, standing and slanted modes [2,19,22]. These various modes are well established due to a unique feature of this configuration in which combustion takes place through a set of laminar matrix flames. This permits the controlled introduction of asymmetries in the fluctuating heat release rate through small changes to the matrix burners without dealing with the additional complexities introduced by turbulence. This combustor has been used to study the dynamics of the pressure field when the system is subject to a spinning mode [19]. One important conclusion is that under the tested conditions there is a phase difference between the pressure field in the combustion chamber and the pressure field in the plenum. Using the same configuration, Prieur et al. [28] establish a stability map as a function of the operating conditions and identify a dual mode region where standing and spinning modes can be observed. It is found that the nature of the limit cycle in which the system saturates is determined by the initial operating conditions. These phenomena are investigated numerically by Orchini et al. [29], by considering the nonlinear interactions between the modes.

This article begins with a brief presentation of the experimental set up. The pressure signals and modes are then expressed in terms of the ansatz proposed by Ghirardo and Bothien [30]. Stability maps corresponding to the different test cases are then described to delimit conditions leading to stable limit cycles. This is followed by a discussion of the effects of symmetry breaking on the location of the pressure antinodal line in the stable limit cycles. A dynamical system is then derived that may suitably represent this configuration by adapting a model proposed by Ghirardo and Gant [31]. This is used to mimic the results of the experimental cases in order to draw some conclusions about the observed behaviour.

## 2. Experimental setup

The MICCA combustor, shown in Fig. 1a, consists of an annular plenum equipped with 16 matrix injectors and an annular combustion chamber open to the atmosphere. The matrix injectors consist of a thick plate with 89 holes of diameter of $d_h = 2$ mm, fea-
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 1.** (a) Photograph of the MICCA combustor equipped with matrix injectors. (b) Schematic of the experimental set up displaying the microphone locations around the plenum. (c) Excerpt of the pressure signals recorded for the experimental case P1A-01 with $\phi = 0.88$, $\bar{u} = 1.16$ m/s (Spinning mode).</figcaption></figure>

as experimentally confirmed by Faure-Beaulieu et al. [10]. One aspect that deserves special attention is that of defining the nature of limit cycles coupled by azimuthal modes and the possible effects of symmetry breaking.

To learn about the onset of a limit cycle it is common to look into the linear regime of the instabilities. This type of analysis can predict whether a specific mode is stable or unstable to acoustically driven perturbations together with the structure of the mode [11,12]. In symmetric combustors, a mode is said to be degenerate [13] if at the same frequency there exist two linearly independent mode structures. For degenerate modes, a linear framework cannot provide information about the amplitude or the nature of the limit cycle, i.e., whether a given azimuthal mode will saturate into a standing mode or a spinning mode. However, if the symmetry is broken, and the degeneracy is resolved, a linear framework can predict the nature of the oscillation, but the amplitude remains undetermined. Nonlinear analysis using the describing function framework [14–16] is better suited to that purpose. One of the key factors in this method is the introduction of a saturation model for the unsteady rate of heat release that allows a finite amplitude response of the perturbations. Such models have been proposed based on phenomenological arguments [14,15], physics informed models such as those based on flame kinematics [17] or experimental observations [16,18,19].

It is shown by Schuermans et al. [15] that the saturation model applied to an annular system leads to a selection between spinning and standing modes. When the amplitude saturates giving rise to a limit cycle, they found that only spinning modes are stable. Using a cubic saturation model for the heat release rate as a function of the acoustic pressure amplitude, Noiray et al. [18] describe the annular combustor dynamics in terms of Van der Pol oscillators and conclude, again, that in rotationally symmetric configurations, the system saturates into a spinning mode. It is also indicated that if the unsteady rate of heat release of azimuthal mode of order $n$ has a nonuniform azimuthal distribution, the magnitude of the $2n$ Fourier component, determines whether this mode saturates into a mixed mode or, if the magnitude is above a critical value, into a standing mode at limit cycle. Addressing the same issues, Ghirardo and Juniper [20] extended the previous formulation and considered only a symmetric distribution of heat release rate, but introduced transverse velocity fluctuations and used the strength of their coupling with the heat release rate as another bifurcation parameter to control whether the system saturates to a standing or a spinning mode. The previous predictions are limited to systems featuring a cubic saturation of the heat release rate as a function of pressure. For a generic time invariant flame describing function (FDF) Ghirardo et al. [21] derive a dynamical system and provide general rules to test the stability of standing and spinning modes. These conditions are used by Laera et al. [22] to test two different operating conditions in the MICCA burner that yield a spinning and a standing mode respectively.

Some of the aforementioned literature concerns cases where the rotational symmetry is broken. As briefly summarized by Poinot [3] breaking the symmetry of a combustor constitutes an appealing way to mitigate azimuthal instabilities. In practice this can be done by staging the injection units in a combustor [18]. However, depending on the way the symmetry is broken, this may either mitigate the instabilities (e.g., [23]) or if a pair of degenerate modes splits, the system may become more unstable [13,24]. In the latter case the system will saturate into a limit cycle [18,25].

The present study concerns features of limit cycles under broken symmetry conditions. Experiments are carried out in the MICCA combustor from EM2C [26,27] equipped with matrix injectors which have been demonstrated to excite well defined spinning, standing and slanted modes [2,19,22]. These various modes are well established due to a unique feature of this configuration in which combustion takes place through a set of laminar matrix flames. This permits the controlled introduction of asymmetries in the fluctuating heat release rate through small changes to the matrix burners without dealing with the additional complexities introduced by turbulence. This combustor has been used to study the dynamics of the pressure field when the system is subject to a spinning mode [19]. One important conclusion is that under the tested conditions there is a phase difference between the pressure field in the combustion chamber and the pressure field in the plenum. Using the same configuration, Prieur et al. [28] establish a stability map as a function of the operating conditions and identify a dual mode region where standing and spinning modes can be observed. It is found that the nature of the limit cycle in which the system saturates is determined by the initial operating conditions. These phenomena are investigated numerically by Orchini et al. [29], by considering the nonlinear interactions between the modes.

This article begins with a brief presentation of the experimental set up. The pressure signals and modes are then expressed in terms of the ansatz proposed by Ghirardo and Bothien [30]. Stability maps corresponding to the different test cases are then described to delimit conditions leading to stable limit cycles. This is followed by a discussion of the effects of symmetry breaking on the location of the pressure antinodal line in the stable limit cycles. A dynamical system is then derived that may suitably represent this configuration by adapting a model proposed by Ghirardo and Gant [31]. This is used to mimic the results of the experimental cases in order to draw some conclusions about the observed behaviour.

## 2. Experimental setup

The MICCA combustor, shown in Fig. 1a, consists of an annular plenum equipped with 16 matrix injectors and an annular combustion chamber open to the atmosphere. The matrix injectors consist of a thick plate with 89 holes of diameter of $d_h = 2$ mm, fea-
<figure><img src="image_7.png" /><figcaption>Fig. 2. The top row shows the two injector geometries used in the experimental campaign. The bottom row shows the three discrete local rotations (A, B, C) that are assigned to the partially blocked injector located at $\theta = 0$.</figcaption></figure>

turing a porosity of $\sigma_1 = 0.33$. A detailed description of the geometry is found in Laera et al. [22]. The configuration is operated in fully premixed conditions using propane/air mixtures at equivalence ratios in the interval $\phi = [0.86, 1.25]$ with bulk velocities measured in each hole of the perforated plate in the interval $\bar{u} = [1.2, 2.5]$ m/s. High speed images of the full annulus (see Fig. 1b for camera position) are acquired using the same set up as Prieur et al. [28]. Pressure measurements $p'_{1-8}$ are taken at 8 different locations $\theta_k = \pi/4(k - 1/4)$, with $k = 1, 2, \dots, 8$ in the plenum of the combustor as shown in Fig. 1b. The signals are sampled at a frequency of 32.768 kHz for 4 s. Typical pressure signals as recorded by the microphones are shown in Fig. 1c. In its baseline configuration the 16 matrix injectors are identical. Efforts were made to ensure that the flow rate was well distributed through these units with no impact on the occurrence of the various types of oscillations.

The present study is focused on the effects of symmetry breaking on the limit cycles of combustion instabilities due to changes in the spatial distribution of the fluctuating heat release rate. The symmetry of the combustor is broken by two means:

1. Changing the injector geometry.

2. Partial blockage of one or more injection units.

The injector geometry is changed by using a different mesh with 37 holes of diameter $d_h = 3$ mm and a porosity $\sigma_2 = 0.31$, as shown in the top row of Fig. 2. The porosities of both meshes are approximately the same ($\sigma_1 \approx \sigma_2$). However, the larger hole diameter produces longer conical flames. Semicircular partial blockages are introduced by custom made matrix injectors which could be discretely rotated about the local burner axis. We consider three different local rotations with respect to the burner axis labelled A, B and C, as shown in the bottom row of Fig. 2.

Each configuration may be used at different operating conditions yielding one of the four types of azimuthal instabilities: spinning, standing, mixed or slanted modes. Figure 3 and Table 1 summarize all the operating conditions and configurations considered in this study. The configurations labelled P1A, P1B and P1C correspond to azimuthal rotations of a single partially blocked injector at 3 different local rotations. The configurations P2 and P3 have 2 or 3 partially blocked injectors respectively in different arrangements. Finally, the configurations D1, DP1 and DP2 use an injector (or injectors) with different mesh characteristics.

### 3. Characterization of pressure signals

Following Ghirardo and Bothien [30], we use the quaternion framework to represent the pressure signals corresponding to the azimuthal instability. The main advantages of this formulation are

<figure><img src="image_8.png" /><figcaption>Fig. 3. Different configurations tested. The semicircles are used to denote partially blocked injectors. Blue and red are used for the injector's mesh with 89 holes and $d_h = 2$ mm. Yellow is used to denote the injector's mesh with 37 holes and $d_h = 3$ mm. The reference angle $\theta = 0$ is the same as in Fig. 1b. Although not shown in Fig. 2 the yellow half semicircles denote partially blocked injectors with $d_h = 3$ mm. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

discussed later on in Section 7.2. For the moment, the pressure signal is assumed to take the form:

$$ p'(t, \theta) = A \cos[n(\theta - \theta_0)] \cos(\chi) \cos(\omega t + \varphi) + A \sin[n(\theta - \theta_0)] \sin(\chi) \sin(\omega t + \varphi), \quad (1) $$

where A is the amplitude, $\chi$ is the nature angle, $\theta_0$ is the orientation angle, $\varphi$ is the phase angle, n is the azimuthal wave number and $\omega$ is the thermoacoustic frequency. The orientation angle spans the interval $[-\pi, \pi]$ and determines the position of the pressure antinodal line. The nature angle spans the interval $[-\pi/4, \pi/4]$ and describes whether the mode is spinning ($\chi = \pm \pi/4$), standing ($\chi = 0$) or mixed. In this study modes with $0.1 \pi/4 < |\chi| < 0.9 \pi/4$, are regarded as mixed modes.

This formulation is valid when the system displays a spinning, standing or a mixed mode. The MICCA combustor can also saturate into a slanted mode, which is a coupled mode having a standing component and an axisymmetric component [32] (see Fig. 4 for a description of the modes). In this case, the pressure fluctuations can be recast as $\bar{p}'(t, \theta) = p'(t, \theta) + p'_L(t)$, with $p'(t, \theta)$ as given in Eq. (1) and $p'_L(t)$ as the axisymmetric component. Since the microphones are equispaced around the annulus, their average corresponds to the value of the axisymmetric component. Accordingly, to analyse the slanted modes, this average is removed and the rest of the signal is treated as a standing mode.

The pressure signals in the plenum can therefore be characterized in terms of the amplitude and the three angles. The power spectral density of all the analysed signals (not shown for brevity) has a peak at a frequency that corresponds to the first azimuthal wave number. This peak is significantly stronger than the rest of the harmonics, therefore, from this point onwards we will consider only the case when $n = 1$. For the configurations P1A-01 and P1B-
<figure><img src="image_8.png" /><figcaption>Fig. 2. The top row shows the two injector geometries used in the experimental campaign. The bottom row shows the three discrete local rotations (A, B, C) that are assigned to the partially blocked injector located at $\theta = 0$.</figcaption></figure>

turing a porosity of $\sigma_1 = 0.33$. A detailed description of the geometry is found in Laera et al. [22]. The configuration is operated in fully premixed conditions using propane/air mixtures at equivalence ratios in the interval $\phi = [0.86, 1.25]$ with bulk velocities measured in each hole of the perforated plate in the interval $\bar{u} = [1.2, 2.5]$ m/s. High speed images of the full annulus (see Fig. 1b for camera position) are acquired using the same set up as Prieur et al. [28]. Pressure measurements $p'_{1-8}$ are taken at 8 different locations $\theta_k = \pi/4(k - 1/4)$, with $k = 1, 2, \dots, 8$ in the plenum of the combustor as shown in Fig. 1b. The signals are sampled at a frequency of 32.768 kHz for 4 s. Typical pressure signals as recorded by the microphones are shown in Fig. 1c. In its baseline configuration the 16 matrix injectors are identical. Efforts were made to ensure that the flow rate was well distributed through these units with no impact on the occurrence of the various types of oscillations.

The present study is focused on the effects of symmetry breaking on the limit cycles of combustion instabilities due to changes in the spatial distribution of the fluctuating heat release rate. The symmetry of the combustor is broken by two means:

1. Changing the injector geometry.

2. Partial blockage of one or more injection units.

The injector geometry is changed by using a different mesh with 37 holes of diameter $d_h = 3$ mm and a porosity $\sigma_2 = 0.31$, as shown in the top row of Fig. 2. The porosities of both meshes are approximately the same ($\sigma_1 \approx \sigma_2$). However, the larger hole diameter produces longer conical flames. Semicircular partial blockages are introduced by custom made matrix injectors which could be discretely rotated about the local burner axis. We consider three different local rotations with respect to the burner axis labelled A, B and C, as shown in the bottom row of Fig. 2.

Each configuration may be used at different operating conditions yielding one of the four types of azimuthal instabilities: spinning, standing, mixed or slanted modes. Figure 3 and Table 1 summarize all the operating conditions and configurations considered in this study. The configurations labelled P1A, P1B and P1C correspond to azimuthal rotations of a single partially blocked injector at 3 different local rotations. The configurations P2 and P3 have 2 or 3 partially blocked injectors respectively in different arrangements. Finally, the configurations D1, DP1 and DP2 use an injector (or injectors) with different mesh characteristics.

### 3. Characterization of pressure signals

Following Ghirardo and Bothien [30], we use the quaternion framework to represent the pressure signals corresponding to the azimuthal instability. The main advantages of this formulation are

<figure><img src="image_9.png" /><figcaption>Fig. 3. Different configurations tested. The semicircles are used to denote partially blocked injectors. Blue and red are used for the injector's mesh with 89 holes and $d_h = 2$ mm. Yellow is used to denote the injector's mesh with 37 holes and $d_h = 3$ mm. The reference angle $\theta = 0$ is the same as in Fig. 1b. Although not shown in Fig. 2 the yellow half semicircles denote partially blocked injectors with $d_h = 3$ mm. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

discussed later on in Section 7.2. For the moment, the pressure signal is assumed to take the form:

$$p'(t, \theta) = A \cos[n(\theta - \theta_0)] \cos(\chi) \cos(\omega t + \varphi) + A \sin[n(\theta - \theta_0)] \sin(\chi) \sin(\omega t + \varphi), \quad (1)$$

where A is the amplitude, $\chi$ is the nature angle, $\theta_0$ is the orientation angle, $\varphi$ is the phase angle, n is the azimuthal wave number and $\omega$ is the thermoacoustic frequency. The orientation angle spans the interval $[-\pi, \pi]$ and determines the position of the pressure antinodal line. The nature angle spans the interval $[-\pi/4, \pi/4]$ and describes whether the mode is spinning ($\chi = \pm \pi/4$), standing ($\chi = 0$) or mixed. In this study modes with $0.1 \pi/4 < |\chi| < 0.9 \pi/4$, are regarded as mixed modes.

This formulation is valid when the system displays a spinning, standing or a mixed mode. The MICCA combustor can also saturate into a slanted mode, which is a coupled mode having a standing component and an axisymmetric component [32] (see Fig. 4 for a description of the modes). In this case, the pressure fluctuations can be recast as $\tilde{p}'(t, \theta) = p'(t, \theta) + p'_L(t)$, with $p'(t, \theta)$ as given in Eq. (1) and $p'_L(t)$ as the axisymmetric component. Since the microphones are equispaced around the annulus, their average corresponds to the value of the axisymmetric component. Accordingly, to analyse the slanted modes, this average is removed and the rest of the signal is treated as a standing mode.

The pressure signals in the plenum can therefore be characterized in terms of the amplitude and the three angles. The power spectral density of all the analysed signals (not shown for brevity) has a peak at a frequency that corresponds to the first azimuthal wave number. This peak is significantly stronger than the rest of the harmonics, therefore, from this point onwards we will consider only the case when $n = 1$. For the configurations P1A-01 and P1B-
<figure><img src="image_9.png" /><figcaption>Fig. 4. Pressure distribution around the annular combustor at 3 different times for the three types of azimuthal instabilities observed with an azimuthal wave number $n = 1$. The red line is the pressure nodal line and $T$ is the period of oscillation. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article)</figcaption></figure>

03 (see Fig. 3) the probability density functions (PDF) of the orientation angle, nature angle and normalized amplitude at two different operating conditions together with high speed snapshots are shown in Fig. 5. The left column corresponds to a spinning mode and the right column corresponds to a standing mode. One key observation of the data in all the cases explored is that the angles' PDFs feature a very small spread. This indicates that the nature of the instability is deterministic, i.e., the background noise does not play an important role in the observed limit cycles, as for instance in some cases reported in Worth and Dawson [8] for turbulent flames. For most of the cases presented in this paper the PDFs do not change significantly if shorter pressure signal lengths are used.

## 4. Stability map

With the characterization of the pressure signals, each of the test cases listed in Table 1 can be categorized in terms of their limit cycles as spinning, standing, mixed or slanted modes. The top plot in Fig. 6 shows the stability map as a function of the operating conditions i.e., equivalence ratio and bulk velocity. The data in the figure has been overlaid on top of the data reported by Prieur et al. [28], which is taken in the same combustor and similar operating conditions but with 16 identical matrix burners (89 holes, $d_h = 2$ mm). The new stability map matches to a large extent the results of Prieur et al. [28]. The slight differences can be attributed to the partial blockage of the injectors in the present study. At a constant mass flow rate, partially blocking one or several injectors accelerates the flow and changes the flame describing function (FDF). Accordingly, we observe the shift mainly in the direction of $\bar{u}$. The bottom plot in Fig. 6 shows the thermoacoustic frequency. This frequency changes depending on the operating conditions and the type of limit cycle in which the system saturates [29]. This is evidenced by noting that standing modes tend to oscillate at a frequency around 500 Hz, while for spinning modes the oscillation frequency shifts to 480 Hz and for slanted modes around 430 Hz.

This is clear by looking at the dual mode region, where stable spinning and standing modes differ in their thermoacoustic frequencies by 10–20 Hz.

The similarity between maps of a rotationally symmetric combustor and the asymmetric ones indicates that breaking the symmetry by partially blocking one or several of the injectors or by changing the mesh, does not change the stability of a predetermined mode in this annular system. The main reason can be attributed to the fact that the partial blockage of one, two or three injectors may correspond only to a weak symmetry breaking mechanism. A stronger symmetry breaking mechanism might be needed or a combination of several weak mechanisms to fully change the stability of a predetermined limit cycle. In the experimental data there might be just one case which enters this category which corresponds to a standing mode in the spinning region pertaining to configuration DP2-01 at $\bar{u} = 1.67$ m/s, $\phi = 1.02$. In this configuration two injectors were blocked and their injection meshes were changed.

While in the weak symmetry breaking scenario the stability of a specific limit cycle may not change, some of the features such as the orientation angle, i.e., the location of the pressure antinode, can be substantially altered.

## 5. Amplitude and nature angle

The peaks of the PDFs of the nature angle $\chi$ and the amplitude $A$ are shown in Fig. 7. Five of the modes are of mixed nature, since $0 < |\chi| < \pi/4$, but still predominantly spinning. These correspond to the first case of configuration P1A-01, configuration DP1-01, and the first three cases of configuration DP2-01 listed in Table 1. Notice that configuration DP2-01 is the same configuration that gives rise to a standing mode in the spinning region (Fig. 6 top). Furthermore, note that the amplitude of the slanted modes is lower than the standing modes. This is because the average of the signal, which is around 120 Pa, has been removed. Thus, the amplitude shown corresponds only to the standing component. Taking the average into account the slanted and the standing clusters overlap. The last observation is that the spinning modes have on average a larger amplitude than the standing or slanted modes (accounting for the average). However, recall that $A$ does not refer to pressure amplitude around the annulus. Instead, at a given azimuthal angle $\theta$ the pressure amplitude of the acoustic field [31,33] is given by:

$$A_p(\theta) = \frac{A}{\sqrt{2}} \left[ 1 + \cos(2\chi) \cos(2n(\theta_0 - \theta)) \right]^{1/2} \quad (2)$$

For a spinning mode, $|\chi| = \pi/4$, we obtain $A_p(\theta) = A/\sqrt{2}$. For a standing mode, $\chi = 0$, we have $A_p(\theta) = A|\cos(n(\theta_0 - \theta))|$. Hence, for the experimental data it is clear that, on average, the maximum pressure amplitude ($A_p$) is larger for standing modes than for spinning modes as observed in Laera et al. [22].

## 6. Orientation angle

It is now interesting to examine the orientation angle $\theta_0$ (azimuthal position of the pressure antinodal line of the standing component) of the pressure field as measured in the plenum during a limit cycle. Before proceeding we first revisit some of the published literature concerning the symmetric combustor. Bourgouin et al. [19] present the histogram of the pressure nodal lines in the plenum and combustion chamber for a spinning mode. These histograms show that there is a preferred orientation for the standing component of the signal. Laera et al. [22] show that two runs at the same operating condition saturate into a standing limit cycle with two possible pressure antinodal line orientations. One of these orientations appears again in the work of Prieur et al. [28].
<figure><img src="image_10.png" /><figcaption>Fig. 4. Pressure distribution around the annular combustor at 3 different times for the three types of azimuthal instabilities observed with an azimuthal wave number n = 1. The red line is the pressure nodal line and T is the period of oscillation. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article)</figcaption></figure>

03 (see Fig. 3) the probability density functions (PDF) of the orientation angle, nature angle and normalized amplitude at two different operating conditions together with high speed snapshots are shown in Fig. 5. The left column corresponds to a spinning mode and the right column corresponds to a standing mode. One key observation of the data in all the cases explored is that the angles' PDFs feature a very small spread. This indicates that the nature of the instability is deterministic, i.e., the background noise does not play an important role in the observed limit cycles, as for instance in some cases reported in Worth and Dawson [8] for turbulent flames. For most of the cases presented in this paper the PDFs do not change significantly if shorter pressure signal lengths are used.

## 4. Stability map

With the characterization of the pressure signals, each of the test cases listed in Table 1 can be categorized in terms of their limit cycles as spinning, standing, mixed or slanted modes. The top plot in Fig. 6 shows the stability map as a function of the operating conditions i.e., equivalence ratio and bulk velocity. The data in the figure has been overlaid on top of the data reported by Prieur et al. [28], which is taken in the same combustor and similar operating conditions but with 16 identical matrix burners (89 holes, $d_h = 2$ mm). The new stability map matches to a large extent the results of Prieur et al. [28]. The slight differences can be attributed to the partial blockage of the injectors in the present study. At a constant mass flow rate, partially blocking one or several injectors accelerates the flow and changes the flame describing function (FDF). Accordingly, we observe the shift mainly in the direction of $\bar{u}$. The bottom plot in Fig. 6 shows the thermoacoustic frequency. This frequency changes depending on the operating conditions and the type of limit cycle in which the system saturates [29]. This is evidenced by noting that standing modes tend to oscillate at a frequency around 500 Hz, while for spinning modes the oscillation frequency shifts to 480 Hz and for slanted modes around 430 Hz.

This is clear by looking at the dual mode region, where stable spinning and standing modes differ in their thermoacoustic frequencies by 10–20 Hz.

The similarity between maps of a rotationally symmetric combustor and the asymmetric ones indicates that breaking the symmetry by partially blocking one or several of the injectors or by changing the mesh, does not change the stability of a predetermined mode in this annular system. The main reason can be attributed to the fact that the partial blockage of one, two or three injectors may correspond only to a weak symmetry breaking mechanism. A stronger symmetry breaking mechanism might be needed or a combination of several weak mechanisms to fully change the stability of a predetermined limit cycle. In the experimental data there might be just one case which enters this category which corresponds to a standing mode in the spinning region pertaining to configuration DP2-01 at $\bar{u} = 1.67$ m/s, $\phi = 1.02$. In this configuration two injectors were blocked and their injection meshes were changed.

While in the weak symmetry breaking scenario the stability of a specific limit cycle may not change, some of the features such as the orientation angle, i.e., the location of the pressure antinode, can be substantially altered.

## 5. Amplitude and nature angle

The peaks of the PDFs of the nature angle $\chi$ and the amplitude $A$ are shown in Fig. 7. Five of the modes are of mixed nature, since $0 < |\chi| < \pi/4$, but still predominantly spinning. These correspond to the first case of configuration P1A-01, configuration DP1-01, and the first three cases of configuration DP2-01 listed in Table 1. Notice that configuration DP2-01 is the same configuration that gives rise to a standing mode in the spinning region (Fig. 6 top). Furthermore, note that the amplitude of the slanted modes is lower than the standing modes. This is because the average of the signal, which is around 120 Pa, has been removed. Thus, the amplitude shown corresponds only to the standing component. Taking the average into account the slanted and the standing clusters overlap. The last observation is that the spinning modes have on average a larger amplitude than the standing or slanted modes (accounting for the average). However, recall that $A$ does not refer to pressure amplitude around the annulus. Instead, at a given azimuthal angle $\theta$ the pressure amplitude of the acoustic field [31,33] is given by:

$$A_p(\theta) = \frac{A}{\sqrt{2}} \left[ 1 + \cos(2\chi) \cos(2n(\theta_0 - \theta)) \right]^{1/2} \quad (2)$$

For a spinning mode, $|\chi| = \pi/4$, we obtain $A_p(\theta) = A/\sqrt{2}$. For a standing mode, $\chi = 0$, we have $A_p(\theta) = A|\cos(n(\theta_0 - \theta))|$. Hence, for the experimental data it is clear that, on average, the maximum pressure amplitude ($A_p$) is larger for standing modes than for spinning modes as observed in Laera et al. [22].

## 6. Orientation angle

It is now interesting to examine the orientation angle $\theta_0$ (azimuthal position of the pressure antinodal line of the standing component) of the pressure field as measured in the plenum during a limit cycle. Before proceeding we first revisit some of the published literature concerning the symmetric combustor. Bourgouin et al. [19] present the histogram of the pressure nodal lines in the plenum and combustion chamber for a spinning mode. These histograms show that there is a preferred orientation for the standing component of the signal. Laera et al. [22] show that two runs at the same operating condition saturate into a standing limit cycle with two possible pressure antinodal line orientations. One of these orientations appears again in the work of Prieur et al. [28].
<figure><img src="image_10.png" /><figcaption>**Fig. 5.** Probability density functions and high-speed snapshots of a spinning mode (left column) and a standing mode (right column). (a) and (b) show PDFs of the orientation angle, nature angle and normalized amplitude; the vertical gray lines denote the span of the burners around the annulus. The black line denotes the location of a partially blocked injector. (c) and (d) show the corresponding high-speed snapshots; the period is $T \approx 2$ ms, the partially blocked injector is indicated with a blue arrow and the pressure node with a dashed line. The spinning mode in the left column corresponds to configuration P1B-03 at $\bar{u} = 1.57$ m/s, $\phi = 1.04$. The nature angle peaks near $2\chi = \pi/2$ indicating a counter-clockwise spinning mode. The standing mode in the right column corresponds to configuration P1A-01 at $\bar{u} = 2.49$ m/s, $\phi = 0.92$. The nature angle peaks near $2\chi = 0$ and the orientation angle shows that the pressure antinodal line lies between two injectors near $\theta_0 = -\pi/16$. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article)</figcaption></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 6.** Stable limit cycles as a function of equivalence ratio and bulk velocity. The markers denote the current experimental cases, see Table 1 for a reference. The shaded areas correspond to data extracted from Prieur et al. [28]. Under asymmetric conditions of the burners the top and bottom plots show the nature and the (thermoacoustic) frequency at which the modes saturate.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 7.** Amplitude vs nature angle for the different configurations. The color code is the same as for Fig. 6. The amplitude of the slanted mode corresponds only to the standing component of the instability.</figcaption></figure>

Since for the spinning case there is a preferred orientation for the antinodal line and for the standing case the antinodal line orientations repeat in different runs, we can conclude that even in the symmetric case there are already some symmetry breaking effects at play. However, the strength of these asymmetries is not nearly as pronounced as partially blocking an injector, given that this effect seems to override any existing symmetry breaking effects as will be discussed next.

Figure 8 shows the PDF of the orientation angle plotted on top of the different configurations. If for the same configuration more than one operating condition yields the same limit cycle, the PDFs shown are the average. It is natural to begin by discussing the orientation angles of standing modes (Fig. 8a). In most of the cases, when a single injector is partially blocked, the pressure antinode in the plenum follows the angular position of the asymmetry, locking onto an azimuthal location between the blocked injector and an adjacent burner. In a few instances (P1A-04 and P1A-05) the nodal line shifts by the span of one or two injectors. Changing the orientation of the partially blocked injector does not appear to have
<figure><img src="image_11.png" /><figcaption>**Fig. 5.** Probability density functions and high-speed snapshots of a spinning mode (left column) and a standing mode (right column). (a) and (b) show PDFs of the orientation angle, nature angle and normalized amplitude; the vertical gray lines denote the span of the burners around the annulus. The black line denotes the location of a partially blocked injector. (c) and (d) show the corresponding high-speed snapshots; the period is $T \approx 2$ ms, the partially blocked injector is indicated with a blue arrow and the pressure node with a dashed line. The spinning mode in the left column corresponds to configuration P1B-03 at $\bar{u} = 1.57$ m/s, $\phi = 1.04$. The nature angle peaks near $2\chi = \pi/2$ indicating a counter-clockwise spinning mode. The standing mode in the right column corresponds to configuration P1A-01 at $\bar{u} = 2.49$ m/s, $\phi = 0.92$. The nature angle peaks near $2\chi = 0$ and the orientation angle shows that the pressure antinodal line lies between two injectors near $\theta_0 = -\pi/16$. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article)</figcaption></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 6.** Stable limit cycles as a function of equivalence ratio and bulk velocity. The markers denote the current experimental cases, see Table 1 for a reference. The shaded areas correspond to data extracted from Prieur et al. [28]. Under asymmetric conditions of the burners the top and bottom plots show the nature and the (thermoacoustic) frequency at which the modes saturate.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 7.** Amplitude vs nature angle for the different configurations. The color code is the same as for Fig. 6. The amplitude of the slanted mode corresponds only to the standing component of the instability.</figcaption></figure>

Since for the spinning case there is a preferred orientation for the antinodal line and for the standing case the antinodal line orientations repeat in different runs, we can conclude that even in the symmetric case there are already some symmetry breaking effects at play. However, the strength of these asymmetries is not nearly as pronounced as partially blocking an injector, given that this effect seems to override any existing symmetry breaking effects as will be discussed next.

Figure 8 shows the PDF of the orientation angle plotted on top of the different configurations. If for the same configuration more than one operating condition yields the same limit cycle, the PDFs shown are the average. It is natural to begin by discussing the orientation angles of standing modes (Fig. 8a). In most of the cases, when a single injector is partially blocked, the pressure antinode in the plenum follows the angular position of the asymmetry, locking onto an azimuthal location between the blocked injector and an adjacent burner. In a few instances (P1A-04 and P1A-05) the nodal line shifts by the span of one or two injectors. Changing the orientation of the partially blocked injector does not appear to have
<table><caption>**Table 1**

Operating conditions tested in each configuration and the limit cycle corresponding to saturation. The type of limit cycle is followed by an asterisk (*) to denote similar operating conditions near $\bar{u} = 2.16$ m/s, $\phi = 1.13$.</caption><thead><tr><th>Id.</th><th>Configuration</th><th>Limit Cycle</th><th>ū (m/s)</th><th>φ</th></tr></thead><tbody><tr><td rowspan="10">○</td><td rowspan="3">P1A-01</td><td>Mixed</td><td>1.16</td><td>0.88</td></tr><tr><td>Standing</td><td>2.49</td><td>0.92</td></tr><tr><td>Standing*</td><td>2.16</td><td>1.13</td></tr><tr><td>P1A-02</td><td>Slanted</td><td>1.40</td><td>1.18</td></tr><tr><td rowspan="2">P1A-03</td><td>Standing*</td><td>2.15</td><td>1.13</td></tr><tr><td>Slanted</td><td>1.43</td><td>1.17</td></tr><tr><td rowspan="4">P1A-04</td><td>Spinning</td><td>2.09</td><td>1.13</td></tr><tr><td>Standing</td><td>2.24</td><td>1.17</td></tr><tr><td>Standing*</td><td>2.18</td><td>1.19</td></tr><tr><td>Slanted</td><td>1.47</td><td>1.22</td></tr><tr><td rowspan="6">P1A-05</td><td>Standing</td><td>2.44</td><td>1.06</td></tr><tr><td>Standing</td><td>2.30</td><td>1.11</td></tr><tr><td>Standing</td><td>2.46</td><td>1.05</td></tr><tr><td>Slanted</td><td>1.49</td><td>1.22</td></tr><tr><td>Standing</td><td>2.39</td><td>1.14</td></tr><tr><td>Slanted</td><td>1.46</td><td>1.22</td></tr><tr><td rowspan="6">P1B-01</td><td>P1B-02</td><td>Spinning</td><td>1.63</td><td>1.01</td></tr><tr><td rowspan="2">P1B-03</td><td>Spinning</td><td>1.57</td><td>1.04</td></tr><tr><td>Standing</td><td>2.37</td><td>1.09</td></tr><tr><td rowspan="2">P1B-04</td><td>Standing</td><td>1.40</td><td>1.22</td></tr><tr><td>Slanted</td><td>2.07</td><td>1.12</td></tr><tr><td rowspan="2">P1C-01</td><td>Standing</td><td>2.07</td><td>1.15</td></tr><tr><td>Standing*</td><td>2.17</td><td>1.15</td></tr><tr><td rowspan="3">P1C-02</td><td>Standing</td><td>1.43</td><td>1.19</td></tr><tr><td>Standing*</td><td>2.18</td><td>1.14</td></tr><tr><td>Spinning*</td><td>2.10</td><td>1.11</td></tr><tr><td rowspan="4">P2-01</td><td rowspan="4">P2-02</td><td>Standing</td><td>2.32</td><td>1.08</td></tr><tr><td>Standing</td><td>2.31</td><td>1.08</td></tr><tr><td>Standing</td><td>2.30</td><td>1.05</td></tr><tr><td>Slanted</td><td>1.47</td><td>1.17</td></tr><tr><td rowspan="3">P3-01</td><td>Spinning*</td><td>2.10</td><td>1.12</td></tr><tr><td>Standing</td><td>2.24</td><td>1.16</td></tr><tr><td>Slanted</td><td>1.60</td><td>1.14</td></tr><tr><td rowspan="4">D1-01</td><td>Spinning*</td><td>2.15</td><td>1.08</td></tr><tr><td>Standing</td><td>1.98</td><td>1.15</td></tr><tr><td>Slanted</td><td>1.60</td><td>1.16</td></tr><tr><td>Slanted</td><td>1.44</td><td>1.19</td></tr><tr><td rowspan="4">D1-02</td><td>Spinning*</td><td>2.19</td><td>1.13</td></tr><tr><td>Spinning</td><td>1.68</td><td>0.95</td></tr><tr><td>Standing</td><td>2.00</td><td>1.10</td></tr><tr><td>Slanted</td><td>1.44</td><td>1.17</td></tr><tr><td rowspan="5">DP1-01</td><td>Mixed</td><td>1.39</td><td>1.07</td></tr><tr><td>Mixed</td><td>1.34</td><td>0.95</td></tr><tr><td>Mixed</td><td>1.15</td><td>0.93</td></tr><tr><td>Mixed</td><td>1.74</td><td>0.96</td></tr><tr><td>Standing</td><td>1.67</td><td>1.02</td></tr><tr><td rowspan="5">DP2-01</td><td>Mixed</td><td>1.39</td><td>1.07</td></tr><tr><td>Mixed</td><td>1.34</td><td>0.95</td></tr><tr><td>Mixed</td><td>1.15</td><td>0.93</td></tr><tr><td>Mixed</td><td>1.74</td><td>0.96</td></tr><tr><td>Standing</td><td>1.67</td><td>1.02</td></tr></tbody></table>

a strong influence on the locking position as shown by comparing P1A-01, P1B-01, and P1C-01. More asymmetries were introduced in two or more azimuthal locations via one or more partially blocked injectors. The azimuthal locations of the partially blocked injectors were chosen with both odd (P2-01, P2-02, P3-01) and even number of injectors between them (P3-01). Odd numbers were chosen to see if the nodal line could be forced into the centre of a burner as opposed to between burners. As shown in the figure the antinodal line has a well defined locking position that follows the introduced asymmetries. It tends to be located directly adjacent to a partially blocked injector (P3-01, DP2-01) as found in the cases of a single partially blocked injector, or in-between (P2-02) which may be defined by the separation between blocked injectors. The antinodal line position seems to be repeatable, as seen for instance in configuration P1A-04, which has two cases with similar operating conditions and the same antinodal line.

The slanted modes (Fig. 8b) show a preferred orientation that tends towards the centre of the modified burner or its neighbour. It also presents a greater variance, which may be due to its nonlinear longitudinal-standing coupled nature.

The spinning modes have a nature angle $|χ|$ close to $π/4$, but not exactly (see Fig. 5a). Hence, these modes can also be deemed as mixed modes with a predominantly spinning component and a weak but definable standing component. Although the latter component is weak, the antinodal line exhibits a preferred orientation angle, as shown in Fig. 8c. Unlike the standing and slanted modes the preferred orientation angle for the spinning modes is, in most cases, located adjacent to injectors which are approximately 90° apart from the partially blocked injector(s). This analysis reveals a preferred orientation for the pressure antinode, which appears to be a function of the distribution of the asymmetry in the heat release rate around the annulus.

## 7. Dynamical system analysis

A dynamical system is now derived in terms of the amplitude and the three angles used to describe the pressure field to provide a simple model that can be used to understand the features observed in the experimental data and perhaps provide insight into the mechanisms at play in mode selection.

One begins by considering the wave equation in an annular combustor with a linear damping term and unsteady heat release rates corresponding to the flames formed by the $N_b = 16$ matrix injectors located at angles $θ_m = (m - 1)π/8$ with $m = 1, 2, ..., N_b$:

$$ \frac{\partial^2 p'}{\partial t^2} + \alpha \frac{\partial p'}{\partial t} - \frac{\bar{c}^2}{R^2} \frac{\partial^2 p'}{\partial \theta^2} = (\gamma - 1) \sum_{m=1}^{N_b} \frac{\partial q'}{\partial t} 2\pi\delta(\theta - \theta_m), \quad (3) $$

where $α$ is the damping coefficient, $\bar{c}$ is the mean speed of sound in the chamber, $R$ is the mean radius of the annular combustor, $γ$ is the ratio of specific heats and $q'$ is the fluctuations in the heat release rate per unit volume. This form of the wave equation and the pressure ansatz given in Eq. (1) are only valid for the study of the dynamics of spinning and standing modes. For slanted modes we require to include extra terms [32] and for that reason their discussion is left for future studies.

Treating this equation as weakly nonlinear, with the nonlinearity arising from the heat release term [14], and using the method of averaging [34], we derive in Appendix A the following dynamical system for the slow flow variables $A, χ, θ_0, φ$ (similar derivations can be found in [25,31,33]):

$$ \begin{bmatrix} \dot{A} \\ \dot{θ}_0 \\ \dot{χ} \\ \dot{φ} \end{bmatrix} = \begin{bmatrix} -\alpha \frac{A}{2} \\ 0 \\ 0 \\ \frac{\omega_n^2}{2\omega} - \frac{\omega}{2} \end{bmatrix} + \sum_{m=1}^{N_b} \frac{F_r^{(m)}}{2} \begin{bmatrix} A + A \cos(2\chi) \cos(2n(\theta_0 - \theta_m)) \\ -\frac{1}{n} \sec(2\chi) \sin(2n(\theta_0 - \theta_m)) \\ -\sin(2\chi) \cos(2n(\theta_0 - \theta_m)) \\ \tan(2\chi) \sin(2n(\theta_0 - \theta_m)) \end{bmatrix} + \sum_{m=1}^{N_b} \frac{F_j^{(m)}}{2} \begin{bmatrix} 0 \\ -\frac{1}{n} \tan(2\chi) \cos(2n(\theta_0 - \theta_m)) \\ \sin(2n(\theta_0 - \theta_m)) \\ 1 + \sec(2\chi) \cos(2n(\theta_0 - \theta_m)) \end{bmatrix}. \quad (4) $$

Here $ω_n = n\bar{c}/R$ is one of the eigenfrequencies of the purely azimuthal modes of the annular combustor and $ω$ is the thermoacoustic frequency of the system. $F_r^{(m)}$ and $F_j^{(m)}$ are the real and imaginary components of the function:

$$ \mathcal{F}^{(m)} = F_r^{(m)} + jF_j^{(m)} = (\gamma - 1) \frac{\bar{q}}{\bar{p}} \text{FDF}(\theta_m), \quad (5) $$

where FDF($θ_m$) is the flame describing function relative to the pressure fluctuations at location $θ_m$:

$$ \text{FDF}(\theta_m) = \frac{\hat{q}(\theta_m)/\bar{q}}{\hat{p}(\theta_m)/\bar{p}}. \quad (6) $$
<table><caption>**Table 1**

Operating conditions tested in each configuration and the limit cycle corresponding to saturation. The type of limit cycle is followed by an asterisk (*) to denote similar operating conditions near $\bar{u} = 2.16$ m/s, $\phi = 1.13$.</caption><thead><tr><th>Id.</th><th>Configuration</th><th>Limit Cycle</th><th>ū (m/s)</th><th>φ</th></tr></thead><tbody><tr><td rowspan="12">○</td><td rowspan="3">P1A-01</td><td>Mixed</td><td>1.16</td><td>0.88</td></tr><tr><td>Standing</td><td>2.49</td><td>0.92</td></tr><tr><td>Standing*</td><td>2.16</td><td>1.13</td></tr><tr><td>P1A-02</td><td>Slanted</td><td>1.40</td><td>1.18</td></tr><tr><td rowspan="2">P1A-03</td><td>Standing*</td><td>2.15</td><td>1.13</td></tr><tr><td>Slanted</td><td>1.43</td><td>1.17</td></tr><tr><td rowspan="4">P1A-04</td><td>Spinning</td><td>2.09</td><td>1.13</td></tr><tr><td>Standing</td><td>2.24</td><td>1.17</td></tr><tr><td>Standing*</td><td>2.18</td><td>1.19</td></tr><tr><td>Slanted</td><td>1.47</td><td>1.22</td></tr><tr><td rowspan="6">P1A-05</td><td>Standing</td><td>2.44</td><td>1.06</td></tr><tr><td>Standing</td><td>2.30</td><td>1.11</td></tr><tr><td>Standing</td><td>2.46</td><td>1.05</td></tr><tr><td>Slanted</td><td>1.49</td><td>1.22</td></tr><tr><td>Standing</td><td>2.39</td><td>1.14</td></tr><tr><td>Slanted</td><td>1.46</td><td>1.22</td></tr><tr><td rowspan="6">P1B-01</td><td>P1B-02</td><td>Spinning</td><td>1.63</td><td>1.01</td></tr><tr><td rowspan="2">P1B-03</td><td>Spinning</td><td>1.57</td><td>1.04</td></tr><tr><td>Standing</td><td>2.37</td><td>1.09</td></tr><tr><td rowspan="2">P1B-04</td><td>Standing</td><td>1.40</td><td>1.22</td></tr><tr><td>Slanted</td><td>2.07</td><td>1.12</td></tr><tr><td rowspan="2">P1C-01</td><td>Standing</td><td>2.07</td><td>1.15</td></tr><tr><td>Standing*</td><td>2.17</td><td>1.15</td></tr><tr><td rowspan="2">P1C-02</td><td>Standing</td><td>1.43</td><td>1.19</td></tr><tr><td>Standing*</td><td>2.18</td><td>1.14</td></tr><tr><td rowspan="4">P2-01</td><td rowspan="4">P2-02</td><td>Spinning*</td><td>2.10</td><td>1.11</td></tr><tr><td>Standing</td><td>2.32</td><td>1.08</td></tr><tr><td>Standing</td><td>2.31</td><td>1.08</td></tr><tr><td>Standing</td><td>2.30</td><td>1.05</td></tr><tr><td rowspan="4">P3-01</td><td rowspan="4">P3-02</td><td>Slanted</td><td>1.47</td><td>1.17</td></tr><tr><td>Spinning*</td><td>2.10</td><td>1.12</td></tr><tr><td>Standing</td><td>2.24</td><td>1.16</td></tr><tr><td>Slanted</td><td>1.60</td><td>1.14</td></tr><tr><td rowspan="4">D1-01</td><td rowspan="4">D1-02</td><td>Spinning*</td><td>2.15</td><td>1.08</td></tr><tr><td>Standing</td><td>1.98</td><td>1.15</td></tr><tr><td>Slanted</td><td>1.60</td><td>1.16</td></tr><tr><td>Slanted</td><td>1.44</td><td>1.19</td></tr><tr><td rowspan="4">D1-03</td><td rowspan="4">D1-04</td><td>Spinning*</td><td>2.19</td><td>1.13</td></tr><tr><td>Spinning</td><td>1.68</td><td>0.95</td></tr><tr><td>Standing</td><td>2.00</td><td>1.10</td></tr><tr><td>Slanted</td><td>1.44</td><td>1.17</td></tr><tr><td rowspan="4">DP1-01</td><td rowspan="4">DP2-01</td><td>Mixed</td><td>1.39</td><td>1.07</td></tr><tr><td>Mixed</td><td>1.34</td><td>0.95</td></tr><tr><td>Mixed</td><td>1.15</td><td>0.93</td></tr><tr><td>Mixed</td><td>1.74</td><td>0.96</td></tr><tr><td rowspan="2">DP2-02</td><td>Standing</td><td>1.67</td><td>1.02</td></tr></tbody></table>

a strong influence on the locking position as shown by comparing P1A-01, P1B-01, and P1C-01. More asymmetries were introduced in two or more azimuthal locations via one or more partially blocked injectors. The azimuthal locations of the partially blocked injectors were chosen with both odd (P2-01, P2-02, P3-01) and even number of injectors between them (P3-01). Odd numbers were chosen to see if the nodal line could be forced into the centre of a burner as opposed to between burners. As shown in the figure the antinodal line has a well defined locking position that follows the introduced asymmetries. It tends to be located directly adjacent to a partially blocked injector (P3-01, DP2-01) as found in the cases of a single partially blocked injector, or in-between (P2-02) which may be defined by the separation between blocked injectors. The antinodal line position seems to be repeatable, as seen for instance in configuration P1A-04, which has two cases with similar operating conditions and the same antinodal line.

The slanted modes (Fig. 8b) show a preferred orientation that tends towards the centre of the modified burner or its neighbour. It also presents a greater variance, which may be due to its nonlinear longitudinal-standing coupled nature.

The spinning modes have a nature angle $|χ|$ close to $π/4$, but not exactly (see Fig. 5a). Hence, these modes can also be deemed as mixed modes with a predominantly spinning component and a weak but definable standing component. Although the latter component is weak, the antinodal line exhibits a preferred orientation angle, as shown in Fig. 8c. Unlike the standing and slanted modes the preferred orientation angle for the spinning modes is, in most cases, located adjacent to injectors which are approximately 90° apart from the partially blocked injector(s). This analysis reveals a preferred orientation for the pressure antinode, which appears to be a function of the distribution of the asymmetry in the heat release rate around the annulus.

## 7. Dynamical system analysis

A dynamical system is now derived in terms of the amplitude and the three angles used to describe the pressure field to provide a simple model that can be used to understand the features observed in the experimental data and perhaps provide insight into the mechanisms at play in mode selection.

One begins by considering the wave equation in an annular combustor with a linear damping term and unsteady heat release rates corresponding to the flames formed by the $N_b = 16$ matrix injectors located at angles $θ_m = (m - 1)π/8$ with $m = 1, 2, ..., N_b$:

$$ \frac{\partial^2 p'}{\partial t^2} + \alpha \frac{\partial p'}{\partial t} - \frac{\bar{c}^2}{R^2} \frac{\partial^2 p'}{\partial \theta^2} = (\gamma - 1) \sum_{m=1}^{N_b} \frac{\partial q'}{\partial t} 2\pi\delta(\theta - \theta_m), \quad (3) $$

where $α$ is the damping coefficient, $\bar{c}$ is the mean speed of sound in the chamber, $R$ is the mean radius of the annular combustor, $γ$ is the ratio of specific heats and $q'$ is the fluctuations in the heat release rate per unit volume. This form of the wave equation and the pressure ansatz given in Eq. (1) are only valid for the study of the dynamics of spinning and standing modes. For slanted modes we require to include extra terms [32] and for that reason their discussion is left for future studies.

Treating this equation as weakly nonlinear, with the nonlinearity arising from the heat release term [14], and using the method of averaging [34], we derive in Appendix A the following dynamical system for the slow flow variables $A, χ, θ_0, φ$ (similar derivations can be found in [25,31,33]):

$$ \begin{bmatrix} \dot{A} \\ \dot{θ}_0 \\ \dot{χ} \\ \dot{φ} \end{bmatrix} = \begin{bmatrix} -\alpha \frac{A}{2} \\ 0 \\ 0 \\ \frac{\omega_n^2}{2\omega} - \frac{\omega}{2} \end{bmatrix} + \sum_{m=1}^{N_b} \frac{F_r^{(m)}}{2} \begin{bmatrix} A + A \cos(2\chi) \cos(2n(\theta_0 - \theta_m)) \\ -\frac{1}{n} \sec(2\chi) \sin(2n(\theta_0 - \theta_m)) \\ -\sin(2\chi) \cos(2n(\theta_0 - \theta_m)) \\ \tan(2\chi) \sin(2n(\theta_0 - \theta_m)) \end{bmatrix} + \sum_{m=1}^{N_b} \frac{F_j^{(m)}}{2} \begin{bmatrix} 0 \\ -\frac{1}{n} \tan(2\chi) \cos(2n(\theta_0 - \theta_m)) \\ \sin(2n(\theta_0 - \theta_m)) \\ 1 + \sec(2\chi) \cos(2n(\theta_0 - \theta_m)) \end{bmatrix}. \quad (4) $$

Here $ω_n = n\bar{c}/R$ is one of the eigenfrequencies of the purely azimuthal modes of the annular combustor and $ω$ is the thermoacoustic frequency of the system. $F_r^{(m)}$ and $F_j^{(m)}$ are the real and imaginary components of the function:

$$ \mathcal{F}^{(m)} \equiv F_r^{(m)} + jF_j^{(m)} = (\gamma - 1) \frac{\bar{q}}{\bar{p}} \text{FDF}(\theta_m), \quad (5) $$

where FDF($θ_m$) is the flame describing function relative to the pressure fluctuations at location $θ_m$:

$$ \text{FDF}(\theta_m) = \frac{\hat{q}(\theta_m)/\bar{q}}{\hat{p}(\theta_m)/\bar{p}}. \quad (6) $$
<figure><img src="image_14.png" /><figcaption>(a) Standing modes.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(b) Slanted modes. Black and green colormaps indicate the antinode angles where the amplitude is largest and smallest respectively (see Fig. 4).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(c) Spinning modes (black colormap) and mixed modes (green colormap).</figcaption></figure>

<caption>**Fig. 8.** PDF of the pressure antinode (orientation angle $\theta_0$) around the annulus. Darker regions are more likely to occur. The configuration is not shown if no operating condition produced the corresponding mode. If for the same configuration more than one operating condition yields the same limit cycle the PDFs shown are the average.</caption>
<figure><img src="image_15.png" /><figcaption>(a) Standing modes.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(b) Slanted modes. Black and green colormaps indicate the antinode angles where the amplitude is largest and smallest respectively (see Fig. 4).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(c) Spinning modes (black colormap) and mixed modes (green colormap).</figcaption></figure>

<caption>**Fig. 8.** PDF of the pressure antinode (orientation angle $\theta_0$) around the annulus. Darker regions are more likely to occur. The configuration is not shown if no operating condition produced the corresponding mode. If for the same configuration more than one operating condition yields the same limit cycle the PDFs shown are the average.</caption>
<figure><img src="image_17.png" /><figcaption>**Fig. 9.** Gain (left) and phase (right) of the experimental FDF<sub>u</sub> that saturates into a standing mode at a frequency of 500 Hz of the injector with 89 holes, d<sub>n</sub> = 2 mm, and σ<sub>1</sub> = 0.33. The operating conditions are φ = 1.11 and ū = 2.01 m/s. FDF<sub>A</sub> × Z is obtained by fitting the experimental data to Eq. (8). The corresponding parameters are β = 3419, κ = 0.2156 1/Pa, ε = 1437 and μ = 0.01712 1/Pa². FDF<sub>B</sub> × Z is obtained by increasing the value of ε to ε = 2155. The impedance Z is obtained using a low order model of the MICCA combustor.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 10.** Phase portraits for the symmetric and asymmetric combustor in the Poincaré sphere. (a) Poincaré sphere of radius A. The red plane corresponds to the meridian plane (A, 2χ, 0) ∪ (A, 2χ, ±π) or (S₁, S₃). The blue plane corresponds to the equatorial plane (A, 0, θ₀) or (S₁, S₂). The round markers show the location of mixed modes in green, standing modes in red and spinning modes in blue. (b) shows the two planes of the Poincaré sphere for a rotationally symmetric combustor and (c) for the combustor with broken symmetry similar to configuration P1A-01. The top row (red color map) shows the meridian plane. The bottom row (blue color map) shows the equatorial plane. The gray dashed lines show the burner angles. The white markers denote repellers or saddles. The black lines (arrows) show trajectories that remain in the selected planes, i.e., the component perpendicular to the plane is 0. The purple arrow in (c) bottom denotes the burner angle of the partially blocked injector in configuration P1A-01. The maximum value for the radius in all figures is A = 30. The color map shows the projected norm of the stream vectors (LHS in Eq. (9)). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article)</figcaption></figure>

The overbar is used to denote time average and the hat to denote the first Fourier component at the frequency of oscillation ω. Accordingly, the describing function is a function of the slow flow variables, i.e., FDF(θ<sub>m</sub>) = FDF(A, θ<sub>0</sub> - θ<sub>m</sub>, χ, ϕ).

By definition the eigenfrequencies ω<sub>n</sub> depend on the geometry and on the speed of sound. On the other hand the damping rate α, the oscillation frequency ω and the FDF depend on the geometry and operating conditions. The experimental determination of the FDF and that of the damping rate raises issues that may be set aside if the model is not meant to be a predictive tool. This

will be the case in the following development. The modelling effort is only intended to provide further insight on the limit cycle behaviour and on mode selection in a combustor under similar symmetry breaking conditions.

### 7.1. Flame describing function model

To derive a model for the FDF, one may begin by analysing one of the experimental FDFs shown in Laera et al. [22]. These FDFs are defined in terms of the velocity perturbations FDF<sub>u</sub> = (ŷ/q)/(û/ū).
<figure><img src="image_18.png" /><figcaption>**Fig. 9.** Gain (left) and phase (right) of the experimental FDF<sub>u</sub> that saturates into a standing mode at a frequency of 500 Hz of the injector with 89 holes, d<sub>n</sub> = 2 mm, and σ<sub>1</sub> = 0.33. The operating conditions are φ = 1.11 and ū = 2.01 m/s. FDF<sub>A</sub> × Z is obtained by fitting the experimental data to Eq. (8). The corresponding parameters are β = 3419, κ = 0.2156 1/Pa, ε = 1437 and μ = 0.01712 1/Pa<sup>2</sup>. FDF<sub>B</sub> × Z is obtained by increasing the value of ε to ε = 2155. The impedance Z is obtained using a low order model of the MICCA combustor.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Fig. 10.** Phase portraits for the symmetric and asymmetric combustor in the Poincaré sphere. (a) Poincaré sphere of radius A. The red plane corresponds to the meridian plane (A, 2χ, 0) ∪ (A, 2χ, ±π) or (S<sub>1</sub>, S<sub>3</sub>). The blue plane corresponds to the equatorial plane (A, 0, θ<sub>0</sub>) or (S<sub>1</sub>, S<sub>2</sub>). The round markers show the location of mixed modes in green, standing modes in red and spinning modes in blue. (b) shows the two planes of the Poincaré sphere for a rotationally symmetric combustor and (c) for the combustor with broken symmetry similar to configuration P1A-01. The top row (red color map) shows the meridian plane. The bottom row (blue color map) shows the equatorial plane. The gray dashed lines show the burner angles. The white markers denote repellers or saddles. The black lines (arrows) show trajectories that remain in the selected planes, i.e., the component perpendicular to the plane is 0. The purple arrow in (c) bottom denotes the burner angle of the partially blocked injector in configuration P1A-01. The maximum value for the radius in all figures is A = 30. The color map shows the projected norm of the stream vectors (LHS in Eq. (9)). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article)</figcaption></figure>

The overbar is used to denote time average and the hat to denote the first Fourier component at the frequency of oscillation ω. Accordingly, the describing function is a function of the slow flow variables, i.e., FDF(θ<sub>m</sub>) = FDF(A, θ<sub>0</sub> - θ<sub>m</sub>, χ, ϕ).

By definition the eigenfrequencies ω<sub>n</sub> depend on the geometry and on the speed of sound. On the other hand the damping rate α, the oscillation frequency ω and the FDF depend on the geometry and operating conditions. The experimental determination of the FDF and that of the damping rate raises issues that may be set aside if the model is not meant to be a predictive tool. This

will be the case in the following development. The modelling effort is only intended to provide further insight on the limit cycle behaviour and on mode selection in a combustor under similar symmetry breaking conditions.

### 7.1. Flame describing function model

To derive a model for the FDF, one may begin by analysing one of the experimental FDFs shown in Laera et al. [22]. These FDFs are defined in terms of the velocity perturbations FDF<sub>u</sub> = (ŷ/q)/(û/ū).
To deduce the FDF with respect to the pressure fluctuations one has to multiply by an impedance $Z = (\hat{p}/\bar{p})/(\hat{u}/\bar{u})$ which describes the injector response to pressure perturbations and is obtained from a low order model. One may then write

$$FDF = \frac{FDF_u}{Z}. \qquad (7)$$

The experimental $FDF_u$ which saturates into a standing mode at a frequency of 500 Hz, is shown with round black markers in Fig. 9. The corresponding operating conditions are $\phi = 1.11$ and $\bar{u} = 2.01$ m/s which from Fig. 6 it is seen to lie in the dual mode region. It can be observed that as $|\hat{u}|/\bar{u}$ increases, the gain decreases, while the phase remains approximately constant until $|\hat{u}|/\bar{u} \ge 0.4$ when it starts diminishing. Taking this observations into account, the following model is used for the describing function:

$$FDF(\theta_m) = \frac{2\beta_m}{1 + \sqrt{1 + \kappa_m^2 A_p(\theta_m)^2}} - \epsilon_m \exp\left(-\frac{\mu_m A_p(\theta_m)^2}{\sqrt{\pi} \ln(2)}\right). \quad (8)$$

In Appendix B the role of the parameters $\beta_m$, $\kappa_m$, $\epsilon_m$, $\mu_m$ is explained and it is also shown that this function corresponds to a model where the unsteady rate of heat release saturates as the pressure amplitude increases. After multiplying Eq. (8) by the impedance $Z$, one can fit the model to the experimental data. The result is shown in Fig. 9 with a blue line and labelled $FDF_A$. The model accurately captures the gain with increasing amplitude, but completely disregards the development of the phase; the phase shown is only due to the impedance $Z$. One can include a function that accounts for the deviation in the phase. However, to keep the model simple and to facilitate the forthcoming analysis Eq. (8) is kept real.

The corresponding plot of $FDF_A$ as a function of $A_p$ is also shown in Appendix B, Fig. B.14b in blue. In this figure one can observe that the FDF in terms of pressure fluctuations is monotonic. Using the tools developed by Ghirardo et al. [21] one finds that, in a symmetric combustor, real FDFs, such as Eq. (8), that decrease monotonically with amplitude only saturate into spinning modes, whereas real FDFs that are non-monotonic can saturate into both standing and spinning modes. Therefore, to ensure that standing modes are observed in the forthcoming analysis the parameter $\epsilon$ is increased until a non-monotonic behaviour is obtained. The result is shown in Figs. 9 and B.14 b in red and labelled $FDF_B$.

Since the model FDF is real, there are two implications on the dynamical system. First, $F_j^{(m)} = 0$ in Eq. (4). Second, the thermoacoustic frequency and the natural frequency<fnref n="1" /> coincide and one may set $\omega = \omega_n$.

## 7.2. Simplified dynamical system

With the previous model assumptions for the FDF and considering the first azimuthal mode corresponding to $n = 1$, the dynamical system reduces to:

$$\begin{bmatrix} \dot{A} \\ \dot{\theta}_0 \\ \dot{\chi} \\ \dot{\varphi} \end{bmatrix} = \begin{bmatrix} -\alpha \frac{A}{2} \\ 0 \\ 0 \\ 0 \end{bmatrix} + \sum_{m=1}^{N_b} \frac{F_r^{(m)}}{2} \begin{bmatrix} A + A \cos(2\chi) \cos(2(\theta_0 - \theta_m)) \\ -\sec(2\chi) \sin(2(\theta_0 - \theta_m)) \\ -\sin(2\chi) \cos(2(\theta_0 - \theta_m)) \\ \tan(2\chi) \sin(2(\theta_0 - \theta_m)) \end{bmatrix}. \quad (9)$$

Before proceeding to the discussion of the results it is important to stress the main advantage of the quaternion formulation of the pressure fluctuations. As can be seen from Eqs. (4) and (9), the variables used to describe the nature of the azimuthal instabilities are state space variables [30]. In other words, from this representation one can directly obtain with no further derivations, the time evolution of (1) the amplitude, (2) pressure antinodal line angular position, (3) spinning, standing or mixed nature of the mode and (4) the phase of the azimuthal oscillation. This is opposed to using derived quantities such as the C-indicator or the spin ratio which require some treatment of the pressure signal as done, for instance, in Worth and Dawson [8].

## 7.3. Phase portraits of the symmetric combustor

To illustrate the evolution of the dynamical system, the phase portraits are shown in two planes of the Poincaré sphere defined by $(A, 2\chi, \theta_0)$ as shown in Fig. 10a. The first plane, shown in red, is the meridian plane $(A, 2\chi, 0) \cup (A, 2\chi, \pm\pi)$ or alternatively $(S_1, S_3)$ in the Cartesian coordinate system. The attractors in this plane can be:

a) Spinning modes if $|2\chi| = \pi/2$. In the figure the two of them are shown with blue round markers.

b) Mixed modes if they lie in the region $0 < |2\chi| < \pi/2$. One of them (in a different plane) is shown with a green round marker.

c) Standing modes if $2\chi = 0$. Two of them are shown with red round markers.

The second plane, shown in blue, is the equatorial plane $(A, 0, \theta_0)$ or alternatively $(S_1, S_2)$ in the Cartesian coordinate system. All the attractors in this plane correspond to standing modes $(2\chi = 0)$.

Given a damping coefficient $\alpha$, the system can display qualitatively different phase portraits. In Fig. 10b we present the phase portraits at $\alpha \bar{p}/((\gamma - 1)\bar{q}) = 2.1 \times 10^4$. This value is chosen:

1. Relative to the magnitude of the FDF. That is to say that we can always re-scale the FDF and $\alpha$, and the phase portraits will remain the same, as long as their relative magnitude is kept the same.

2. To reproduce as much as possible some of the conditions observed in the experimental campaign.

For a symmetric combustor, the top and bottom rows of Fig. 10b show the meridian and the equatorial planes respectively. The dynamical system has stable standing modes and stable spinning modes. On the meridian plane (top) there are stable attractors at the origin and at $(25.57, 2 \times (\pm\pi/4), 0)$ corresponding to spinning modes. On the equatorial plane (bottom) there are stable attractors at the origin and at $(19.81, 0, (2k - 17)\pi/16)$ with $k = 1, 2, \dots, N_b$ corresponding to standing modes. These positions coincide with angles in between burners. As in the experimental cases (see Fig. 7) the spinning modes have larger amplitude than the standing modes.<fnref n="2" /> Furthermore, at small amplitudes the system is linearly stable. Observe that all the trajectories in this region point towards the origin in both planes. However, at moderate amplitudes the system can saturate into standing or spinning modes.<fnref n="3" /> The initial conditions of the system will determine the type of limit cycle at which the system saturates. This is reminiscent of the dual mode region in Fig. 6 and thus the hysteresis region analysed by Prieur et al. [28].

<footnote n="1">¹ The thermoacoustic frequency and the natural frequency of a system are usually close to each other. The acoustic or natural frequency of the MICCA combustor has been numerically computed to be around 450 Hz [35], while the thermoacoustic frequencies lie between 420 and 510 Hz according to Fig. 6.</footnote>

<footnote n="2">² While the amplitude $A$ has $A^{sp} = 25.57 > A^{st} = 19.81$, the pressure distribution around the annulus $A_p$ has $A_p^{sp} = 25.57/\sqrt{2} < \max(A_p^{st}) = 19.81$.</footnote>

<footnote n="3">³ Under these conditions the system is capable of triggering.</footnote>
To deduce the FDF with respect to the pressure fluctuations one has to multiply by an impedance $Z = (\hat{p}/\bar{p})/(\hat{u}/\bar{u})$ which describes the injector response to pressure perturbations and is obtained from a low order model. One may then write

$$FDF = \frac{FDF_u}{Z}. \qquad (7)$$

The experimental $FDF_u$ which saturates into a standing mode at a frequency of 500 Hz, is shown with round black markers in Fig. 9. The corresponding operating conditions are $\phi = 1.11$ and $\bar{u} = 2.01$ m/s which from Fig. 6 it is seen to lie in the dual mode region. It can be observed that as $|\hat{u}|/\bar{u}$ increases, the gain decreases, while the phase remains approximately constant until $|\hat{u}|/\bar{u} \ge 0.4$ when it starts diminishing. Taking this observations into account, the following model is used for the describing function:

$$FDF(\theta_m) = \frac{2\beta_m}{1 + \sqrt{1 + \kappa_m^2 A_p (\theta_m)^2}} - \epsilon_m \exp\left(-\frac{\mu_m A_p (\theta_m)^2}{\sqrt{\pi} \ln(2)}\right). \quad (8)$$

In Appendix B the role of the parameters $\beta_m$, $\kappa_m$, $\epsilon_m$, $\mu_m$ is explained and it is also shown that this function corresponds to a model where the unsteady rate of heat release saturates as the pressure amplitude increases. After multiplying Eq. (8) by the impedance $Z$, one can fit the model to the experimental data. The result is shown in Fig. 9 with a blue line and labelled $FDF_A$. The model accurately captures the gain with increasing amplitude, but completely disregards the development of the phase; the phase shown is only due to the impedance $Z$. One can include a function that accounts for the deviation in the phase. However, to keep the model simple and to facilitate the forthcoming analysis Eq. (8) is kept real.

The corresponding plot of $FDF_A$ as a function of $A_p$ is also shown in Appendix B, Fig. B.14b in blue. In this figure one can observe that the FDF in terms of pressure fluctuations is monotonic. Using the tools developed by Ghirardo et al. [21] one finds that, in a symmetric combustor, real FDFs, such as Eq. (8), that decrease monotonically with amplitude only saturate into spinning modes, whereas real FDFs that are non-monotonic can saturate into both standing and spinning modes. Therefore, to ensure that standing modes are observed in the forthcoming analysis the parameter $\epsilon$ is increased until a non-monotonic behaviour is obtained. The result is shown in Figs. 9 and B.14 b in red and labelled $FDF_B$.

Since the model FDF is real, there are two implications on the dynamical system. First, $F_j^{(m)} = 0$ in Eq. (4). Second, the thermoacoustic frequency and the natural frequency<fnref n="1" /> coincide and one may set $\omega = \omega_n$.

## 7.2. Simplified dynamical system

With the previous model assumptions for the FDF and considering the first azimuthal mode corresponding to $n = 1$, the dynamical system reduces to:

$$\begin{bmatrix} \dot{A} \\ \dot{\theta}_0 \\ \dot{\chi} \\ \dot{\varphi} \end{bmatrix} = \begin{bmatrix} -\alpha \frac{A}{2} \\ 0 \\ 0 \\ 0 \end{bmatrix} + \sum_{m=1}^{N_b} \frac{F_r^{(m)}}{2} \begin{bmatrix} A + A \cos(2\chi) \cos(2(\theta_0 - \theta_m)) \\ -\sec(2\chi) \sin(2(\theta_0 - \theta_m)) \\ -\sin(2\chi) \cos(2(\theta_0 - \theta_m)) \\ \tan(2\chi) \sin(2(\theta_0 - \theta_m)) \end{bmatrix}. \quad (9)$$

Before proceeding to the discussion of the results it is important to stress the main advantage of the quaternion formulation of the pressure fluctuations. As can be seen from Eqs. (4) and (9), the variables used to describe the nature of the azimuthal instabilities are state space variables [30]. In other words, from this representation one can directly obtain with no further derivations, the time evolution of (1) the amplitude, (2) pressure antinodal line angular position, (3) spinning, standing or mixed nature of the mode and (4) the phase of the azimuthal oscillation. This is opposed to using derived quantities such as the C-indicator or the spin ratio which require some treatment of the pressure signal as done, for instance, in Worth and Dawson [8].

## 7.3. Phase portraits of the symmetric combustor

To illustrate the evolution of the dynamical system, the phase portraits are shown in two planes of the Poincaré sphere defined by $(A, 2\chi, \theta_0)$ as shown in Fig. 10a. The first plane, shown in red, is the meridian plane $(A, 2\chi, 0) \cup (A, 2\chi, \pm\pi)$ or alternatively $(S_1, S_3)$ in the Cartesian coordinate system. The attractors in this plane can be:

a) Spinning modes if $|2\chi| = \pi/2$. In the figure the two of them are shown with blue round markers.

b) Mixed modes if they lie in the region $0 < |2\chi| < \pi/2$. One of them (in a different plane) is shown with a green round marker.

c) Standing modes if $2\chi = 0$. Two of them are shown with red round markers.

The second plane, shown in blue, is the equatorial plane $(A, 0, \theta_0)$ or alternatively $(S_1, S_2)$ in the Cartesian coordinate system. All the attractors in this plane correspond to standing modes $(2\chi = 0)$.

Given a damping coefficient $\alpha$, the system can display qualitatively different phase portraits. In Fig. 10b we present the phase portraits at $\alpha \bar{p}/((\gamma - 1)\bar{q}) = 2.1 \times 10^4$. This value is chosen:

1. Relative to the magnitude of the FDF. That is to say that we can always re-scale the FDF and $\alpha$, and the phase portraits will remain the same, as long as their relative magnitude is kept the same.

2. To reproduce as much as possible some of the conditions observed in the experimental campaign.

For a symmetric combustor, the top and bottom rows of Fig. 10b show the meridian and the equatorial planes respectively. The dynamical system has stable standing modes and stable spinning modes. On the meridian plane (top) there are stable attractors at the origin and at $(25.57, 2 \times (\pm\pi/4), 0)$ corresponding to spinning modes. On the equatorial plane (bottom) there are stable attractors at the origin and at $(19.81, 0, (2k - 17)\pi/16)$ with $k = 1, 2, \dots, N_b$ corresponding to standing modes. These positions coincide with angles in between burners. As in the experimental cases (see Fig. 7) the spinning modes have larger amplitude than the standing modes.<fnref n="2" /> Furthermore, at small amplitudes the system is linearly stable. Observe that all the trajectories in this region point towards the origin in both planes. However, at moderate amplitudes the system can saturate into standing or spinning modes.<fnref n="3" /> The initial conditions of the system will determine the type of limit cycle at which the system saturates. This is reminiscent of the dual mode region in Fig. 6 and thus the hysteresis region analysed by Prieur et al. [28].

<footnote n="1">¹ The thermoacoustic frequency and the natural frequency of a system are usually close to each other. The acoustic or natural frequency of the MICCA combustor has been numerically computed to be around 450 Hz [35], while the thermoacoustic frequencies lie between 420 and 510 Hz according to Fig. 6.</footnote>

<footnote n="2">² While the amplitude $A$ has $A^{sp} = 25.57 > A^{st} = 19.81$, the pressure distribution around the annulus $A_p$ has $A_p^{sp} = 25.57/\sqrt{2} < \max(A_p^{st}) = 19.81$.</footnote>

<footnote n="3">³ Under these conditions the system is capable of triggering.</footnote>
<figure><img src="image_19.png" /><figcaption>Fig. 11. Region in the vicinity of the mixed modes.</figcaption></figure>

## 7.4. Phase portraits of the combustor under weak broken symmetry conditions

In the experimental cases P1A, P1B and P1C, one of the injectors was partially blocked. To represent this in the model we assume that blocking half of the injector corresponds to a reduction of the parameter $\epsilon_m$ in the FDF. As explained in Appendix B this translates into a small change in the rate of saturation. The FDF with a reduction of approximately 33% of the original value of $\epsilon_m$ is plotted in blue in Fig. 9 and labelled $FDF_A$. The main difference is that at low amplitudes the magnitude of the FDF is slightly higher. To analyse the consequences on the phase portraits, we consider configuration P1A-01. This configuration features the partially blocked injector at an angle $\theta_1 = 0$ and therefore its effects can be analysed in the same planes used for the symmetric combustor.

For the same damping coefficient as in the symmetric case, the new phase portraits are shown in Fig. 10c. For reference, the partially blocked injector is indicated with a purple arrow in the equatorial plane. The new phase portraits show mixed modes and stable standing modes. In the meridian plane<fnref n="4" /> (top) the spinning modes split into two and the new modes are slightly shifted, see Fig. 11. The new modes are mixed modes with a predominantly spinning component. They are located at $(25.60, 2 \times \pm 0.782, 0)$, with a preferred pressure antinodal line orientation in the direction $\theta_0 = 0$, i.e., the direction of the partially blocked injector. This differs from the experimental data where the antinode angle of the spinning modes is roughly 90 degrees apart from the partially blocked injector (see Fig. 8c).

In addition, notice that as $|\chi| \to \pi/4$ the second and fourth elements of Eq. (9) become infinitely large. As mentioned at the end of Appendix A, there cannot exist fixed points with $|\chi| = \pi/4$ if the symmetry is broken, except at $A = 0$. Therefore, due to the huge gradients in the vicinity of this region, the saddle node at the origin extends its field of influence to all amplitudes, thus explaining the dynamics of the trajectories seen in Fig. 11.

The equatorial plane presents the most notable change. The attractor at the origin splits, and two standing modes emerge at $(2.00, 0, 0)$. Moreover, the attractors on either side of the injectors located at $\theta = \pm \pi/2$ disappear. Furthermore, the rest of the stable standing modes have an amplitude of $A = 19.93$ and are located in the following positions:

* $\approx \pi/16$ radians away from the perturbed burner, i.e., in between the perturbed injector and the adjacent injectors. This is similar to what is observed in cases P1A-01, P1B-01 and P1C-01 and with the appropriate rotation of the frame of reference cases P1A-03, P1A-05, P1B-03, D1-01, D1-02 from Fig. 8a.

<figure><img src="image_20.png" /><figcaption>Fig. 12. Comparison of pressure antinode (orientation angle $\theta_0$) from experiments and from the model. The inner part of the annulus shows the experimental results from Fig. 8a. The outer part of the annulus shows in green the pressure antinodes predicted by the model. Note: For configuration P2-02 the experimental results show that there is a pressure antinode line shown in very light gray (barely visible) just below the burner located at $\theta = 0$. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article)</figcaption></figure>

* $\approx 3\pi/16$ radians away from the perturbed injector, i.e., in between two injectors at a position one injector away from the perturbed unit. With the appropriate rotation of the frame of reference, this is similar to case P1A-04.

In the previous analysis we perturbed a single FDF to simulate one partially blocked injector. One may extend the analysis and perturb several FDFs to emulate configurations P2-01, P2-02 and P3-01. Disregarding the stable standing modes at very low amplitudes, the pressure antinodes predicted by the model compared to the experiments are shown in Fig. 12. As explained above, for the configuration P1A-01 the model predicts four stable pressure antinodal line orientations out of which one matches the experimental results. For configuration P2-01 the model predicts four stable pressure antinodal line positions and in the experiment one observes two antinode line positions. The most probable one (shown in black) is not captured by the model, but the least probable (shown in light gray) is. For configuration P2-02 the model predicts four antinodal line orientations and in the experiments one observes two. Both of them are correctly predicted by the model. Finally for configuration P3-01 the model predicts two antinodal lines out of which one is observed in experiments.

## 7.5. Phase portraits of the combustor under strong broken symmetry conditions

The previous analysis reveals that the orientation angle $\theta_0$ of the stable standing modes is very sensitive to small changes in the distribution of unsteady heat release rate, if the original distribution is azimuthally homogeneous. As will be shown next, if one starts from a distribution that is not azimuthally homogeneous the sensitivity of the orientation angle is greatly reduced. However, before one can assess the sensitivity it is natural to define a measure of the azimuthal inhomogeneity of the system.

The strength of the variation in the distribution of the unsteady heat release rate may be estimated by means of the $C_{2n}$ indicator<fnref n="5" />

<footnote n="4">⁴ The meridian plane is symmetric along the vertical line. Therefore all the comments referring to the right side ($\theta = 0$) are also applicable to the left side ($\theta = |\pi|$).</footnote>

<footnote n="5">⁵ The $C_{2n}$ indicator is chosen because all the changes are introduced in the linear parameters of the FDF ($\epsilon$ and $\beta$) if the variations were introduced in the nonlinear parameters, such as $\kappa$ or $\mu$, the $N_{2n}$ indicator [21] would be more appropriate.</footnote>
<figure><img src="image_20.png" /><figcaption>**Fig. 11.** Region in the vicinity of the mixed modes.</figcaption></figure>

### 7.4. Phase portraits of the combustor under weak broken symmetry conditions

In the experimental cases P1A, P1B and P1C, one of the injectors was partially blocked. To represent this in the model we assume that blocking half of the injector corresponds to a reduction of the parameter $\epsilon_m$ in the FDF. As explained in Appendix B this translates into a small change in the rate of saturation. The FDF with a reduction of approximately 33% of the original value of $\epsilon_m$ is plotted in blue in Fig. 9 and labelled $FDF_A$. The main difference is that at low amplitudes the magnitude of the FDF is slightly higher. To analyse the consequences on the phase portraits, we consider configuration P1A-01. This configuration features the partially blocked injector at an angle $\theta_1 = 0$ and therefore its effects can be analysed in the same planes used for the symmetric combustor.

For the same damping coefficient as in the symmetric case, the new phase portraits are shown in Fig. 10c. For reference, the partially blocked injector is indicated with a purple arrow in the equatorial plane. The new phase portraits show mixed modes and stable standing modes. In the meridian plane<fnref n="4" /> (top) the spinning modes split into two and the new modes are slightly shifted, see Fig. 11. The new modes are mixed modes with a predominantly spinning component. They are located at $(25.60, 2 \times \pm 0.782, 0)$, with a preferred pressure antinodal line orientation in the direction $\theta_0 = 0$, i.e., the direction of the partially blocked injector. This differs from the experimental data where the antinode angle of the spinning modes is roughly 90 degrees apart from the partially blocked injector (see Fig. 8c).

In addition, notice that as $|\chi| \to \pi/4$ the second and fourth elements of Eq. (9) become infinitely large. As mentioned at the end of Appendix A, there cannot exist fixed points with $|\chi| = \pi/4$ if the symmetry is broken, except at $A = 0$. Therefore, due to the huge gradients in the vicinity of this region, the saddle node at the origin extends its field of influence to all amplitudes, thus explaining the dynamics of the trajectories seen in Fig. 11.

The equatorial plane presents the most notable change. The attractor at the origin splits, and two standing modes emerge at $(2.00, 0, 0)$. Moreover, the attractors on either side of the injectors located at $\theta = \pm \pi/2$ disappear. Furthermore, the rest of the stable standing modes have an amplitude of $A = 19.93$ and are located in the following positions:

* $\approx \pi/16$ radians away from the perturbed burner, i.e., in between the perturbed injector and the adjacent injectors. This is similar to what is observed in cases P1A-01, P1B-01 and P1C-01 and with the appropriate rotation of the frame of reference cases P1A-03, P1A-05, P1B-03, D1-01, D1-02 from Fig. 8a.

<figure><img src="image_21.png" /><figcaption>**Fig. 12.** Comparison of pressure antinode (orientation angle $\theta_0$) from experiments and from the model. The inner part of the annulus shows the experimental results from Fig. 8a. The outer part of the annulus shows in green the pressure antinodes predicted by the model. Note: For configuration P2-02 the experimental results show that there is a pressure antinode line shown in very light gray (barely visible) just below the burner located at $\theta = 0$. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article)</figcaption></figure>

* $\approx 3\pi/16$ radians away from the perturbed injector, i.e., in between two injectors at a position one injector away from the perturbed unit. With the appropriate rotation of the frame of reference, this is similar to case P1A-04.

In the previous analysis we perturbed a single FDF to simulate one partially blocked injector. One may extend the analysis and perturb several FDFs to emulate configurations P2-01, P2-02 and P3-01. Disregarding the stable standing modes at very low amplitudes, the pressure antinodes predicted by the model compared to the experiments are shown in Fig. 12. As explained above, for the configuration P1A-01 the model predicts four stable pressure antinodal line orientations out of which one matches the experimental results. For configuration P2-01 the model predicts four stable pressure antinodal line positions and in the experiment one observes two antinode line positions. The most probable one (shown in black) is not captured by the model, but the least probable (shown in light gray) is. For configuration P2-02 the model predicts four antinodal line orientations and in the experiments one observes two. Both of them are correctly predicted by the model. Finally for configuration P3-01 the model predicts two antinodal lines out of which one is observed in experiments.

### 7.5. Phase portraits of the combustor under strong broken symmetry conditions

The previous analysis reveals that the orientation angle $\theta_0$ of the stable standing modes is very sensitive to small changes in the distribution of unsteady heat release rate, if the original distribution is azimuthally homogeneous. As will be shown next, if one starts from a distribution that is not azimuthally homogeneous the sensitivity of the orientation angle is greatly reduced. However, before one can assess the sensitivity it is natural to define a measure of the azimuthal inhomogeneity of the system.

The strength of the variation in the distribution of the unsteady heat release rate may be estimated by means of the $C_{2n}$ indicator<fnref n="5" />

<footnote n="4">⁴ The meridian plane is symmetric along the vertical line. Therefore all the comments referring to the right side ($\theta = 0$) are also applicable to the left side ($\theta = |\pi|$).</footnote>

<footnote n="5">⁵ The $C_{2n}$ indicator is chosen because all the changes are introduced in the linear parameters of the FDF ($\epsilon$ and $\beta$) if the variations were introduced in the nonlinear parameters, such as $\kappa$ or $\mu$, the $N_{2n}$ indicator [21] would be more appropriate.</footnote>
<figure><img src="image_21.png" /><figcaption>**Fig. 13.** Reference for (a) see Fig. 10. (b) and (c) show trajectories of the fixed points as the FDF is varied. × denotes the initial location of the fixed points and ○ denotes the final position. The thick lines are used to trace the trajectories of the fixed points as one of the parameters is varied. In red we show the stable fixed points (attractors) and in black the unstable fixed points (saddles and repellers). Starting from a configuration with strong azimuthal variation of the FDF, in (b) the fixed points barely move when $ε_1$ is reduced, in (c) the fixed points change drastically when B is lowered and the strong azimuthal variation is removed. Starting from 7 fixed points, decreasing the value of B leads to 23 fixed points due to the occurrence of multiple bifurcations. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article)</figcaption></figure>

<table><caption>**Table 2**
Summary of the strength of the azimuthal variations in the heat release rate as given by the $C_{2n}/C_0$ index.</caption><thead><tr><th>Configuration</th><th>Condition</th><th>C<sub>2n</sub>/C<sub>0</sub></th></tr></thead><tbody><tr><td>ε<sub>m</sub> = 2155</td><td>Symmetric</td><td>0</td></tr><tr><td>ε<sub>1</sub> = 1437</td><td>Weak broken symmetry</td><td>0.034</td></tr><tr><td>ε<sub>m</sub> = 2155(1 + sin(2θ<sub>m</sub>))</td><td>Strong broken symmetry</td><td>0.472</td></tr></tbody></table>

introduced by Noiray et al. [18]. Since the changes also affect the mean flame response, it is more instructive to look at the ratio $C_{2n}/C_0$. For a discrete number of burners with a real FDF this can be calculated as:

$$C_{2n} = \sum_{m=1}^{N_b} \text{FDF}(\theta_m) \sin(2n(\theta_m - \theta_0)), \quad C_0 = \sum_{m=1}^{N_b} \text{FDF}(\theta_m) \quad (10)$$

The FDF is evaluated at $A_p = 0$, the azimuthal wave number is set to $n = 1$, and $θ_0$ is chosen to maximize the value of $C_{2n}$. For the symmetric combustor this index vanishes $C_{2n}/C_0 = 0$. For the combustor where the FDF of the burner located at $θ_1$ was modified by changing the parameter $ε_1$ from 2155 to 1437, one finds $C_{2n}/C_0 = 0.034$, i.e., weak broken symmetry conditions. Next consider a case with strong broken symmetry conditions, which is obtained by introducing an azimuthal variation in the FDF which increases the magnitude of its $2n$ Fourier component [18,25]. This can be achieved by modifying the parameter $ε_m$ such that $ε_m = 2155(1 + B \sin(2θ_m))$. Setting $B = 1$ one finds that $C_{2n}/C_0 = 0.472$. The results for the different configurations are summarized in Table 2.

To assess the sensitivity of the pressure antinodal line to small changes in the distribution of unsteady heat release rate, one may begin with the combustor under strong broken symmetry conditions. For this configuration the equatorial plane is shown in Fig. 13a. There are two stable standing modes located at $(18.48, 0, -π/4)$ and $(18.48, 0, 3π/4)$. From this configuration one may proceed by slightly decreasing the parameter $ε_1$ from 2155 to 1437 and trace the trajectories of the fixed points as shown in Fig. 13b. The effects are clearly negligible. In contrast, consider now this configuration but with a decrease of the parameter B to zero. The trajectories, traced in Fig. 13c, indicate that the changes to the fixed points are very prominent and, as expected, one recovers the

plane shown in Fig. 10c (bottom). This analysis confirms that in configurations with a homogeneous distribution of unsteady heat release rate small changes to the distribution can translate to big changes in the phase portraits due to the relocation of the fixed points.

Extending these conclusions to the experimental cases, one can argue that the MICCA combustor does not have a strong azimuthal variation of the unsteady heat release rate, because the position of the pressure antinodes change when the injection pattern changes. Furthermore, as shown in the last example, where a strong azimuthal variation is introduced in the FDF which increases the magnitude of its $2n$ Fourier component, strong asymmetries are required to completely change the phase portrait of the system. This reinforces the possibility that an assortment of weak symmetry breakers are required to change the stability of a limit cycle as seen in the standing case of the configuration DP2-01.

## 8. Conclusions

This study comprises an experimental investigation of the effects of symmetry breaking on the stability of limit cycles. The desired effect is achieved by partially blocking one or several injectors or by changing their geometry in the MICCA combustor. Results show that breaking the symmetry by changing the distribution of the unsteady heat release rate does not change the stability of a predetermined limit cycle when the symmetry breaking mechanism is weak. Under this scenario the main feature that responds to these changes is the pressure antinodal line (orientation angle) which locks in at specific positions that depend on the pattern of the blocked injector(s).

By using a dynamical system and a modelled flame describing function to represent the combustion response under similar conditions, one can predict some features of the limit cycles that match the experimental cases. For the cases where a standing limit cycle develops, the model is capable of predicting several pressure antinodal line positions, some of which coincide with those found experimentally. Using the model we further conclude that for a standing azimuthal mode of order $n$, the strength of the $2n$ Fourier component of the azimuthal variation of the unsteady rate of heat release, determines the sensitivity of the pressure antinodal line position. When the azimuthal variation is weak, the pressure antinodal line orientation is very sensitive. When the variation is
<figure><img src="image_22.png" /><figcaption>**Fig. 13.** Reference for (a) see Fig. 10. (b) and (c) show trajectories of the fixed points as the FDF is varied. × denotes the initial location of the fixed points and ○ denotes the final position. The thick lines are used to trace the trajectories of the fixed points as one of the parameters is varied. In red we show the stable fixed points (attractors) and in black the unstable fixed points (saddles and repellers). Starting from a configuration with strong azimuthal variation of the FDF, in (b) the fixed points barely move when $ε_1$ is reduced, in (c) the fixed points change drastically when B is lowered and the strong azimuthal variation is removed. Starting from 7 fixed points, decreasing the value of B leads to 23 fixed points due to the occurrence of multiple bifurcations. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article)</figcaption></figure>

<table><caption>**Table 2**
Summary of the strength of the azimuthal variations in the heat release rate as given by the $C_{2n}/C_0$ index.</caption><thead><tr><th>Configuration</th><th>Condition</th><th>C<sub>2n</sub>/C<sub>0</sub></th></tr></thead><tbody><tr><td>ε<sub>m</sub> = 2155</td><td>Symmetric</td><td>0</td></tr><tr><td>ε<sub>1</sub> = 1437</td><td>Weak broken symmetry</td><td>0.034</td></tr><tr><td>ε<sub>m</sub> = 2155(1 + sin(2θ<sub>m</sub>))</td><td>Strong broken symmetry</td><td>0.472</td></tr></tbody></table>

introduced by Noiray et al. [18]. Since the changes also affect the mean flame response, it is more instructive to look at the ratio $C_{2n}/C_0$. For a discrete number of burners with a real FDF this can be calculated as:

$$C_{2n} = \sum_{m=1}^{N_b} \text{FDF}(\theta_m) \sin(2n(\theta_m - \theta_0)), \quad C_0 = \sum_{m=1}^{N_b} \text{FDF}(\theta_m) \quad (10)$$

The FDF is evaluated at $A_p = 0$, the azimuthal wave number is set to $n = 1$, and $\theta_0$ is chosen to maximize the value of $C_{2n}$. For the symmetric combustor this index vanishes $C_{2n}/C_0 = 0$. For the combustor where the FDF of the burner located at $\theta_1$ was modified by changing the parameter $ε_1$ from 2155 to 1437, one finds $C_{2n}/C_0 = 0.034$, i.e., weak broken symmetry conditions. Next consider a case with strong broken symmetry conditions, which is obtained by introducing an azimuthal variation in the FDF which increases the magnitude of its $2n$ Fourier component [18,25]. This can be achieved by modifying the parameter $ε_m$ such that $ε_m = 2155(1 + B \sin(2\theta_m))$. Setting $B = 1$ one finds that $C_{2n}/C_0 = 0.472$. The results for the different configurations are summarized in Table 2.

To assess the sensitivity of the pressure antinodal line to small changes in the distribution of unsteady heat release rate, one may begin with the combustor under strong broken symmetry conditions. For this configuration the equatorial plane is shown in Fig. 13a. There are two stable standing modes located at $(18.48, 0, -\pi/4)$ and $(18.48, 0, 3\pi/4)$. From this configuration one may proceed by slightly decreasing the parameter $ε_1$ from 2155 to 1437 and trace the trajectories of the fixed points as shown in Fig. 13b. The effects are clearly negligible. In contrast, consider now this configuration but with a decrease of the parameter B to zero. The trajectories, traced in Fig. 13c, indicate that the changes to the fixed points are very prominent and, as expected, one recovers the

plane shown in Fig. 10c (bottom). This analysis confirms that in configurations with a homogeneous distribution of unsteady heat release rate small changes to the distribution can translate to big changes in the phase portraits due to the relocation of the fixed points.

Extending these conclusions to the experimental cases, one can argue that the MICCA combustor does not have a strong azimuthal variation of the unsteady heat release rate, because the position of the pressure antinodes change when the injection pattern changes. Furthermore, as shown in the last example, where a strong azimuthal variation is introduced in the FDF which increases the magnitude of its $2n$ Fourier component, strong asymmetries are required to completely change the phase portrait of the system. This reinforces the possibility that an assortment of weak symmetry breakers are required to change the stability of a limit cycle as seen in the standing case of the configuration DP2-01.

## 8. Conclusions

This study comprises an experimental investigation of the effects of symmetry breaking on the stability of limit cycles. The desired effect is achieved by partially blocking one or several injectors or by changing their geometry in the MICCA combustor. Results show that breaking the symmetry by changing the distribution of the unsteady heat release rate does not change the stability of a predetermined limit cycle when the symmetry breaking mechanism is weak. Under this scenario the main feature that responds to these changes is the pressure antinodal line (orientation angle) which locks in at specific positions that depend on the pattern of the blocked injector(s).

By using a dynamical system and a modelled flame describing function to represent the combustion response under similar conditions, one can predict some features of the limit cycles that match the experimental cases. For the cases where a standing limit cycle develops, the model is capable of predicting several pressure antinodal line positions, some of which coincide with those found experimentally. Using the model we further conclude that for a standing azimuthal mode of order $n$, the strength of the $2n$ Fourier component of the azimuthal variation of the unsteady rate of heat release, determines the sensitivity of the pressure antinodal line position. When the azimuthal variation is weak, the pressure antinodal line orientation is very sensitive. When the variation is
strong the pressure antinodal line position barely changes. These results confirm that the unsteady heat release distribution in the MICCA combustor is rather homogeneous given that the pressure antinodal line position changes according to the position of the blocked injectors.

## Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgments

The authors would like to thank Davide Laera for providing the experimental FDF data.

## Appendix A. Derivation of the dynamical system

Previous derivations of the dynamical system [25,31,33] have been carried out using quaternion algebra, which adds a layer of complexity for those who are not familiar with this formalism. In this appendix we avoid these complexities and sketch the derivation of the dynamical system for the case with a discrete number of burners using multiple scales theory. We begin by considering the wave equation:

$$ \frac{\partial^2 p'}{\partial t^2} + \alpha \frac{\partial p'}{\partial t} - \frac{\bar{c}^2}{R^2} \frac{\partial^2 p'}{\partial \theta^2} = (\gamma - 1) \sum_{m=1}^{N_b} \frac{\partial q'}{\partial t} 2\pi \delta(\theta - \theta_m), \quad (A.1) $$

The nonlinearity in the equation arises from the unsteady heat release term. As such, this equation can be treated as a weakly nonlinear oscillator. This problem is dominated by two time scales, a fast time scale $t_f = 1/\omega$ which is proportional to the period of oscillation and a slow time scale $t_s = 2/(\alpha - \eta)$ associated with the growth rate and reflected in the rate of change of the state variables. Here $\eta$ is the linear coefficient of the r.h.s. of Eq. (A.1) proportional to $\partial p'/\partial t$. Using the following shorthand notation:

$$ \begin{aligned} c_{\theta_0} &= \cos(n(\theta - \theta_0)), & s_{\theta_0} &= \sin(n(\theta - \theta_0)), \\ c_{\chi} &= \cos(\chi), & s_{\chi} &= \sin(\chi), \\ c_{\varphi} &= \cos(\omega t + \varphi), & s_{\varphi} &= \sin(\omega t + \varphi), \end{aligned} $$

One may now introduce the pressure ansatz and write its first derivative with respect to time:

$$ p'(t, \theta) = A(c_{\theta_0} c_{\chi} c_{\varphi} + s_{\theta_0} s_{\chi} s_{\varphi}), \quad (A.2) $$

$$ \frac{\partial p'}{\partial t} = A\omega(s_{\theta_0} s_{\chi} c_{\varphi} - c_{\theta_0} c_{\chi} s_{\varphi}). \quad (A.3) $$

<figure><img src="image_22.png" /><figcaption>Fig. B.14. (a) Heat release rate as a function of pressure amplitude for the saturation models given by Eqs. (B.3) and (B.5). (b) Flame describing function of each of the saturation models. The values for the parameters are: $\beta = 3419$, $\kappa = 0.2156$ 1/Pa, $\epsilon_A = 1437$, $\epsilon_B = 2155$ and $\mu = 0.01712$ 1/Pa$^2$. The dashed lines are the analytic approximation to FDF$_{A,B}$ given by Eq. (8).</figcaption></figure>

These two expressions may be regarded as defining the slow flow variables $A, \theta_0, \chi, \varphi$. Eq. (A.3) requires the following equation to hold:

$$ \begin{aligned} 0 = \dot{A}(c_{\theta_0} c_{\chi} c_{\varphi} + s_{\theta_0} s_{\chi} s_{\varphi}) + An\dot{\theta}_0(s_{\theta_0} c_{\chi} c_{\varphi} - c_{\theta_0} s_{\chi} s_{\varphi}) \\ + A\dot{\chi}(s_{\theta_0} c_{\chi} s_{\varphi} - c_{\theta_0} s_{\chi} c_{\varphi}) + A\dot{\varphi}(s_{\theta_0} s_{\chi} c_{\varphi} - c_{\theta_0} c_{\chi} s_{\varphi}). \end{aligned} \quad (A.4) $$

From Eq. (A.3) it follows that:

$$ \begin{aligned} \frac{\partial^2 p'}{\partial t^2} ={}& \omega(\dot{A}(s_{\theta_0} s_{\chi} c_{\varphi} - c_{\theta_0} c_{\chi} s_{\varphi}) - An\dot{\theta}_0(c_{\theta_0} s_{\chi} c_{\varphi} + s_{\theta_0} c_{\chi} s_{\varphi}) \\ & + A\dot{\chi}(s_{\theta_0} c_{\chi} c_{\varphi} + c_{\theta_0} s_{\chi} s_{\varphi}) - A(\dot{\varphi} + \omega)(c_{\theta_0} c_{\chi} c_{\varphi} + s_{\theta_0} s_{\chi} s_{\varphi})). \end{aligned} \quad (A.5) $$

Before averaging the equations, we first treat the unsteady heat release term by considering only the first harmonic:

$$ q'(\theta, t) = \frac{1}{2} (\hat{q} e^{j\omega t} + \text{c.c.}). \quad (A.6) $$

Using the definition of the flame describing function Eq. (6), we substitute $\hat{q}$ in Eq. (A.6) to obtain:

$$ \frac{\partial q'}{\partial t} = \frac{1}{2} \left( \bar{q} \text{ FDF} \frac{\hat{p}}{\bar{p}} e^{j\omega t} j\omega + \text{c.c.} \right). \quad (A.7) $$

Here $\hat{p}$ is defined via $p' = \text{Re}\{\hat{p}e^{j\omega t}\}$ and equals:

$$ \hat{p} = A(c_{\theta_0} c_{\chi} + j s_{\theta_0} s_{\chi})(\cos(\varphi) + j \sin(\varphi)). \quad (A.8) $$

Proceed by substituting the above definitions into (A.1). To remove the dependency on the azimuthal angle $\theta$ we multiply Eqs. (A.1) and (A.4) by $\cos(n\theta)$ and average over the circumference. We repeat the same process with $\sin(n\theta)$, providing a total of four equations for the four slow flow variables. Since the derivatives of the slow flow variables appear linearly in the equations, the set of equations can be recast in vector form as:

$$ M_1(\mathbf{x}, t) \dot{\mathbf{x}} = \mathbf{b}_1(\mathbf{x}, t). \quad (A.9) $$

Where $\mathbf{x} = [A, \theta_0, \chi, \varphi]^T$, $M_1$ is a matrix and $\mathbf{b}_1$ is a vector both of which depend on time and $\mathbf{x}$. We proceed by splitting the matrix $M_1$ into a matrix that contains the temporal terms ($c_\varphi$ and $s_\varphi$), say $B_2(\mathbf{x}, t)$, and a matrix that contains the rest of the terms, say $M_2(\mathbf{x})$, such that $M_1(\mathbf{x}, t) = B_2(\mathbf{x}, t)M_2(\mathbf{x})$. The matrix $B_2$ is invertible in the entire domain of the slow flow variables, therefore we can multiply Eq. (A.9) by $B_2^{-1}$ from the left to obtain:

$$ M_2(\mathbf{x}) \dot{\mathbf{x}} = (B_2(\mathbf{x}, t))^{-1} \mathbf{b}_1(\mathbf{x}, t). \quad (A.10) $$

We perform temporal averaging over one cycle on all the terms in the right:

$$ \mathbf{b}_2(\mathbf{x}) = \frac{\omega}{2\pi} \int_0^{2\pi} (B_2(\mathbf{x}, t))^{-1} \mathbf{b}_1(\mathbf{x}, t) \, dt. \quad (A.11) $$

<figure><img src="image_23.png" /></figure>
strong the pressure antinodal line position barely changes. These results confirm that the unsteady heat release distribution in the MICCA combustor is rather homogeneous given that the pressure antinodal line position changes according to the position of the blocked injectors.

## Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgments

The authors would like to thank Davide Laera for providing the experimental FDF data.

## Appendix A. Derivation of the dynamical system

Previous derivations of the dynamical system [25,31,33] have been carried out using quaternion algebra, which adds a layer of complexity for those who are not familiar with this formalism. In this appendix we avoid these complexities and sketch the derivation of the dynamical system for the case with a discrete number of burners using multiple scales theory. We begin by considering the wave equation:

$$ \frac{\partial^2 p'}{\partial t^2} + \alpha \frac{\partial p'}{\partial t} - \frac{\bar{c}^2}{R^2} \frac{\partial^2 p'}{\partial \theta^2} = (\gamma - 1) \sum_{m=1}^{N_b} \frac{\partial q'}{\partial t} 2\pi \delta(\theta - \theta_m), \quad (A.1) $$

The nonlinearity in the equation arises from the unsteady heat release term. As such, this equation can be treated as a weakly nonlinear oscillator. This problem is dominated by two time scales, a fast time scale $t_f = 1/\omega$ which is proportional to the period of oscillation and a slow time scale $t_s = 2/(\alpha - \eta)$ associated with the growth rate and reflected in the rate of change of the state variables. Here $\eta$ is the linear coefficient of the r.h.s. of Eq. (A.1) proportional to $\partial p'/\partial t$. Using the following shorthand notation:

$$ \begin{aligned} c_{\theta_0} &= \cos(n(\theta - \theta_0)), & s_{\theta_0} &= \sin(n(\theta - \theta_0)), \\ c_{\chi} &= \cos(\chi), & s_{\chi} &= \sin(\chi), \\ c_{\varphi} &= \cos(\omega t + \varphi), & s_{\varphi} &= \sin(\omega t + \varphi), \end{aligned} $$

One may now introduce the pressure ansatz and write its first derivative with respect to time:

$$ p'(t, \theta) = A(c_{\theta_0} c_{\chi} c_{\varphi} + s_{\theta_0} s_{\chi} s_{\varphi}), \quad (A.2) $$

$$ \frac{\partial p'}{\partial t} = A\omega(s_{\theta_0} s_{\chi} c_{\varphi} - c_{\theta_0} c_{\chi} s_{\varphi}). \quad (A.3) $$

<figure><img src="image_23.png" /><figcaption>Fig. B.14. (a) Heat release rate as a function of pressure amplitude for the saturation models given by Eqs. (B.3) and (B.5). (b) Flame describing function of each of the saturation models. The values for the parameters are: $\beta = 3419$, $\kappa = 0.2156$ 1/Pa, $\epsilon_A = 1437$, $\epsilon_B = 2155$ and $\mu = 0.01712$ 1/Pa$^2$. The dashed lines are the analytic approximation to FDF$_{A,B}$ given by Eq. (8).</figcaption></figure>

These two expressions may be regarded as defining the slow flow variables $A, \theta_0, \chi, \varphi$. Eq. (A.3) requires the following equation to hold:

$$ \begin{aligned} 0 = \dot{A}(c_{\theta_0} c_{\chi} c_{\varphi} + s_{\theta_0} s_{\chi} s_{\varphi}) + An\dot{\theta}_0(s_{\theta_0} c_{\chi} c_{\varphi} - c_{\theta_0} s_{\chi} s_{\varphi}) \\ + A\dot{\chi}(s_{\theta_0} c_{\chi} s_{\varphi} - c_{\theta_0} s_{\chi} c_{\varphi}) + A\dot{\varphi}(s_{\theta_0} s_{\chi} c_{\varphi} - c_{\theta_0} c_{\chi} s_{\varphi}). \end{aligned} \quad (A.4) $$

From Eq. (A.3) it follows that:

$$ \frac{\partial^2 p'}{\partial t^2} = \omega(\dot{A}(s_{\theta_0} s_{\chi} c_{\varphi} - c_{\theta_0} c_{\chi} s_{\varphi}) - An\dot{\theta}_0(c_{\theta_0} s_{\chi} c_{\varphi} + s_{\theta_0} c_{\chi} s_{\varphi}) \\ + A\dot{\chi}(s_{\theta_0} c_{\chi} c_{\varphi} + c_{\theta_0} s_{\chi} s_{\varphi}) - A(\dot{\varphi} + \omega)(c_{\theta_0} c_{\chi} c_{\varphi} + s_{\theta_0} s_{\chi} s_{\varphi})). \quad (A.5) $$

Before averaging the equations, we first treat the unsteady heat release term by considering only the first harmonic:

$$ q'(\theta, t) = \frac{1}{2}(\hat{q} e^{j\omega t} + \text{c.c.}). \quad (A.6) $$

Using the definition of the flame describing function Eq. (6), we substitute $\hat{q}$ in Eq. (A.6) to obtain:

$$ \frac{\partial q'}{\partial t} = \frac{1}{2} \left( \bar{q} \text{ FDF} \frac{\hat{p}}{p} e^{j\omega t} j\omega + \text{c.c.} \right). \quad (A.7) $$

Here $\hat{p}$ is defined via $p' = \text{Re}\{\hat{p}e^{j\omega t}\}$ and equals:

$$ \hat{p} = A(c_{\theta_0} c_{\chi} + j s_{\theta_0} s_{\chi})(\cos(\varphi) + j \sin(\varphi)). \quad (A.8) $$

Proceed by substituting the above definitions into (A.1). To remove the dependency on the azimuthal angle $\theta$ we multiply Eqs. (A.1) and (A.4) by $\cos(n\theta)$ and average over the circumference. We repeat the same process with $\sin(n\theta)$, providing a total of four equations for the four slow flow variables. Since the derivatives of the slow flow variables appear linearly in the equations, the set of equations can be recast in vector form as:

$$ M_1(\mathbf{x}, t) \dot{\mathbf{x}} = \mathbf{b}_1(\mathbf{x}, t). \quad (A.9) $$

Where $\mathbf{x} = [A, \theta_0, \chi, \varphi]^T$, $M_1$ is a matrix and $\mathbf{b}_1$ is a vector both of which depend on time and $\mathbf{x}$. We proceed by splitting the matrix $M_1$ into a matrix that contains the temporal terms ($c_\varphi$ and $s_\varphi$), say $B_2(\mathbf{x}, t)$, and a matrix that contains the rest of the terms, say $M_2(\mathbf{x})$, such that $M_1(\mathbf{x}, t) = B_2(\mathbf{x}, t)M_2(\mathbf{x})$. The matrix $B_2$ is invertible in the entire domain of the slow flow variables, therefore we can multiply Eq. (A.9) by $B_2^{-1}$ from the left to obtain:

$$ M_2(\mathbf{x}) \dot{\mathbf{x}} = (B_2(\mathbf{x}, t))^{-1} \mathbf{b}_1(\mathbf{x}, t). \quad (A.10) $$

We perform temporal averaging over one cycle on all the terms in the right:

$$ \mathbf{b}_2(\mathbf{x}) = \frac{\omega}{2\pi} \int_0^{2\pi} (B_2(\mathbf{x}, t))^{-1} \mathbf{b}_1(\mathbf{x}, t) \, dt. \quad (A.11) $$

<figure><img src="image_24.png" /></figure>
Matrix $M_2$ is not always invertible. Its determinant is $(A^3 n \omega^2 \cos(2\chi))/16$, which implies that the matrix cannot be inverted when $A=0$ or $|\chi|=\pi/4$. Therefore, we split the matrix into the product of an invertible matrix, say $B_3$, and a noninvertible matrix, say $M_3$, which contains all the terms that make the determinant zero or undefined, such that $M_2=B_3 M_3$. Multiplying from the left by $B_3^{-1}$ produces the following system of equations:

$$M_3(\mathbf{x})\dot{\mathbf{x}} = (B_3(\mathbf{x}))^{-1}\mathbf{b}_2(\mathbf{x}). \qquad (A.12)$$

If we consider a meridian plane $\theta_0$ = constant or the equatorial plane as in Fig. 10, the difficulties raised by the non-invertibility of matrix $M_3$ do not arise. Then, solving for $\dot{\mathbf{x}}$ produces the dynamical system shown in Eq. (4). For the case when $A=0$, Eq. (A.12) reduces to $\dot{A}=0$. For the limiting case $|\chi|=\pi/4$, $M_3$ is not invertible, then, after dividing Eq. (A.12) by A we obtain:

$$\begin{bmatrix} \frac{\dot{A}}{A} \\ \frac{\dot{\chi}}{\chi} \end{bmatrix} \begin{bmatrix} n\dot{\theta}_0 + \sin(2\chi)\dot{\varphi} \\ \cos(2\chi)\dot{\varphi} \end{bmatrix} = \begin{bmatrix} \left(\frac{\omega_n^2}{2\omega} - \frac{\omega}{2}\right)\sin(2\chi) \\ 0 \\ \left(\frac{\omega_n^2}{2\omega} - \frac{\omega}{2}\right)\cos(2\chi) \end{bmatrix} + \sum_{m=1}^{N_b} \frac{F_r^{(m)}}{2} \begin{bmatrix} 1 + \cos(2\chi)\cos(2n(\theta_0 - \theta_m)) \\ -\cos(2\chi)\sin(2n(\theta_0 - \theta_m)) \\ -\sin(2\chi)\cos(2n(\theta_0 - \theta_m)) \\ \sin(2\chi)\sin(2n(\theta_0 - \theta_m)) \\ 0 \\ \sin(2\chi) \\ \sin(2n(\theta_0 - \theta_m)) \\ \cos(2\chi) + \cos(2n(\theta_0 - \theta_m)) \end{bmatrix} + \sum_{m=1}^{N_b} \frac{F_r^{(m)}}{2} \begin{bmatrix} 1 \\ 0 \\ \mp\cos(2n(\theta_0 - \theta_m)) \\ \pm\sin(2n(\theta_0 - \theta_m)) \\ 0 \\ \pm 1 \\ \sin(2n(\theta_0 - \theta_m)) \\ \cos(2n(\theta_0 - \theta_m)) \end{bmatrix}, \quad (A.13)$$

with $F_r$ and $F_j$ defined in Eq. (5). This equation is the analogue of Eq. (19) in Ghirardo and Gant [31] for a discrete number of burners. Substituting $\chi = \pm\pi/4$ gives:

$$\begin{bmatrix} \frac{\dot{A}}{A} \\ n\dot{\theta}_0 \pm \dot{\varphi} \\ \dot{\chi} \\ 0 \end{bmatrix} = \begin{bmatrix} -\frac{\alpha}{2} \\ \pm\left(\frac{\omega_n^2}{2\omega} - \frac{\omega}{2}\right) \\ 0 \\ 0 \end{bmatrix} + \sum_{m=1}^{N_b} \frac{F_r^{(m)}}{2} \begin{bmatrix} 1 \\ 0 \\ \mp\cos(2n(\theta_0 - \theta_m)) \\ \pm\sin(2n(\theta_0 - \theta_m)) \end{bmatrix} + \sum_{m=1}^{N_b} \frac{F_j^{(m)}}{2} \begin{bmatrix} 0 \\ \pm 1 \\ \sin(2n(\theta_0 - \theta_m)) \\ \cos(2n(\theta_0 - \theta_m)) \end{bmatrix}, \quad (A.14)$$

The equation in the bottom row only holds in the case of a symmetric combustor. When the symmetry is broken there can only exist mixed modes in accordance to the findings in Noiray et al. [18], Faure-Beaulieu and Noiray [25] and the results shown in this study.

## Appendix B. Flame describing function

In this appendix two models for the saturation of the heat release rate as a function of pressure are examined. It is here convenient to write pressure fluctuations from Eq. (1) using the definition of the amplitude of the pressure fluctuations around the annulus Eq. (2):

$$p'(t, \theta) = A_p(\theta) \cos(\omega t + \varphi^+(\theta)). \quad (B.1)$$

where $\varphi^+$ is a phase angle (which differs from $\varphi$). One may then proceed by recasting Eq. (6) using the definition for the first Fourier components:

$$\mathrm{FDF}(A_p, \omega) = \frac{\bar{p}}{\bar{q}} \frac{1}{A_p e^{j\varphi^+}} \frac{\omega}{\pi} \int_0^{2\pi/\omega} q'(t) e^{-j\omega t} \, dt. \quad (B.2)$$

The nonlinear saturation of the heat release rate as a function of pressure may be expressed in the form:

$$q'(t) = \frac{\bar{q}}{\bar{p}} \frac{\beta}{\kappa} \arctan(\kappa p'(t)). \quad (B.3)$$

where $\beta$ is a constant of proportionality and $\kappa$ is a nonlinear saturation term. For this model the corresponding FDF is readily calculated to be:

$$\mathrm{FDF} = \frac{2\beta}{1 + \sqrt{1 + \kappa^2 A_p^2}}. \quad (B.4)$$

This FDF is monotonically decreasing and has the following limiting behaviour: $FDF \to 0$ as $A_p \to \infty$ and $FDF \to \beta$ as $A_p \to 0$. The heat release rate as a function of the pressure amplitude and the corresponding FDF are shown in green in Fig. B.14. Next we alter the slope of the saturation function at low and moderate amplitudes by introducing a Gaussian term:

$$q'_{A,B}(t) = \frac{\bar{q}}{\bar{p}} \left( \frac{\beta}{\kappa} \arctan(\kappa p'(t)) - \epsilon p'(t) \exp(-\mu p'(t)^2) \right). \quad (B.5)$$

Here $\epsilon$ and $\mu$ are constants that control the behaviour of the Gaussian, modifying the rate of saturation. The FDF for this function is computed by numerical integration. Depending on the magnitude of the parameter $\epsilon$, the FDF can change from being monotonic to being non-monotonic. In Fig. B.14 the monotonic behaviour is shown in blue, corresponding to $q'_A(t)$ and $FDF_A$. An increase in $\epsilon$ renders the FDF non-monotonic as shown in the red curves corresponding to $q'_B(t)$ and $FDF_B$. Observe that the change in slope of saturation of $q'(t)$ substantially alters the FDF at small and moderate amplitudes. To avoid numerical integration in the dynamical system we use the analytic approximation given by Eq. (8) (obtained by least squares fitting of the Gaussian integral), which is shown in dashed lines and is in good agreement with the one found numerically.

## References

[1] W. Krebs, P. Flohr, B. Prade, S. Hoffmann, Thermoacoustic stability chart for high-intensity gas turbine combustion systems, Combust. Sci. Technol. 174 (7) (2002) 99-128, doi:10.1080/00102200208984089.

[2] J. Bourgouin, D. Durox, J. Moeck, T. Schuller, S. Candel, A new pattern of instability observed in an annular combustor: the slanted mode, Proc. Combust. Inst. 35 (3) (2015) 3237-3244, doi:10.1016/j.proci.2014.06.029.

[3] T. Poinsot, Prediction and control of combustion instabilities in real engines, Proc. Combust. Inst. 36 (1) (2017) 1-28, doi:10.1016/j.proci.2016.05.007.

[4] A.P. Dowling, Y. Mahmoudi, Combustion noise, Proc. Combust. Inst. 35 (1) (2015) 65-100, doi:10.1016/j.proci.2014.08.016.

[5] J. O'Connor, V. Acharya, T. Lieuwen, Transverse combustion instabilities: acoustic, fluid mechanic, and flame processes, Prog. Energy Combust. Sci. 49 (2015) 1-39, doi:10.1016/j.pecs.2015.01.001.

[6] P. Wolf, G. Staffelbach, L.Y.M. Gicquel, J.-D. Müller, T. Poinsot, Acoustic and large eddy simulation studies of azimuthal modes in annular combustion chambers, Combust. Flame 159 (11) (2012) 3398-3413, doi:10.1016/j.combustflame.2012.06.016.

[7] N.A. Worth, J.R. Dawson, Modal dynamics of self-excited azimuthal instabilities in an annular combustion chamber, Combust. Flame 160 (11) (2013) 2476-2489, doi:10.1016/j.combustflame.2013.04.031.

[8] N.A. Worth, J.R. Dawson, Effect of equivalence ratio on the modal dynamics of azimuthal combustion instabilities, Proc. Combust. Inst. 36 (3) (2017) 3743-3751, doi:10.1016/j.proci.2016.06.115.

[9] N. Noiray, B. Schuermans, On the dynamic nature of azimuthal thermoacoustic modes in annular gas turbine combustion chambers, Proc. Royal Soc. A: Math. Phys. Eng. Sci. 469 (2151) (2013) 20120535, doi:10.1098/rspa.2012.0535.

[10] A. Faure-Beaulieu, T. Indlekofer, J.R. Dawson, N. Noiray, Experiments and low-order modelling of intermittent transitions between clockwise and anticlockwise spinning thermoacoustic modes in annular combustors, Proc. Combust. Inst. 38 (4) (2021) 5943-5951, doi:10.1016/j.proci.2020.05.008.
Matrix $M_2$ is not always invertible. Its determinant is $(A^3 n \omega^2 \cos(2\chi))/16$, which implies that the matrix cannot be inverted when $A=0$ or $|\chi|=\pi/4$. Therefore, we split the matrix into the product of an invertible matrix, say $B_3$, and a noninvertible matrix, say $M_3$, which contains all the terms that make the determinant zero or undefined, such that $M_2=B_3 M_3$. Multiplying from the left by $B_3^{-1}$ produces the following system of equations:

$$M_3(\mathbf{x})\dot{\mathbf{x}} = (B_3(\mathbf{x}))^{-1}\mathbf{b}_2(\mathbf{x}). \qquad (\text{A.12})$$

If we consider a meridian plane $\theta_0 = \text{constant}$ or the equatorial plane as in Fig. 10, the difficulties raised by the non-invertibility of matrix $M_3$ do not arise. Then, solving for $\dot{\mathbf{x}}$ produces the dynamical system shown in Eq. (4). For the case when $A=0$, Eq. (A.12) reduces to $\dot{A}=0$. For the limiting case $|\chi|=\pi/4$, $M_3$ is not invertible, then, after dividing Eq. (A.12) by A we obtain:

$$\begin{bmatrix} \frac{\dot{A}}{A} \\ \dot{\chi} \\ \cos(2\chi)\dot{\phi} \end{bmatrix} = \begin{bmatrix} \left(\frac{\omega_n^2}{2\omega} - \frac{\omega}{2}\right) \sin(2\chi) \\ 0 \\ \left(\frac{\omega_n^2}{2\omega} - \frac{\omega}{2}\right) \cos(2\chi) \end{bmatrix} + \sum_{m=1}^{N_b} \frac{F_r^{(m)}}{2} \begin{bmatrix} 1 + \cos(2\chi) \cos(2n(\theta_0 - \theta_m)) \\ -\cos(2\chi) \sin(2n(\theta_0 - \theta_m)) \\ -\sin(2\chi) \cos(2n(\theta_0 - \theta_m)) \\ \sin(2\chi) \sin(2n(\theta_0 - \theta_m)) \\ 0 \\ \sin(2\chi) \\ \sin(2n(\theta_0 - \theta_m)) \\ \cos(2\chi) + \cos(2n(\theta_0 - \theta_m)) \end{bmatrix} + \sum_{m=1}^{N_b} \frac{F_r^{(m)}}{2} \begin{bmatrix} 1 \\ 0 \\ \mp \cos(2n(\theta_0 - \theta_m)) \\ \mp \sin(2n(\theta_0 - \theta_m)) \\ 0 \\ \pm 1 \\ \sin(2n(\theta_0 - \theta_m)) \\ \cos(2n(\theta_0 - \theta_m)) \end{bmatrix}, \qquad (\text{A.13})$$

with $F_r$ and $F_j$ defined in Eq. (5). This equation is the analogue of Eq. (19) in Ghirardo and Gant [31] for a discrete number of burners. Substituting $\chi = \pm \pi/4$ gives:

$$\begin{bmatrix} \frac{\dot{A}}{A} \\ n\dot{\theta}_0 \pm \dot{\phi} \\ \dot{\chi} \\ 0 \end{bmatrix} = \begin{bmatrix} -\frac{\alpha}{2} \\ \pm\left(\frac{\omega_n^2}{2\omega} - \frac{\omega}{2}\right) \\ 0 \\ 0 \end{bmatrix} + \sum_{m=1}^{N_b} \frac{F_r^{(m)}}{2} \begin{bmatrix} 1 \\ 0 \\ \mp\cos(2n(\theta_0 - \theta_m)) \\ \pm\sin(2n(\theta_0 - \theta_m)) \end{bmatrix} + \sum_{m=1}^{N_b} \frac{F_j^{(m)}}{2} \begin{bmatrix} 0 \\ \pm 1 \\ \sin(2n(\theta_0 - \theta_m)) \\ \cos(2n(\theta_0 - \theta_m)) \end{bmatrix}, \qquad (\text{A.14})$$

The equation in the bottom row only holds in the case of a symmetric combustor. When the symmetry is broken there can only exist mixed modes in accordance to the findings in Noiray et al. [18], Faure-Beaulieu and Noiray [25] and the results shown in this study.

## Appendix B. Flame describing function

In this appendix two models for the saturation of the heat release rate as a function of pressure are examined. It is here convenient to write pressure fluctuations from Eq. (1) using the definition of the amplitude of the pressure fluctuations around the annulus Eq. (2):

$$p'(t, \theta) = A_p(\theta) \cos(\omega t + \varphi^+(\theta)). \qquad (\text{B.1})$$

where $\varphi^+$ is a phase angle (which differs from $\varphi$). One may then proceed by recasting Eq. (6) using the definition for the first Fourier components:

$$\mathrm{FDF}(A_p, \omega) = \frac{\bar{p}}{\bar{q}} \frac{1}{A_p e^{j\varphi^+}} \frac{\omega}{\pi} \int_0^{2\pi/\omega} q'(t) e^{-j\omega t} \, dt. \qquad (\text{B.2})$$

The nonlinear saturation of the heat release rate as a function of pressure may be expressed in the form:

$$q'(t) = \frac{\bar{q}}{\bar{p}} \frac{\beta}{\kappa} \arctan(\kappa p'(t)). \qquad (\text{B.3})$$

where $\beta$ is a constant of proportionality and $\kappa$ is a nonlinear saturation term. For this model the corresponding FDF is readily calculated to be:

$$\mathrm{FDF} = \frac{2\beta}{1 + \sqrt{1 + \kappa^2 A_p^2}}. \qquad (\text{B.4})$$

This FDF is monotonically decreasing and has the following limiting behaviour: $FDF \to 0$ as $A_p \to \infty$ and $FDF \to \beta$ as $A_p \to 0$. The heat release rate as a function of the pressure amplitude and the corresponding FDF are shown in green in Fig. B.14. Next we alter the slope of the saturation function at low and moderate amplitudes by introducing a Gaussian term:

$$q'_{A,B}(t) = \frac{\bar{q}}{\bar{p}} \left( \frac{\beta}{\kappa} \arctan(\kappa p'(t)) - \epsilon p'(t) \exp(-\mu p'(t)^2) \right). \qquad (\text{B.5})$$

Here $\epsilon$ and $\mu$ are constants that control the behaviour of the Gaussian, modifying the rate of saturation. The FDF for this function is computed by numerical integration. Depending on the magnitude of the parameter $\epsilon$, the FDF can change from being monotonic to being non-monotonic. In Fig. B.14 the monotonic behaviour is shown in blue, corresponding to $q'_A(t)$ and $FDF_A$. An increase in $\epsilon$ renders the FDF non-monotonic as shown in the red curves corresponding to $q'_B(t)$ and $FDF_B$. Observe that the change in slope of saturation of $q'(t)$ substantially alters the FDF at small and moderate amplitudes. To avoid numerical integration in the dynamical system we use the analytic approximation given by Eq. (8) (obtained by least squares fitting of the Gaussian integral), which is shown in dashed lines and is in good agreement with the one found numerically.

## References

[1] W. Krebs, P. Flohr, B. Prade, S. Hoffmann, Thermoacoustic stability chart for high-intensity gas turbine combustion systems, Combust. Sci. Technol. 174 (7) (2002) 99–128, doi:10.1080/00102200208984089.

[2] J. Bourgouin, D. Durox, J. Moeck, T. Schuller, S. Candel, A new pattern of instability observed in an annular combustor: the slanted mode, Proc. Combust. Inst. 35 (3) (2015) 3237–3244, doi:10.1016/j.proci.2014.06.029.

[3] T. Poinot, Prediction and control of combustion instabilities in real engines, Proc. Combust. Inst. 36 (1) (2017) 1–28, doi:10.1016/j.proci.2016.05.007.

[4] A.P. Dowling, Y. Mahmoudi, Combustion noise, Proc. Combust. Inst. 35 (1) (2015) 65–100, doi:10.1016/j.proci.2014.08.016.

[5] J. O'Connor, V. Acharya, T. Lieuwen, Transverse combustion instabilities: acoustic, fluid mechanic, and flame processes, Prog. Energy Combust. Sci. 49 (2015) 1–39, doi:10.1016/j.pecs.2015.01.001.

[6] P. Wolf, G. Staffelbach, L.Y.M. Gicquel, J.-D. Müller, T. Poinot, Acoustic and large eddy simulation studies of azimuthal modes in annular combustion chambers, Combust. Flame 159 (11) (2012) 3398–3413, doi:10.1016/j.combustflame.2012.06.016.

[7] N.A. Worth, J.R. Dawson, Modal dynamics of self-excited azimuthal instabilities in an annular combustion chamber, Combust. Flame 160 (11) (2013) 2476–2489, doi:10.1016/j.combustflame.2013.04.031.

[8] N.A. Worth, J.R. Dawson, Effect of equivalence ratio on the modal dynamics of azimuthal combustion instabilities, Proc. Combust. Inst. 36 (3) (2017) 3743–3751, doi:10.1016/j.proci.2016.06.115.

[9] N. Noiray, B. Schuermans, On the dynamic nature of azimuthal thermoacoustic modes in annular gas turbine combustion chambers, Proc. Royal Soc. A: Math. Phys. Eng. Sci. 469 (2151) (2013) 20120535, doi:10.1098/rspa.2012.0535.

[10] A. Faure-Beaulieu, T. Indlekofer, J.R. Dawson, N. Noiray, Experiments and low-order modelling of intermittent transitions between clockwise and anticlockwise spinning thermoacoustic modes in annular combustors, Proc. Combust. Inst. 38 (4) (2021) 5943–5951, doi:10.1016/j.proci.2020.05.008.
[11] S.R. Stow, A.P. Dowling, Thermoacoustic oscillations in an annular combustor, Volume 2: Coal, Biomass and Alternative Fuels; Combustion and Fuels; Oil and Gas Applications; Cycle Innovations, American Society of Mechanical Engineers, 2001, doi:10.1115/2001-gt-0037.

[12] S. Evesque, W. Polifke, Low-order acoustic modelling for annular combustors: validation and inclusion of modal coupling, Volume 1: Turbo Expo 2002, ASMEDC, 2002, doi:10.1115/gt2002-30064.

[13] M. Bauerheim, P. Salas, F. Nicoud, T. Poinsot, Symmetry breaking of azimuthal thermo-acoustic modes in annular cavities: a theoretical study, J. Fluid Mech. 760 (2014) 431–465, doi:10.1017/jfm.2014.578.

[14] A.P. Dowling, Nonlinear self-excited oscillations of a ducted flame, J. Fluid Mech. 346 (1997) 271–290, doi:10.1017/s0022112097006484.

[15] B. Schuermans, C. Paschereit, P. Monkewitz, Non-linear combustion instabilities in annular gas-turbine combustors, 44th AIAA Aerospace Sciences Meeting and Exhibit, American Institute of Aeronautics and Astronautics, 2006, doi:10.2514/6.2006-549.

[16] N. Noiray, D. Durox, T. Schuller, S. Candel, A unified framework for nonlinear combustion instability analysis based on the flame describing function, J. Fluid Mech. 615 (2008) 139–167, doi:10.1017/s0022112008003613.

[17] A.P. Dowling, A kinematic model of a ducted flame, J. Fluid Mech. 394 (1999) 51–72, doi:10.1017/s0022112099005686.

[18] N. Noiray, M. Bothien, B. Schuermans, Investigation of azimuthal staging concepts in annular gas turbines, Combust. Theor. Model. 15 (5) (2011) 585–606, doi:10.1080/13647830.2011.552636.

[19] J.-F. Bourgouin, D. Durox, J.P. Moeck, T. Schuller, S. Candel, Characterization and modeling of a spinning thermoacoustic instability in an annular combustor equipped with multiple matrix injectors, J. Eng. Gas Turbine Power 137 (2) (2014), doi:10.1115/1.4028257.

[20] G. Ghirardo, M.P. Juniper, Azimuthal instabilities in annular combustors: standing and spinning modes, Proc. Royal Soc. A: Math. Phys. Eng. Sci. 469 (2157) (2013) 20130232, doi:10.1098/rspa.2013.0232.

[21] G. Ghirardo, M.P. Juniper, J.P. Moeck, Weakly nonlinear analysis of thermoacoustic instabilities in annular combustors, J. Fluid Mech. 805 (2016) 52–87, doi:10.1017/jfm.2016.494.

[22] D. Laera, T. Schuller, K. Prieur, D. Durox, S.M. Camporeale, S. Candel, Flame describing function analysis of spinning and standing modes in an annular combustor and comparison with experiments, Combust. Flame 184 (2017) 136–152, doi:10.1016/j.combustflame.2017.05.021.

[23] J.P. Moeck, M. Paul, C.O. Paschereit, Thermoacoustic instabilities in an annular rijke tube, Volume 2: Combustion, Fuels and Emissions, Parts A and B, ASMEDC, 2010, doi:10.1115/gt2010-23577.

[24] G.A. Mensah, L. Magri, A. Orchini, J.P. Moeck, Effects of asymmetry on thermoacoustic modes in annular combustors: a higher-order perturbation study, J. Eng Gas Turbine Power 141 (4) (2018), doi:10.1115/1.4041007.

[25] A. Faure-Beaulieu, N. Noiray, Symmetry breaking of azimuthal waves: slow-flow dynamics on the bloch sphere, Phys. Rev. Fluids 5 (2) (2020), doi:10.1103/physrevfluids.5.023201.

[26] J.-F. Bourgouin, D. Durox, J.P. Moeck, T. Schuller, S. Candel, Self-sustained instabilities in an annular combustor coupled by azimuthal and longitudinal acoustic modes, Volume 1B: Combustion, Fuels and Emissions, American Society of Mechanical Engineers, 2013, doi:10.1115/gt2013-95010.

[27] J.-F. Bourgouin, D. Durox, T. Schuller, J. Beaunier, S. Candel, Ignition dynamics of an annular combustor equipped with multiple swirling injectors, Combust. Flame 160 (8) (2013) 1398–1413, doi:10.1016/j.combustflame.2013.02.014.

[28] K. Prieur, D. Durox, T. Schuller, S. Candel, A hysteresis phenomenon leading to spinning or standing azimuthal instabilities in an annular combustor, Combust. Flame 175 (2017) 283–291, doi:10.1016/j.combustflame.2016.05.021.

[29] A. Orchini, G.A. Mensah, J.P. Moeck, Effects of nonlinear modal interactions on the thermoacoustic stability of annular combustors, J. Eng. Gas Turbine Power 141 (2) (2018), doi:10.1115/1.4040768.

[30] G. Ghirardo, M.R. Bothien, Quaternion structure of azimuthal instabilities, Phys. Rev. Fluids 3 (11) (2018), doi:10.1103/physrevfluids.3.113202.

[31] G. Ghirardo, F. Gant, Averaging of thermoacoustic azimuthal instabilities, J. Sound Vib. (2020) 115732, doi:10.1016/j.jsv.2020.115732.

[32] J.P. Moeck, D. Durox, T. Schuller, S. Candel, Nonlinear thermoacoustic mode synchronization in annular combustors, Proc. Combust. Inst. 37 (4) (2019) 5343–5350, doi:10.1016/j.proci.2018.05.107.

[33] G. Ghirardo, F. Gant, F. Boudy, M. Bothien, Protection and identification of thermoacoustic azimuthal modes, Proceedings of ASME Turbo Expo 2020 Turbo machinery Technical Conference and Exposition (2020).

[34] R. Grimshaw, Nonlinear Ordinary Differential Equations, Blackwell Scientific Publications, 1990.

[35] P.E. Buschmann, G.A. Mensah, J.P. Moeck, Intrinsic thermoacoustic modes in an annular combustion chamber, Combust. Flame 214 (2020) 251–262, doi:10.1016/j.combustflame.2019.11.006.
[11] S.R. Stow, A.P. Dowling, Thermoacoustic oscillations in an annular combustor, Volume 2: Coal, Biomass and Alternative Fuels; Combustion and Fuels; Oil and Gas Applications; Cycle Innovations, American Society of Mechanical Engineers, 2001, doi:10.1115/2001-gt-0037.

[12] S. Evesque, W. Polifke, Low-order acoustic modelling for annular combustors: validation and inclusion of modal coupling, Volume 1: Turbo Expo 2002, ASMEDC, 2002, doi:10.1115/gt2002-30064.

[13] M. Bauerheim, P. Salas, F. Nicoud, T. Poinsot, Symmetry breaking of azimuthal thermo-acoustic modes in annular cavities: a theoretical study, J. Fluid Mech. 760 (2014) 431–465, doi:10.1017/jfm.2014.578.

[14] A.P. Dowling, Nonlinear self-excited oscillations of a ducted flame, J. Fluid Mech. 346 (1997) 271–290, doi:10.1017/s0022112097006484.

[15] B. Schuermans, C. Paschereit, P. Monkewitz, Non-linear combustion instabilities in annular gas-turbine combustors, 44th AIAA Aerospace Sciences Meeting and Exhibit, American Institute of Aeronautics and Astronautics, 2006, doi:10.2514/6.2006-549.

[16] N. Noiray, D. Durox, T. Schuller, S. Candel, A unified framework for nonlinear combustion instability analysis based on the flame describing function, J. Fluid Mech. 615 (2008) 139–167, doi:10.1017/s0022112008003613.

[17] A.P. Dowling, A kinematic model of a ducted flame, J. Fluid Mech. 394 (1999) 51–72, doi:10.1017/s0022112099005686.

[18] N. Noiray, M. Bothien, B. Schuermans, Investigation of azimuthal staging concepts in annular gas turbines, Combust. Theor. Model. 15 (5) (2011) 585–606, doi:10.1080/13647830.2011.552636.

[19] J.-F. Bourgouin, D. Durox, J.P. Moeck, T. Schuller, S. Candel, Characterization and modeling of a spinning thermoacoustic instability in an annular combustor equipped with multiple matrix injectors, J. Eng. Gas Turbine Power 137 (2) (2014), doi:10.1115/1.4028257.

[20] G. Ghirardo, M.P. Juniper, Azimuthal instabilities in annular combustors: standing and spinning modes, Proc. Royal Soc. A: Math. Phys. Eng. Sci. 469 (2157) (2013) 20130232, doi:10.1098/rspa.2013.0232.

[21] G. Ghirardo, M.P. Juniper, J.P. Moeck, Weakly nonlinear analysis of thermoacoustic instabilities in annular combustors, J. Fluid Mech. 805 (2016) 52–87, doi:10.1017/jfm.2016.494.

[22] D. Laera, T. Schuller, K. Prieur, D. Durox, S.M. Camporeale, S. Candel, Flame describing function analysis of spinning and standing modes in an annular combustor and comparison with experiments, Combust. Flame 184 (2017) 136–152, doi:10.1016/j.combustflame.2017.05.021.

[23] J.P. Moeck, M. Paul, C.O. Paschereit, Thermoacoustic instabilities in an annular rijke tube, Volume 2: Combustion, Fuels and Emissions, Parts A and B, ASMEDC, 2010, doi:10.1115/gt2010-23577.

[24] G.A. Mensah, L. Magri, A. Orchini, J.P. Moeck, Effects of asymmetry on thermoacoustic modes in annular combustors: a higher-order perturbation study, J. Eng Gas Turbine Power 141 (4) (2018), doi:10.1115/1.4041007.

[25] A. Faure-Beaulieu, N. Noiray, Symmetry breaking of azimuthal waves: slow-flow dynamics on the bloch sphere, Phys. Rev. Fluids 5 (2) (2020), doi:10.1103/physrevfluids.5.023201.

[26] J.-F. Bourgouin, D. Durox, J.P. Moeck, T. Schuller, S. Candel, Self-sustained instabilities in an annular combustor coupled by azimuthal and longitudinal acoustic modes, Volume 1B: Combustion, Fuels and Emissions, American Society of Mechanical Engineers, 2013, doi:10.1115/gt2013-95010.

[27] J.-F. Bourgouin, D. Durox, T. Schuller, J. Beaunier, S. Candel, Ignition dynamics of an annular combustor equipped with multiple swirling injectors, Combust. Flame 160 (8) (2013) 1398–1413, doi:10.1016/j.combustflame.2013.02.014.

[28] K. Prieur, D. Durox, T. Schuller, S. Candel, A hysteresis phenomenon leading to spinning or standing azimuthal instabilities in an annular combustor, Combust. Flame 175 (2017) 283–291, doi:10.1016/j.combustflame.2016.05.021.

[29] A. Orchini, G.A. Mensah, J.P. Moeck, Effects of nonlinear modal interactions on the thermoacoustic stability of annular combustors, J. Eng. Gas Turbine Power 141 (2) (2018), doi:10.1115/1.4040768.

[30] G. Ghirardo, M.R. Bothien, Quaternion structure of azimuthal instabilities, Phys. Rev. Fluids 3 (11) (2018), doi:10.1103/physrevfluids.3.113202.

[31] G. Ghirardo, F. Gant, Averaging of thermoacoustic azimuthal instabilities, J. Sound Vib. (2020) 115732, doi:10.1016/j.jsv.2020.115732.

[32] J.P. Moeck, D. Durox, T. Schuller, S. Candel, Nonlinear thermoacoustic mode synchronization in annular combustors, Proc. Combust. Inst. 37 (4) (2019) 5343–5350, doi:10.1016/j.proci.2018.05.107.

[33] G. Ghirardo, F. Gant, F. Boudy, M. Bothien, Protection and identification of thermoacoustic azimuthal modes, Proceedings of ASME Turbo Expo 2020 Turbo machinery Technical Conference and Exposition (2020).

[34] R. Grimshaw, Nonlinear Ordinary Differential Equations, Blackwell Scientific Publications, 1990.

[35] P.E. Buschmann, G.A. Mensah, J.P. Moeck, Intrinsic thermoacoustic modes in an annular combustion chamber, Combust. Flame 214 (2020) 251–262, doi:10.1016/j.combustflame.2019.11.006.