2.1 Problem geometry

We adopt cylindrical coordinates $z,r,\unicode[STIX]{x1D703}$ , with the $z$ -axis corresponding to the axis of the combustion chamber and $\unicode[STIX]{x1D703}$ in $[0\,,\,2\unicode[STIX]{x03C0})$ . We assume that a number $N_{b}$ of equal burners are equispaced along the annulus and that each of the $N_{b}$ sectors has the same geometry. The problem is then invariant to the group of $N_{b}$ -fold rotational symmetry $C_{N_{b}}$ , with the fundamental domain being a sector spanning the angle $\unicode[STIX]{x0394}\unicode[STIX]{x1D703}\equiv 2\unicode[STIX]{x03C0}/N_{b}$ .

We assume that the flames are acoustically compact, so that the fluctuating heat release rate is concentrated at the locations of the burners:

where $\unicode[STIX]{x1D6FF}$ is the Dirac delta in three dimensions, $\overline{r}$ is the radial position of the burners and the plane $z=0$ is at the interface between the combustion chamber and the burners, which are located at the azimuthal positions $\unicode[STIX]{x1D703}_{j}$ , equispaced by $\unicode[STIX]{x0394}\unicode[STIX]{x1D703}$ . We number the burners in anticlockwise direction, and we choose a frame of reference so that the first burner is positioned at

The addition of the coefficient $\unicode[STIX]{x1D701}$ is arbitrary and corresponds to a simple rotation of the frame of reference, which will be useful later. This is exemplified in figure 1 for an odd $N_{b}=5$ number of burners.

Figure 1. Example of the position of the frame of reference for different values of $\unicode[STIX]{x1D701}$ for a number $N_{b}=5$ of burners. The burners are represented with black discs, and the large circles are the internal and external walls of the combustion chamber. In general the position at the angle $\unicode[STIX]{x1D703}=\unicode[STIX]{x03C0}/4$ : $\bullet$ for $\unicode[STIX]{x1D701}=0$ is occupied by a burner; $\bullet$ for $\unicode[STIX]{x1D701}=2$ is equispaced between two adjacent burners; $\bullet$ for $\unicode[STIX]{x1D701}=1$ is $3\unicode[STIX]{x0394}\unicode[STIX]{x1D703}/4$ far from the preceding burner and $\unicode[STIX]{x0394}\unicode[STIX]{x1D703}/4$ far from the next.

Some designs feature a mean velocity $U_{\unicode[STIX]{x1D703}}$ in the azimuthal direction along the combustion chamber annulus. A non-zero azimuthal mean flow $U_{\unicode[STIX]{x1D703}}$ has two effects. Firstly, a non-zero $U_{\unicode[STIX]{x1D703}}$ makes one of the two spinning modes rotate faster and the other slower, and makes standing modes slowly rotate with pressure and velocity nodes moving at the speed of the mean azimuthal flow. See for example Wolf et al. (Reference Wolf, Staffelbach, Balakrishnan, Roux and Poinsot2010) for numerical evidence and a discussion, and refer to Bauerheim et al. (Reference Bauerheim, Salas, Nicoud and Poinsot2014), Bauerheim, Cazalens & Poinsot (Reference Bauerheim, Cazalens and Poinsot2015) for a detailed analysis of this first effect of $U_{\unicode[STIX]{x1D703}}$ in a linear framework. Secondly, a non-zero $U_{\unicode[STIX]{x1D703}}$ bends the flames in the azimuthal direction, orthogonally to the burner’s axis. This leads to a loss of axisymmetry of the mean flame shape, and this loss of axisymmetry is in turn a necessary condition for the flame to have a non-zero linear response to azimuthal velocity perturbations, as proven by Acharya & Lieuwen (Reference Acharya and Lieuwen2014). In most cases however $U_{\unicode[STIX]{x1D703}}$ is very small compared to the speed of sound, and is fixed to zero in the following. This introduces $N_{b}$ reflection planes that are parallel to the combustor axis and pass through the middle of one segment, so that the problem is invariant to the pyramidal group of symmetry  $C_{N_{b}v}$ .

2.2 Flame response

In this study, we assume that the flame response to the acoustic field is known, both in the linear and nonlinear regime. A common modelling approach consists of expressing the fluctuating heat release rate $q_{j}$ of the $j$ th flame in terms of only the acoustic axial velocity $v_{j}$ just upstream of the burner. Doing so, we assume that the azimuthal, acoustic velocity $u$ does not affect the response. This last point is proved theoretically in the linear limit for axisymmetric premixed flames in Acharya et al. (Reference Acharya, Shreekrishna and Lieuwen2012). This influence is experimentally verified to be small at low amplitudes of transverse forcing for the cases of a burner positioned at pressure/velocity nodes and for the case where it is swept by a spinning wave, where both $u$ and $v$ are excited at the same time (Saurabh et al. Reference Saurabh, Steinert, Moeck and Paschereit2014). This effect is usually not taken into account because little is known in the nonlinear regime, i.e. at amplitudes of oscillation typical of self-excited thermoacoustic oscillations. In this paper we make the same assumption, but point out that the nonlinear effect of the transverse azimuthal velocity $u$ on each flame has been investigated by Ghirardo & Juniper (Reference Ghirardo and Juniper2013). It does not affect the linear stability properties of the system, but it does affect the nonlinear dynamics, and can explain stable standing solutions in axisymmetric annular chambers.

