2.1. The MICCA annular combustor and injectors’ characteristics

The laboratory-scale annular combustor MICCA, shown in figure 1(a), is used to investigate thermoacoustic instabilities coupled by azimuthal modes. The combustion chamber is formed by two cylindrical quartz walls of height $l=400\,\textrm{mm}$ , of outer diameter 300 mm for the inner quartz and inner diameter 400 mm for the outer quartz. The backplane comprises 16 regularly spaced injection units delivering liquid heptane in the form of a hollow cone spray of droplets. The air flow rate is controlled with two Bronkhorst EL-FLOW flow meters and the fuel flow rate with a Bronkhorst CORI-FLOW controller. Eight Brüel & Kjær microphones are mounted on waveguides and plugged on the chamber backplane to record the pressure fluctuations at a sampling rate $f_s=32\,768\,\textrm{Hz}$ . An array of 8 photomultipliers equipped with an OH* filter centred at 310 nm records the light emitted by eight adjacent flames (see figure 1 b). A mask is placed in front of each flame to ensure that each PM only records the light emitted by one flame. A cylindrical mask is also placed inside the inner cylindrical quartz to hide the flames in the background, as can be seen in figure 1(c), showing the MICCA combustor under operation.

Table 1. Injectors’ characteristics: swirl number ( $S_N$ ), pressure drop ( $\Delta p$ ) and pressure drop coefficient ( $\sigma$ ), obtained from measurements under cold flow conditions for an air mass flow rate of 2.3 g s $^{-1}$ (Vignat et al. Reference Vignat, Rajendram Soundararajan, Durox, Vié, Renaud and Candel2021).

Figure 1. (a) The MICCA annular combustor with an array of eight photomultipliers. (b) Microphones (labelled ‘MX’) and photomultiplier positions (labelled ‘PMX’). (c) View of the MICCA combustor under operation.

Figure 2. Exploded view of an injection unit, showing its main components (a). A top view of the swirler appears in (b).

The injectors, for which an exploded view is presented in figure 2(a), comprise four main elements: an air distributor, an atomizer, a swirler and a terminal plate. Changing the swirler (tangential channels’ radius and orientation figure 2 b) modifies the pressure drop and the swirl number of the units. Two types of injectors are used in this work: a low-swirl low-pressure drop and a high-swirl high-pressure drop, designated respectively as ‘injector S’ and ‘injector U’. The characteristics of the two injectors are gathered in table 1. Additional information on the dynamics of the flames formed by these two injectors is given by Rajendram Soundararajan et al. (Reference Rajendram Soundararajan, Durox, Renaud, Vignat and Candel2022b ) and Latour et al. (Reference Latour, Durox, Renaud and Candel2024a ).

In the present study, the annular combustor is operated at a thermal power $\mathcal {P}=118\,\textrm{kW}$ and a global equivalence ratio $\phi =0.9$ . Under these operating conditions, MICCA is stable when equipped with S-injectors and unstable when 16 U-injectors are mounted on the backplane.

2.2. Modal structure and oscillation frequency

An acoustic analysis of the MICCA combustor has to be carried out to identify the eigenmodes susceptible to being involved in the combustion/acoustics coupling. As discussed in Latour et al. (Reference Latour, Durox, Renaud and Candel2024a ), the combustion chamber in MICCA is decoupled from the plenum by the injectors that are weakly transparent to acoustic waves and introduce an important area change between the chamber and the plenum. The coupling of the injector ports’ acoustics with the combustion chamber acoustics can also be neglected, as shown in the acoustic analysis of MICCA presented in the supplementary material avilable at https://doi.org/10.1017/jfm.2025.10. One can hence assume that the modes that need to be considered are those of the chamber, where the backplane can be assimilated to a rigid wall and the outlet can be modelled as being open to the atmosphere. Experiments indicate that the first azimuthal–first longitudinal (1A1L) mode is involved in the combustion/acoustics coupling in MICCA (Latour et al. Reference Latour, Durox, Renaud and Candel2024a ). Neglecting the radial dependence, the pressure field can then be cast in the form

(2.1) \begin{equation} p(x, \theta ,t)_{1A1L} = [A_+\exp (i \theta -i \omega t +i\phi _+) +A_- \exp (-i \theta -i \omega t+i\phi _-)] \psi _{1L}(x) ,\end{equation}

where $A_+ \exp (i\phi _+)$ and $A_- \exp (i\phi _-)$ correspond to the complex amplitudes of the CCW and CW spinning waves respectively, $\theta$ designates the azimuthal coordinate considered positive in the CCW direction and $\omega$ is the angular frequency. In the previous expression, $\psi _{1L}(x)=\cos [\pi x/(2l^{\prime})]$ is the axial wave function satisfying the boundary conditions on the chamber backplane and at its exhaust, with $l^{\prime}=l+\delta _a$ , $l$ being the chamber length and $\delta _a$ the end correction. For MICCA, pressure measurements near the chamber exhaust indicate that $\delta _a \simeq 0.044\,\textrm{m}$ so that $l^{\prime}\simeq 0.44\,\textrm{m}$ (Laera et al. Reference Laera, Prieur, Durox, Schuller, Camporeale and Candel2017).

The eigenfrequency corresponding to the 1A1L mode is

(2.2) \begin{equation} f_{1A1L}= \left[ \left(\frac {c}{\mathcal {P}_a}\right)^2+ \left(\frac {c}{4l^{\prime}}\right)^2 \right]^{1/2} ,\end{equation}

where $\mathcal {P}_a= 2 \pi R$ is the mean perimeter of the system and $c$ the speed of sound. In the MICCA experiments, the perimeter is $\mathcal {P}_a \simeq 1.1\,\textrm{m}$ . Assuming an average temperature in the chamber of $T\simeq 1400\,\textrm{K}$ (estimated from exhaust gas temperature measurements carried out by Vignat (Reference Vignat2020) with a thermocouple in the single-sector counterpart of the MICCA annular set-up), defining a speed of sound $c=754\,\textrm{m}\,\textrm{s}^-{^1}$ , the eigenfrequency of the 1A1L mode is $f_{1A1L}\simeq 808\,\textrm{Hz}$ .

The complex wave amplitudes $A_+ \exp (i\phi _+)$ and $A_- \exp (i\phi _-)$ may be retrieved from the eight microphone signals by solving an over-determined system of linear equations using a least squares algorithm. One may then deduce the instability amplitude, defined as

(2.3) \begin{equation} A = \big[{A_+^2+A_-^2 }\big]^{1/2} .\end{equation}

It is worth noting that $A$ is equal to the quaternion formulation amplitude divided by $\sqrt 2$ . Finally, the frequency of the instability, $f$ , is experimentally determined from the power spectral density of the pressure signals, calculated using Welch’s method applied to 63 blocks of 8192 samples and a Hamming window weighting with a 50 % overlap. The frequency resolution for these conditions is $\Delta f = 4\,\textrm{Hz}$ .

At this point it is worth recalling that the frequency of oscillation of a thermoacoustic system (‘closed-loop’ case) is close to the eigenfrequency of the mode (the ‘open-loop’ frequency) but does not coincide with it because of the shift introduced by the presence of the flames. This shift depends on the flame dynamical characteristics (FDF) and on the set of parameters that also govern the growth rate (Schuller et al. Reference Schuller, Poinsot and Candel2020). One may show (see Schuller et al. Reference Schuller, Poinsot and Candel2020 or § 4) that the relative frequency shift, $\Delta \omega / \omega _0$ , is linked to the growth rate, $\omega _i$ , by $\Delta \omega / \omega _0 \simeq - (\omega _i / \omega _0) \tan \varphi _p$ , where $\varphi _p$ represents the phase between the pressure and the HRR signals. Since $\omega _i / \omega _0\lt 1$ and $ \varphi _p \simeq 0$ in an unstable situation, one deduces that $|\Delta \omega / \omega _0| \lt \lt 1$ . The previous argument indicates why the shift in frequency is in most cases below 5 % of the modal eigenfrequency, as can be seen in the data compiled by Ghirardo et al. (Reference Ghirardo, Juniper and Bothien2018) or in figures showing the evolution of thermoacoustic systems in the frequency/growth rate plane reported by Noiray et al. (Reference Noiray, Durox, Schuller and Candel2008), Paliès et al. (Reference Paliès, Durox, Schuller and Candel2011) and Laera et al. (Reference Laera, Prieur, Durox, Schuller, Camporeale and Candel2017) and confirmed by experimental results obtained in MICCA (Rajendram Soundararajan et al. Reference Rajendram Soundararajan, Durox, Renaud, Vignat and Candel2022b ; Latour et al. Reference Latour, Durox, Renaud and Candel2024a ).

