Weakly nonlinear analysis of thermoacoustic oscillations in can-annular combustors

Abstract

Can-annular combustors feature clusters of thermoacoustic eigenvalues, which originate from the weak acoustic coupling between $N$ identical cans at the downstream end. When instabilities occur, one needs to consider the nonlinear interaction between all $N$ modes in the unstable cluster in order to predict the steady-state behaviour. A nonlinear reduced-order model for the analysis of this phenomenon is developed, based on the balance equations for acoustic mass, momentum and energy. Its linearisation yields explicit expressions for the $N$ complex-valued eigenfrequencies that form a cluster. To treat the nonlinear equations semianalytically, a Galerkin projection is performed, resulting in a nonlinear system of $N$ coupled oscillators. Each oscillator represents the dynamics of a global mode that oscillates in the whole can-annular combustor. The analytical expressions of the equations reveal how the geometrical and thermofluid parameters affect the thermoacoustic response of the system. To gain further insights, the method of averaging is applied to obtain equations for the slow-time dynamics of the amplitude and phase of each mode. The averaged system, whose solutions compare very well with those of the full oscillator equations, is shown to be able to predict complex transient dynamics. A variety of dynamical states are identified in the steady-state oscillatory regime, including push–push (in-phase) and spinning oscillations. Notably, the averaged equations are able to predict the existence of synchronised states. These states occur when the frequencies of two (or more) unstable modes with nominally different frequencies lock onto a common frequency as a result of nonlinear interactions.

1. Introduction

Thermoacoustic oscillations arise from the coupling of an unsteady heating process with the acoustic modes of an enclosure (Rayleigh 1878). This phenomenon can be found in a variety of technical systems, and it has been a particularly plaguing challenge for the development of modern low-emission gas turbine combustors (Keller 1995; Dowling 1997; Candel 2002; Lieuwen & Yang 2005; Poinsot 2017). The vast majority of gas turbine combustors are equipped with annular or can-annular combustor architectures. These systems nominally feature a high degree of spatial symmetry, which leads to interesting properties in the thermoacoustic mode structure, most notably, the occurrence of degenerate mode pairs (Noiray & Schuermans 2013; Bauerheim, Nicoud & Poinsot 2016; Magri, Bauerheim & Juniper 2016; Mensah et al. 2018). The linear modal structure of annular (Noiray, Bothien & Schuermans 2011) and can-annular (Ghirardo et al. 2019) thermoacoustic systems has been scrutinised to some extent in the past, and can, in fact, be considered as relatively well understood.

In the present work, we explore nonlinear thermoacoustic phenomena in can-annular systems. For this purpose, we conduct a weakly nonlinear analysis of a generic can-annular configuration using the method of averaging. A thorough understanding of the linear modal structure of can-annular systems is a prerequisite; therefore, a concise summary of the relevant elements is given next.

1.1. Linear structure of thermoacoustic modes in can-annular combustors

The spectrum of a can-annular system is characterised by the existence of clusters of eigenvalues that form in the vicinity of the thermoacoustic eigenfrequencies of a single isolated can, i.e. in the absence of coupling between the cans (Bethke et al. 2002; Panek, Farisco & Huth 2017; Ghirardo et al. 2019; von Saldern, Orchini & Moeck 2021c), as illustrated in figure 1. These mode clusters are linked to the weak coupling between nominally identical cans – the weaker the coupling, the closer the modes are in frequency and growth rate. On a conceptual level, it is straightforward to see that a system consisting of weakly coupled identical subsystems features mode clusters. Let us consider a system of $N$ identical combustor cans, initially uncoupled (figure 1a,c). In the composite system, each eigenvalue associated with an isolated can appears $N$ times, and because the cans are all identical, these eigenvalues are the same. They are, therefore, (at least) $N$-fold degenerate. In general, this degeneracy will be semisimple since the states in the individual cans are independent and, hence, span the full eigenspace. If we now introduce a weak coupling between adjacent cans (figure 1b,d), this $N$-fold degeneracy is unfold. Since the coupling is weak and the sensitivity of the eigenvalues is bounded, all of the $N$ eigenvalues will be located in the vicinity of the initial $N$-fold degenerate eigenvalue associated with an isolated can. In other words, they form a cluster.

Figure 1. Spectrum of can-annular systems. Qualitative depiction of the emergence of eigenvalue clusters through the weak coupling of identical subsystems. (a,c) Set of $N$ identical, uncoupled cans. Each eigenvalue of the composite system is $N$-fold degenerate (two are shown, one in the stable and one in the unstable half-plane). (b,d) Set of $N$ identical, weakly coupled cans. The eigenvalues form clusters around the eigenvalues of the uncoupled cans, which leads to the existence of many closely spaced unstable eigenvalues, a feature of can-annular combustors.

Not all of the $N$ unfold eigenvalues in the weakly coupled system will be distinct. Because of the discrete rotational symmetry of the system, most of the eigenvalues will be two-fold degenerate, more precisely, all except for those with azimuthal orders $0$ and $N/2$ (Ghirardo et al. 2019; von Saldern et al. 2021c). Here, we assume that $N$ is even, which is the technically relevant case. When $N$ is odd, all modes except those with azimuthal orders being 0 or multiples of $N$ are degenerate. In summary, the composite system consisting of $N$ identical uncoupled cans features $N$-fold degenerate eigenvalues, with algebraic and geometric multiplicity $N$. When a weak coupling between the cans is introduced, the $N$-fold degenerate eigenvalues split into a cluster consisting of two simple eigenvalues and $N/2-1$ degenerate eigenvalues with algebraic and geometric multiplicity two. The existence of eigenvalue clusters, featuring modes with similar frequencies and growth rates, is, hence, a general feature of can-annular thermoacoustic systems. This is in contrast to most of the literature in thermoacoustics, which considers a single unstable mode or a pair of degenerate unstable modes. When many modes are simultaneously unstable, it is not straightforward to predict the final oscillatory state; sophisticated nonlinear analysis techniques have to be utilised. Moreover, such a scenario provides ample opportunity for non-trivial nonlinear phenomena, in particular, synchronisation between modes that are close in frequency and growth rate (von Saldern, Moeck & Orchini 2021a). This is further elaborated on in the next section.

1.2. Nonlinear effects and synchronisation

Synchronisation generally refers to the emergence of phase coherence between oscillatory states. This phenomenon has been identified in an abundance of physical, chemical and biological systems since Huygens’ observations of the synchronisation of pendulum clocks. More recently, it has been reported on in thermoacoustic systems in a number of studies. To put the present work into proper context, it is important to acknowledge that a crucial aspect is what states are considered to synchronise. Various forms of synchronisation have been considered in thermoacoustic systems, including phase-locking of an oscillatory state to an external forcing (e.g. Kashinath, Li & Juniper 2018; Guan et al. 2019), or the synchronous oscillation of fields (e.g. Pawar et al. 2017; Mondal, Pawar & Sujith 2018). However, these studies are not directly relevant for the synchronisation of self-excited modes that we consider in the present work. To this extent, we would like to recall some basic facts about synchronisation from the seminal book by Pikovsky, Rosenblum & Kurths (2001).

The relevant setting in the present context is the mutual synchronisation of self-sustained oscillators (Pikovsky et al. 2001, see Chap. 4). As highlighted by Pikovsky et al. (2001), synchronisation is only a meaningful concept when the oscillators are self-sustained (possibly self-excited) and if they are weakly coupled. To understand why the first condition is required for synchronisation, consider two oscillators, only one of which is self-excited, with the other being stable. When these are now coupled, even if only weakly, the stable oscillator will always follow the self-excited one. This is not an example of synchronisation, but of response to a forcing. To understand the second condition, consider a scenario in which the coupling is strong. Then, one oscillatory state will always leave a significant trace in the other states it is coupled to. Again, these states share some phase coherence, however, not through a process of synchronisation, merely by virtue of strong coupling. Generally, synchronisation between oscillatory states is favoured through proximity of resonance frequencies. As the resonance frequencies drift apart, an increased level (nonetheless weak) of coupling is required to achieve synchronisation, a phenomenon often represented in the form of Arnold tongues (Pikovsky et al. 2001; Balanov et al. 2009). Also note that synchronisation can occur when the oscillation frequency of one state is close to a multiple (harmonic) or fractional multiple (subharmonic) of the resonance frequency of another state.

Now considering a thermoacoustic can-annular system, there are two possible choices for what to consider as states/oscillators when studying synchronisation.

1. (I) Since the cans are typically weakly coupled, one can consider the thermoacoustic modes of the individual uncoupled cans as oscillatory states that may synchronise under the effect of the coupling. The coupling between adjacent cans is of acoustic or aeroacoustic nature and has reactive and diffusive components; it is linear at leading order but may exhibit finite-amplitude effects (Pedergnana et al. 2021). Chen et al. (2021) and Guan et al. (2021) reported mutual synchronisation of two coupled, nominally identical thermoacoustic oscillators. In the former study, synchronisation was assessed in a numerical model under variation of the resonator length (frequency detuning) and resonator outlet proximity (coupling strength). Phase-coupled solutions were found to be more likely when frequency detuning is small and the coupling strength large, consistent with basic synchronisation theory. Guan et al. (2021) observed similar phenomena experimentally and were able to produce some of them with a low-order model based on coupled van-der-Pol oscillators. Moon et al. (2020a) investigated synchronisation between two nominally identical turbulent combustors, acoustically communicating via an adjustable cross-talk gap. When the combustors were operated at identical conditions, in-phase and phase-opposition synchronised states occurred sporadically in a competitive manner, with asynchronous states existing at low coupling strength. Guan et al. (2022) consider a system composed of four thermoacoustic oscillators and identify synchronisation and other nonlinear phenomena in symmetric and asymmetric settings.

2. (II) Alternatively, one can consider the thermoacoustic modes of the full can-annular system as global oscillatory states that may synchronise. In a symmetric setting, these modes will be of Bloch type (Mensah, Campa & Moeck 2016; von Saldern et al. 2021c). Furthermore, since the thermoacoustic modes of the full system are biorthogonal, they are linearly uncoupled. They are coupled only through the nonlinear terms in the flame response (and possibly in the can-to-can coupling). As such, this coupling can also be considered weak.

