Experiments and low-order modelling of intermittent transitions between clockwise and anticlockwise spinning thermoacoustic modes in annular combustors

Annular combustion chambers of gas turbines and aircraft engines are subject to unstable azimuthal thermoacoustic modes leading to high amplitude acoustic waves propagating in the azimuthal direction. For certain operating conditions, the propagating direction of the wave switches randomly. The strong turbulent noise prevailing in gas turbine combustors is a source of random excitation for the thermoacoustic modes and can be the cause of these switching events. A low-order model is proposed to describe qualitatively this property of the dynamics of thermoacoustic azimuthal modes. This model is based on the acoustic wave equation with a destabilizing thermoacoustic source term to account for the ﬂame’s response and a stochastic term to account for the turbulent combustion noise. Slow-ﬂow averaging is applied to describe the modal dynamics on times scales that are slower than the acoustic pulsation. Under certain conditions, the model reduces formally to a Fokker-Planck equation describing a stochastic di ﬀ usion process in a potential landscape with two symmetric wells: One well corresponds to a mode propagating in the clockwise direction, the other well corresponds to a mode propagating in the anticlockwise direction. When the level of turbulent noise is su ﬃ cient, the stochastic force makes the mode jump from one well to the other at random times, reproducing the phenomenon of direction switching. Experiments were conducted on a laboratory scale annular combustor fea-turing 12 hydrogen-methan ﬂames. System identiﬁcation techniques were used to ﬁt the model on the experimental data, allowing to extract the potential shape and the intensity of the stochastic excitation. The statistical predictions obtained from the Fokker-Planck equation on the mode’s behaviour and the direction switching time are in good agreement with the experiments.