The longitudinal fluctuating velocity $v_{j}$ oscillating in the $j$ th burner can be expressed as a linear time-invariant operator of the acoustic pressure difference $\unicode[STIX]{x0394}p_{j}$ between the two sides of the burner, which are the chamber and the plenum. The burners are assumed to be acoustically compact (Blimbaum et al. Reference Blimbaum, Zanchetta, Akin, Acharya, O’Connor, Noble and Lieuwen2012), which allows them to be modelled as lumped elements. However, if we consider one mode oscillating harmonically in time, and we assume that the burner transfer function of the lumped element does not depend on the amplitude of oscillation (as validated for example in Ćosić, Moeck & Paschereit (Reference Ćosić, Moeck and Paschereit2014)), then $\unicode[STIX]{x0394}p_{j}$ is proportional to $p_{j}$ , and one can express $v_{j}$ as a linear operator of the local value of the pressure in the chamber $p_{j}$ . The same reasoning applies also to two degenerate azimuthal modes oscillating at the same frequency, as will be the case in the following.

In particular, we model the fluctuating heat release rate as a time-invariant operator  ${\mathcal{Q}}$ :

The operator ${\mathcal{Q}}$ contains all the complexity of the relation between $p_{j}$ and $q_{j}$ , and is nonlinear. We restrict our study to operators ${\mathcal{Q}}$ that, excited with a harmonic input $p=A\cos (\unicode[STIX]{x1D714}t)$ , respond strongly at the same input frequency $\unicode[STIX]{x1D714}$ and weakly at higher harmonics $2\unicode[STIX]{x1D714},3\unicode[STIX]{x1D714},\ldots$  . This is a feature of flames, acting like a low-pass filter on the acoustic input (Schuller, Durox & Candel Reference Schuller, Durox and Candel2003). This, together with the acoustics being a narrow-band filter, is one of the reasons why frequency-domain calculations truncated at the first harmonic have proven successful in thermoacoustics even for limit-cycle calculations. We also assume that ${\mathcal{Q}}$ is a stable operator, i.e. the fluctuations of the heat release rate are present only if an external acoustic wave perturbs the flame. This means for example that we do not consider the flame response to its own acoustic field (Assier & Wu Reference Assier and Wu2014) if it leads to a linearly unstable flame.

We will study the problem both in the time and frequency domains. We refer with the calligraphic symbol ${\mathcal{Q}}$ to the time-domain operator mapping pressure perturbations to fluctuations in the heat release rate. In the frequency domain, we can calculate the corresponding describing function, which we label with the uppercase $Q$ , defined as (Gelb & Vander Velde Reference Gelb and Vander Velde1968):

The real and imaginary parts of $Q(A,\unicode[STIX]{x1D714})$ express the amplitudes of the components of the fluctuating heat release rate, i.e. the output of the operator, respectively, in phase and in quadrature with the sinusoidal pressure input. In particular it is $\text{Re}[Q(A,\unicode[STIX]{x1D714})]$ that leads to a contribution to the energy production term $q(t)p(t)$ in the Rayleigh criterion: if positive, the fluctuating heat release rate is partially in phase with respect to the pressure input and the energy production term ${\mathcal{Q}}[p(t)]p(t)$ in the Rayleigh criterion is positive over a limit cycle. One can then define the gain $G$ and the phase lag $\unicode[STIX]{x1D719}$ of the flame response as:

Notice that for a model with a constant time delay $\unicode[STIX]{x1D70F}$ between the pressure and the fluctuating heat release rate we have $\unicode[STIX]{x1D719}(A,\unicode[STIX]{x1D714})=+\unicode[STIX]{x1D714}\unicode[STIX]{x1D70F}$ . The sign convention of $+\text{i}\unicode[STIX]{x1D714}t$ in the exponential in (2.4) is historical, and we point out that it is the opposite of the Fourier transform that we will use later.

The response of the flame is always bounded, i.e. the gain is always between $0$ and $G_{max}$ . We also assume that the describing function is a continuous function of the amplitude $A$ and of the frequency $\unicode[STIX]{x1D714}$ . This is usually an observed property of the experimental data (see for example Palies et al. (Reference Palies, Durox, Schuller and Candel2011)), although it is possible that the flame will abruptly extinguish above a certain amplitude of forcing, typically because of blow-off or flash-back events.

We also observe that the level of acoustic damping is typically constant or decreases with the amplitude of oscillation (Ćosić, Reichel & Paschereit Reference Ćosić, Reichel and Paschereit2012). This means that the system arrives at a limit cycle because the flame response decreases with amplitude, not because the damping increases with amplitude. Since for convenience we prefer to not set a lower bound for the level of acoustic damping, we characterize the existence of an amplitude at which the energy balance is obtained by assuming that $\lim _{A\rightarrow \infty }|Q(A,\unicode[STIX]{x1D714})|=0$ .

2.3 Governing equations

Making a zero Mach number assumption, the inhomogeneous wave equation governing the problem is, as derived for example by Nicoud et al. (Reference Nicoud, Benoit, Sensiau and Poinsot2007):