As pointed out by one reviewer, other aspects of the system, like the temperature distribution linked to the presence of the flames, may also intervene, and affect the modal structure of the azimuthal mode, a phenomenon not accounted for in (2.2), used to determine the frequency of the 1A1L mode of the combustion chamber in MICCA. The temperature distribution can, for example, be easily taken into account by making use of a Helmholtz solver to obtain the modal distributions and eigenfrequencies, as done in previous works (Bourgouin et al. Reference Bourgouin, Durox, Moeck, Schuller and Candel2013; Laera et al. Reference Laera, Prieur, Durox, Schuller, Camporeale and Candel2017). It was found that the obtained eigenfrequencies were not very different from those calculated using theoretical expressions if the mean speed of sound $c$ , and hence the temperature in the combustion chamber, are suitably determined, for example, from exhaust gas temperature measurements, as proposed in this work from results reported by Vignat (Reference Vignat2020). The shift in frequency being estimated to be less than a few per cent of the open-loop frequency (obtained from self-sustained oscillation measurements in MICCA), and the calculated value of the open-loop frequency (808 Hz) from the burnt gas temperature estimate falling in the range of frequencies corresponding to observed oscillations (around 800 Hz), one can have confidence in the estimated value. It is, however, interesting to determine the effect of an error in the estimation of the temperature on the calculated frequency of oscillation. To that end, one may consider, for instance, an error $\Delta T= 100\,\textrm{K}$ . The variation in modal frequency is then given by $\Delta f/f \simeq \Delta c/c = (1/2) \Delta T/ T$ . Taking $T= 1400\,\textrm{K}$ , one finds that $\Delta f/ f \simeq 3.5\, \%$ and for a frequency $f=800\,\textrm{Hz}$ this would induce a variation $\Delta f \simeq 28\,\textrm{Hz}$ .

4.1. Polynomial heat release rate fluctuation formulation: ${\dot Q}^{\prime}= \beta ^* p- \kappa _1^* p^3$

The starting point is the non-homogeneous wave equation, derived from the combination of the mass, momentum and energy equations for low-Mach flows and perfect gases, to which a damping term is added in the form $\alpha ({\partial p}/{\partial t})$ to account for the losses in the system

(4.1) \begin{equation} \frac {\partial ^2 p}{\partial t^2}-\rho c^2 \boldsymbol {\nabla } \boldsymbol {\cdot } \frac {1}{\rho }\boldsymbol {\nabla } p =(\gamma -1) \frac {\partial {\dot q}^{\prime}}{\partial t} -\alpha \frac {\partial p}{\partial t} .\end{equation}

In this expression, ${\dot q}^{\prime}$ , $\rho$ , $c$ and $\gamma$ respectively designate the volumetric HRR fluctuations, the mean density, the mean speed of sound and the specific heat ratio. The steps leading to the derivation of (4.1) are recalled, for instance, in Nicoud et al. (Reference Nicoud, Benoit, Sensiau and Poinsot2007, Reference Noiray, Bothien and Schuermans2011) or Ghirardo et al. (Reference Ghirardo, Boudy and Bothien2018).

At this stage, it is interesting to discuss the zero-Mach-number assumption used to write (4.1), which is common to most thermoacoustic investigations of annular systems (Nicoud et al. Reference Nicoud, Benoit, Sensiau and Poinsot2007, Reference Noiray, Bothien and Schuermans2011; Ghirardo et al. Reference Ghirardo, Juniper and Moeck2016). In annular combustors that are typically used in gas turbines, the Mach number is low to ensure that the residence time of reactants is sufficient for flame stabilization, to allow complete conversion into products and minimize pressure losses associated with heat addition. The analysis by Nicoud et al. (Reference Nicoud, Benoit, Sensiau and Poinsot2007) indicates that the mean flow terms can be neglected if the characteristic Mach number is small compared with the ratio of the flame dimension to the typical acoustic wavelength and this condition is generally fulfilled. There are, however, some more subtle effects of the presence of a mean swirling flow on the acoustic modes in an annular system. If, for example, the azimuthal velocity is high enough, the global rotation of the flow will suppress the degeneracy of the azimuthal modes and the CW and CCW modes will feature different eigenfrequencies and growth rates (Bauerheim, Cazalens & Poinsot Reference Bauerheim, Cazalens and Poinsot2015). It is also indicated by Faure-Beaulieu et al. (Reference Faure-Beaulieu, Pedergnana and Noiray2023) that, when a swirl is imposed in a certain direction, a statistical preference is observed toward mixed states propagating against the swirl. However, in the situation at hand, the global rotation velocity is so small that the two eigenfrequencies cannot be distinguished in the spectral analysis of the microphone signals and the statistical preference is expected to be small compared with the symmetry-breaking effects induced by injector staging of the kind investigated in the present work.

A second assumption made in writing (4.1) is that the acoustics may be treated as a linear process. This is also widely adopted because the levels of relative pressure fluctuations in typical gas turbine combustion systems remain below a few per cent of the chamber pressure, in contrast with rocket thrust chambers where relative fluctuation levels may reach up to 40 % of the chamber pressure. The reader can find further details on the nonlinear acoustics in rocket engines in Culick (Reference Culick1994). In the case of gas turbines, the level of oscillation does not exceed a few per cent. In the MICCA experimental set-up, even at a pressure fluctuations level of 5000 Pa (5 % of the chamber pressure), the microphone signals remain sinusoidal while the light intensity signal from OH $^*$ radicals detected by photomultipliers facing the flames become highly nonlinear (Prieur et al. Reference Prieur, Durox, Schuller and Candel2018). One may then safely assume that the acoustics is linear and that the main source of nonlinearity in the system is linked to the unsteady HRRs and that this nonlinearity can be represented with a FDF (Dowling Reference Dowling1997; Noiray et al. Reference Noiray, Durox, Schuller and Candel2008).

Finally, it is also worth noting that the $\alpha$ appearing in this expression corresponds to the damping rate of acoustic energy and is equal to twice the damping rate of pressure. The value of this damping rate can be obtained in various ways using, for instance, system identification methods, as proposed in Boujo et al. (Reference Boujo, Denisov, Schuermans and Noiray2016). Another possibility, used in this study (see § 5), relies on source term measurements at the limit cycle, as exemplified in Durox et al. (Reference Durox, Schuller, Noiray, Birbaud and Candel2009) or Latour et al. (Reference Latour, Durox, Renaud and Candel2024b ), which will be discussed, along with the modelling hypothesis, in § 5.

In flames that are nearly isobaric, the product $\rho c^2$ is nearly constant and may be introduced in the divergence operator in (4.1). One may also assume that the flame is compact with respect to the acoustic wavelength and consider that the HRR fluctuations are concentrated at one point, $\boldsymbol {x}_f$ . One may then write ${\dot q}^{\prime}= \delta (\boldsymbol {x}-\boldsymbol {x}_f){\dot Q}^{\prime}$ , where $\delta$ is the Dirac function and ${\dot Q}^{\prime}$ corresponds to the HRR fluctuations integrated over the volume of the flame. Hence, the wave equation may now be replaced by

(4.2) \begin{equation} \frac {\partial ^2 p}{\partial t^2}-\boldsymbol {\nabla } \boldsymbol {\cdot } c^2 \boldsymbol {\nabla } p =(\gamma -1) \delta (\boldsymbol {x}-\boldsymbol {x}_f) \frac {\partial {\dot Q}^{\prime}}{\partial t} -\alpha \frac {\partial p}{\partial t} .\end{equation}

The normal modes of this equation, in the absence of heat release and damping, are such that

(4.3) \begin{equation} \boldsymbol {\nabla } \boldsymbol {\cdot } c^2 \boldsymbol {\nabla } \psi _n + \omega _n^2 \psi _n=0 ,\end{equation}

where $\omega _n$ and $\psi _n$ represent the modal eigenvalues and eigenfunctions, respectively. The pressure field can then be expanded on the orthogonal basis formed by the eigenmodes $\psi _n$ , solutions of the homogeneous wave equation

(4.4) \begin{equation} p=\sum _n \eta _n(t) \psi _n(\boldsymbol {x}) ,\end{equation}

where $\eta _n$ represents the amplitude of the nth mode. For homogeneous boundary conditions, the normal modes $\psi _n$ are orthogonal (see for example Nicoud et al. Reference Nicoud, Benoit, Sensiau and Poinsot2007) and one may write

(4.5) \begin{equation} \int _V \psi _n \psi _m \textrm{d}V= \Lambda _n \delta _{mn} \, \, \textrm {with} \,\, \delta _{nn}= 1 \,\, \textrm { and } \,\, \delta _{mn}= 0 \,\, \textrm {when} \,\, m \ne n .\end{equation}

Introducing the modal expansion (4.4) in the wave equation (4.2) and projecting the result on the normal modes, one obtains a set of differential equations for the modal amplitudes $\eta _n(t)$

(4.6) \begin{equation} {\ddot \eta _n}+ \alpha {\dot \eta _n} + \omega _n^2 \eta _n = \frac {\gamma -1}{\Lambda _n} \frac {\partial {\dot Q}^{\prime}}{\partial t} \psi _n(x_f) \,\, \textrm {for} \,\, n=1,2\ldots .\end{equation}

For the sake of simplicity, it will be assumed in what follows that a single mode, characterized by an eigenfrequency $\omega _0$ and a modal amplitude $\eta$ , is involved in the combustion/acoustics coupling. This situation physically corresponds, for instance, to a coupling by a longitudinal mode in a single-injector set-up. This assumption is made to analyse the effects of the HRR formulation in the simplest possible framework. In this case, the pressure field reads

(4.7) \begin{equation} p=\eta (t) \psi (\boldsymbol {x}) .\end{equation}

One possibility, introduced by Noiray et al. (Reference Noiray, Bothien and Schuermans2011), consists in writing ${\dot Q}^{\prime}$ in the form of a third-order polynomial of the pressure ${\dot Q}^{\prime}=\beta ^* p-\kappa ^*p^3$ , which, using (4.7), reads