Introduction
Thermoacoustic instabilities are a major issue in the development of new aeroengines combustion chambers aligned with the current ecological challenges because they induce vibrations that can severely damage the combustor and the turbine. The seriousness of this problem lies in the fact that despite decades of intense research, these instabilities are still very difficult to predict, which calls for further advances in their modelling and understanding [1]. The present study falls within this context and deals with the modelling of an intriguing phenomenon occuring in annular aeroengine combustors: the intermittent transitions between clockwise and counterclockwise spinning thermoacoustic modes. Indeed, in modern aeroengines, the combustor is annular, which has advantages in terms of weight, maintenance and manufacturing, and the instabilities frequently involve the first azimuthal eigenmodes of the annulus. Based on acoustic measurements in the combustors of practical gas turbines [2,3], in academic annular chambers [4][5][6] and on high fidelity numerical simulations [7], these self-sustained azimuthal thermoacoustic modes have been classified as spinning modes that propagate at the speed of sound in the clockwise or the anticlockwise direction, standing modes whose nodal line stays at a constant azimuthal position or drifts slowly compared to the speed of sound, or mixed modes that result from the sum of a standing mode and a spinning mode. Many studies dealing with the dynamics of these modes have been published over the last ten years, covering several aspects of the underlying physics such as the effects of spatial asymmetries, turbulent combustion noise, mean azimuthal flow or specific types of nonlinearities in the flame response, e.g. [8][9][10][11][12][13]. However, no detailed studies have been published on the intermittent transition events between mixed and standing modes and between clockwise-spinning and anticlockwise-spinning mixed modes. This paper aims at filling this gap with new experimental data and a loworder model, which is derived from first principles and describes the stochastic dynamics of these transitions in the same way as thermally activated barrier crossing in reaction rate theory [14]. The experimental setup and data are presented in section 2. Section 3 presents the low-order model. In section 4, the validity of this model is assessed using the experimental data, and its parameters are identified. Section 5 introduces the potential representation. The calibrated model is then used in section 6 to shed light on the mechanisms associated with the intermittent transitions of the spinning direction.  Figure 1a shows the experimental setup used in this study. It consists in an annular combustion chamber at atmospheric pressure, with 12 regularly spaced flames anchored on axisymmetric bluff bodies. The fuel is a gaseous mixture of CH 4 (contributing to 30% of the thermal power) and H 2 (70% of the thermal power) that are premixed in a plenum. The inner and outer walls of the chamber are cooled with water. The setup is described more in detail in [15]. Acoustic measurements are performed with 6 pressure transducers Kulite XCS-093-05D flush mounted in the burner's pipes at equidistant azimuthal positions. All the experiments are run with the same thermal power of 12×6 = 72 kW and equivalence ratios Φ comprised between 0.4 and 0.62. This is achieved by keeping the fuel mass flow constant at 320 sl/min and by varying the air mass flow between 2600 sl/min and 1670 sl/min. The mass flows are controlled with several Alicat mass flow controllers. Figure 1b shows the evolution of the acoustic amplitude for different equivalence ratios: The signature of a supercritical Hopf bifurcation is detected from the acoustic time traces when Φ is increased: below 0.4, the system is stable. The Hopf point is located at Φ = 0.51 : above this value, a thermoacoustic instability is detected. To ensure stationarity of the thermoacoustic dynamics, acoustic records are performed at least 60 s after ignition at each operating points so that thermal transient behavior of the combustor has ended. At this time, the measured inner wall temperature and acoustic root-mean-square amplitudes are stabilized. In this paper, we focus on three specific operating points to illustrate our conclusions: Φ = 0.525, Φ = 0.55 and Φ = 0.60. The equivalence ratios at these conditions are beyond the Hopf point and the associated acoustic spec- tra are dominated by one single peak close to 1.1 kHz, which corresponds to an acoustic eigenmode with one wavelength spanning the circumference of the annular chamber. Figure 2 shows portions of the corresponding acoustic pressure time series. One can see in Fig.  2b, 2c, 2d that the acoustic oscillations are almost sinusoidal even before applying any spectral filtering. Figure 2b corresponds to Φ=0.525, and one can conclude from the time traces that around t = 40 s, the acoustic oscillation corresponds to a standing mode. This is because all the microphone signals are in phase or exactly in phase opposition and exhibit different amplitudes. In Fig. 2c, the signals are not in phase and the amplitude depends on the angular position, which corresponds, at that instant, to a mixed mode. In Fig. 2d, the signals exhibit phase differences of about 120 • with similar amplitudes which also corresponds to a mixed mode, but with a stronger spinning component. At a given axial position, the pressure field associated to the first azimuthal mode is represented with the quaternion ansatz introduced in [16]:

Experiments
where i, j and k are the quaternion imaginary units. The real-valued azimuthal mode state variables ϕ, θ, A, and χ are slowly varying in time with respect to the acoustic period 2π/ω. They describe the slow dynamics of the mode [16,17]: ϕ is a slow temporal phase drift, θ is the direction of the antinodal line of the pressure oscillations and A is the oscillation amplitude. The nature angle χ characterizes the mode type: χ = 0 for a pure standing mode, χ = π/4 for a pure anticlockwise spinning mode, χ = −π/4 for a pure clockwise spinning mode, −π/4 > χ > 0 and 0 > χ > −π/4 for a anticlockwise and for a clockwise mixed mode respectively. The real pressure p(Θ, t) is obtained by taking the real part of the quaternion ansatzp. This quaternion representation describes the modal dynamics in a non-ambiguous way [16]. To obtain the time traces of the slow variables A, χ, θ and ϕ, we first write the pressure field associated to the first azimuthal mode with the ansatz p(Θ, t) = η 1 (t) cos(Θ) + η 2 (t) sin(Θ). Since the number of microphones is greater than 2, the signals η 1 and η 2 are obtained from the filtered microphones time traces by inverting an overdetermined system with a least-square method, with an error smaller than 5% for the present experiments. Then, the procedure to obtain the slow variables from the bivariate signal (η 1 , η 2 ) is described in [16]. time in the same direction. For Φ = 0.525 and Φ = 0.60, θ shows some jumps between two values distant of π, which are in fact describing the same antinodal line [16]. In Fig. 4a, one can see that ϕ does not show a preferential value, its PDF is flat. The time derivative of ϕ gives information on the frequency drift. The PDF oḟ ϕ shown in Fig. 4b suggests that this drift remains close to 0 and is small compared to the instability frequency. The nature angle χ behaves completely differently between the three cases: for Φ = 0.525, it fluctuates randomly around 0, which corresponds to a predominant standing mode state. For Φ = 0.55, χ switches randomly between a positive value and its negative, which indicates a statistically dominant mixed mode whose propagating direction changes randomly. For Φ = 0.60, χ is locked on a negative value close to −π/4, corresponding to an almost pure spinning mode propagating in the clockwise direction.