In the equation, subscript 0 refers to mean quantities, which depend on $\boldsymbol{x}$ only, while subscript 1 refers to fluctuating quantities, which depend on $\boldsymbol{x}$ and $t$ . In this paper we assume that the density $\unicode[STIX]{x1D70C}_{0}$ and the isentropic coefficient $\unicode[STIX]{x1D6FE}$ are uniform. This hypothesis can possibly be lifted, but the equations become complicated without adding more insight. We also non-dimensionalize the equations, with respect to a spatial length scale $D$ (say the radius of the annular chamber) and the acoustic time scale $D/c$ , with $c$ being the mean speed of sound in the chamber. We assume an ideal gas, so that $\unicode[STIX]{x1D70C}_{0}c^{2}=\unicode[STIX]{x1D6FE}p_{0}$ . We introduce the non-dimensional fluctuating pressure $p$ and fluctuating heat release rate $q$ as

In the new non-dimensional coordinates, equation (2.7) simplifies to

where we substituted the expression (2.1) and (2.3) for the fluctuating heat release rate $q$ . We adopt the following convention for the definition of the Fourier transform:

By multiplying all terms of (2.9) by $\text{e}^{-\text{i}\unicode[STIX]{x1D714}t}/\unicode[STIX]{x03C0}$ and integrating over the time $t$ we obtain the inhomogeneous Helmholtz equation:

Each of the integrals in the sum on the right-hand side can be rewritten as

Notice that we assume that the response at the frequency $\unicode[STIX]{x1D714}$ of ${\mathcal{Q}}[p(\boldsymbol{x},t)]$ only depends on the amplitude $|\hat{p}(\boldsymbol{x})|$ of the solution at the same frequency $\unicode[STIX]{x1D714}$ . This is correct as long as all other harmonics are negligible, i.e. the filtering hypothesis holds (Gelb & Vander Velde Reference Gelb and Vander Velde1968). We also point out that in the last passage of (2.12) the complex conjugation denoted by the asterisk appears because of the different sign convention in the exponential in the definitions (2.4) and (2.10). Substituting (2.12) in (2.11) we obtain:

Equation (2.13) must be accompanied by suitable boundary conditions. At the combustor walls these will be zero normal gradient conditions for the pressure. At the axial extremes of the domain, the combustor inlet and outlet, the boundary conditions will in general be of impedance type, $\hat{p}=Z(\unicode[STIX]{x1D714})\hat{u}$ , with $Z(\unicode[STIX]{x1D714})$ a complex-valued function.

2.4 Eigenmodes’ degeneracy

We linearise equation (2.13) with respect to the amplitude $|\hat{p}(\boldsymbol{x})|$ of the solution and obtain:

where $L(\unicode[STIX]{x1D714})$ is the transfer function of ${\mathcal{Q}}[p(t)]$ at infinitesimal amplitudes. The set of solutions of the eigenvalue problem (2.14) is $\{(\hat{\unicode[STIX]{x1D713}}_{n}(\boldsymbol{x}),\unicode[STIX]{x1D70E}_{n}+\text{i}\hat{\unicode[STIX]{x1D714}}_{n})$ with $\,\,\unicode[STIX]{x1D70E}_{n},\hat{\unicode[STIX]{x1D714}}_{n}\in \mathbb{R}\,\,,n=1,2,\ldots \}$ where $\hat{\unicode[STIX]{x1D713}}_{n}(\boldsymbol{x})$ is the complex-valued eigenvector describing the shape of the mode and $\unicode[STIX]{x1D70E}_{n}+\text{i}\hat{\unicode[STIX]{x1D714}}_{n}$ is the corresponding eigenvalue. The modes and their eigenvalues can be calculated using a Helmholtz solver (Nicoud et al. Reference Nicoud, Benoit, Sensiau and Poinsot2007) or a thermoacoustic network model of the problem (Stow & Dowling Reference Stow and Dowling2001; Schuermans, Bellucci & Paschereit Reference Schuermans, Bellucci and Paschereit2003).

We are particularly interested in azimuthal modes, i.e. solutions that are $n$ -periodic in $\unicode[STIX]{x1D703}$ in $[0\,,\,2\unicode[STIX]{x03C0}]$ , with $n$ called the azimuthal wavenumber of the mode. As discussed by Moeck, Paul & Paschereit (Reference Moeck, Paul and Paschereit2010), Bauerheim et al. (Reference Bauerheim, Salas, Nicoud and Poinsot2014), an azimuthal mode of wavenumber $n$ belongs to an eigenspace of dimension two because of the rotational symmetry of the problem. There are however exceptions, when $n$ is a multiple of $N_{b}/2$ in the case of an even number of burners $N_{b}$ , and when $n$ is a multiple of $N_{b}$ in the case of $N_{b}$ odd. We refer the reader also to § 5.4 of Salas (Reference Salas2013) for a concise summary of these two cases. These exceptions are non-degenerate cases, i.e. their modes belong to an eigenspace of dimension one, and occur because the rotational symmetry is discrete. We focus the analysis on the degenerate case where the dimensionality is two because: (i)  $N_{b}$ is usually large and the excited modes typically have a low azimuthal wavenumber $n$ (up to $n=4$ in Seume et al. (Reference Seume, Vortmeyer, Krause, Hermann, Hantschk, Gleis, Vortmeyer and Orthmann1998)); (ii) the non-degenerate case does not give rise to the rich dynamics that can be observed in the degenerate case. We study thermoacoustic oscillations of azimuthal modes with $n=1$ in the following, but the same analysis can be generalized to higher $n$ , as long as the case stays degenerate.