(4.8) \begin{equation} {\dot Q}^{\prime}= \beta ^*\eta \psi _f- \kappa _1^* \eta ^3 \psi _f^3 ,\end{equation}

where $\psi _f= \psi ({\boldsymbol x}_f)$ . It is then convenient to define

(4.9) \begin{equation} \beta = [(\gamma -1)/\Lambda ] \beta ^* \psi _f^2 \quad \textrm {and} \quad \kappa _1= [(\gamma -1)/\Lambda ] \kappa _1^* \psi _f^4 .\end{equation}

With these notations, the dynamical equation governing $\eta$ becomes

(4.10) \begin{equation} {\ddot \eta }+ \alpha {\dot \eta } + \omega _0^2 \eta =\beta {\dot \eta } - 3{\kappa _1} \eta ^2 {\dot \eta }.\end{equation}

In this expression, $\beta$ is the linear rate of growth corresponding to small pressure perturbations and $\kappa _1$ is positive and governs the saturation process taking place at large oscillation amplitudes. The damping coefficient $\alpha$ is the rate at which $\eta ^2$ decays as a function of time when the right-hand side of (4.10) is equal to zero.

Solutions of this nonlinear differential equation may be obtained by making use of the method of averaging introduced by Krylov & Bogoliubov (Reference Krylov and Bogoliubov1950). This standard method is used here, as in most previous studies on azimuthal combustion instabilities referenced in the introduction. As pointed out by a reviewer, this method may not always yield the best approximation of the slow-flow variables, and other approaches, such as the multiple time scales method, can be used, as exemplified in Sirignano & Krieg (Reference Sirignano and Krieg2016). The reader is also referred to Cole (Reference Cole1968), Nayfeh & Mook (Reference Nayfeh and Mook1979) and Verhulst (Reference Verhulst1996) for additional information and the comparison of different approaches for the obtention of slow-flow variable equations. However, the advantage of the averaging method is that it gives access to the slow-flow variables with a reasonable amount of calculations and is sufficient to reveal the effects of the different HRR formulations.

The averaging method is quite standard, but the main steps are nevertheless provided here to facilitate the understanding of the derivation. The solution is first written in terms of a slowly varying amplitude $a$ and a phase $\varphi$

(4.11) \begin{equation} \eta (t)= a(t) \cos [\omega _0 t + \varphi (t)]. \end{equation}

To introduce this expression in (4.10), one needs to calculate the first derivative of $\eta$

(4.12) \begin{equation} {\dot \eta }= {\dot a} \cos (\omega _0 t + \varphi )-a \omega _0 \sin (\omega _0 t + \varphi ) - a {\dot \varphi }\sin (\omega _0 t + \varphi ). \end{equation}

In the method of averaging, it is standard to impose at this stage

(4.13) \begin{equation} {\dot a} \cos (\omega _0 t + \varphi )- a {\dot \varphi }\sin (\omega _0 t + \varphi )=0 ,\end{equation}

so that ${\dot \eta }= - a \omega _0\sin (\omega _0 t + \varphi )$ , and one may then proceed to calculate the second derivative of $\eta$

(4.14) \begin{equation} {\ddot \eta }=-a \omega _0^2 \cos (\omega _0 t + \varphi )-{\dot a}\omega _0\sin (\omega _0 t + \varphi ) -a\omega _0 {\dot \varphi }\cos (\omega _0 t + \varphi ) .\end{equation}

Inserting the previous expression in (4.10) and noting $\phi = \omega _0t +\varphi$ , one obtains, after some straightforward calculations,

(4.15) \begin{equation} {\dot a}\sin \phi +a {\dot \varphi }\cos \phi =-\alpha a \sin \phi + \beta a \sin \phi - 3 \kappa _1 a^3 \cos ^2 \phi \sin \phi. \end{equation}

This last equation, together with (4.13), form a linear system. The determinant of this system is $\Delta =1$ , and it is a simple matter to obtain

(4.16) \begin{eqnarray} {\dot a}&= -\alpha a \sin ^2 \phi + \beta a \sin ^2 \phi - 3 \kappa _1 a^3 \cos ^2 \phi \sin ^2\phi ,\end{eqnarray}

(4.17) \begin{eqnarray} a{\dot \varphi }&=-\alpha a \sin \phi \cos \phi + \beta a \sin \phi \cos \phi - 3 \kappa _1 a^3 \cos ^3 \phi \sin \phi .\end{eqnarray}

These equations may then be averaged over a period of oscillation $T= 2 \pi / \omega _0$ by taking into account that $a$ and $\varphi$ on the right-hand side do not vary over that period. The averages on the left-hand side are $[a(t+T)-a(T)]/T$ and $[\varphi (t+T)-\varphi (T)]/T$ , which represent the slow variable derivatives with respect to time

(4.18) \begin{align}\frac {\textrm{d}\it a}{\textrm{d}\it t} = [-\alpha /2 + \beta /2 -3 \kappa _1 a^2/ 8]a ,\end{align}

(4.19) \begin{align} a \frac {\textrm{d}\varphi }{\textrm{d}\it t} = 0. \end{align}

According to this model, the phase $\varphi$ remains constant, while the rate of change of the amplitude $a$ is given by the differential equation (4.18), which may also be written

(4.20) \begin{equation} \frac {1}{a}\frac {\textrm{d}\it a}{\textrm{d}\it t}= \frac { (\beta - \alpha )}{2} -\frac {3}{8} \kappa _1 a^2 .\end{equation}

When $\beta \lt \alpha$ , the right-hand side of this equation is negative and the amplitude decays from its initial value at a rate that is always greater than $(\beta - \alpha )/{2}$ . When $\beta \gt \alpha$ , the previous equation has a stationary solution with a finite amplitude corresponding to a limit cycle

(4.21) \begin{equation} a_s= [(4/3)(\beta -\alpha )/(\kappa _1)]^{1/2} .\end{equation}

This behaviour is typical of a Van der Pol oscillator and one may ask whether the stationary solution is stable. This question can be settled by considering the time evolution of a small perturbation in the amplitude, $a=a_s+\epsilon$ . Inserting this expression in the dynamical equation for the amplitude and only retaining first-order terms in $\epsilon$ , one obtains, after some straightforward calculations,

(4.22) \begin{equation} \frac {\textrm{d}\epsilon }{\textrm{d}\it t}=-\frac {3}{4}\kappa _1 a_s^2 \epsilon .\end{equation}

The perturbation in amplitude diminishes exponentially if $\kappa _1$ is positive and one concludes that the stationary solution corresponds to a stable limit cycle.

If now the flame response is delayed by a time lag $\tau _p$ , such that ${\dot Q}^{\prime}=\beta ^*p(t-\tau _p)-\kappa ^*p(t-\tau _p)^3$ , one gets the following system for the slow-flow variable equations:

(4.23) \begin{align} \frac {1}{a}\frac {\textrm{d}\it a}{\textrm{d}\it t} = -\frac {1}{2}[\alpha - (\beta -3 \kappa _1 a^2/ 4) \cos (\omega _0 \tau _p)] ,\end{align}

(4.24) \begin{align} \frac {\textrm{d}\varphi }{\textrm{d}\it t} = -\frac {1}{2}(\beta - 3 \kappa _1 a^2/4) \sin (\omega _0 \tau _p) .\end{align}

Compared with (4.18) and (4.19), one can see that the introduction of a time delay in the flame response induces a drift in the slow-flow variable $\varphi$ and modifies the growth rate and the amplitude of the limit cycle.

However, the problem is that the model used in this section to represent HRR fluctuations as a function of pressure does not allow for an easy and flexible adaptation to experimental flame dynamics data, such as those reported in § 3, figure 4, commonly expressed in terms of FDFs, as a function of slow-flow variables such as the oscillation level. In addition, the amplitude of the stationary solution (4.21) is proportional to the square root of $\beta -\alpha$ , and we will see later on (§ 6) that experiments do not comply with this law. For these reason, an alternative HRR model is examined in § 4.2.

4.2. Heat release rate expressed as a function of the oscillation amplitude: ${\dot Q}^{\prime}= g^*(a)p$

One may now consider a second model that is more flexible and easily adaptable to the experimental data reported in § 3 (figure 4). The choice of this formulation is motivated by the fact that experimental FDF data available in the literature (and used as input in the modelling framework of the kind considered here) are commonly presented as a function of oscillation amplitude level or slow-flow variables (see for instance Ghirardo et al. (Reference Ghirardo, Juniper and Moeck2016) who use experimental FDF data reported as a function of the level of oscillation or Ghirardo et al. (Reference Ghirardo, Nygård, Cuquel and Worth2021) who consider flame response functions expressed in terms of slow-flow variables, like the nature angle, $\chi$ ). An expression for the HRR as a function of slow-flow variables is better suited for practical applications and corresponds to the natural way of presenting flame dynamics data found in the literature. Investigating this kind of formulation hence enables a simpler comparison between different flame dynamics data.

The central idea is that the ratio between the HRR and pressure fluctuations is a function of the oscillation amplitude level $g^*(a)$ , so that ${\dot Q}^{\prime}=g^*(a)p$ . Considering, as in § 4.1, a coupling by a single mode of modal amplitude $\eta$ ( $p=\eta \psi$ ), the HRR term reads