Theoretical model
We introduce now a low-order model that describes the dynamics of azimuthal thermoacoustic eigenmodes in an idealized one dimensional annular combustor. The four parameters of this model will be afterward identified from the experimental data presented in the previous section. It will be shown that the calibrated model quantitatively reproduce the dynamics of the slow-flow variables, which shows that it can used as a minimal model for describing the state of eigenmodes exhibiting an azimuthal component in real 3 dimensional combustors. It is important to stress that this model is not intended to be used for prediction of thermoacoustic stability in real systems. The chamber is modelled as a 1D annular waveguide of radius R. The azimuthal angular coordinate is referred to as Θ. The wave equation for the acoustic pressure in the chamber in presence of a thermoacoustic source term is where c is the speed of sound, α is the acoustic damping, γ is the heat capacity ratio, Ξ is a random forcing representing the turbulent heat release rate fluctuations, anḋ Q is the coherent component of the heat release rate of the flames. The latter is responsible of the thermoacoustic instability phenomenon, and it is modelled with a 3 rd order nonlinearity: where β[1 + c 2 cos(2Θ)]p influences the linear stability of the thermoacoustic modes and κp 3 governs the saturation to limit cycles [3,8,18,19]. To account for spatial asymmetries in the coherent heat release rate fluctuations, only the second term c 2 of the spatial Fourier expansion of the source term distribution was kept. Indeed, it was shown in previous works that it is the only term having an effect on the dynamics of the first azimuthal mode [8]. In [19], the quaternion formalism of Eq. (1) is introduced in the wave equation (2), and methods of spatial and slow time averaging are applied to obtain a first order system of coupled Langevin equations for A, χ, θ and ϕ. Based on the experimental data presented in Fig. 3, we consider that the fluctuations of the antinodal line direction are weak enough to set θ as a constant in the equations. Under these assumptions, the equations for A and χ are: where ν = (β − α)/2. The constant Γ is the intensity of the noise coming from the projection of the random fluctuation Ξ on the first azimuthal mode shape [19]. ζ A and ζ χ are uncorrelated gaussian noises of intensities Γ/2ω 2 . In the considered experimental time traces, A is almost constant and will be replaced by its mean value A 0 in the equation for χ. In the equation for A, χ appears in the terms cos(2χ) and cos(4χ). However, for the considered experimental conditions, the standard deviation of cos(2χ) never exceeds 0.15. Under the assumption that c 2 β/4 is of the same order as ν or smaller, the fluctuations of cos(2χ) are neglected and χ is replaced by the average value of its absolute value χ 0 = |χ| . The standard deviation of cos(4χ) in the experiments is always under 0.5, which has small effect in the sum (5 + cos(4χ)). Therefore, χ is also replaced by χ 0 in this This is a pre-print version. Published in Proceedings of the combustion institute (2020) Vol. 38, p. 5943-5951, DOI: doi.org/10.1016/j.proci.2020.05.008 term. The two equations are in this way decoupled: In a more compact formȦ = F A (A) + ζ A andχ = F χ (χ) + ζ χ /A 0 , where all the deterministic terms are grouped in the functions F A and F χ , which are the drift coefficients of the Fokker-Planck (FP) equations governing the evolution of the probability density P(A) and P(χ). The FP equations have also a diffusion coefficient. The latter is equal to Γ/2ω 2 in the FP equation for P(A).