Let us consider the differences between the two settings introduced above. We assume here that the can-annular system features discrete rotational symmetry of the order of the number of cans, $N$. This implies that all cans, as well as the coupling between them, are identical. In this case, all the states in approach (I) have nominally the same frequencies and growth rates. Because of the direct (linear) coupling and identical eigenfrequencies, synchronised (in the sense of phase-locked) oscillations among the cans are strongly favoured. We can therefore consider that this approach is perhaps more appropriate in an asymmetric case, when the uncoupled cans have slightly different resonance frequencies. On the other hand, the oscillatory states considered in approach (II) will generally have different resonance frequencies. Yet, when the can-to-can coupling is weak, these differences will be small, as explained in § 1.1. Synchronisation is then a possible scenario, but it is challenging to predict under what conditions it occurs and which modes will be involved.

Another point to consider is that approach (II) can be adopted for studying synchronisation phenomena in annular or essentially longitudinal thermoacoustic systems too. Approach (I), on the other hand, would be unsuitable in this context. This is because the individual flames in an annular or sequential combustor cannot be considered as oscillators (assuming that there are no unstable intrinsic thermoacoustic modes). Note, in any case, that synchronisation of thermoacoustic modes in an annular or longitudinal combustor is generally much less likely because the modes are typically not close in frequency. However, it can occur in special circumstances, one example being the manifestation of the slanted mode (Bourgouin et al. 2015; Moeck et al. 2019), which results from the synchronisation of an azimuthal and an axisymmetric thermoacoustic mode. Synchronisation between two longitudinal modes with one having approximately twice the frequency of the other was observed in a sequential combustion system and successfully modelled based on two coupled stochastic oscillators (Bonciolini & Noiray 2019). Acharya, Bothien & Lieuwen (2018) investigated the nonlinear interaction of two thermoacoustic modes with close eigenfrequencies based on an oscillator model derived from the Euler equations using a Galerkin expansion. In the present work we will also adopt approach (II) and consider synchronisation between thermoacoustic modes of the full system.

The evolution of the dynamics of a combustor with multiple linearly unstable eigenvalues is governed by the nonlinear modal coupling between the unstable modes. There exist in the literature studies that have focused on assessing this interplay for specific scenarios, with few (two or three) linearly unstable modes. In particular, by using the method of averaging, Noiray et al. (2011) showed how a single pair of degenerate azimuthal modes in an annular configuration saturated by a cubic nonlinearity always converges to spinning solutions; asymmetries are required to stabilise standing azimuthal wave solutions (Noiray & Schuermans 2013). In more recent work, these results have been extended to account also for the effect of noise on the dynamics, showing that (combustion) noise tends to promote the existence of mixed modes (Faure-Beaulieu et al. 2021). Instead, Moeck & Paschereit (2012), Orchini & Juniper (2016) and Acharya et al. (2018) considered the coexistence of two non-degenerate unstable modes, and derived equations for the (slow) evolution of the oscillation amplitudes, accounting for the modal interaction, using the describing function method or a two time scales approach. They discussed how these equations are able to predict the existence of periodic solutions, in which one of the modes dominates over the other, or quasiperiodic solutions, in which both modes participate in the steady-state oscillations. One significant challenge in these types of analysis is that predicting the existence and stability of multimode solutions requires rather detailed knowledge of the nonlinear terms in the system – a single-input describing function is not sufficient. Moeck et al. (2019) considered the nonlinear interaction of three simultaneously unstable modes in an annular combustor, specifically an axisymmetric mode and a pair of degenerate azimuthal modes. This picture becomes significantly more complicated for can-annular systems, because (at least) $N$ modes – some of which are degenerate – need to be retained in a nonlinear analysis, and their coupling accounted for. Unfortunately, the large number of nonlinear terms that arise in the study of a can-annular system inhibits an analytical analysis of the solutions and their stability.

1.3. Scope

This work aims at developing a weakly nonlinear framework for describing the nonlinear dynamics of can-annular thermoacoustic systems, which are characterised by the presence of clusters of linearly unstable modes. In § 2 we derive from first principles the nonlinear equations that govern the coupled dynamics in the configuration considered. The linear properties of the equations are discussed in § 3, and explicit expressions for the growth rates and frequencies of the eigenvalues as a function of the geometrical and coupling parameters are derived. In § 4 we employ the method of averaging to derive equations for the slow-dynamics of the system. These nonlinear, coupled equations govern the evolution of the amplitudes and frequencies of all the modes that contribute to the oscillations. Lastly, in § 5 we numerically verify the validity of the derived equations, and we use them to discuss what dynamical solutions are supported by can-annular systems, with a focus on the existence of synchronised states in these configurations.

2. Nonlinear time-domain equations for the dynamics of thermoacoustic modes in can-annular systems

We consider a can-annular system composed of $N$ identical combustor cans coupled to their neighbours at the downstream end by small apertures, which allow for acoustic communication between neighbouring cans. Each can hosts an acoustically compact heat source. The cans form a periodic array, with the $N$th can connected to the first. The geometry is shown schematically in figure 2. The thermoacoustic response of each can, labelled by the index $j$, is described by the conservation of mass, momentum and energy. For the purpose of this study, mean flow effects in the cans are neglected by invoking a low Mach-number assumption (Culick 2006). The linearised equations for the conservation of mass, momentum and energy, respectively, read