(4.25) \begin{equation} {\dot Q}^{\prime}= g^*(a) \eta \psi _f .\end{equation}

To simplify notations, it is convenient to define

(4.26) \begin{equation} g(a)= \frac {\gamma -1}{\Lambda } g^*(a) \psi _f^2 ,\end{equation}

and the dynamical equation for $\eta$ now becomes

(4.27) \begin{equation} {\ddot \eta }+ \alpha {\dot \eta } + \omega _0^2 \eta =g^{\prime}(a){\dot a}\eta + g(a) {\dot \eta } .\end{equation}

Introducing $\phi = \omega _0 t + \varphi$ and noting that $\dot \eta = -a\omega _0 \sin (\omega _0 t + \varphi )$ , a few calculations yield

(4.28) \begin{equation} {\dot a}\left[\sin \phi + \frac {g^{\prime}(a)a}{\omega _0}\cos \phi \right] +a {\dot \varphi }\cos \phi = -\alpha a \sin \phi + g(a)a \sin \phi . \end{equation}

The system of equations now reads

(4.29) \begin{equation} \begin{bmatrix} \cos \phi & &-\sin \phi \\ \sin \phi +\frac {g^{\prime}(a)a}{\omega _0}\cos \phi & &\cos \phi \end{bmatrix} \begin{bmatrix} {\dot a}\\ a{\dot \varphi } \end{bmatrix} = \begin{bmatrix} 0 \\ -\alpha a \sin \phi + g(a)a \sin \phi \end{bmatrix} ,\end{equation}

and the determinant is

(4.30) \begin{equation} \Delta =1 + \frac {g^{\prime}(a)a}{\omega _0} \sin \phi \cos \phi .\end{equation}

System 4.29 may then be solved, yielding

(4.31) \begin{align} {\dot a} = \frac {1}{\Delta }[-\alpha a + g(a)a] \sin ^2 \phi ,\end{align}

(4.32) \begin{align} a{\dot \varphi } = \frac {1}{\Delta }[-\alpha a + g(a)a]\sin \phi \cos \phi .\end{align}

Averaging over a period is now complicated because the determinant appears in the denominator but it is possible to make an approximation by considering that $g^{\prime}(a)a/\omega _0\lt \lt 1$ (see Appendix B, where the order of magnitude of this term is estimated) and write

(4.33) \begin{align} {1}/{\Delta } \simeq 1- \frac{g^{\prime}(a)a}{\omega _0}\sin \phi \cos \phi. \end{align}

Introducing this expression in (4.31) and (4.32) and integrating over a period yields

(4.34) \begin{align} \frac {\textrm{d}\it a}{\textrm{d}\it t} = \frac {1}{2} [-\alpha +g(a)]a ,\end{align}

(4.35) \begin{align} \frac {\textrm{d}\varphi }{\textrm{d}\it t} = [-\alpha + g(a)]\frac {[-g^{\prime}(a)a]}{8\omega _0} .\end{align}

These equations indicate that both $a$ and $\varphi$ slowly vary with time. An example of resolution of (4.27) when $g(a) = \beta -\kappa _2 a$ is shown in figure 5(a). The envelope is obtained by integrating (4.34) for the slow variable $a$ . A steady solution exists if $g(a)\gt \alpha$ in an interval of amplitudes and if $g(a)-\alpha$ vanishes for an amplitude $a_s$ . Here again, one may ask whether the stationary solution corresponds to a stable limit cycle. For this, one can consider a small perturbation in amplitude, such that $a=a_s+ \epsilon$ in the vicinity of the stationary solution (or solutions), defined by $g(a_s)-\alpha =0$ . Only retaining first-order terms with respect to $\epsilon$ , one obtains

(4.36) \begin{equation} \frac {\textrm{d}\epsilon }{\textrm{d}\it t}=\frac {1}{2} g^{\prime}(a_s) a_s \epsilon .\end{equation}

The perturbation vanishes exponentially if $g^{\prime}(a_s)\lt 0$ and, in that case, the stationary solution corresponds to a stable limit cycle. However, if $g^{\prime}(a_s)\gt 0$ , the small perturbation grows exponentially and the stationary solution is unstable.

To fix the ideas, one may consider the case where $g(a)$ is a linearly decreasing function of the amplitude, $f(a)=\beta -\kappa _2 a$ . In that case

(4.37) \begin{align} \frac {\textrm{d}\it a}{\textrm{d}\it t} = \frac {1}{2} ( \beta -\alpha - \kappa _2 a)a ,\end{align}

(4.38) \begin{align} \frac {\textrm{d}\varphi }{\textrm{d}\it t} = ( \beta -\alpha - \kappa _2 a)\frac {\kappa _2 a}{8\omega _0} .\end{align}

If $\beta \gt \alpha$ , the stationary solution has an amplitude $a_s=(\beta -\alpha )/\kappa _2$ and the rate of change of $\varphi$ is positive when the amplitude $a$ is below $a_s$ . This situation corresponds to a positive shift in the frequency of oscillation. We will see later on that the linear relation between $\beta -\alpha$ and the limit-cycle amplitude deduced in this case, is much closer to what is observed experimentally (see the end of § 6 and Appendix E).

Figure 5. (a) Typical solution of the second-order oscillator (4.27) when $g(a) = \beta -\kappa _2 a$ . The envelope is obtained by integrating (4.34) for the slow variable $a$ . (b) Limit-cycle oscillation amplitude $a_s$ as a function of $\beta$ for the polynomial HRR formulation (black curve) and the new formulation of the form ${\dot Q}^{\prime}= g^*(a)p$ (red curve) for different values of the saturation coefficients $\kappa _1$ and $\kappa _2$ . Here, $\alpha$ designates the damping rate of the system.

It is also instructive to model the situation where the flame response is delayed by a time lag $\tau _p$ , and in this case, ${\dot Q}^{\prime}= g^*(a(t-\tau _p)) p(t-\tau _p)$ . Since the order of magnitude of the delay is well below that of a period of oscillation, one may assume that $a(t-\tau _p) \approx a(t)$ ( $a$ is a slow-flow variable assumed constant over a period of oscillation) and write ${\dot Q}^{\prime}= g^*(a)p(t-\tau _p)$ . For the pressure field of (4.7), one has

(4.39) \begin{equation} {\dot Q}^{\prime}= g^*(a) \eta (t-\tau _p) \psi _f .\end{equation}

Introducing $g(a)= [{(\gamma -1)}/{\Lambda }]g^*(a) \psi _f^2$ and applying the method of averaging, one obtains the following slow-flow variable equations:

(4.40) \begin{align} \frac {1}{a} {\dot a} &= -\frac {1}{2} [\alpha - g(a)\cos \omega _0 \tau _p] + \frac {g^{\prime}(a) a}{\omega _0}\left[\frac {3}{8} \alpha \sin \omega _0 \tau _p - \frac {1}{8} g(a) \sin 2 \omega _0 \tau _p\right] ,\end{align}

(4.41) \begin{align} {\dot \varphi } &= -\frac {1}{2} g(a) \sin \omega _0 \tau _p +\frac {1}{8} \frac {g^{\prime}(a) a}{\omega _0}[\alpha \cos \omega _0 \tau _p - g(a) \cos 2\omega _0 \tau _p] .\end{align}

One checks that (4.34) and (4.35) are retrieved when $\tau _p=0$ . Noting, in addition, that $g^{\prime}(a)a/\omega _0 \lt \lt 1$ , one finds that a non-zero delay induces a change in the gain of the form $+(1/2) g(a) \cos \omega _0 \tau _p$ and a shift in frequency $\Delta \omega$ with respect to $\omega _0$ such that $\Delta \omega \simeq -(1/2) g(a) \sin \omega _0 \tau _p$ . Using these two expressions, one retrieves the result, verified experimentally, and quoted in § 2, that the relative shift in frequency is generally quite small.

At this stage, it is also worth plotting the limit-cycle amplitudes $a_s$ as a function of $\beta$ for the two HRR formulations ( ${\dot Q}^{\prime}= \beta p- \kappa _1 p^3$ and ${\dot Q}^{\prime}= (\beta -\kappa _2 a) p$ ) and for different values of the saturation constants $\kappa _1$ and $\kappa _2$ , as shown in figure 5(b). One may note that, in both cases, the effective growth rate is such that $\omega^{\prime}_i(0)={1}/{2} (\beta -\alpha )$ . The diagram shown in figure 5 is also reminiscent of the one plotted in Latour et al. (Reference Latour, Durox, Renaud and Candel2024a ), (figure 16 of that reference), where the limit-cycle amplitudes obtained were found to be in a quasi-linear relation with the effective linear growth rate. We will see later on that this behaviour can be theoretically explained with the model derived in what follows.

Finally, it is interesting to note that, by injecting a quadratic expression of the form $g(a)=\beta -\kappa a^2$ in (4.40) and (4.41) and remembering that the term associated with $g^{\prime}(a)a/\omega _0$ is negligible, one retrieves (4.23) and (4.24), obtained for the HRR model discussed in § 4.1. However, the HRR formulation of the form ${\dot Q}^{\prime}= g^*(a)p$ has the advantage of being flexible and easily adaptable to pressure-based FDF data reported in terms of slow-flow variables (or velocity based FDFs if an impedance is used).

5.1. Governing equations