We assume that these modes are close to their Hopf bifurcation. In other words, we assume that all other modes are stable, and only azimuthal modes of wavenumber $n=1$ take part in the oscillation. The Floquet–Bloch theorem (Chap. VIII pp. 139–140 Brillouin Reference Brillouin1953; Mensah & Moeck Reference Mensah and Moeck2015) ensures that one of the $n=1$ solutions can generally be written in the form $\widetilde{\unicode[STIX]{x1D712}}(z,r)\text{e}^{\text{i}\unicode[STIX]{x1D703}}$ , where $\widetilde{\unicode[STIX]{x1D712}}(z,r)$ is periodic in $\unicode[STIX]{x1D703}$ with period $2\unicode[STIX]{x03C0}/N_{b}$ , i.e. one burner segment. The dependence of $\unicode[STIX]{x1D712}$ on $\unicode[STIX]{x1D703}$ can be in principle be taken into account. Since it is of secondary importance when compared to the long-wave component $\text{e}^{\text{i}\unicode[STIX]{x1D703}}$ , it is neglected in the following in favour of a clearer exposition. Because of the reflection symmetry of the problem, there exists a second solution of the eigenspace that is symmetric with respect to the first, with structure $\widetilde{\unicode[STIX]{x1D712}}(z,r)\text{e}^{-\text{i}\unicode[STIX]{x1D703}}$ . We refer to these two modes in the following as spinning modes, because their phases linearly increase or decrease in the azimuthal direction.

By linearly combining these two spinning modes we can obtain two solutions $\unicode[STIX]{x1D713}_{1}$ and $\unicode[STIX]{x1D713}_{2}$ that have a constant phase as a function of $\unicode[STIX]{x1D703}$ :

We refer to these modes in the following as standing modes, because if observed only at the burners’ location they have pressure nodes and pressure antinodes fixed in time and in space. By direct substitution one can prove that they are orthogonal:

where $\unicode[STIX]{x1D6FA}$ is the volume of the combustion chamber. One proves by direct substitution also that

We normalise the standing modes by fixing the value of $\widetilde{\unicode[STIX]{x1D712}}(0,\overline{r})$ to 1 at the burners’ positions $(z,r,\unicode[STIX]{x1D703})=(0,\overline{r},\unicode[STIX]{x1D703}_{j})$ , so that

2.5 Weakly nonlinear analysis

We study the solution of the nonlinear problem as a perturbation of its linear, neutrally stable counterpart. We obtain the latter by changing the problem (2.14), by substituting $\widetilde{\unicode[STIX]{x1D709}}L(\unicode[STIX]{x1D714})$ in place of $L(\unicode[STIX]{x1D714})$ , with $\widetilde{\unicode[STIX]{x1D709}}$ a real, non-negative coefficient, so that for $\widetilde{\unicode[STIX]{x1D709}}=1$ we recover the original equations. We then look for the value $\unicode[STIX]{x1D709}$ of the coefficient $\widetilde{\unicode[STIX]{x1D709}}$ such that the growth rate of the first two modes $\unicode[STIX]{x1D70E}_{1,2}(\widetilde{\unicode[STIX]{x1D709}})$ is zero. In other words, we linearly increase/decrease the gain of the flame response to the level that makes the system neutrally stable. Notice that by looking at (2.14), one may guess that this may happen only for $\widetilde{\unicode[STIX]{x1D709}}=0$ . This is not generally the case, due to the presence of partially transmitting boundary conditions or sources of damping, such as Helmholtz resonators and/or acoustic liners that remove fluctuation energy from the system. In this study, we consider only linear damping. Nonlinear acoustic damping effects at the boundaries can be characterised with a describing function (Schuller et al. Reference Schuller, Tran, Noiray, Durox, Ducruix and Candel2009) in the frequency domain, and its time-domain realization (Ghirardo et al. Reference Ghirardo, Ćosić, Juniper and Moeck2015) in the time domain. We refer to quantities evaluated for $\widetilde{\unicode[STIX]{x1D709}}=\unicode[STIX]{x1D709}$ by dropping the tilde, so that the eigenmodes are indicated with $\unicode[STIX]{x1D713}_{1}(\boldsymbol{x})$ and $\unicode[STIX]{x1D713}_{2}(\boldsymbol{x})$ , and their real-valued eigenfrequency is $\unicode[STIX]{x1D714}_{1}=\unicode[STIX]{x1D714}_{2}$ .

For later use, we write (2.14) for the first two modes to obtain

We can multiply all terms of (2.18) by $\unicode[STIX]{x1D713}_{k}^{\ast }$ and integrate over the domain $\unicode[STIX]{x1D6FA}$ :

where, by exploiting the fact that $\sum _{j=1}^{N_{b}}c_{j}^{2}=\sum _{j=1}^{N_{b}}s_{j}^{2}=N_{b}/2$ , we introduced the quantities

and $\unicode[STIX]{x1D707}=\unicode[STIX]{x1D707}_{1}=\unicode[STIX]{x1D707}_{2}$ is a real valued normalisation factor that is the same for the two modes, defined as:

In both right-hand sides of (2.20) the frequency $\unicode[STIX]{x1D714}_{k}=\unicode[STIX]{x1D714}_{1}=\unicode[STIX]{x1D714}_{2}$ is much larger than the other terms, so that $\unicode[STIX]{x1D714}_{0}\approx \unicode[STIX]{x1D714}_{k}$ and $\unicode[STIX]{x1D6FC}\ll |\unicode[STIX]{x1D714}_{k}|$ . This follows from the weakly nonlinear nature of thermoacoustic problems. Equations (2.20) define the equivalent acoustic damping $\unicode[STIX]{x1D6FC}$ and natural frequency $\unicode[STIX]{x1D714}_{0}$ of the system when the flame response is uniformly reduced in the annulus to the point of making the system neutrally stable, i.e. at $\widetilde{\unicode[STIX]{x1D709}}=\unicode[STIX]{x1D709}$ . This paragraph led to (2.19), which will be used in the following.

At the value $\widetilde{\unicode[STIX]{x1D709}}=\unicode[STIX]{x1D709}$ no dissipation/gain of energy in a limit cycle occurs in the linearised system for the dominant mode, and the exact solution of the problem is

where $X_{1}$ and $X_{2}$ are two complex-valued constants and the decaying terms depend on the initial condition and they will be neglected in the following because they converge to zero in time after the initial transient. We now study the original nonlinear problem (2.9) as a perturbation of this, as $\widetilde{\unicode[STIX]{x1D709}}$ changes from $\unicode[STIX]{x1D709}$ to 1. We choose as ansatz in the frequency domain

In (2.23), $\unicode[STIX]{x1D700}=1-\unicode[STIX]{x1D709}$ is the perturbation parameter and expresses the deviation of the exact nonlinear solution from the first term at the right-hand side, which is the solution of the linear problem, with coefficients $\hat{\unicode[STIX]{x1D702}}_{k}$ that will be calculated next. This deviation occurs because of the onset of higher harmonics in time and in space, and because of the structural change in the equations that can affect slightly the shape of the modes. Using (2.17), the expression for the pressure field at the burners’ location reads

We then substitute this ansatz into (2.13), multiply all terms by $-\unicode[STIX]{x1D713}_{1}(\boldsymbol{x})^{\ast }/(\text{i}\unicode[STIX]{x1D714})$ , and integrate over the domain $\unicode[STIX]{x1D6FA}$ . We first exploit the orthogonality properties (2.16) between the two degenerate modes, and then substitute (2.19) and obtain:

Notice how on the left-hand side, the frequency $\unicode[STIX]{x1D714}_{1}$ is the frequency of the mode $\unicode[STIX]{x1D713}_{1}$ , so that in principle $\text{i}\unicode[STIX]{x1D714}_{1}\hat{\unicode[STIX]{x1D702}}_{1}$ is not the Fourier transform of $\unicode[STIX]{x2202}\unicode[STIX]{x1D702}_{1}/\unicode[STIX]{x2202}t$ . However, the frequency of the nonlinear system $\unicode[STIX]{x1D714}$ is close to $\unicode[STIX]{x1D714}_{1}$ , and we can make this approximation. We take the inverse Fourier transform of all terms and obtain

where the second equation was obtained analogously by symmetry, the coefficients $c_{j}$ and $s_{j}$ were introduced in (2.17) and $\unicode[STIX]{x1D707}$ was defined in (2.21). The equations (2.26) are investigated in the rest of the paper and describe the temporal evolution of the amplitudes $\unicode[STIX]{x1D702}_{1}(t)$ and $\unicode[STIX]{x1D702}_{2}(t)$ of the two standing modes $\unicode[STIX]{x1D713}_{1}$ and $\unicode[STIX]{x1D713}_{2}$ . The left-hand side of equation (2.26a ) describes a damped oscillator with natural frequency $\unicode[STIX]{x1D714}_{0}$ , defined in (2.20a ), and the damping $\unicode[STIX]{x1D6FC}$ , defined in (2.20b ). The right-hand side of equation (2.26a ) is the ensemble response of the flames. The argument of the operator ${\mathcal{Q}}$ is the local value of the pressure at the $j$ th burner truncated at zero order, as can be observed in (2.24). Notice how the local contribution of the $j$ th burner in the sum on the right-hand side is weighted by the term $\unicode[STIX]{x1D707}c_{j}$ , which is a measure of how large the mode is at that burner.

2.6 The final oscillator model

We now fix $\unicode[STIX]{x1D707}=1$ in (2.26) without loss of generality, because the same effect could be obtained by rescaling ${\mathcal{Q}}$ . In the following we denote $\unicode[STIX]{x1D714}_{0}$ simply as $\unicode[STIX]{x1D714}$ because we remain in the time domain. By neglecting the correction in $O(\unicode[STIX]{x1D700})$ and by denoting the time derivative with a prime, the system (2.26) is:

where the functions $f_{k}$ on the right-hand side are: This refactoring of the equations allows a more concise discussion in § 3.

This section rigorously derived a set of governing equations for the instantaneous amplitude of oscillations of two standing modes. These equations are consistent with existing low-order models, which are based instead on the projection of the equations on a Galerkin basis $\{\cos \unicode[STIX]{x1D703},\sin \unicode[STIX]{x1D703}\}$ which was chosen because it matches with experiments. A great advantage of this approach is that we can calculate the coefficients with a Helmholtz solver or a network model so that the model can be used as a predictive tool.

3.1 Temporal averaging

In this section we apply the method of time averaging to the equations (2.33) by assuming that the terms $f_{k}(\unicode[STIX]{x1D702}_{1}^{\prime },\unicode[STIX]{x1D702}_{2}^{\prime })-\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D702}_{k}^{\prime }$ are small. From a physical standpoint, this means that we are assuming that the dynamics of the acoustic waves (the left-hand side in (2.33)) is only slightly influenced by the net contribution of flame response and acoustic damping (the right-hand side in (2.33)). This can be observed experimentally in the time traces of pressure signals: usually thermoacoustic systems take many periods of oscillation to stabilize to a periodic solution, meaning that the amount of net energy contributed to the acoustics in one oscillation cycle is small.