(2.1a)$$\begin{gather} \frac{\partial \rho'_j}{\partial t} + \bar{\rho}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u}'_j = 0, \end{gather}$$

(2.1b)$$\begin{gather}\frac{\partial \boldsymbol{u}'_j}{\partial t} + \frac{1}{\bar{\rho}} \boldsymbol{\nabla} p'_j = 0, \end{gather}$$

(2.1c)$$\begin{gather}\frac{\partial p'_j}{\partial t} - \bar{c}^2 \frac{\partial \rho'_j}{\partial t} = (\gamma-1)\dot{q}'_j, \end{gather}$$

where $\boldsymbol {u}' \equiv [u'_x,u'_y,u'_z]$ is the acoustic velocity field, $\bar {\rho }$ the mean density, $p'$ the unsteady pressure fluctuation, $\bar {c}$ the speed of sound (assumed constant), $\gamma$ the heat capacity ratio, and $\dot {q}'$ the unsteady heat release rate due to combustion in a small, acoustically compact volume. These equations may be combined to obtain a second-order partial differential equation for the acoustic pressure fluctuations

(2.2)\begin{equation} \frac{\partial^2 p'_j}{\partial t^2} - \bar c^2 \nabla^2 p'_j = (\gamma-1)\frac{\partial \dot{q}'_j}{\partial t}, \end{equation}

which is a wave equation with a source term on the right-hand side, due to the generation of acoustic waves by the unsteady heat release rate. Once a model for $\dot {q}'$ and boundary conditions are provided, the thermoacoustic equation (2.2) can be solved.

Figure 2. Schematic of the can-annular geometry considered in this study. (a) Top view of the $N$ connected cans, with depth $H$ and width $B$ in the azimuthal direction $\theta$. (b) Unwrapped side view of the geometry. The cans have height $L$, contain unsteady heat release rate sources $q'$ and are connected by a downstream aperture of axial extension $l_g$. This gap allows for acoustic communication with the neighbouring cans. Note that the azimuthal angle $\theta$ corresponds to the Cartesian $y$ coordinate in the unwrapped representation.

2.1. Coupling between cans

When the cans are connected through small apertures, one needs to account for the acoustic communication from a can to the neighbouring ones. From a can perspective, this can be taken into account as an acoustic mass source/sink along the extension of the aperture (Fournier et al. 2021; von Saldern et al. 2021c). The aperture is assumed to have surface area $A_g \equiv l_g H$, where $l_g$ is the axial extension of the apertures and $H$ their depth. By denoting with $\varOmega _{j,j+1}$ the surface of the aperture that allows for acoustic communication between cans $j$ and can $j+1$, a uniform azimuthal transverse velocity is assumed to be present in the aperture

(2.3)\begin{equation} u'_{y,j}(\boldsymbol{x},t) = \left\{ \begin{array}{@{}ll} u'_{j,j+1}(t) & \mbox{if } \boldsymbol{x} \in \varOmega_{j,j+1}, \\ 0 & \mbox{otherwise.} \end{array} \right. \end{equation}

The characteristic frequency of the cluster we are interested in modelling is assumed to be below the cut-on frequency of transverse oscillations in a can, which depends on the speed of sound and the transverse dimensions of the can (Rienstra & Hirschberg 2004). Typical cut-on frequencies for transverse modes in realistic configurations are of the order of a few kilohertz. When transverse oscillations are cut off, the acoustics in each can is well-approximated by planar wave propagation along the axial direction $x$ (Rienstra & Hirschberg 2004). While it is true that near the connection between the cans and the gap a non-planar near field exists, as highlighted in Ghirardo et al. (2019), non-planar waves rapidly decay for frequencies below the lowest transverse mode cut-on frequency. The depiction of the acoustics established here is analogous to that of conventional T-junctions used, e.g. to represent the acoustics at the interface of a duct and a transversally mounted Helmholtz resonator (Bothien & Wassmer 2015). By considering frequencies for which all transverse modes are cut off, we can integrate the conservation equations (2.1) over the cross-section of the can, $A_c\equiv BH$, where $B$ is the width of a can. For the conservation of mass (2.1a) this reads

(2.4)\begin{equation} \int_0^B {{\rm d}y} \int_0^H {\rm d}z \left[\frac{\partial \rho'_j}{\partial t} + \bar{\rho} \frac{\partial u'_{x,j}}{\partial x} + \bar{\rho} \frac{\partial u'_{y,j}}{\partial y} + \bar{\rho} \frac{\partial u'_{z,j}}{\partial z}\right] = 0. \end{equation}

Furthermore, under the planar wave assumption the acoustic variables are only functions of $x$, the axial coordinate, and (2.4) simplifies to

(2.5)\begin{align} &A_c\left[ \frac{\partial \rho'_j}{\partial t} + \bar{\rho} \frac{\partial u'_{x,j}}{\partial x}\right] + \bar{\rho} H \left[u'_{y,j}(\kern0.7pt y=B) - u'_{y,j}(\kern0.7pt y=0)\right]\varTheta_g(x) \nonumber\\ &\quad + \bar{\rho}B \left[ u'_{z,j}(z=H) - u'_{z,j}(z=0) \right]= 0. \end{align}

In (2.5), $u'_{y,j}(\kern0.7pt y=B)$ is the transverse acoustic velocity on the boundary separating cans $j$ and $j+1$, and $u'_{y,j}(\kern0.7pt y=0)$ that on the boundary separating cans $j-1$ and $j$. These are non-zero only along the extension of the aperture, as per (2.3), which is accounted for by introducing the Heaviside function $\varTheta (x)$ and by defining

(2.6)\begin{equation} \varTheta_g(x) = \varTheta(x-x_g)\varTheta(x_g+l_g-x), \end{equation}

where $x_g$ is the upstream end of the aperture. Furthermore, $u'_{z,j} = 0$ at $z=0$ and $z=H$ because of the presence of physical walls, at which the velocity vanishes. To ease the notation, from now onward we will drop the subscript $_x$ from the axial component of the velocity in the cans. Also note that the gap between two combustor cans occupies zero volume in our model, which essentially means that it is assumed acoustically compact.

By integrating also the momentum and energy equations, the conservation laws read

(2.7a)$$\begin{gather} \frac{\partial \rho'_j}{\partial t} + \bar{\rho} \frac{\partial u'_j}{\partial x} = \bar{\rho} \frac{H}{A_c} \left[u'_{j-1,j} - u'_{j,j+1} \right] \varTheta_g(x), \end{gather}$$

(2.7b)$$\begin{gather}\frac{\partial u'_j}{\partial t} + \frac{1}{\bar{\rho}}\frac{\partial p'_j}{\partial x} = 0, \end{gather}$$

(2.7c)$$\begin{gather}\frac{\partial p'_j}{\partial t} - \bar{c}^2 \frac{\partial \rho'_j}{\partial t} = (\gamma-1)\dot{q}'_j. \end{gather}$$

In the absence of transverse acoustic fluctuations ($u'_{j,j+1} = u'_{j-1,j} = 0$), these equations reduce to the classic one-dimensional thermoacoustic equations in a duct. To obtain a Helmholtz-like equation analogous to (2.2) that accounts for the presence of transverse acoustic fluctuations between the cans, we substitute the conservation of mass (2.7a) into the conservation of energy (2.7c),

(2.8)\begin{equation} \frac{\partial p'_j}{\partial t} + \bar{\rho} \bar{c}^2 \frac{\partial u'_j}{\partial x} = \bar{\rho} \bar{c}^2 \frac{H}{A_c} \left[u'_{j-1,j} - u'_{j,j+1}\right] \varTheta_g(x) + (\gamma-1)\dot{q}'_j, \end{equation}

and then combine the latter with the conservation of momentum to obtain a second-order partial differential equation for the acoustic pressure, in which a linear damping term with damping coefficient $\alpha$ is included (Noiray et al. 2011):

(2.9)\begin{equation} \frac{\partial^2 p'_j}{\partial t^2} + \alpha\frac{\partial p'_j}{\partial t} - \bar{c}^2 \frac{\partial^2 p'_j}{\partial x^2} = \bar{\rho} \bar{c}^2 \frac{H}{A_c} \left[\frac{\partial u'_{j-1,j}}{\partial t} - \frac{\partial u'_{j,j+1}}{\partial t}\right] \varTheta_g(x) + (\gamma-1)\frac{\partial \dot{q}'_j}{\partial t}. \end{equation}

The structure of the thermoacoustic equation (2.9) for can-annular systems makes it clear how the dynamics in each can is driven by two sources: the classic volumetric energy addition given by the heat release rate $\dot {q}'$ in that can, and the can-annular-specific acoustic fluctuations at the connections with the neighbouring cans.

2.2. Closure of the equations

To close (2.9) we need to express the source terms, viz. heat release rate and transverse acoustic velocity fluctuations, as functions of the pressure variables. For the heat release rate, we shall employ the simple but widely used cubic saturation relation between the unsteady heat release and the pressure fluctuations, by defining (Noiray et al. 2011; Ghirardo, Juniper & Bothien 2018)

(2.10)\begin{equation} \dot{q}_j'(x,t) \equiv \left[\beta p_j'(x,t) - \kappa p_j'(x,t)^3 \right]\delta(x_f), \end{equation}

where $\delta (x)$ denotes Dirac's delta. Flames utilised in gas turbine combustors are typically velocity sensitive; however, the pressure and velocity fluctuations at the flame are directly related through the upstream can impedance. In principle, the combustor cans may communicate with their neighbours through the upstream plenum, but this is usually considered of less importance compared with the cross-talk gap at the turbine inlet (Ghirardo et al. 2019).

The acoustic communication between the cans is present only when the pressure difference between neighbouring cans is non-zero, which in turn causes transverse acoustic velocity fluctuations through the gaps. These can then be linked to the pressure difference between two cans by an impedance $\zeta$, sometimes expressed in terms of the Rayleigh conductivity $K_R$ (Howe 1997):

(2.11a,b)\begin{equation} \frac{\hat{p}_{j}-\hat{p}_{j+1}}{\hat{u}_{j,j+1}} = \bar{\rho} \bar{c} \zeta, \quad \zeta = \frac{s A_g}{\bar{c} K_R}. \end{equation}

The latter equation is expressed in the frequency domain, for which we adopt the convention $p'(x,t) \rightarrow \hat {p}(x) {\rm e}^{st}$, with $s\equiv \sigma + \mathrm {i} \omega$, $\sigma$ denoting the growth rate and $\omega$ the angular frequency.

By considering (2.11a,b) also for the pressure difference between can $j$ and $j-1$, one can express the transverse acoustic velocity fluctuations as a function of the acoustic pressure fluctuations in neighbouring cans as (von Saldern et al. 2021c)

(2.12)\begin{equation} \hat{u}_{j-1,j} - \hat{u}_{j,j+1} = \frac{1}{{\bar{\rho} \bar{c}}} \frac{1}{\zeta} \left( \hat{p}_{j-1}-2\hat{p}_{j}+\hat{p}_{j+1} \right). \end{equation}

The latter relation is valid at the connection between the cans and the annular section, $x=x_g$, and should be interpreted as an integral average over the extension of the gap, $l_g$, for an acoustically non-compact connection (von Saldern, Orchini & Moeck 2021b).

In the present study, we consider the Rayleigh conductivity to be a complex-valued constant. This is a fair approximation when the analysis is restricted to the vicinity of a frequency $s = \mathrm {i}\omega _0$ of interest, which corresponds to the frequency at which the cluster of eigenvalues investigated in our analysis is found. Following Howe (1997) (note that due to a different convention in the definition of the Laplace transform, the Rayleigh conductivity (2.13) is the complex conjugate of that reported by Howe (1997)), we define

(2.13)\begin{equation} K_R\equiv 2 r_g(\varGamma+\mathrm{i}\varDelta), \end{equation}

where $r_g \equiv \sqrt {A_g/{\rm \pi} }$ is a characteristic length (equivalent radius) of the aperture. The real part of the Rayleigh conductivity represents a stiffness term, and, when $s=\mathrm {i}\omega _0$, corresponds to the imaginary part of $\zeta$, the can-to-can reactance. Similarly, the imaginary part of $K_R$ represents a damping term, and it corresponds to the real part of the impedance $\zeta$, the can-to-can resistance.

By substituting the expression for the Rayleigh conductivity in (2.12), and recalling its relation (2.11a,b) with the can-to-can impedance $\zeta$, we obtain

(2.14)\begin{equation} s\left[\hat{u}_{j-1,j} - \hat{u}_{j,j+1}\right] = 2r_g\frac{\varGamma+\mathrm{i}\varDelta} {\rho A_g} \left( \hat{p}_{j-1}-2\hat{p}_{j}+\hat{p}_{j+1} \right), \end{equation}

which, in the vicinity of the frequency of interest $s = \mathrm {i}\omega _0$, can be written in the time domain as

(2.15)\begin{equation} \frac{\partial}{\partial t}\left[u'_{j-1,j} - u'_{j,j+1}\right] = \frac{2r_g} {\rho A_g} \left[\varGamma + \frac{\varDelta} {\omega_0} \frac{\partial}{\partial t}\right] \left( p'_{j-1}-2p'_{j}+p'_{j+1} \right). \end{equation}

The latter can be substituted in (2.9), leading to a system of coupled wave equations describing the dynamics of the can-annular configuration

(2.16)\begin{align} \frac{\partial^2 p'_j}{\partial t^2} + \alpha\frac{\partial p'_j}{\partial t} - \bar{c}^2 \frac{\partial^2 p'_j}{\partial x^2} &= \bar{c}^2 \frac{H2r_g}{A_c A_g} \left[\varGamma + \frac{\varDelta} {\omega_0} \frac{\partial}{\partial t}\right] \left( p'_{j-1}-2p'_{j}+p'_{j+1} \right) \varTheta_g(x)\nonumber\\ &\quad + (\gamma-1)\beta\frac{\partial p'_j}{\partial t}\delta(x_f) - (\gamma-1)\kappa\frac{\partial p'^3_j}{\partial t}\delta(x_f). \end{align}

Equation (2.16) reveals that there exist two coupling mechanisms between neighbouring cans: one proportional to the pressure difference (reactive coupling, proportional to $\varGamma$) and one proportional to difference in the rate of change of the pressure (diffusive coupling, proportional to $\varDelta$).

2.3. Galerkin projection

To study (2.16), we employ the Galerkin method to project the dynamics onto a set of prescribed can mode shapes. Because we are interested in studying the dynamics in the vicinity of a thermoacoustic eigenvalue, it is reasonable to assume that, in each can, the dynamics is dominated by a single Galerkin mode

(2.17)\begin{equation} p_j'(x,t) = \eta_{j}(t)\psi_{j}(x). \end{equation}

The mode shapes $\psi _{j}(x)$ are prescribed as the passive thermoacoustic modes with frequency $\omega _0$, identified in each single, uncoupled can, and are subject to the normalisation condition

(2.18)\begin{equation} \int_0^L |\psi_{j}(x)|^2 \, {{\rm d}x} = 1. \end{equation}

In this study, we shall consider a can-annular set-up with nominal rotational symmetry, so that the Galerkin mode is identical in all cans, $\psi _{j}(x) = \psi (x)$, as guaranteed by Bloch's theorem (von Saldern et al. 2021c).

By substituting the Galerkin ansatz in (2.16), multiplying the equation by $\psi ^*$, and integrating over $x$, we obtain a set of coupled equations for the amplitudes of the modes:

(2.19)\begin{align} \ddot{\eta}_{j} + \omega_0^2 \eta_{j} &= \left((\gamma-1) |\psi(x_f)|^2 \beta - \alpha \right)\dot{\eta}_{j} \nonumber\\ &\quad +\bar{c}^2\frac{H2r_g}{A_g A_c} \int_{x_g}^{x_g+l_g}|\psi(x)|^2 \,\mathrm{d}x \, \varGamma \left[ \left( \eta_{j-1}-\eta_{j}\right)+\left(\eta_{j+1}-\eta_{j} \right) \right]\nonumber\\ &\quad + \bar{c}^2 \frac{H2r_g}{A_c A_g} \int_{x_g}^{x_g+l_g}|\psi(x)|^2 \,\mathrm{d}x \, \frac{\varDelta} {\omega_0} \left[ \left( \dot{\eta}_{j-1}-\dot{\eta}_{j}\right)+\left(\dot{\eta}_{j+1}-\dot{\eta}_{j} \right) \right] +\mathrm{NL}(\eta_j,\dot{\eta}_j). \end{align}

On the left-hand side, we have isolated the oscillatory dynamics of the mode of interest in can $j$ in the absence of (i) damping, (ii) coupling with the other cans and (iii) heat release rate dynamics. On the right-hand side we distinguish three contributions to the dynamics: (i) a linear damping/driving term, related to the linear flame response, $\beta$, and the acoustic damping, $\alpha$; (ii) two linear coupling terms, one driven by the reactive coupling $\varGamma$, the other driven by the diffusive coupling $\varDelta$; (iii) nonlinear terms ($\mathrm {NL}$) that couple the dynamics of the modes and, in our model, originate from the cubic saturation of the heat release rate. For the nonlinear cubic saturation of the heat release rate employed in our analysis, the nonlinear term reads

(2.20)\begin{equation} \mathrm{NL}(\eta_j,\dot{\eta}_j) ={-} 3 (\gamma-1) |\psi(x_f)|^4 \kappa \eta_j^2 \dot{\eta}_j. \end{equation}

The coefficients $\varDelta$ and $\varGamma$ of the coupling terms between the oscillators in (2.19) are all identical because we have assumed a perfect rotational symmetry, implying that the coupling type/strength remains constant (uniform) throughout the entire network of oscillators. This symmetry assumption could be nonetheless relaxed, resulting in can-specific coupling terms $\varDelta _j$ and $\varGamma _j$. This asymmetric scenario will, however, not be considered in this study.

We emphasise that, even though the Galerkin modes are identical, the dynamics of the can-annular system is nonetheless described by a total of $N$ second-order equations, one for the acoustic oscillations in each can. The collective behaviour of the oscillation amplitudes in each can and their phase patterns are, however, related to the dynamics of the global modes of the entire can-annular combustor, as shown in the next section.

3. Linear response

By neglecting the nonlinear terms in (2.19), we investigate the linear stability of the rotationally symmetric, coupled can-annular system. We introduce the scaled parameters

(3.1)\begin{equation} \left.\begin{gathered} \tilde{\varGamma}\equiv \bar{c}^2 \frac{H2r_g}{A_cA_g} \int_{x_g}^{x_g+l_g}|\psi(x)|^2 \,\mathrm{d}x \, \frac{\varGamma}{\omega_0},\quad \tilde{\varDelta}\equiv \bar{c}^2 \frac{H2r_g}{A_cA_g} \int_{x_g}^{x_g+l_g}|\psi(x)|^2 \,\mathrm{d}x \, \frac{\varDelta}{\omega_0}, \\ \tilde{\beta}\equiv (\gamma-1) |\psi(x_f)|^2 \beta, \end{gathered}\right\} \end{equation}

so that (2.19) can be written in the more compact form

(3.2)\begin{equation} \ddot{\eta}_j + \omega_0^2\eta_j = \left[\tilde{\beta} - \alpha -2\tilde{\varDelta}\right ] \dot{\eta}_j -2 \omega_0\tilde{\varGamma}\eta_j + \omega_0\tilde{\varGamma} \eta_{j-1}+ \tilde{\Delta} \dot{\eta}_{j-1} + \omega_0\tilde{\varGamma}\eta_{j+1} +\tilde{\Delta}\dot{\eta}_{j+1}. \end{equation}

The set of $N$ equations for $j=1,\ldots,N$ is cast in state-space form by introducing the state variables $\zeta _j \equiv \dot {\eta }_j$, so that $\dot {\zeta }_j = \ddot {\eta }_j$. By collecting all the state variables in the state vector $\boldsymbol {z} \equiv [\eta _1,\zeta _1,\ldots,\eta _N,\zeta _N]^{\rm T}$, (3.2) can be rewritten in matrix form as

(3.3)\begin{equation} \frac{{\rm d}\boldsymbol{z}}{{\rm d}t} =\mathcal{M} \boldsymbol{z}, \end{equation}

where the matrix $\mathcal {M}$ is block-circulant and has the form

(3.4)\begin{equation} \mathcal{M} = \begin{bmatrix} \mathcal{A} & \mathcal{B} & 0 & \ldots & 0 & \mathcal{B} \\ \mathcal{B} & \mathcal{A} & \mathcal{B} & \ldots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \ldots & \mathcal{A} & \mathcal{B} \\ \mathcal{B} & 0 & 0 & \ldots & \mathcal{B} & \mathcal{A} \end{bmatrix} \end{equation}

with

(3.5a,b)\begin{equation} \mathcal{A} \equiv \begin{bmatrix} 0 & 1 \\ -\omega_0^2-2\omega_0\tilde{\varGamma} & \tilde{\beta}-\alpha-2\tilde{\varDelta} \end{bmatrix}, \quad \mathcal{B} \equiv \begin{bmatrix} 0 & 0 \\ \omega_0\tilde{\varGamma} & \tilde{\varDelta} \end{bmatrix}. \end{equation}

In the absence of can-to-can coupling ($\tilde {\varGamma } = 0$ and $\tilde {\varDelta } = 0$ when $r_g=0$), the matrix $\mathcal {B}$ vanishes, and $\mathcal {M}$ becomes block-diagonal. Each block represents the acoustic dynamics in a decoupled can, with natural (undamped) frequency $\omega _0$ and damping coefficient $\alpha -\tilde {\beta }$. In this limit, the matrix has $N$-fold degenerate eigenvalues. This degeneracy manifests in the ambiguity between the phase of the oscillations in the cans: indeed, with no acoustic communication, any arbitrary phase pattern is possible. As we shall see, the coupling between the cans resolves (some of) the degeneracy: only specific phase patterns are allowed, and different oscillation patterns have slightly different frequencies, which leads to the creation of a cluster of eigenvalues.

To calculate the eigenvalues of (3.3), we recall that block-circulant matrices can always be block-diagonalised (Olson et al. 2014). This is accomplished by multiplying $\mathcal {M}$ from the left and from the right by the unitary matrix $\mathcal {E}_N \otimes \mathcal {I}_2$, where $\otimes$ denotes the Kronecker product, $\mathcal {I}_2$ is the two-dimensional identity matrix, and $\mathcal {E}_N$ is the Fourier matrix of size $N$, defined by

(3.6) \begin{equation} \mathcal{E}_N \equiv \frac{1}{\sqrt{N}}\begin{bmatrix} 1 & 1 & 1 & \ldots & 1 \\ 1 & w_N & w_N^2 & \ldots & w_N^{N-1} \\ 1 & w_N^2 & w_N^4 & \ldots & w_N^{2(N-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & w_N^{N-1} & w_N^{2(N-1)} & \ldots & w_N^{(N-1)(N-1)} \end{bmatrix}, \end{equation}

with $w_N \equiv {\rm e}^{\mathrm {i}({2{\rm \pi} }/{N})}$. This yields

(3.7)\begin{equation} (\mathcal{E}_N \otimes \mathcal{I}_2)^H \mathcal{M} (\mathcal{E}_N \otimes \mathcal{I}_2) = \mathcal{D} = \begin{bmatrix} \mathcal{D}_0 & & & \boldsymbol{0} \\ & \mathcal{D}_1 & & \\ & & \ddots & \\ \boldsymbol{0} & & & \mathcal{D}_{N-1} \end{bmatrix}, \end{equation}

where $\mathcal {D}$ is a $2N\times 2N$ block-diagonal matrix, whose blocks $\mathcal {D}_{k-1}$ have dimensions $2\times 2$, and are defined by (Olson et al. 2014)

(3.8)\begin{equation} \mathcal{D}_{k-1} \equiv \mathcal{A} + \mathcal{B}\omega_N^{(k-1)} + \mathcal{B}\omega_N^{(N-1)(k-1)} = \mathcal{A} + 2\mathcal{B}\cos\left(\frac{2{\rm \pi}}{N}(k-1)\right), \quad k = 1,\ldots,N. \end{equation}

By recognising $k-1$ as the Bloch number $b$, the block matrix $\mathcal {D}_b$ reads

(3.9)\begin{align} \mathcal{D}_b = \begin{bmatrix} 0 & 1 \\ -\omega_0^2-4\omega_0\tilde{\varGamma}\sin^2\left(\dfrac{{\rm \pi} b}{N}\right) & \tilde{\beta}-\alpha-4\tilde{\varDelta}\sin^2\left(\dfrac{{\rm \pi} b}{N}\right) \end{bmatrix}, \quad \mbox{for } b = 0,\ldots,N-1, \end{align}

and the dynamics of the $b$th mode is described by the linear eigenvalue problem

(3.10)\begin{equation} \mathcal{D}_b \boldsymbol{v}_b = s_b \boldsymbol{v}_b. \end{equation}

The growth rate and angular frequency of the eigenvalues $s_b \equiv \sigma _b + \mathrm {i} \omega _b$ of the matrix $\mathcal {D}_b$, defined in (3.9), are

(3.11a)$$\begin{gather} \sigma_b = \frac{\tilde{\beta}-\alpha}{2}-2\tilde{\varDelta}\sin^2\left(\frac{{\rm \pi} b}{N}\right), \end{gather}$$

(3.11b)$$\begin{gather}\omega_b = \sqrt{\omega_0^2 -\left(\frac{\tilde{\beta}-\alpha}{2}\right)^2 + 2\left[\left(\tilde{\beta}-\alpha\right)\tilde{\varDelta}+2\omega_0\tilde{\varGamma} \right]\sin^2\left(\frac{{\rm \pi} b}{N}\right) - 4\tilde{\varDelta}^2 \sin^4\left(\frac{{\rm \pi} b}{N}\right)}. \end{gather}$$

Equations (3.11) are analytical expressions that describe the distribution of all the can-annular eigenvalues in a cluster associated with a single-can mode as a function of the system parameters (nominal frequency, nominal damping and real and imaginary parts of the coupling impedance). The Bloch number $b$ takes value from 0 to $N-1$, or equivalently from $-N/2+1$ to $N/2$ if $N$ is even, or from $-(N-1)/2$ to $(N-1)/2$ if $N$ is odd. Because $\sin ^2({\rm \pi} b/ N)$ is an even function, it holds that $s_{b} = s_{-b}$, leading to two-fold degeneracy of all modes with $b\neq 0$ and $b\neq N/2$. For $b=0$ the expressions retrieve the eigenvalue of a single, isolated can, given by

(3.12)\begin{equation} s_0 = \frac{\tilde{\beta}-\alpha}{2} \pm \mathrm{i}\omega_0\sqrt{1 -\left(\frac{\tilde{\beta}-\alpha}{2\omega_0}\right)^2 }. \end{equation}

The eigenvectors $\boldsymbol {z}_b$ of the linear system (3.3) are given by $(\mathcal {E}_N \otimes \mathcal {I}_2) \boldsymbol {v}_b$, where $\boldsymbol {v}_b$ is the $2\times 1$ eigenvector in the block-diagonal formulation. Because of the nature of $(\mathcal {E}_N \otimes \mathcal {I}_2)$, it can be shown (Olson et al. 2014) that both the states contained in the eigenvectors $\boldsymbol {z}_b$ follow a discrete Fourier mode pattern of order $b$.

4. Weakly nonlinear dynamics

When performing a weakly nonlinear analysis, all the modes that are linearly unstable, and thus can contribute to the dynamics, need to be retained. Due to the formation of eigenvalue clusters in can-annular systems, if the uncoupled eigenvalue $\omega _0$ is marginally stable, so are all the eigenvalues of the weakly coupled can-annular system – as per (3.11a) when $\tilde {\varDelta }$ is small compared with $(\tilde {\beta } - \alpha )$. All the $N$ modes in this cluster must therefore be retained. In our analysis, it is assumed that a single cluster is (marginally) unstable.

Weakly nonlinear analysis of thermoacoustic oscillations in longitudinal and annular configurations has been conducted previously (Ghirardo, Juniper & Moeck 2016; Orchini, Rigas & Juniper 2016). In these works, however, the interaction of multiple linearly unstable modes with close but initially unequal frequencies and growth rates, which is the generic setting for can-annular systems, was not considered. Rather than describing the amplitude of the oscillations in each can, it is more appropriate for a nonlinear analysis to describe the amplitude of the global modes. This is because, from the results of the linear analysis in § 3, we know that the can-annular pressure mode shapes in a cluster have the form of Fourier modes with azimuthal order up to $N/2$. We wish to track the evolution of the amplitude of these physically motivated mode shapes, which exponentially grow in time in the linear regime, and lead to nonlinear modal interactions when the amplitudes reach non-negligible values. Thus, by using the results of the linear analysis, and by recalling that the $y$ direction in figure 2 is equivalent to the azimuthal direction $\theta$, the global pressure dynamics can be described as

(4.1)\begin{align} p_j(x,t) = p(x,\theta_j,t) &= \eta_{0c}(t)\psi(x) + \eta_{N/2c}(t)\psi(x)\cos(N/2 \, \theta_j) \nonumber\\ &\quad +\sum_{m=1}^{N/2-1} \left[ \eta_{mc}(t)\psi(x)\cos(m\theta_j) + \eta_{ms}(t)\psi(x)\sin(m\theta_j) \right]. \end{align}

The subscript $_{mc}$ indicates the mode with structure $\cos (m\theta )$ for $m=0,\ldots,N/2$, whereas $_{ms}$ the mode with structure $\sin (m\theta )$ for $m=1,\ldots,N/2-1$. Note that, from the linear stability results of § 3, we have exploited the fact that the $m=0$ (push–push) and $m=N/2$ (push–pull) are non-degenerate.

Using the ansatz (4.1) in the can-annular Helmholtz equation (2.16), and defining the scaled saturation constant $\tilde {\kappa } \equiv (\gamma -1) |\psi (x_f)|^4 \kappa$, together with the coefficients (3.1), yields a set of second-order, coupled, nonlinear ordinary differential equations governing the dynamics of each mode:

(4.2)\begin{equation} \ddot{\eta}_{mc/s} + \omega_0^2 \eta_{mc/s} = (\tilde{\beta} - \alpha) \dot{\eta}_{mc/s} - 4 \sin^2\left(\frac{m{\rm \pi}}{N}\right) \left(\omega_0\tilde{\varGamma}\eta_{mc/s} + \tilde{\Delta} \dot{\eta}_{mc/s}\right) - 3\tilde{\kappa}\mathcal{N}_{mc/s}, \end{equation}

with $\mathcal {N}_{mc/s}$ defined by

(4.3a,b)\begin{align} \mathcal{N}_{mc} = \frac{1}{(1+\delta_{m0}){\rm \pi}}\int_0^{2{\rm \pi}} \cos(m\theta) \frac{\partial p'}{\partial t} p'^2 \,{\rm d}\theta, \quad \mathcal{N}_{ms} = \frac{1}{\rm \pi}\int_0^{2{\rm \pi}} \sin(m\theta) \frac{\partial p'}{\partial t} p'^2 \,{\rm d}\theta. \end{align}

The Kronecker delta $\delta _{m0}$ in the normalisation accounts for the fact that all modes with $m\neq 0$ average to ${\rm \pi}$ when integrated over the circumference, whereas the mode with $m=0$ averages to $2{\rm \pi}$. The expressions for $\mathcal {N}$ can be made explicit by means of trigonometric identities, but become very large for increasing values of $N$. For $N=8$ cans, each $\mathcal {N}$ contains over 100 terms. Their explicit expressions are thus not fully reported in the manuscript, but an example calculation of one of these terms is provided in the Appendix A.

Note that the left-hand side is identical for all the modes. This reflects the fact that, in the absence of coupling, the cans are free to oscillate at the same frequency with any circumferential phase pattern. The linear corrections due to the can-to-can coupling to the frequencies (given by the terms proportional to $\tilde {\varGamma }$) and to the growth rates (given by the terms proportional to $\tilde {\varDelta }$) are considered as perturbations of the undamped oscillators, all sharing the same natural frequency $\omega _0$.

To perform a weakly nonlinear analysis, we write

(4.4a,b)\begin{equation} \eta_{mc/s} = A_{mc/s} \cos\left(\omega_0 t + \phi_{mc/s}\right), \quad \dot{\eta}_{mc/s} ={-} \omega_0 A_{mc/s} \sin\left(\omega_0 t + \phi_{mc/s}\right), \end{equation}

where $A$ and $\phi$ are slowly varying functions of time. Their dynamics can be retrieved by a moving average of the equations over a period $T=2{\rm \pi} /\omega _0$, i.e. an oscillation period of the fast time scale defined by the uncoupled oscillators. The averaged quantities are denoted by an overbar, and are defined by

(4.5)\begin{equation} \bar{g}(t) = \langle g\rangle(t) \equiv \frac{1}{T} \int_{t}^{t+T} g(s) \,{\rm d}s. \end{equation}

After rewriting the original equation in standard form, $\ddot {\eta } + \omega _0^2\eta + h(\eta,\dot {\eta })=0$, first-order, coupled nonlinear ordinary differential equations that govern the (slow) evolution of the amplitudes and phases can be obtained (Strogatz 1994):

(4.6a,b)\begin{equation} \dot{\bar{A}} = \left\langle\frac{ h(A,\phi)}{\omega_0}\sin(\omega_0t+\phi)\right\rangle \quad \mbox{and} \quad \dot{\bar{\phi}} =\left\langle \frac{ h(A,\phi)}{\omega_0 A}\cos(\omega_0t+\phi)\right\rangle. \end{equation}

By approximating the instantaneous values on the right-hand side of (4.6a,b) with the averaged values, the evolution of the amplitudes and phases of (4.2) is given by

(4.7a)$$\begin{gather} \dot{\bar{A}}_{mc/s} = \frac{\tilde{\beta} - \alpha}{2} \bar{A}_{mc/s} - 2\tilde{\varDelta} \sin^2\left(\frac{m{\rm \pi}}{N}\right) \bar{A}_{mc/s} + 3\tilde{\kappa} \left\langle \frac{\mathcal{N}_{mc/s}}{\omega_0} \sin(\omega_0t + \bar{\phi}_{mc/s}) \right\rangle, \end{gather}$$

(4.7b)$$\begin{gather}\dot{\bar{\phi}}_{mc/s} = 2\tilde{\varGamma}\sin^2\left(\frac{m {\rm \pi}}{N}\right) + 3\tilde{\kappa}\left\langle \frac{\mathcal{N}_{mc/s}}{\omega_0A_{mc/s}} \cos(\omega_0t + \bar{\phi}_{mc/s}) \right\rangle. \end{gather}$$

In the absence of the nonlinear terms, the amplitudes of the modes grow exponentially, according to the linear term in (4.7a), with the same growth rates as those predicted by (3.11a). The phases of all modes with $m\neq 0$ drift at a constant rate, indicating a shift in frequency of $2\tilde {\varGamma }\sin ^2(m{\rm \pi} /N)$, as per the constant term in (4.7b). This frequency correction is not identical to that predicted by the linear analysis, as in (3.11b). This result is, however, correct in the limits considered for the method of averaging, in which even the linear deviations from the unperturbed oscillator, proportional to $(\tilde {\beta }-\alpha )$, $\tilde {\varGamma }$ and $\tilde {\varDelta }$, are considered small and expanded to first order. Indeed, a Taylor expansion of (3.11b) around $\omega _0$ retrieves, at first order, the same frequency shift as the method of averaging. Higher-order corrections to the frequency can be retrieved by the method of averaging by considering higher-order terms when averaging the equations. However, since the frequency drift is generally small, these higher-order effects will not be considered in this study.

5. Results and discussion

In this section we present and discuss results obtained by integrating the slow dynamics equations for the evolution of the amplitude and the phase of an $N=8$ can-annular configuration. Although the equations derived in §§ 3 and 4 are valid for an arbitrary number $N$ of cans, in order to present quantitative results, we need to fix $N$ to a specific value. In general, it should be expected that the results will depend on the specific value of $N$, since, as $N$ increases, qualitatively different types of solutions exist. The specific value $N=8$ was chosen to balance the complexity of the dynamics to the manageability of the nonlinear terms in the equations. The configuration considered here models the atmospheric can-annular experimental set-up assembled at the Norwegian University of Science and Technology and described by Buschmann, Worth & Moeck (2023). We refer to the latter for details about the geometry. For all the calculations in the following, we have set $(\tilde {\beta } - \alpha )/2 \approx 23\ {\rm s}^{-1}$, $\omega _0 \approx 2{\rm \pi} \times 180\ {\rm rad}\ {\rm s}^{-1}$ and $\tilde {\kappa } \approx 50$, and varied the strength of the (aero)acoustic coupling between the cans by varying $\tilde {\varDelta }$ and $\tilde {\varGamma }$. The value of the driving $(\tilde {\beta }- \alpha )$ term was chosen to be approximately 2 % that of $\omega _0$ because (i) the method of averaging requires $(\tilde {\beta } - \alpha ) \ll \omega _0$ to yield accurate results at first order and (ii) it was observed that the magnitude of this value is consistent with typical positive growth rates found in Orchini et al. (2022), where the same experimental set-up of Buschmann et al. (2023) was modelled. The nonlinear saturation coefficient $\kappa$ was instead tuned in order to have oscillation amplitudes with a maximum value close to unity.

5.1. Fast and slow dynamics

To validate the method of averaging presented in § 4, we integrate numerically both the fast evolution of the coupled oscillators, (4.2), and the slow, averaged evolution of the amplitudes and phases of all the linearly unstable modes, (4.7).

Figure 3 shows the evolution of the coefficients $\eta _{jm}$ with thinner lines; a fast oscillation is apparent. To emphasise the oscillating nature of these time series, the last 20 ms, containing approximately three full cycles, are stretched in each panel. With the provided initial conditions, the oscillations converge to a push–push mode ($b=0$). The decay of the amplitudes of the modes with $m=1,2$ is relatively slow and not visible in the reported data, which show only the evolution of the dynamics for a total time of 1 s. Figure 3 shows also the evolution of the modal amplitudes as predicted by the method of averaging with thicker, darker lines. The initial conditions of the averaged equations are chosen to be identical to those of the oscillator dynamics, by exploiting the relation between the oscillator states and the amplitudes and phases of the modes provided by (4.3a,b) at $t=0$ s. Despite the highly non-monotonic behaviour of the oscillation amplitudes, the evolution of the amplitudes, as approximated by the method of averaging, compares very favourably with the real amplitudes of the oscillators. The averaged solutions can, in fact, also predict the modulation of the oscillation amplitudes, which is due to beating and nonlinear modal interaction, both typical of the transient dynamics of systems with multiple, closely spaced eigenvalues.

Figure 3. Time integration of the fast (thin, lighter colours) and slow (thick, darker colours) dynamics of the considered weakly nonlinear can-annular system. Each panel shows the evolution of the amplitude of one of the $N=8$ modes that contribute to the dynamics; similar colours denote degenerate modes. Results obtained for $\tilde {\varGamma }=54$ and $\tilde {\varDelta }=0$.

To assess the description of the phase dynamics from the method of averaging, we need to extract phase information from the dynamics of the oscillators. For each real-valued time trace $\eta _{mj}(t)$, we perform a Hilbert transform to obtain the imaginary part of the complex-valued analytic signal $\mathcal {H}[\eta _{mj}](t)$, which satisfies ${\rm Re}[\mathcal {H}[\eta _{mj}](t)] = \eta _{mj}(t)$. The magnitude, $|\mathcal {H}[\eta _{mj}]|$, and the phase angle, $\varphi _{mj} \equiv \angle (\mathcal {H}[\eta _{mj}])$ of these complex oscillatory signals correspond to the instantaneous amplitude and phase of the dynamics of the oscillator $\eta _{mj}$. The phase $\varphi _{mj}$ contains information on both the natural phase of the uncoupled oscillators, defined by $\phi _0 = \omega _0 t$, and on the deviation from the latter of the coupled, nonlinear system. The method of averaging approximates the latter, and to allow for a comparison, we define

(5.1)\begin{equation} \phi_{mj} \equiv \varphi_{mj} - \phi_0 \equiv \angle(\mathcal{H}[\eta_{mj}]) - \omega_0 t. \end{equation}

Figure 4 shows a comparison between the instantaneous phases of the oscillators – solid lines, lighter colours – obtained with the Hilbert transform, and the averaged phases obtained by integrating directly the slow dynamics equations (4.7) – dashed lines, darker colours. Except for some inaccuracy introduced by the finite-length Hilbert transform at the beginning and end of the time series, the comparison between the phases is strikingly good.

Figure 4. Time integration of the fast (solid, lighter colours) and slow (dashed, darker colours) dynamics of the considered weakly nonlinear can-annular system. The time interval shown corresponds to 4 seconds. The panels show the instantaneous phase of the oscillations, defined as a deviation from the natural phase of the uncoupled oscillators, $\phi _0 = \omega _0 t$. Similar colours refer to degenerate modes. Results obtained for $\tilde {\varGamma }=54$ and $\tilde {\varDelta }=0$.

Although not shown, analogous comparisons have been performed for various values of the coupling parameters $\tilde {\varDelta }$ and $\tilde {\varGamma }$, and a good match between the averaged and exact oscillator dynamics has always been observed. These results validate the derivation and implementation of the averaged equations, and for this reason in the following sections we will only show and discuss results obtained from the averaged equations.

5.2. Initial conditions and supported solutions

To study which types of solutions the system supports, we rely on the numerical integration of the averaged equations from a variety of initial conditions. In this way, we assess whether specific solutions are supported by the equations. By fixing $\tilde {\varDelta } = 0$ and $\tilde {\varGamma } = 21$, we alternately set the initial amplitude of one of the modes to one, while setting the amplitudes of the remaining modes to $10^{-5}$. All phases are arbitrarily chosen to start at $\bar {\phi }_{mj}=0$. This particular set of initial conditions is chosen to assess the stability of standing solutions. Because we have chosen a standing basis for decomposing the pressure dynamics in the cans, a pure standing solution is found if, for example, the amplitude of one mode with $0< m< N/2$ is non-zero, and all other amplitudes are zero.

Figure 5 shows the evolution of the amplitudes and phases of all the considered initial conditions. Only the evolution of the phase difference between the modes that participate in the steady-state, stable oscillation are shown for each case. Starting from figure 5(a), initialising the oscillations to a push–push mode results in a stable push–push state (mode 0c). In the steady-state, the phase of this mode is flat, denoting that no frequency drift occurs, and the final steady-state oscillations have the same frequency of the uncoupled system, $\omega _0$, which is a general characteristic of push–push modes (von Saldern et al. 2021a,c). When initialising the oscillations to a standing mode with $m=1$, figure 5(c,d), instead, the nonlinear modal interaction triggers the participation of other modes. After a transient, the amplitudes of all modes drop to zero, except those of the two standing modes with $m=1$, which converge to the same value. In the steady state, their phases exhibit a constant drift, and are shifted apart by ${\rm \pi} /2$. In conjunction with the amplitudes reaching the same level, this indicates that the system has converged to a spinning solution since it holds that

(5.2)\begin{align} A\cos(\omega_0 t +\phi(t)) \cos(\theta) + A\cos(\omega_0 t+\phi(t)+{\rm \pi}/2) \sin(\theta) = A \cos(\omega_0 t + \phi(t) + \theta). \end{align}

The slope of the phase in the steady state reaches a value of $6.15\ {\rm rad}\ {\rm s}^{-1}$ for both modes. This compares very well with the linear correction to the unperturbed angular frequency predicted by the method of averaging, $\Delta \omega \equiv 2\tilde {\varGamma }\sin ^2(m{\rm \pi} /N)$, as per (4.7b).

Figure 5. Evolution of the amplitudes and phases of the coupled dynamics of an $N=8$ can-annular system when the oscillations are initialised to (a) a push–push mode, (b) a push–pull mode, (c,d) a standing mode with $m=1$, (e,f) a standing mode with $m=2$, (g,h) a standing mode with $m=3$. The coupling parameters are set to $\tilde {\varGamma }=21$ and $\tilde {\varDelta }=0$.

Similar results hold for cases in which standing solutions with $m=2$, figure 5(e,f), and $m=3$, figure 5(g,h), were initialised. In all these cases, after a transient, the solutions diverge from the initial standing state and converge to spinning modes with order $m$ corresponding to that of the initial standing mode. This result is consistent with the findings of Noiray et al. (2011), where it was shown that standing oscillations are not stable solutions in symmetric annular configurations saturated by a cubic nonlinearity. The angular frequency drifts calculated from the slopes of the phases are $\Delta \omega =21\ {\rm rad}\ {\rm s}^{-1}$ and $\Delta \omega = 35.8\ {\rm rad}\ {\rm s}^{-1}$, for oscillations of spinning modes with $m=2$ and $m=3$, respectively. Once again, these values coincide with those predicted by the linear term in (4.7b).

5.3. Synchronisation

A different behaviour is observed when the solutions are initialised to a push–pull mode (figure 5b). In this case, multiple modes with different values of $m$ are triggered by the nonlinear modal interaction, and survive in the final, steady-state oscillations. In particular, the $m=4$ push–pull mode, both modes with $m=3$ and one of the modes with $m=2$ participate in the finite-amplitude oscillation. Figure 5 shows only a total integration time of 4 s; at larger times the amplitudes of both standing modes with $m=3$ converge to the same amplitude, $A_{3j} = 0.67$, but with a phase difference of zero. This corresponds to a standing mode, but rotated by ${\rm \pi} /4$ with respect to those used as a basis, and with an amplitude of $\sqrt {2}A_3$. The amplitudes of the modes with $m=2$ and $m=4$ converge to $A_{2s} = 0.24$, $A_{4c} = 0.88$, respectively, and their phases have a constant phase difference with respect to the modes having $m=3$. This is a signature of synchronisation. By synchronisation we mean the adjustment of rhythms of oscillators due to their weak interaction, as defined by Pikovsky et al. (2001). As discussed in § 1.2, two conditions must be satisfied for synchronisation to be possible: (i) the oscillators must be self-sustained when decoupled; and (ii) a weak coupling must be established between the oscillators. By using the results derived in § 3, we can verify that these conditions are easily satisfied for can-annular geometries. When the cans are decoupled, $\tilde {\varGamma } = \tilde {\varDelta } = 0$, the eigenvalues of all the $N$ modes converge to (3.12). Thus, provided that the isolated can system is thermoacoustically unstable, $(\tilde {\beta } - \alpha )>0$, condition (i) is satisfied because there exist a set of $N$ self-sustained oscillation patterns in the system, one in each can. A (weak) coupling between the cans is present when $\tilde {\varGamma } \neq 0$ and/or $\tilde {\varDelta } \neq 0$, satisfying condition (ii). Furthermore, the linear effect of the weak coupling can be quantitatively seen from (3.11). All the growth rates $\sigma _b$ and angular frequencies $\omega _b$ are close to the frequency of an isolated can. If the growth rate of an uncoupled thermoacoustic mode is positive, $(\tilde {\beta }-\alpha )/2 > 0$, then one should expect all the eigenvalues of the coupled can-annular system to have a positive growth rate, because they form a cluster. All these modes grow exponentially in the linear regime and interact with each other during the transient response. The interaction between modes may suppress the growth of some modes and enhance the growth of others. Furthermore, the nonlinear interaction causes a drift of the oscillation frequencies as the amplitude of oscillation increases. Under appropriate conditions, these frequencies may lock to the same value: if modes with different linear oscillation frequencies persist in the saturated oscillation state with the same frequency, synchronisation has occurred.

Identifying explicit conditions for synchronisation is a non-trivial task, because one needs to find synchronised states that are supported by the governing equations (4.7), and verify that these solutions are stable. Already for $N=8$, this problem is complex to treat analytically because the phase space is 16-dimensional, and there exists a very large number of possible synchronised states. It is, however, analytically treatable for $N=3$. The expressions of the governing equations (4.7) including all the nonlinear terms can be explicitly reported for this case:

(5.3a)\begin{align} \dot{\bar{A}}_{0c} &= \frac{\tilde{\beta} - \alpha}{2} \bar{A}_{0c} \nonumber\\ &\quad - \frac{3\tilde{\kappa}}{16} \bar{A}_{c0}\left[\bar{A}_{c1}^2 \cos(2\bar{\phi}_{c0}-2\bar{\phi}_{c1})+\bar{A}_{s1}^2 \cos(2\bar{\phi}_{c0}-2\bar{\phi}_{s1})+2\bar{A}_{c0}^2+2\bar{A}_{c1}^2+2\bar{A}_{s1}^2\right], \end{align}

(5.3b)\begin{align} \dot{\bar{A}}_{1c} &= \frac{\tilde{\beta} - \alpha}{2} \bar{A}_{1c} - 3\tilde{\Delta} \bar{A}_{1c}\nonumber\\ & \quad- \frac{3\tilde{\kappa}}{32}\bar{A}_{c1}\left[4\bar{A}_{c0}^2\cos(2\bar{\phi}_{c0} - 2\bar{\phi}_{c1}) + \bar{A}_{s1}^2\cos(2\bar{\phi}_{c1} - 2\bar{\phi}_{s1}) + 8\bar{A}_{c0}^2 + 3\bar{A}_{c1}^2 + 2\bar{A}_{s1}^2\right], \end{align}

(5.3c)\begin{align} \dot{\bar{A}}_{1s} &= \frac{\tilde{\beta} - \alpha}{2} \bar{A}_{1s} - 3\tilde{\Delta}\bar{A}_{1s} \nonumber\\ &\quad - \frac{3\tilde{\kappa}}{32} \bar{A}_{s1}\left[4\bar{A}_{c0}^2\cos(2\bar{\phi}_{c0} - 2\bar{\phi}_{s1}) + \bar{A}_{c1}^2\cos(2\bar{\phi}_{c1} - 2\bar{\phi}_{s1}) + 8 \bar{A}_{c0}^2 + 2 \bar{A}_{c1}^2 + 3 \bar{A}_{s1}^2\right], \end{align}

(5.3d)$$\begin{gather} \dot{\bar{\phi}}_{0c} ={+} \frac{3\tilde{\kappa}}{16} \left[\bar{A}_{c1}^2\sin(2\bar{\phi}_{c0} - 2 \bar{\phi}_{c1}) + \bar{A}_{s1}^2 \sin(2\bar{\phi}_{c0} - 2\bar{\phi}_{s1}) \right], \end{gather}$$

(5.3e)$$\begin{gather}\dot{\bar{\phi}}_{1c} = 3\tilde{\varGamma}+ \frac{3\tilde{\kappa}}{8} \left[ \frac{\bar{A}_{s1}^2}{4} \sin(2\bar{\phi}_{c1} - 2\bar{\phi}_{s1}) - \bar{A}_{c0}^2\sin(2\bar{\phi}_{c0} - 2\bar{\phi}_{c1})\right], \end{gather}$$

(5.3f)$$\begin{gather}\dot{\bar{\phi}}_{1s} = 3\tilde{\varGamma}-\frac{ 3\tilde{\kappa}}{8} \left[ \bar{A}_{c0}^2\sin(2\bar{\phi}_{c0} - 2\bar{\phi}_{s1}) + \frac{\bar{A}_{c1}^2}{4}\sin(2\bar{\phi}_{c1} - 2\bar{\phi}_{s1})\right]. \end{gather}$$

One can verify that these equations are equivalent to equation (5) presented by Moeck et al. (2019), which were derived to model synchronisation between an axisymmetric and an azimuthal mode in an annular combustor. We note that the physical meaning of the coupling terms is different: the reactive ($\hat {\varDelta }$) and resistive ($1-\hat {\sigma }$) coupling terms in the work of Moeck et al. (2019) are directly related to differences in the frequencies/growth rates of longitudinal and azimuthal modes in an annular combustion chamber. On the other hand, the reactive ($\tilde {\varGamma }$) and resistive ($\tilde {\varDelta }$) terms in the can-annular model presented in this work stem from the Rayleigh conductivity that provides a weak coupling between otherwise identical cans. Although the geometry and the detailed interpretation of the parameters are different, the physics underlying the thermoacoustic phenomenon is the same, and it should not come as a surprise that equations with the same structure appear here.

Given the equivalent representation of the dynamics, we can directly apply the results of Moeck et al. (2019) about synchronisation, shown in their figure 3, to the can-annular system considered here: for $N=3$, synchronisation is possible only for specific combinations of the reactive and resistive coefficients. In particular, we can exploit their result to conclude that, for $N=3$, $\tilde {\varDelta }$ must be negative for synchronisation to occur – the resistance $\tilde {\varDelta }$ is proportional to the parameter ($1-\hat {\sigma }$) of Moeck et al. (2019). An interesting question is whether a negative resistance is a strict requirement for synchronisation in can-annular systems with a larger number of cans, say $N=8$. The extension of this result to a larger number of cans could not formally be proven. Instead, we relied on extensive numerical calculations for the $N=8$ case to investigate when synchronisation occurs. The results suggest that our model supports synchronised solutions only for negative values of the resistivity. It was reported in Pedergnana et al. (2021) that such a negative resistance may be induced in a frequency band by aeroacoustic effects that excite the shear layer formed by the grazing flow along the can-to-can apertures.

For our calculations, we consider coupling parameters in the range $\tilde {\varDelta }/(\tilde {\beta }- \alpha ) \in [-0.5,0.5]$ for the resistance and $\tilde {\varGamma }/(\tilde {\beta } - \alpha ) \in [-0.5, 0.25]$ for the reactance. Starting from the smallest values of both coupling parameters, we initialise (4.7) close to the unstable fixed points, with amplitudes equal to $10^{-5}$ and phases randomly distributed in the range $[0,2{\rm \pi} ]$. We integrate the equations until $t=200$ s, to ensure convergence of the results. We then assess if the final state is synchronised by determining how many modes contribute to the final steady-state oscillations (with amplitudes larger than 0.05), and by comparing the slopes of the phases – in other words, the frequencies. Each subsequent calculation is then initialised from the end state of a neighbouring state, following a continuation-like approach. Because the first identified state is synchronised, the approach we follow is able to identify when the stability of synchronised solutions is lost. Figure 6 shows the results by using on the left a colour scheme that assigns a unique colour to states with $n$ synchronised modes, where $n=0$ means that no synchronisation occurs, and on the right the same colour scheme as in figure 3 to show which mode has the largest amplitude in the final, possibly synchronised state. It is observed that synchronised states are stable when $\tilde {\varDelta }<0$, consistent with the analytical findings of Moeck et al. (2019) for $N=3$. However, in contrast to the latter study, the window of synchronisation is much larger. This is likely because in that study only one synchronisation pattern was possible, whereas in the can-annular configuration considered in the present work, a large number of possible synchronised states exist. Furthermore, in the identified synchronised states the push–pull (${4c}$) mode has the largest amplitude – an example can be seen in figure 5(b). On the other hand, when synchronisation does not happen, for most $\tilde {\varDelta }>0$, the push–push ($0c$) mode dominates the dynamics. We note that analogous results were reported in Pedergnana & Noiray (2022), where dominant push–push oscillations for positive values of the resistive coupling and dominant push–pull oscillations for negative values of the resistive coupling were reported. We, however, also observe that there exist regions of positive resistance where synchronisation persists. We, however, emphasise that these solutions have been obtained by tracking synchronised solutions from negative resistance to positive resistance and that figure 6 provides no information on the basins of attraction of these synchronised solutions, nor on the possible existence of other stable oscillatory states. We also highlight that the synchronised map shown in figure 6 exhibits a clear symmetry around the $\tilde {\varGamma } = 0$ axis and that, provided that $\tilde {\varDelta }<0$, the larger is the magnitude of $\tilde {\varGamma }$ the fewer number of modes participate in a synchronised state oscillation.

Figure 6. Synchronisation map for a can-annular system with $N=8$ cans. (a) Number of modes that participate in the synchronised final state. (b) Label of the mode that dominates the oscillations. The results were obtained by time-marching equations (4.7), and by using a continuation algorithm while slowly increasing $\tilde {\varDelta }$ and $\tilde {\varGamma }$. Numerical results indicate that the synchronisation of states with different azimuthal orders is likely when $\tilde {\varDelta }<0$.

The results shown here correspond to an eight-can system. As pointed out previously, existence and stability of oscillatory solutions will depend on the number of cans $N$. It is difficult to make any general statement on how the dynamical picture changes when $N$ is increased or decreased. One obvious outcome, however, is that for odd $N$, the push–pull mode does not exist. Since the push–pull mode was often observed to be dominant in the synchronised states for our $N=8$ computations, one can expect that the modal dynamics would be quite different for $N=7$ or 9.

To provide more insight into what a synchronised state is and what it is not, we report in figure 7 details about three specific calculations for $\tilde {\varGamma } = 21$, for which the phases were initialised to zero. When $\tilde {\varDelta } = 0$, case (a), after a transient only the amplitudes of the $\cos (\theta )$ and $\sin (\theta )$ structures remain non-zero, whereas the others decay. The $\cos (\theta )$ and $\sin (\theta )$ structures oscillate with the same frequency; however, this is not a synchronised state because these modes are degenerate by construction, and have the same oscillation frequency and growth rates already in the linear regime – see (3.11) for positive/negative $b$. When $\tilde {\varDelta }= 5$, the steady-state dynamics is dominated by the $b=0$ mode only. Because a single global mode participates in the oscillations, synchronisation cannot occur. We note, however, that a $b=0$ mode corresponds to a push–push oscillation. Some authors refer to this state also as a synchronised state, because the phases of the oscillations in all the cans are aligned (Moon et al. 2020b; Pedergnana & Noiray 2022). This is not consistent with our classification of synchronisation, because the alignment of the oscillation phases between the cans for the $b=0$ mode is a linear effect, as we discussed in § 3, whereas synchronisation between modes is an inherently nonlinear phenomenon (Pikovsky et al. 2001; Strogatz 2004). Lastly, when $\tilde {\varDelta } = -4$, four modes participate in the oscillation, $\sin (2\theta )$, $\cos (3\theta )$, $\sin (3\theta )$ and $\cos (4\theta )$. Despite the fact that the linear frequencies of these modes are different, nonlinear effects drift these frequencies and synchronise them to the same value. That in the steady state all modes participating in the dynamics have synchronised to the same oscillation frequency can be seen from the fact that the slopes of the averaged phase drifts are identical for all modes.

Figure 7. Evolution of the averaged amplitudes and phases for $\tilde {\varGamma } = 21$: (a) $\tilde {\varDelta }=0$, converging to a spinning solution; (b) $\tilde {\varDelta }=5$, converging to a push–push solution; (c) $\tilde {\varDelta }=-4$, converging to a synchronised solution. Animated versions of these graphs, which make it visually easy to understand the dynamics, are provided as supplementary material available at https://doi.org/10.1017/jfm.2024.4.

It is interesting to consider these results in view of available experimental observations. However, a detailed quantitative comparison is not possible for several reasons. Most of the experimental studies to date investigate two-can or four-can systems (Moon et al. 2020a,b; Chen et al. 2021; Guan et al. 2022), which host a smaller number of Bloch modes. The laboratory-scale can-annular system investigated by Buschmann et al. (2023) corresponds to the configuration considered in the present work, but a quantitative comparison would require detailed information about the linear and nonlinear flame responses, which is not available from the experiment. Lastly, most of the can-annular systems studied experimentally are turbulent, which adds significant noise to the thermoacoustic dynamics and partly obscures the attractor. Consequently, only general, qualitative aspects of our results can be compared with the various experimental observations. Eigenvalue clusters comprising all Bloch numbers and organised around the longitudinal can modes are a general feature of can-annular systems, as we have argued in § 3. Traces of these eigenvalue clusters can be seen in pressure spectra as multiple close acoustic resonances (Moon et al. 2020b; Buschmann et al. 2023). Furthermore, while in single-burner and annular multiburner systems, combustion instabilities most commonly manifest as single-mode oscillations, our calculations show that multimode oscillations routinely occur in a can-annular system, with a push–push or push–pull mode often being dominant. This is consistent with experimental observations by Moon et al. (2020b), Guan et al. (2022) and Buschmann et al. (2023).

6. Conclusions

In this work we have studied linear and nonlinear characteristics of can-annular combustors. Can-annular (and annular) combustors are the most common combustor architectures in stationary gas turbines for power generation. It is therefore important to understand the complexity of thermoacoustic phenomena that occur in these systems. In contrast to single-flame combustors, which are studied predominantly in laboratory settings, multiflame systems exhibit more complex dynamics. And while annular and can-annular systems share some qualitative similarities, foremost a high degree of discrete rotational symmetry, there are also some marked differences. In the latter type of combustor, each flame is situated in an almost isolated acoustic environment and communicates with neighbouring cans only weakly through a small gap at the turbine inlet. The weak coupling in conjunction with the combustor cans all being nominally identical gives rise to mode clusters, which spawn from the modes of an isolated can and spread farther with increasing coupling.

Clusters of unstable thermoacoustic modes, hence, generically occur in can-annular combustors. With a group of eigenvalues having close frequencies and similar growth rates, enhanced modal interaction is expected and synchronisation is indeed possible. Since these are nonlinear phenomena, predictive modelling is challenging and cannot be accomplished by relying, for example, on the describing function technique. While the latter has been proven to be a successful tool for single-mode limit cycles in thermoacoustic (and, of course, other) systems (Dowling 1997; Noiray et al. 2008), it is not able to properly capture nonlinear multimodal interaction (Moeck & Paschereit 2012; Orchini & Juniper 2016). In the presence of multiple linearly unstable modes with comparable growth rates, a proper nonlinear framework has to be employed to successfully predict the final oscillatory state the system evolves to.

Starting from the linearised Euler equations with a nonlinear flame model, we established a weakly nonlinear framework capable of representing the aforementioned nonlinear phenomena in generic can-annular combustors. Wave-like equations for the acoustic pressure field in the individual cans were derived, taking into account the communication with the neighbouring cans through explicit coupling terms, and the feedback from the flame response. Galerkin projection was then used to reduce the infinite-dimensional problem to a set of coupled nonlinear oscillator equations, which represent the temporal evolution of the thermoacoustic system modes. Application of the method of averaging further reduced this set of equations to a slow-time state space spanned by the oscillator amplitudes and phases.

We applied this framework to represent a laboratory can-annular configuration with eight combustor cans coupled at the downstream end. The aeroacoustic coupling between neighbouring cans takes into account reactive and resistive effects. With a set of flame model and coupling parameters, a number of stable single-mode limit cycles exist and can be accessed in numerical simulations of the model by providing different (single-mode) initial conditions. When initialising the model with a push–pull mode, in other words, a mode with azimuthal order equal to half the number of combustor cans, a multimode periodic solution is established that involves several modes at finite amplitude. All the participating modes oscillate at the same frequency and feature fixed phase relations. This is, hence, an occurrence of synchronisation, which has previously been observed in other scenarios and contexts in thermoacoustic systems (Bonciolini & Noiray 2019; Moeck et al. 2019; Moon et al. 2020a; Chen et al. 2021; Guan et al. 2021). In the present case, synchronisation takes place between thermoacoustic system modes and is, consequently, purely nonlinear. Considering the effect of the can-to-can coupling, we found that the resistance essentially has to be negative for synchronised solutions to exist (i.e. the coupling has to be amplifying). Depending on both reactive and resistive coupling components, a different number of modes can synchronise.

The solutions of the averaged system with those of the original oscillator equations agree remarkably well, as evident in figure 3. Even non-trivial amplitude development during initial transients is captured with high accuracy. This highlights the suitability of this method even in the presence of nonlinear multimode interaction. The method of averaging has been used before for weakly nonlinear analysis of thermoacoustic systems, but, to the best of our knowledge, not in settings as complex as in the present case, with a large number of active modes, all having similar frequencies and growth rates. We utilised a simple cubic nonlinearity in this work to model the saturation of the flame response. While this model agrees with experimentally observed behaviour to some extent (Noiray et al. 2011), the sensitivity of our observations to changes in the flame nonlinearity would be interesting to assess.