The starting point is the non-homogeneous wave equation derived in § 4.1

(5.1) \begin{equation} \frac {\partial ^2 p}{\partial t^2}-\boldsymbol {\nabla } \boldsymbol { \cdot } c^2\boldsymbol {\nabla }p = S -\alpha \frac {\partial p}{\partial t}. \end{equation}

The source term $S$ is now expressed as the sum of the contributions of $N$ acoustically compact flames behaving like point sources

(5.2) \begin{equation} S=\sum _{j=1}^N (\gamma -1) \frac {\partial {\dot q}^{\prime}_j}{\partial t} ,\end{equation}

with ${\dot q}^{\prime}_j$ the HRR fluctuations at flame $j$ , expressed as

(5.3) \begin{equation} {\dot q}^{\prime}_j={{\dot Q}^{\prime}}_j \delta (\boldsymbol {x}-\boldsymbol {x}_j) ,\end{equation}

where $\delta$ is the Dirac function centred at the position $\boldsymbol {x}_j$ of the $j$ th flame and ${{\dot Q}^{\prime}}_j$ the HRR fluctuation integrated over the volume of that flame.

The modal eigenfunctions pertaining to the MICCA annular combustor are discussed in Appendix C. In what follows, only two eigenmodes will be retained in the expression of the pressure field of (4.4), both corresponding to the same eigenfrequency

(5.4) \begin{equation} p = \eta _1(t) \psi _1(\boldsymbol {x}) + \eta _2(t) \psi _2(\boldsymbol {x}) ,\end{equation}

with $\psi _1$ and $\psi _2$ , two orthogonal 1A1L modes, expressed as $\psi _1(\boldsymbol {x})=\cos (\theta ) \psi _f$ and $\psi _2(\boldsymbol {x})=\sin (\theta ) \psi _f$ , $\psi _f$ representing the 1L axial wave function evaluated at the flame barycentre position, $x_f$ ( $\psi _{1L}(x_f)=\cos [\pi x_f/(2l^{\prime})]$ , presented in § 2), and $\eta _1$ and $\eta _2$ corresponding to the associated modal amplitudes.

As in § 4 and discussed in the supplementary material, the coupling with the injector ports is not taken into account: the changes in section induced by the injection units decouple the chamber and the plenum cavity and there is no plenum mode matching the eigenfrequency of the 1A1L chamber mode, eliminating the possibility of veering.

Using the orthogonality of the eigenmodes basis and integrating over the volume of the combustor, one gets the following equation for the amplitude $\eta _n$ of the $n$ th mode:

(5.5) \begin{equation} {\ddot \eta }_n + \alpha {\dot \eta }_n + \omega _n^2 \eta _n=\frac {1}{\Lambda _n} \int _{V} S \psi _n \textrm{d}V,\end{equation}

where $n=1,2$ and $\Lambda _n=\int _{V} \psi _n \psi _n^* \textrm{d}V$ is the normalization factor. It is shown in Appendix C that, for azimuthal modes having a non-zero azimuthal number $n$ , $\Lambda _n=V/4$ , where $V$ is the combustor volume. Injecting (5.2) and (5.3) in the right-hand side of (5.5), one finally obtains

(5.6) \begin{equation} \frac {1}{\Lambda _n} \int _{V} S \psi _n \textrm{d}V= \frac {(\gamma -1)}{\Lambda _n} \sum _{j=1}^N \frac {\partial {{\dot Q}^{\prime}}_j}{\partial t} \psi _n(x_j). \end{equation}

5.2. Heat release rate model

One now needs an analytical expression for the HRR response of the flame to acoustic disturbances. As discussed in § 4.2, the HRR ${{\dot Q}^{\prime}}_j$ at the $j$ th flame is now expressed as a function $\mathcal{G}$ of the local amplitude, at the angular frequency of interest $\omega _0$ and local pressure with a time delay

(5.7) \begin{equation} {{\dot Q}^{\prime}}_j= \mathcal {G}(\Pi _j, \omega _0) p_j(t-\tau _p), \end{equation}

where the local amplitude at flame $j$ is the slow-flow variable, $\Pi _j=(p_{rms})_j/(\rho U_b^2)$ , which will be expressed later in terms of the variables of interest. It is also shown in § 5.3 how $\mathcal {G}$ is linked to the pressure based FDF at the frequency $\omega _0$ of interest.

The local amplitude $\Pi _j$ being a slow-flow variable assumed constant over a period of oscillation and, in addition, as the terms associated with the ${\dot {\mathcal {G}}}$ term were shown to be negligible in the SFVE in § 4, one may write

(5.8) \begin{equation} \frac {\partial {{\dot Q}^{\prime}}_j}{\partial t}= \mathcal {G}(\Pi _j,\omega _0) {\dot p}_j(t-\tau _p) ,\end{equation}

and the governing equations projected on the modes $\psi _1$ and $\psi _2$ read

(5.9) \begin{align} {\ddot \eta }_1 & + \alpha {\dot \eta }_1 + \omega _0^2 \eta _1 =\frac {\gamma -1}{\Lambda }\sum _{j=1}^N \mathcal { G}(\Pi _j,\omega _0) [{\dot \eta }_1(t-\tau _p) \psi _1(x_j)+{\dot \eta }_2 (t-\tau _p) \psi _2(x_j)] \psi _1(x_j) \nonumber\\ {\ddot \eta }_2 & + \alpha {\dot \eta }_2 + \omega _0^2 \eta _2 =\frac {\gamma -1}{\Lambda }\sum _{j=1}^N \mathcal { G}(\Pi _j,\omega _0) [{\dot \eta }_1(t-\tau _p) \psi _1(x_j)+{\dot \eta }_2 (t-\tau _p) \psi _2(x_j)] \psi _2(x_j) .\end{align}

5.3. Expression of the HRR using the pressure-based FDF

One now seeks to link the time domain expression involving the function $\mathcal{G}$ to the pressure-based FDF, $\mathcal { F}_p=G_{p_j}(\Pi _j,\omega )e^{i\varphi _{p_j}}$ , as defined in (3.1), where $\Pi _j$ is the reduced pressure oscillation root-mean-square (r.m.s.) value and $G_{p_j}$ and $\varphi _{p_j}$ are the FDF gain and phase of the $j$ th flame. It is important to remember at this point that the dynamics of the system only takes place around the eigenfrequency of the 1A1L mode, identified to be involved in the combustion/acoustics coupling in MICCA. One hence looks for a solution of the pressure field in the form

(5.10) \begin{equation} p_j= A_1 \cos (\omega _0 t+\varphi _1) \psi _1(x_j) + A_2 \cos (\omega _0 t+\varphi _2) \psi _2(x_j),\end{equation}

where $A_1, A_2, \varphi _1$ and $\varphi _2$ are slowly varying compared with the oscillation period $T=2\pi /\omega _0$ . One can also identify from (5.4), $\eta _1(t)= A_1 \cos (\omega _0 t+\varphi _1)$ and $\eta _2=A_2 \cos (\omega _0 t+\varphi _2)$ .

It is convenient to introduce at this stage the analytic signals ${\widehat p}_j$ and ${\widehat {\dot Q}}_j$ , such that the pressure and HRR signals are the real parts of these complex quantities: $p_j= \textrm {Re}({\widehat p}_j)$ and ${{\dot Q}^{\prime}}_j= \textrm {Re}( {\widehat {\dot Q}}_j)$ . Now, the complex pressure signal may be written in terms of the eigenmodes $\psi _1$ and $\psi _2$ as

(5.11) \begin{equation} {\widehat p}_j = {{\widehat \eta }_1} \psi _1(x_j) + {{\widehat \eta }_2} \psi _2(x_j) , \end{equation}

where the complex amplitudes $\widehat {\eta }_1$ and ${\widehat \eta }_2$ are expressed as ${\widehat \eta }_1=A_1 \exp {(-i\omega _0 t-i\varphi _1)}$ and ${\widehat \eta }_2=A_2 \exp {(-i\omega _0 t-i\varphi _2)}$ .

The complex HRR signal ${\widehat {\dot Q}}_j$ at flame $j$ may now be linked to the complex pressure signal by the describing function $\mathcal {F}_p$ , expressed in terms of a gain and a phase and considered at the angular frequency $\omega _0$ :

(5.12) \begin{equation} {\widehat {\dot Q}}_j=G_{p_j}(\Pi _j,\omega _0)e^{i\varphi _{p_j}} \frac {{\dot Q}_0}{\rho U_b^2} \widehat {p}_j .\end{equation}

This is just the describing function extension of a standard result of linear system theory which indicates that, when the input to the system is a complex sinusoidal signal at the frequency $\omega _0$ , its complex signal output features the same frequency and is equal to the product of the transfer function at that frequency by the complex signal input. Injecting (5.11) in (5.12), one obtains

(5.13) \begin{equation} {\widehat {\dot Q}}_j=G_{p_j}(\Pi _j,\omega _0) \frac {{\dot Q}_0}{\rho U_b^2} [{\widehat \eta }_1 e^{i\varphi _{p_j}} \psi _1(x_j)+{\widehat \eta }_2 e^{i\varphi _{p_j}} \psi _2(x_j) ]. \end{equation}

One may now deduce the HRR signal by taking the real part of the complex HRR signal in (5.13)