We express each oscillator as a harmonic oscillation (Sanders & Verhulst Reference Sanders and Verhulst2007):

with $A_{k}(t),\unicode[STIX]{x1D711}_{k}(t)$ the envelope and the instantaneous phase of the $k$ th oscillator, slowly varying with respect to the fast time variable $t$ , usually referred to as slow variables of the problem.

The expression of the pressure field at the burners’ position can be rewritten in terms of $A_{k},\unicode[STIX]{x1D711}_{k}$ by substituting (3.1) and (2.17) into (2.24) and neglecting terms of order $O(\unicode[STIX]{x1D700})$ :

Applying the method of temporal averaging for a delayed system (Wahi & Chatterjee Reference Wahi and Chatterjee2004) we obtain the equations governing the temporal evolution of the slow variables:

where $\unicode[STIX]{x1D711}\equiv \unicode[STIX]{x1D711}_{1}-\unicode[STIX]{x1D711}_{2}$ is the phase difference between the two oscillators and the angled brackets denote averaging over a limit cycle:

The period of averaging $2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D714}$ is a constant, where $\unicode[STIX]{x1D714}$ is the frequency of oscillation appearing in (2.26). Keeping $\unicode[STIX]{x1D714}$ constant has a slight effect on the accuracy of the method and will be discussed later.

We discuss in § A of the supplementary material available at http://dx.doi.org/ 10.1017/jfm.2016.494 how to find an analytical solution of the terms (3.4). We present here only the key physical features of the solution. In the integral (3.4) the function $f_{k}$ consists of the sum of the contributions of each burner. Each burner responds to the local value of the pressure field as described by (3.2), which depends on the two modes $A_{1},A_{2}$ . However both modes $A_{1}$ and $A_{2}$ oscillate at the same frequency $\unicode[STIX]{x1D714}$ , so that each burner responds to a single harmonic input. The averaged response of one burner to a harmonic signal is by definition (with some minor multiplicative coefficient adjustments) the describing function of the flame. Then the integrals (3.4) can be rewritten in terms of the gain $G$ and of the phase lag $\unicode[STIX]{x1D719}$ of a single flame.

By substituting the solution of (3.4) into (3.3), we obtain the slow flow equations:

where $R_{j}$ in (3.5) is defined as:

and the terms $A_{1}c_{j}$ and $A_{2}s_{j}$ are the amplitudes of the pressure of the two modes $\unicode[STIX]{x1D702}_{1},\unicode[STIX]{x1D702}_{2}$ at the $j$ th burner. The amplitude $R_{j}$ is then the slowly varying amplitude of oscillation of the pressure at the $j$ th burner. Notice that both the gain $G$ and the phase lag $\unicode[STIX]{x1D719}$ of each burner depend on the local amplitude of oscillation $R_{j}$ , and that $R_{j}$ plays a role in the equations only as their forcing amplitude. The spatial structure of $R_{j}$ for standing and spinning waves is examined in § 3.2.