Model parameters calibration
It is now interesting to make use of the experimental data to identify the model parameters, and therefore calibrate it for each of the stationary operating conditions investigated in this work. This can be achieved by following the framework proposed in [20,21], which has been already applied in a thermoacoustic context, e.g. [22]. It is based on the fact that one can compute from time series the drift and diffusion coefficients of a FP equation by using their statistical definition [20,21,23]: The noise intensity Γ is obtained with the second order moment for the time traces of A In Eq. 6, the transition probability P(a, t + τ|A, t) is the probability of the amplitude a at a time t + τ knowing that the amplitude at the time t was A. These probabilities are extracted for τ 1/∆ f with ∆ f the bandwidth of the filter applied on the microphones timetraces [24]. The limits for τ → 0 are obtained by extrapolation, as explained in [24]. The next step is to fit the theoretical model on the extracted drift and diffusion coefficients to identify the parameters Γ/ω 2 , κ and c 2 β. Although the noise intensity does not depend on the amplitude in the model, Figure 5a shows that the extracted diffusion coefficient presents a parabolic shape. This is due to limitations in the extrapolation method for the extreme values of A. However, a validation on numerically generated data has shown that a reliable estimation of the noise intensity Γ/2ω 2 can be obtained by weighting these extrapolated points with the experimental probability density function of A for the fit of Γ/2ω 2 . Once Γ/2ω 2 is known, ν and κ are identified using the extracted drift function F A (Fig 5b). With the knowledge of Γ/2ω 2 and κ, c 2 β is identified using F χ (Fig. 5c). As a final remark, it is important to note that modeling the nonlinear saturation with a cubic term only is too simplistic for the considered operating conditions, which are far away from the bifurcation point [25]. Consequently, although the identified values for ν, c 2 and κ allows to capture the topology of the phase space, they shall only be considered as indicative values.