(5.14) \begin{align} {{\dot Q}^{\prime}}_j & = G_{p_j}(\Pi _j,\omega _0) \frac {{\dot Q}_0}{\rho U_b^2} [A_1 \cos (\omega _0 t+\varphi _1-\varphi _{p_j})\psi _1(x_j)\nonumber\\& \quad +A_2 \cos (\omega _0 t+\varphi _2-\varphi _{p_j})\psi _2(x_j)] .\end{align}

It is here convenient to define $\tau _{p_j}=\varphi _{p_j}/\omega _0$ and write the previous expression as

(5.15) \begin{align} {{\dot Q}^{\prime}}_j & =G_{p_j}(\Pi _j,\omega _0) \frac {{\dot Q}_0}{\rho U_b^2} [A_1 \cos (\omega _0 (t-\tau _{p_j})+\varphi _1)\psi _1(x_j)\nonumber\\& \quad +A_2 \cos (\omega _0 (t-\tau _{p_j})+\varphi _2)\psi _2(x_j)] .\end{align}

One identifies in this expression the delayed pressure signal $p_j(t-\tau _{p_j})$ so that one may write

(5.16) \begin{equation} {{\dot Q}^{\prime}}_j=G_{p_j}(\Pi _j,\omega _0) \frac {{\dot Q}_0}{\rho U_b^2} p_j(t-\tau _{p_j}) .\end{equation}

Slow-flow variables $A_1, A_2, \varphi _1, \varphi _2$ and $\Pi _j$ being assumed constant over a period of oscillation, the time derivative of the HRR simply involves the rate of change of the pressure delayed by $\tau _{p_j}$

(5.17) \begin{equation} \frac {\partial {{\dot Q}^{\prime}}_j}{\partial t} =G_{p_j}(\Pi _j,\omega _0) \frac {{\dot Q}_0}{\rho U_b^2} {\dot p}_j(t-\tau _{p_j}). \end{equation}

Comparing (5.8) and (5.17), one finds that

(5.18) \begin{equation} \mathcal {G}(\Pi _j, \omega _0) = G_{p_j}(\Pi _j,\omega _0)\frac {{\dot Q}_0}{\rho U_b^2} \end{equation}

and the dynamical system of (5.5) finally becomes

(5.19) \begin{align} {\ddot \eta }_1 & + \alpha {\dot \eta }_1 + \omega _0^2 \eta _1 =\frac {\gamma -1}{\Lambda }\sum _{j=1}^N G_{p_j} \frac {{\dot Q}_0}{\rho U_b^2} [{\dot \eta }_1(t-\tau _{p_j}) \psi _1(x_j)+{\dot \eta }_2 (t-\tau _{p_j}) \psi _2(x_j)] \psi _1(x_j) \nonumber\\ {\ddot \eta }_2 & + \alpha {\dot \eta }_2 + \omega _0^2 \eta _2 =\frac {\gamma -1}{\Lambda }\sum _{j=1}^N G_{p_j} \frac {{\dot Q}_0}{\rho U_b^2} [{\dot \eta }_1(t-\tau _{p_j}) \psi _1(x_j)+{\dot \eta }_2 (t-\tau _{p_j}) \psi _2(x_j)] \psi _2(x_j), \end{align}

where $G_{p_j}=G_{p_j}(\Pi _j,\omega _0)$ and the local oscillation amplitude $\Pi _j$ , expressed in terms of the slow-flow variables, reads

(5.20) \begin{equation} \Pi _j=\frac {1}{\rho U_b^2 \sqrt {2}} \big[A_1^2\psi _1^2+A_2^2\psi _2^2+2A_1A_2\psi _1\psi _2\cos (\varphi _2-\varphi _1)\big]^{1/2}. \end{equation}

Using, in addition, the standard hypothesis of the method of averaging, the time derivative of the amplitude $\eta$ is

(5.21) \begin{align} {\dot \eta }_k&=-A_k(t) \omega _0 \sin (\omega _0 t+\varphi _k) \, \, \textrm {for}\,\, k=1,2,\end{align}

and one can express the slow-flow variables as a function of $\eta _k$ and ${\dot \eta }_k$

(5.22) \begin{equation} A_k= \left[{\eta _k^2+\left(\frac {{\dot \eta }_k}{\omega _0}\right)^2}\right]^{1/2}, \,\, \textrm {and} \,\, A_j A_k \cos (\varphi _j-\varphi _k)= \eta _j\eta _k+\frac {1}{\omega _0^2} {\dot \eta }_j {\dot \eta }_k .\end{equation}

The dynamics can hence be described by two coupled second-order oscillators. The coupling between the two eigenmodes appears on the right-hand side of these two equations, corresponding to the projection of the source term $S$ on each eigenmode. Equations (5.19) can be solved directly, as done, for instance, in Noiray et al. (Reference Noiray, Bothien and Schuermans2011). The integration is here carried out by making use of the FDF representation obtained in § 3. Typical examples of the direct integration of 5.19 are provided for staging patterns $\mathcal {C}_0$ , $\mathcal {C}_4$ , $\mathcal {L}_4$ and $\mathcal {A}_4$ and shown in figure 7 in the form of plots in the $(\eta _1, \eta _2)$ plane. Circular trajectories in this plane found for $\mathcal {C}_0$ or for $\mathcal {A}_4$ indicate that the oscillation at the limit cycle takes the form of a spinning mode. The elongated trajectories characterizing $\mathcal {C}_4$ or $\mathcal {L}_4$ are typical of standing mode oscillations reflecting the symmetry breaking induced by the placement of four S injectors.

Figure 7. Direct integration of the ordinary differential system with delay (5.19), for an initial condition $[\eta _1 =75, \dot {\eta }_1=50, \eta _2 =15, \dot {\eta }_2=10]$ , for configurations $\mathcal {C}_0$ , $\mathcal {C}_4$ , $\mathcal {L}_4$ and $\mathcal {A}_4$ .

5.4. Slow-flow equations

A further understanding of the system’s dynamics may be gained by deriving a set of equations for the slow-flow variables $A_1, A_2, \varphi _1$ and $\varphi _2$ using the method of averaging. Solutions of (5.19) are now sought in the standard form

(5.23) \begin{align}\eta _k & =A_k(t) \cos (\omega _0 t+\varphi _k) \nonumber\\ {\dot \eta }_k & =-A_k(t) \omega _0 \sin (\omega t+\varphi _k) \, \, \textrm {for}\,\, k=1,2. \end{align}

Using the same procedure as that described in § 4 and integrating over a period of oscillation $T=2\pi /\omega _0$ , this system yields the following equations for the slow-flow variables $A_1$ , $A_2$ , $\varphi _1$ and $\varphi _2$ :

(5.24) \begin{align} {\dot A}_1 & =-\frac {\alpha }{2} A_1 +\frac {\gamma -1}{2\Lambda } \frac {{\dot Q}_0}{\rho U_b^2}\sum _{j=1}^N G_{p_j} \big[A_1 \cos (\varphi _{p_j}) \psi _1(x_j)+A_2 \cos (\varphi _{p_j}\nonumber\\& \quad +\varphi _{12})\psi _2(x_j)\big]\psi _1(x_j)\nonumber\\ {\dot A}_2& =-\frac {\alpha }{2} A_2 +\frac {\gamma -1}{2\Lambda } \frac {{\dot Q}_0}{\rho U_b^2}\sum _{j=1}^N G_{p_j} \big[A_2 \cos (\varphi _{p_j}) \psi _2(x_j)+A_1 \cos (\varphi _{p_j}\nonumber\\& \quad +\varphi _{21})\psi _1(x_j)\big]\psi _2(x_j)\nonumber\\ A_1 {\dot \varphi }_1& = \frac {\gamma -1}{\Lambda }\frac {{\dot Q}_0}{\rho U_b^2} \sum _{j=1}^N G_{p_j} \big[-A_1 \sin (\varphi _{p_j}) \psi _1^2(x_j) +A_2 \sin (\varphi _{21}-\varphi _{p_j})\psi _1(x_j)\psi _2(x_j)\big]\nonumber\\ A_2 {\dot \varphi }_2& = \frac {\gamma -1}{\Lambda } \frac {{\dot Q}_0}{\rho U_b^2}\sum _{j=1}^N G_{p_j} \big[-A_2 \sin (\varphi _{p_j}) \psi _2^2(x_j) +A_1 \sin (\varphi _{12}-\varphi _{p_j})\psi _1(x_j)\psi _2(x_j)\big].\end{align}

In these equations $\varphi _{12}=\varphi _1-\varphi _2$ , $\varphi _{21}=\varphi _2-\varphi _1$ and the gains and phases, $G_{p_j} =G_j (\Pi _j,\omega _0)$ and $\varphi _{p_j}=\varphi _{p_j}(\Pi _j,\omega _0)$ , correspond to the flames established by the different injectors, with $\Pi _j=(p_{rms})_j/(\rho U_b^2)$ , the dimensionless pressure amplitude at flame $j$ , which is expressed in terms of the slow-flow variables in (5.20).

It is also interesting to deduce an evolution equation for the difference $\varphi _1-\varphi _2$ . This can be done by multiplying the third equation in the set (5.24) by $A_2$ and subtracting the fourth equation in this set after multiplication by $A_1$ . This yields