The effect of the acoustic damping coefficient $\unicode[STIX]{x1D6FC}$ (first term on the right-hand side of (3.5a ), (3.5b ) is to push individually the two standing modes to smaller amplitudes, and it opposes the effect of the flame response (second term on the right-hand side of (3.5a ), (3.5b ). In order to reach a limit cycle a balance between the two terms must be reached in both equations.

The synchronization of the two modes is described by (3.5c ), and determines the standing/spinning nature of the solutions. The synchronisation does not depend on the linear properties of the system, but only on its nonlinear saturation features. In fact, (3.5c ) depends only indirectly on the amplitudes through the dependence on $R_{j}$ if we fix a certain ratio $A_{1}/A_{2}$ .

We do not discuss the linear stability of the fixed point at the origin $(A_{1},A_{2},\unicode[STIX]{x1D711})=(0,0,\unicode[STIX]{x1D711})$ , because it leads to the same results discussed earlier in § 2.5.1.

In summary, we applied the method of averaging to the dynamic equations of the two oscillators (2.33), which were in terms of the four variables $\{\unicode[STIX]{x1D702}_{1},\unicode[STIX]{x1D702}_{1}^{\prime },\unicode[STIX]{x1D702}_{2},\unicode[STIX]{x1D702}_{2}^{\prime }\}$ . The original equations may exhibit limit-cycle solutions, oscillating at a fast, acoustic time scale with frequency $\unicode[STIX]{x1D714}$ . The resulting equations (3.5) can be rewritten as

They describe the dynamic evolution of three variables, which are oblivious of the fast acoustic time scale: the two amplitudes of the standing modes $A_{1}$ and $A_{2}$ and their phase difference $\unicode[STIX]{x1D711}$ . The limit-cycle solutions of the oscillators are fixed points $(A_{1},A_{2},\unicode[STIX]{x1D711})$ of the new set of equations. The time scale of this process depends in the linear regime on the relative strengths of the linear flame response and the acoustic damping, and in the nonlinear regime on the nonlinear saturation of the gain $G$ and of the phase $\unicode[STIX]{x1D719}$ of the flames.

3.2 Standing and spinning waves

In this section we discuss the structure of standing and spinning waves. Waves are considered as possible initial conditions of the problem at a certain instant of time, and the system can drift away from this initial state as time evolves. This differs from standing and spinning solutions, which are waves that are also periodic solutions of the problem.

Some results presented here are well known in the literature, and are presented only for reference. In particular we prove in this section that a point in the state space of the averaged system with coordinates $(A_{1},A_{2},\unicode[STIX]{x1D711})=(A,A,k\unicode[STIX]{x03C0}/2)$ is always a standing or a spinning wave:

This follows from the structure of the pressure field (3.2), and is not a property of the dynamical equations (3.5). Some other results, regarding the structure of the slowly varying, local envelope of pressure oscillation $R_{j}$ , are new and have implications in the dynamic equations (3.5).

3.3 Amplitudes of standing and spinning solutions

In this section we discuss the existence and the amplitude $A$ of standing and spinning waves as limit cycles of (2.33), i.e. fixed points of (3.5). We then discuss separately the two cases of spinning (3.8b ) and standing waves (3.8a ).

The following implication holds for all values of $\unicode[STIX]{x1D701}$ considered in (2.2):

The proof, discussed in the supplementary materials in § B, exploits the symmetries of the equations and does not add physical insight to the problem.

We were not able to prove that the solutions with coordinates $(A,A,k\unicode[STIX]{x03C0}/2)$ , which can be calculated with (3.15), plus the solutions due to the rotational symmetry of the system are all the possible fixed points of (3.5). It is in general difficult to determine all the fixed points of a nonlinear dynamical system, in this case in three dimensions. In the following we will assume that there are no other solutions. In all the cases studied we could not numerically find any other fixed points, and the system did not converge to other solutions in the time domain.

3.4 Stability of standing and spinning solutions

The previous section investigated the amplitudes of standing and spinning solutions. This and the next section discusses the stability of these solutions, and presents implications for experiments/simulations and industrial applications. We are here discussing the stability of periodic solutions of the equations (2.33), which are fixed points of (3.5). By evaluating the eigenvalues of the Jacobian of these fixed points we can establish necessary and sufficient conditions for the stability of the periodic solutions. This can be done analytically, as discussed in § C of the supplementary materials. Consistently with the system having three variables, we find three eigenvalues and require all three growth rates to be negative in order for the solution to be an attractor.

5.1 The describing function of the experiment

In this subsection we discuss the describing function $Q$ , which in this manuscript expresses the heat release rate response of one flame as a function of the local value of the pressure in the combustion chamber. This will be crucial later in § 5.2.

The describing function $\tilde{Q}(A_{u},\unicode[STIX]{x1D714})$ is measured in a single burner test rig as a function of the forced harmonic axial acoustic velocity just upstream of the flame. We assume good transferability of the measurements in the single burner test rig to the annular configuration, as discussed by Bourgouin et al. (Reference Bourgouin, Durox, Moeck, Schuller and Candel2015). Five values of the root mean square amplitude $A_{u}=\hat{u} _{rms}/\overline{u}$ of this velocity were tested, respectively $0.1,0.2,0.3,0.4,0.5$ . This describing function $\tilde{Q}$ is a function of the longitudinal acoustic velocity at the same location. We can write that

where $Z=\hat{p}/(\unicode[STIX]{x1D70C}c\hat{u} )$ is the impedance of the whole part of the system upstream of the flames for an azimuthal mode of order $n=1$ , $A_{p}$ is the amplitude in terms of the pressure and the complex-conjugate is required because of the sign convention discussed after (2.6). Notice that the ratio of the amplitudes of pressure and velocity, i.e. the gain of $Z$ in (5.1), does not lead to a change of sign in the definition (3.25) of $N_{2}$ , since it appears only as a linear multiplication factor. On the other hand, the phase between pressure and velocities, i.e. $\arg [Z]$ , plays a crucial role in the expression of $N_{2n}$ since it can change the sign of $\text{Re}[Q]$ .

We first model the whole experiment as a network of acoustic elements (Evesque & Polifke Reference Evesque and Polifke2002; Schuermans et al. Reference Schuermans, Bellucci and Paschereit2003; Hirschberg & Rienstra Reference Hirschberg and Rienstra2015) and calculate the eigenvalues of the system. The model consists of the plenum and the combustion chamber and of the transition ducts and burners between plenum and combustion chamber, based on the detailed description of the modelling that also Bourgouin et al. (Reference Bourgouin, Durox, Moeck, Schuller and Candel2015) carry out. We do not discuss the details of the linear acoustic modelling because it does not introduce any element of novelty, but report that when the unsteady flame response is switched off we find three modes that are first-order azimuthal with frequencies at 501, 1171, 2457 Hz, which are within $4$ % of the values calculated by Bourgouin et al. (Reference Bourgouin, Durox, Moeck, Schuller and Candel2015), giving confidence in the validity of the model. We then evaluate (5.1) at $\unicode[STIX]{x1D714}_{s}$ , where $Z(\unicode[STIX]{x1D714}_{s})=0.035\text{e}^{1.879\text{i}}$ was extracted from the model, and the values of $\tilde{Q}$ come partly from Figure 11 of Bourgouin et al. (Reference Bourgouin, Durox, Moeck, Schuller and Candel2015) at 500 Hz and from values reported in the same paper. The resulting nonlinear saturation is reported in figure 8. The predicted phase response for a spinning mode at $A_{u}^{sp}\approx 0.45$ matches the phase difference between chemiluminescence and pressure measured in the rig for the spinning mode, showing good agreement with the assumption of $A_{u}^{sp}=0.5$ made by Bourgouin et al. (Reference Bourgouin, Durox, Moeck, Schuller and Candel2015).

Figure 8. Best reconstruction of the describing function $Q$ of one burner at $\unicode[STIX]{x1D714}_{s}$ as a function of the pressure oscillation in the annulus, scaled however on the horizontal axis with the local amplitude of the acoustic velocity, for ease of comparison. Available data are dotted, while the lines are obtained with inter-extrapolation. (a) The phase $\arg [Q]$ at a value $A_{u}\approx 0.45$ is close to the two values of phase difference between heat release rate and pressure measured experimentally at two azimuthal locations in the annulus. (b) The part of the heat release rate $\text{Re}[Q]$ contributing to the Rayleigh criterion differs significantly from the gain $|Q|$ at large amplitudes. This is because the phase in (a) departs from close to zero as the amplitude increases.

In principle, one would calculate the amplitudes $A_{u}^{sp}$ and $A_{u}^{st}$ by applying respectively equations (3.16) and (3.20), which balance acoustic energy sinks and sources. This however requires a characterization of the acoustic damping of the first azimuthal mode in the experiment, which is not available. We here assume that the standing solution exists with an amplitude of oscillation not larger than $A_{u}=0.5$ , because we are limited by the available experimental data.

5.2 Stability of the standing mode in the experiment

In this subsection we check that the theory predicts that the standing solutions are unstable as found in the experiment. We assume that the number of burners $N_{b}=16$ is large enough to use the simplified condition (3.25), and verify that $N_{2}$ is negative in this case (we do so especially to present a clearer physical interpretation of the condition in terms of integrals. The discussion in terms of discrete sums over the burners using (3.24c ) is harder to interpret.) Since it was not possible in the previous section to predict the amplitude of the standing mode, we study the coefficient $N_{2}$ for the available experimental amplitudes $A_{u}$ . We calculate numerically $N_{2}$ for these amplitudes, obtaining for increasing velocities the values $-16,-17,-26,-58,-174$ , all multiplied by $10^{-4}$ . All values of $N_{2n}$ are negative, so standing solutions at these amplitudes are unstable, and the theory is consistent with the experiment. We can conclude that $N_{2}$ is negative also as described in the first paragraph after equation (3.27), by observing that $\text{Re}[Q]$ in figure 8(b) stays positive and is decreasing with amplitude.

For completeness we describe quantitatively in figure 9 the structure of standing modes in the experiment at different azimuthal locations. All five panels span half an acoustic wavelength, from one pressure node to the next, as exemplified in figure 9(a). We show in figure 9(b) that the flames at the pressure antinode are exposed to a larger amplitude of oscillation and their response is reduced as presented in figure 8(b). The integrand appearing in the definition (3.25) of $N_{2}$ is the product of two functions, reported in figure 9(b,c). We present in figure 9(d) their product, so that the integral of this product is $N_{2}$ .

Before commenting on figure 9(d), we introduce in figure 9(e) the shape of the integrand of $N_{2}$ for the stable standing solution of the example discussed in § 4 and pinpointed in figure 6, so that the reader can visualize side by side two possible cases of an unstable (left) and stable (right) standing mode. In both figures 9(d) and 9(e), the integrand $\text{Re}[Q]\sin (2\unicode[STIX]{x1D703})$ resembles a sombrero, and $N_{2}$ is positive if the central, positive part of the sombrero is larger than the two negative sides. In figure 9(d) we observe that at large amplitudes the tip of the sombrero lowers, due to the drop of the response of the flame at the pressure antinode presented in figure 9(b). Instead, in figure 9(e) the tip of the sombrero is broad, due to the response of the flame that has a gain that is larger at non-zero amplitudes. This confirms the theoretical conclusions reached at the end of § 3.4.2.3 for a set of not very restrictive systems: flames that respond strongly in the linear regime (at pressure nodes) and weakly in the nonlinear regime (at pressure antinodes) lead to unstable standing solutions, with the amplitude of the saturated nonlinear regime fixed by the balance of acoustic energy gains and losses.

Figure 9. (a) View of the amplitude of the standing mode on one half of the annular chamber, between two pressure nodes at $-\unicode[STIX]{x03C0}/4$ and $3\unicode[STIX]{x03C0}/4$ , with a pressure antinode in the middle. In (a), (b) and (d) we colour results for modes at a larger amplitude with a lighter tone of grey. In (b) we show how on a given line the flames at the pressure antinode respond with a smaller $\text{Re}[Q]$ , consistently with figure 8(b). This term $\text{Re}[Q]$ is the first factor in the definition (3.25) of $N_{2}$ that determines the stability of standing modes. The second factor $\sin (2\unicode[STIX]{x1D703})$ is reported in (c) and is negative/positive at pressure nodes/antinodes where flames respond respectively in the linear/nonlinear regime. In (d) we report the integrand that defines $N_{2}$ . By eye one observes that the positive part is smaller than the negative part, especially at large amplitudes, so that the overall integral $N_{2}$ is always negative, and a standing solution at each of the amplitudes would be unstable. In (e) we present for comparison with (d) the integrand defining $N_{2}$ for the standing solution presented in figure 6, which leads instead to a stable standing solution.