Potential description
The equation for χ is now written in a potential form: where the potential U is defined as: Equation (9) is formally similar to the equation of a particle of negligible mass evolving in a potential landscape and undergoing a stochastic force ζ χ /A 0 and a strong fluid friction [14]. It has to be noted that in eq. (9), the noise intensity Γ affects not only the amplitude of the stochastic noise ζ χ , but also the shape of the potential landscape U. The expression (10) of the potential has three terms: the first term, proportional to A 2 0 cos(4χ), grows with the limit cycle amplitude A 0 and promotes spinning modes. The second term, proportional to c 2 β cos(2χ), grows with the asymmetry coefficient c 2 and promotes standing modes. The third term increases with the ratio Γ/A 2 0 , which depends on the relative importance of the noise compared to the limit cycle amplitude, it also promotes standing modes and creates potential barriers in ±π/4, preventing the occurence of pure spinning modes. These terms and their contribution to the potential shape are represented on Fig. 6. Depending on the relative importance of these three terms, two situations can be encountered: the one-well potential leading to a dynamic centered around standing modes, and the two-well potential leading to a dynamic centered around two mixed modes spinning in opposite directions. The first situation occurs when (3κ/16)A 2 0 ≤ c 2 β/2 + Γ 2 /(ω 2 A 2 0 ). Then, the contribution of the asymmetry and the deterministic contribution of the turbulent noise in U predominate over the term coming from the nonlinear saturation, and F χ has only one stable equilibrium point in χ = 0, corresponding to a one-well potential around the standing mode, corresponding to the blue line in Fig. 7b. The noise excitation will make χ wander randomly around its equilibrium point, but the dynamics will be centered around a standing mode. The second situation occurs when (3κ/16)A 2 0 > c 2 β/2 + Γ 2 /(ω 2 A 2 0 ). Then, the term coming from the nonlinear saturation in F χ (χ) predominates over the terms coming from the asymmetry and the deterministic contribution of the turbulent noise, and F χ has two symmetric stable equilibrium point ±χ eq corresponding to symmetric mixed modes propagating in opposite directions, and one unstable equilibrium point in χ = 0. This corresponds to a potential landscape featuring two wells in ±χ eq separated by a potential barrier in χ = 0. Figure 7b shows examples of potentials presenting this two-well shape (red and green lines). The noise excitation will make χ wander randomly around one of the mixed modes and eventually trigger random jumps from one well to the other, making the propagating direction of the wave switch randomly. For the cases showing two potential wells, the rate of jumps from one well to the other depends on the noise intensity. When χ is in one of the wells, Kramers theory will allow us to compute the probability that χ crosses the barrier after a certain time [14]. The probability density function of χ is governed by the Fokker-Planck equation which is equivalent to the Langevin equation (9). We assume that χ is initially in one of the wells, e.g the positive one, χ = χ eq . This is modelled with an initial condition P(χ, 0) = δ(χ−χ eq ) for the probability density. An absorbing boundary condition is placed in χ = 0. We then consider P c (t) = 1 − P(χ > 0, t), which is the probability that χ crossed the barrier in 0 at a time τ ≤ t. P c is initially 0 and increases with t because of the absorbing condition. P c (t) = t 0 ρ(τ)dτ with ρ the probability density function of the first passage time though the barrier. When the potential and the noise intensity are known, the Fokker-Planck equation can be solved numerically, giving access to ρ by differentiation of the cumulative distribution function P c . the nature angle between the three considered cases. Fig. 7b shows the potential landscapes associated with the calibrated model, whose parameters where identified from experimental data for each of the three cases, as explained in the previous section. The Φ = 0.525 case (blue line) shows a single potential well in χ = 0, so that the stochastic forcing leads to a predominantly standing mode. The Φ = 0.55 (red line) and Φ = 0.60 cases (green line) both present two symmetric potential The turbulent noise is therefore strong enough to drive χ from one well to the other, causing the observed phenomenon of intermittent transitions between clockwise and counterclockwise mixed modes. Conversely, the Φ = 0.60 case (green) presents a slightly smaller level of noise, due to the stronger amplitude, and a higher potential barrier. The noise level is small compared to the height of the potential barrier, making any switching event very unlikely: χ will remain confined in one well, corresponding to a constant spinning direction. The choice of the well depends only on the initial conditions. Fig. 7c compares the theoretical stationary PDFs against the experimental ones. It has to be noted that for the Φ = 0.60 case (green), the solution of the stationary Fokker-Planck equation gives a bimodal distribution while the experiments show a single peak. This is because the solution of the stationary Fokker-Planck equation gives the repartition obtained for an infinite time. The duration of the experiment is too short to observe any switching event. For Φ = 0.55, the probability density function of switching time ρ has been obtained numerically and compared to experimental statistics, showing good agreement (Fig. 8).

Conclusion
A low order model was used to describe the dynamics of the first azimuthal mode in an idealized annular combustion chamber subject to turbulent combustion noise. Considering a wave equation with coherent and stochastic source terms, we obtain an equation describing the modal dynamics as a stochastic diffusion process in a potential landscape. This equation gives a model for the spontaneous symmetry breaking that is observed on a real combustor, and which leads to either a clockwise spinning or a counter-clockwise spinning mode for sufficiently large values of the bifurcation parameter. The decomposition of the model's potential landscape allows us 1) to attribute the predominance of a spinning to the nonlinear saturation of the flame, while turbulent forcing and spatial asymmetries favor predominance of standing modes, and 2) to explain the intermittent transitions induced by the turbulent forcing when the potential barrier separating the two counter-spinning attractors has a height that is similar to the normalized stochastic forcing intensity. The model was able to reproduce quantitatively the statistics of these intermittent transitions, showing that it is a suited minimal model for describing the complex topology of the phase space.