(5.25) \begin{align} A_1A_2 (\dot {\varphi }_1- \dot {\varphi }_2) & = \frac {(\gamma -1)}{\Lambda } \frac {{\dot Q}_0}{\rho U_b^2} \sum _{j=1}^N G_{p_j} \Big \{A_1A_2 \sin (\varphi _{p_j})[\psi _2^2(x_j)- \psi _1^2(x_j)]\nonumber\\& \quad -A_1^2 \sin (\varphi _1-\varphi _2-\varphi _{p_j})\psi _1(x_j)\psi _2(x_j)\nonumber\\& \quad + A_2^2 \sin (\varphi _2-\varphi _1-\varphi _{p_j})\psi _1(x_j)\psi _2(x_j) \Big \} .\end{align}

A comprehensive discussion on the theoretical conditions for the stability of fixed standing and spinning solutions of the system formed by (5.24) can be found in Ghirardo et al. (Reference Ghirardo, Juniper and Moeck2016). Additional terms due to turbulence-induced acoustic forcing were not considered in the present work, because, as will be seen in § 5.5, the dynamical equations, in the form proposed in (5.24), enable us to reasonably retrieve the features observed experimentally, which are mainly dominated by symmetry-breaking effects induced by injectors’ staging. There is, however, one exception, in the special case were all injectors are of the same kind: we will see that the dynamical equations predict a limit cycle in the form of a purely spinning wave, while experiments indicate that the system continuously switches between spinning modes and mixed modes of various types. In that case, the stochastic term representing turbulent disturbances might probably be necessary to obtain a suitable match with observations (on that point see Ghirardo et al. Reference Ghirardo, Boudy and Bothien2018).

5.5. Application to staging patterns in MICCA-spray

The system of ordinary differential equations (ODEs) obtained in § 5.4 is now solved for different configurations mixing U- and S-injectors in MICCA, shown in figure 6. For all the staging patterns tested, the instability frequencies lie between 775 and 820 Hz and the oscillation amplitudes between 100 and 1400 Pa. Further details on the experimental results may be found in Latour et al. (Reference Latour, Durox, Renaud and Candel2024a ).

Expressions for $G_{p_j}(\Pi _j,\omega _0)$ and $\varphi _{p_j}(\Pi _j,\omega _0)$ are those deduced from experimental data for the two injection units U and S, collected for a range of frequencies close to the eigenfrequency of the 1A1L mode of the MICCA combustor and reported in § 3. One can reasonably model the gain of the pressure-based FDF as $G_{p_j}=\beta -\kappa \Pi _j$ , as discussed in § 3. The coefficients $\beta$ and $\kappa$ for injectors U and S are deduced from a linear regression of the $G_p=G_p(\Pi )$ data points. The phases, $\varphi _S$ and $\varphi _U$ , are assumed to be constant and equal to the mean value, as shown by the experimental results reported in figure 4, where the phases for the two injectors take nearly constant values throughout the range of pressure amplitudes tested. Noting that the damping rate $\alpha$ describes the rate of change in acoustic energy (i.e. of pressure squared), the damping rate used for the numerical integration is assigned a value that is twice that obtained by Latour et al. (Reference Latour, Durox, Renaud and Candel2024b ), where the damping rate for the pressure is extracted from direct measurements of the Rayleigh source term while operating the MICCA combustor at limit cycle. It is also assumed that the damping rate value does not change with the staging pattern and may be represented by a constant value.

As pointed out by a reviewer, there are cases where the damping depends nonlinearly on the level of oscillation. For example, damping by Helmholtz resonators is a function of the level of oscillation, as shown by Zinn & Lores (Reference Zinn and Lores1972) or Ćosić, Reichel & Paschereit (Reference Ćosić, Reichel and Paschereit2012). The nonlinearity is, in this case, linked to the vortex shedding from the resonator outlet. The use of perforated liners or quarter wave cavities also introduces nonlinearities as indicated for example by Schuller et al. (Reference Schuller, Tran, Noiray, Durox, Ducruix and Candel2009). However, MICCA is not equipped with such devices and the damping nonlinearity is less probable. In addition, the damping rates of acoustic energy estimated in MICCA from energy balance considerations at various limit-cycle amplitudes for a range of limit-cycle oscillation levels (670–1400 Pa) by Latour et al. (Reference Latour, Durox, Renaud and Candel2024b ) show no amplitude dependence, as can be seen in figure 8. There are only minor variations in the damping rate for the different oscillation levels and one may safely conclude that nonlinearities need not be taken into consideration.

Using the parameters gathered in table 2, the differential equations are integrated with the MATLAB ode45 solver. Typical results shown in figure 9 pertain to the four staging patterns $\mathcal {C}_0, \mathcal {C}_4, \mathcal {L}_4$ and $\mathcal {A}_4$ .

Table 2. Pressure-based FDFs and damping rate used in the numerical simulations.

Figure 8. Damping rate values as a function of the oscillation level ${A}$ obtained from source term and acoustic energy estimates in MICCA during self-sustained limit-cycle oscillations for different staging configurations Latour et al. (Reference Latour, Durox, Renaud and Candel2024b ). The annular combustor is operated at a thermal power $\mathcal {P}=118\,\textrm{kW}$ and an equivalence ratio $\phi =0.9$ .

In the $\mathcal {C}_0$ case, the amplitudes $A_1$ and $A_2$ grow to the limit cycle in approximately $25$ ms and take the same value, while the phase difference $\varphi _1-\varphi _2$ tends to $\pi /2$ , indicating that the mode is a spinning wave. For configurations $\mathcal {C}_4$ and $\mathcal {L}_4$ , the amplitudes reach constant levels while the phase difference $\varphi _1-\varphi _2$ tends to zero. The oscillation now corresponds to a standing mode. The anti-nodal line angle $\theta _0$ may be deduced from the values of $A_1$ and $A_2$ at limit cycle, using the expression $\theta _0 = \arccos [A_1/ (A_1^2+ A_2^2)^{1/2}]$ . The nodal line found analytically is aligned with the diameter joining the two groups of S-injectors in diametrically opposed locations for $\mathcal {L}_4$ and passes at equal distance of the four S injectors in configuration $\mathcal {C}_4$ . These analytical results correspond to those found experimentally (Latour et al. Reference Latour, Durox, Renaud and Candel2024a ). Finally, in case $\mathcal {A}_4$ , the growth in amplitude is less rapid than in the $\mathcal {C}_0$ case and the levels at limit cycle are also lower. Similarly to configuration $\mathcal {C}_0$ , $A_1$ is equal to $A_2$ while $\varphi _1-\varphi _2$ tends to $\pi /2$ , which correspond to a spinning wave.

Figure 9. Slow-flow variables $A_1, A_2$ and $\varphi _1-\varphi _2$ time evolutions for configurations $\mathcal {C}_0$ , $\mathcal { C}_4$ , $\mathcal {L}_4$ and $\mathcal {A}_4$ , obtained from the resolution of the ordinary differential equations for the initial conditions $A_1=50\,\textrm{Pa}$ , $A_2=75\,\textrm{Pa}$ , $\varphi _1-\varphi _2=0.1\,\textrm{rad}$ .

It is next interesting to examine the predicted limit-cycle amplitudes and compare them with the values measured experimentally. This comparison may be carried out by making use of the link between the amplitude $A=(A_+^2+A_-^2)^{1/2}$ that is used in the experimental determination of the limit-cycle amplitudes (see § 2, (2.3) and Latour et al. Reference Latour, Durox, Renaud and Candel2024a ) and the amplitudes $A_1$ and $A_2$ of the eigenmodes $\psi _1$ and $\psi _2$ , intervening in (5.24). It is shown in Appendix D that

(5.26) \begin{equation} A=\big(A_+^2+A_-^2\big)^{1/2}=\frac {1}{{2^{1/2}}} \big(A_1^2+A_2^2\big)^{1/2} .\end{equation}

Figure 10. (a) Limit-cycle oscillation amplitudes obtained by integrating the system of ODEs for the slow-flow variables for the three types of staging configurations investigated (type $\mathcal{C}$ , $\mathcal{L}$ and $\mathcal{A}$ ). (b,c,d) Comparison of the limit-cycle oscillation amplitudes measured experimentally (Latour et al. Reference Latour, Durox, Renaud and Candel2024a ) and obtained by solving the system of ODEs for the slow-flow variables for type $\mathcal{C}$ (b), $\mathcal{L}$ (c) and $\mathcal{A}$ (d) configurations. Here $N_S$ corresponds to the number of S-injectors.

The predicted limit-cycle oscillation amplitudes $A$ , deduced by solving the slow-flow equations for $A_1$ and $A_2$ and making use of 5.26, are displayed in figure 10. The amplitude $A$ decreases when the number of S-injectors is increased and the trends corresponding to type $\mathcal{C}$ arrangements clearly differ from those pertaining to $\mathcal{L}$ configurations. A comparison with the experimental data from Latour et al. (Reference Latour, Durox, Renaud and Candel2024a ) is carried out in figure 10(b,c,d) where the $\mathcal{C}$ , $\mathcal{L}$ and $\mathcal{A}$ configurations appear in three separate graphs. The match obtained between experimental data and model prediction is, despite some differences, quite satisfactory. The limit-cycle amplitudes cover the same range with a maximum value of slightly less than 1500 Pa. Configurations $\mathcal {C}_8$ and $\mathcal {A}_8$ , which are experimentally stable, are suitably predicted by the model as having a low level of oscillation. The decrease in amplitude observed when the staging patterns is evolving from $\mathcal {C}_0$ to $\mathcal {C}_8$ is well retrieved but the calculated amplitudes are slightly higher than those measured experimentally.

The trends observed experimentally for the $\mathcal{L}$ -type staging patterns, where the limit-cycle amplitude levels decrease more slowly than for type $\mathcal{C}$ configurations as the number of S-injectors is augmented, are also found with the model, although larger discrepancies are observed between experimental data and simulation for these configurations. The greatest differences in limit-cycle amplitudes are found for $\mathcal {C}_4$ , $\mathcal {L}_4$ and $\mathcal {A}_4$ . This may be due to a different damping rate value for these staging arrangements (Latour et al. Reference Latour, Durox, Renaud and Candel2024b ) that is not accounted for in the modelling. It is also possible that some adjustment of the gain $G_p^S$ corresponding to the S-injectors might have improved the prediction, but no attempt was made to go in that direction. As indicated in § 3, the gain and phase values for S-injectors are less well determined because these units are located in regions where the oscillation level is low and this influences the measurement precision. Another possibility could be that the gain $G_p$ is influenced to some extent by the nature of the mode as inferred by Nygård, Ghirardo & Worth (2021) and Ghirardo et al. (Reference Ghirardo, Nygård, Cuquel and Worth2021), and that this might have to be accounted for in the pressure-based FDF formulation. However, although some discrepancies exist between experimental data and model predictions, this is perhaps the first time where calculations are shown to be capable of predicting limit-cycle levels of oscillation for a large set of experiments.

It is next interesting to examine the distribution of the r.m.s. of the pressure field as a function of the azimuthal position deduced from the numerical integration of the slow-flow equations. This distribution is shown in figure 11 for all the staging configurations. The first row in this figure corresponds to $\mathcal{C}$ configurations. For $\mathcal {C}_0$ , where all injectors are of the U-type, the r.m.s. pressure distribution is uniform, as expected for a spinning mode. As S-injectors are being added to replace U-injectors, the r.m.s. pressure distribution features two maximum and two minimum values. The minimum is reached at the location of the S-injector (case $\mathcal {C}_1$ ) or on the diameter passing through the centre of the $\mathcal{C}$ configuration as exemplified by $\mathcal { C}_2$ , $\mathcal {C}_4$ and $\mathcal {C}_6$ . In the last two cases, the minimum r.m.s. pressure is close to zero, corresponding to a standing mode. In case $\mathcal {C}_8$ , the calculated level of oscillation is very low, indicating that the system is stable. The second row shows the r.m.s. pressure distributions for $\mathcal{L}$ and $\mathcal{A}$ cases. The minimum values are observed on the diameter that passes through the S-injectors (case $\mathcal {L}_{2}$ ) or on the median diameter of the S-injector configurations ( $\mathcal {L}_{4}$ , $\mathcal {L}_{6}$ , $\mathcal {L}_{8}$ ). In the last three cases, the minimum value is close to zero, a clear signature of a standing mode. The cases $\mathcal {A}_4$ and $\mathcal {A}_8$ feature a uniform distribution of r.m.s. pressure as the symmetry is not broken by the presence of the 4 or 8 S-injectors. In case $\mathcal {A}_4$ , the amplitude is reduced but not to the point observed experimentally (this configuration is marginally unstable). In $\mathcal {A}_8$ , the level is quite low and the system is stable.

Figure 11. Root mean square of the pressure distribution as a function of the azimuthal position $\theta$ obtained from the numerical integration of the ODEs for the slow-flow variables. The red vertical lines correspond to S-injectors positions.

The previous calculations already show typical consequences of the breaking of symmetry in various configurations. These features generally retrieve those found in experiments. For example, it is known that the presence of a single S-injector in $\mathcal {C}_1$ is insufficient to break the symmetry and the pressure distribution only shows weak undulations. In cases $\mathcal {C}_4$ and $\mathcal {C}_6$ and in arrangements $\mathcal {L}_{4}$ , $\mathcal {L}_{6}$ and $\mathcal {L}_{8}$ , symmetry is broken, giving rise to a standing mode and this agrees well with what is found experimentally (see Latour et al. Reference Latour, Durox, Renaud and Candel2024b ). It is, however, interesting to pursue the examination of the nature of the unstable mode by computing the spin ratio, initially defined by Bourgouin et al. (Reference Bourgouin, Durox, Moeck, Schuller and Candel2013) in terms of azimuthal wave amplitudes $A_+, A_-$ as

(5.27) \begin{equation} s_r= \frac {A_+ - A_-}{A_+ + A_-} .\end{equation}

It is worth remembering that the nature angle $\chi$ , that is used in many recent studies of azimuthal instabilities (see for instance Ghirardo & Bothien Reference Ghirardo and Bothien2018; Faure-Beaulieu & Noiray Reference Faure-Beaulieu and Noiray2020; Indlekofer et al. Reference Indlekofer, Faure-Beaulieu, Dawson and Noiray2022), is linked to the the spin ratio by $\chi = \tan ^{-1}s_r$ .

Now, $s_r$ may be expressed in terms of the amplitudes $A_1$ and $A_2$ of the eigenmodes $\psi _1$ and $\psi _2$ using relations derived in Appendix D. One gets

(5.28) \begin{equation} s_r= \frac {A_1A_2 \sin (\varphi _1-\varphi _2)\cos 2\theta _0}{A_1^2\cos ^2 \theta _0-A_2^2\sin ^2 \theta _0} ,\end{equation}

where $\theta _0$ is the anti-nodal line position, which is also expressed in terms of $A_1$ and $A_2$ in Appendix D.

Statistics for the experimental nodal line position, spin ratio time series and joint probability density functions for $A_+$ and $A_-$ can be found in Latour et al. (Reference Latour, Durox, Renaud and Candel2024a ). The nature of the mode at limit cycle is here represented by placing the limit-cycle characteristics of the unstable configurations, obtained from the resolution of the ODEs, in the ( $|s_r|$ , $A$ ) plane, as shown in figure 12. In cases $\mathcal {C}_0$ and $\mathcal {A}_4$ , one finds $|s_r|=1$ , corresponding to a purely spinning mode. While calculations for symmetric unstable configurations lead to a spinning mode, experiments indicate that the spin ratio continuously changes between -1 and +1. This difference in behaviour between calculations and experiments confirms findings by Noiray et al. (Reference Noiray, Bothien and Schuermans2011) and Faure-Beaulieu & Noiray (Reference Faure-Beaulieu and Noiray2020) and is due, according to these authors, to the absence of a stochastic forcing. Such a forcing, which idealizes the effect of turbulence, is needed to obtain a variable spin ratio taking values over the interval $[-1,+1]$ .

It is also concluded from figure 12 that:

1. (i) Weak asymmetry in the staging pattern leads to a mixed mode ( $\mathcal {C}_1$ , $\mathcal {C}_2$ and $\mathcal {L}_2$ ).

2. (ii) Strong asymmetry in the staging pattern gives rise to a standing mode ( $\mathcal {C}_4$ , $\mathcal {C}_6$ , $\mathcal { L}_4$ , $\mathcal {L}_6$ and $\mathcal {L}_8$ ), with the nodal line aligned with the S-injectors median diameter. When standing modes prevail, the nodal line position predicted analytically corresponds to that observed experimentally.

Figure 12. Unstable configurations in the ( $A$ , $s_r$ ) plane.

5.6. Sensitivity analysis

It is natural at this stage to examine the sensitivity of the limit-cycle amplitudes to the choice of the FDF gain formulation and fitting coefficients. This is done here for the linear and quadratic fits shown in figure 13 for injectors U and S. The linear fit used in § 5.4 is first compared with a quadratic FDF gain formulation of the form $G_p=\beta -\kappa \Pi ^2$ . The predicted limit-cycle amplitudes are shown for the linear and quadratic expressions in figure 14(a) for type $\mathcal{C}$ configurations. One can see that both formulations lead to similar trends and oscillation amplitude values.

The sensitivity of the calculated amplitudes to the regression coefficients is then tested by changing by $\pm 5\, \%$ the fitted $\beta$ and $\kappa$ values for injector U. Results of this analysis are plotted in figure 14(b,c) for the linear fit, for type $\mathcal{C}$ configurations. One can see that the calculated results are most sensitive to $\beta$ as a $\pm 5\, \%$ change in the value of this parameter induces a $\pm 200\,\textrm{Pa}$ variation in the limit-cycle amplitude. A $\pm 5\, \%$ change in $\kappa$ has a lesser impact on the predicted limit-cycle oscillation amplitudes. A similar sensitivity analysis for the quadratic fit, not reproduced here, leads to the same conclusions.

Figure 13. Linear (a,b) and quadratic (c,d) fits for the pressure-based FDF gain experimental data for injector U (a,c) and injector S (b,d).

Figure 14. Comparison between a linear and a quadratic expression for the FDF gain on the limit-cycle oscillation amplitude (a). Effects of a $\pm 5\%$ change in the $\beta$ (b) and $\kappa$ (c) coefficients values on the analytical limit-cycle oscillation amplitudes for type $\mathcal{C}$ configurations. Here $N_S$ corresponds to the number of S-injectors.