Symmetry breaking modelling for azimuthal combustion dynamics

Annular combustors can exhibit azimuthal thermoacoustic instabilities, which can rotate as a spinning wave at the speed of sound in the azimuthal direction, oscillate as a standing wave with pressure nodes ﬁxed in space, or be a linear combination of these. These oscillations happen if a positive feedback loop between acoustics and the response of the ﬂames to the acoustics in the annulus occurs. This paper discusses how two di ﬀ erent explicit symmetry breaking mechanisms a ﬀ ect the dynamics of these waves. We ﬁrst show how small di ﬀ erences between the ﬂame responses lead to one strong topological change in the dynamical system phase space, making the system prefer orientation angles at two azimuthal locations, one opposite of the other in the annulus, as found in the experiments. This symmetry breaking is modelled by directly perturbing the ﬂame responses around the annulus with some scatter, to represent the e ﬀ ect of manufacturing tolerances of the burners. We then consider recent experimental evidence that the heat release rate of the ﬂames depends on the spinning direction (clockwise or anticlockwise) when the system is spinning. In particular we model one experiment in which the ﬂame response is found to be stronger when the wave rotates in the anticlockwise direction. We show that the statistics of the resulting model are qualitatively very similar to the experimental results showing a preference for spinning states in the anticlockwise direction.


Introduction
Azimuthal instabilities are typical of symmetrical combustion systems, either as unwanted acoustic oscillations, like in annular combustors [1], or as part of the combustion concept, as in rotating detonating combustion [2] .A common approach in the modelling of these instabilities assumes that the flame responds, by means of the acoustic admittance of the system upstream of the flames, to the local value of acoustic pressure 1 at the flame position [1].Under this assumption, in the fully symmetric case one proves in absence of noise that spinning solutions always exist as stable periodic attractors, and that when standing periodic solutions exist the N 2n criterion can be used to ascertain their stability [5].This is in agreement with the experiments for both spinning [5] and standing solutions [6].In presence of noise, as (Giulio Ghirardo) 1 See [3,4] for works where it is instead the azimuthal velocity that drives the flames is typical of many industrial applications and of turbulent flames, one proves that the system state is driven away from spinning solutions, also in agreement with the experiments [7].
The first part of this paper lifts the assumption of full symmetry, and investigates how a slight loss of discrete symmetry affects the system.This has been previously studied theoretically in the deterministic case, for a simple cubic flame response, and for a homogeneous heat release distribution along the annulus [8].We will lift all these assumptions except fixing a certain flame response, and show how the pressure antinodes of the azimuthal instability are pushed towards two azimuthal orientations in the annulus, as found in the experiment.
The second part of the paper considers recent advances in the response of the flame to also the azimuthal acoustic velocity.Theoretical works have shown that the nonlinear response of the flame to azimuthal velocity fluctuations can affect the dynamics of the solution [9] and push the system towards standing states.Later experimental works have confirmed this effect, which originates either by a loss of mirror symmetry of the  flame [10,11] or by the complex interaction of neigbouring flames [12] or by both.This paper models this effect and shows how we can recover a statistical preference of the acoustic field for anticlockwise solutions as in the experiment.
We review from the literature the governing equations and the experiment in §2.We then discuss the baseline, reference axisymmetric case in §3.We consider the first symmetry breaking mechanism, related to manufacturing imperfections, in §4, and the second mechanism in §5.We draw the conclusions in §6. .The radius A describes the amplitude of acoustic oscillation, the nature angle 2χ describes whether the system is spinning (at the poles) or standing (on the equator) and the angle nθ 0 describes the location of the pressure antinode of the standing component of the instability [13].

The model
In this section we briefly recall the features of interest of the experiment and recapitulate the governing equations from the literature.We study a laboratory annular combustor with bluff-body swirling flames that has been discussed in detail in [14][15][16].In the following we will consider a variant of this setup with M = 18 burners, as sketched in Fig. 1 and described recently in [12], operated at an equivalence ratio of 0.75 and a bulk flow velocity U b =20m/s at the burner exit.At this condition, a top view of the measured mean chemiluminescence is presented in Fig. 2, at which the system exhibits thermoacoustic pulsations of azimuthal order n = 1 at a frequency of approximately 1.7 kHz.The acoustic pressure field in the annular combustor, at the cross-section just downstream of the burners can be written as [13]: where θ is the azimuthal coordinate defined in Fig. 1, the three variables (A, nθ 0 , χ) depend on the time t and are described in Fig. 3, and ϕ also depends on time and is the temporal phase of the instability.The resulting governing equations are obtained by substituting the ansatz (1) into the fluctuating mass and momentum conserva-tion equations [7]: 2π 0 e i2n(θ−θ 0 ) e kχ + e −kχ Q θ (A p (θ))dθ e kχ (2) µ and are described in the rest of this section.Equation ( 2) is quaternion-valued, so that i, j, k are distinct imaginary units [7,17].The left hand side of (2) describes the evolution of the four variables of interest {A, nθ 0 , χ, ϕ} as a function of time.The real part describes the the evolution of the amplitude A, the imaginary i and j units encode spatial and temporal phase, while k encodes the nature angle χ.The right hand side describes the vector field that the four variables follow, which will be presented in the figures later.In particular the second line in (2) is the contribution over the annulus of Q θ , which is the projection of the describing functions of the flames and of the acoustic losses on the azimuthal mode [18]: where the last term describes the acoustic damping of the system, assumed linear, by means of the coefficient α .In the simplest cases the flame responds to the axial acoustic velocity measured just upstream of the flame position.It is possible to express such velocity as function of the acoustic pressure in the combustion chamber by mean of the acoustic impedance of the whole system just upstream of the flame [5,19].As a consequence, the response of the flames Q θ,flames depends on the local amplitude A p (θ) of acoustic pressure at the azimuthal location θ [1], which can be calculated as Exploiting the fact that each flame is acoustically compact, we can rewrite Q θ,flames as: where δ(x) is the Dirac distribution, and θ m is the location of the m-th burner as presented in Fig. 1, and the projected flame responses Q m = Q m,r + jQ m, j may be different from one another, and are assumed here to be real-valued for simplicity, as in many other studies [1,18,20].Under these assumptions the thermoacoustic oscillation frequency ω matches the natural oscillation frequency ω 0 and the first term in the last row of (2) cancels out.Ending the analysis of (2), the last two terms on the third line depend on the intensity σ of the background noise.These terms arise because of the random HRR fluctuations of the turbulent flame [21,22] and are modelled here as stochastic [23], with µ being a quaternion-valued additive white gaussian noise process.A detailed description of the model ( 2) and its derivation can be found in [7].

Baseline: the axisymmetric case with a simple nonlinear saturation in absence of noise
Because the describing functions of the flames are not available, we choose a simple saturation model with monotonic response for the axisymmetric case.In particular we constrain the flame response such that in the linear regime it is equal to a linear effective gain β m : where for the axisymmetric case the flames have the same constant effective gain β m = β sym .Since the flame transfer functions of the flames are not known, we fix first the linear growth rate of the system azimuthal instability to the value of ν 0 /ω 0 = 0.16 in absence of acoustic damping, on the high end of common values in the literature discussed in [19].The resulting value of the effective gains β m = β sym , equal for all the flames in this symmetric case, can then be calculated as β sym /ω 0 = (2ν 0 /ω 0 )/(πM).We then add a level of acoustic damping α equal to β m /5.The flame strength and the acoustic damping are kept fixed in the following, are common tunable parameters in low-order models for azimuthal instabilities and are typically identified in applications.Most of the following results have been tested also for different values, and should apply as long as β sym > α to make the system linearly unstable.The three variables {A, nθ 0 , 2χ} are the spherical coordinate of a point in Fig. 3, which describes the system trajectory as function of time t, driven by the vector field described by the right hand side of (2).This threedimensional vector field is hard to present as a whole, therefore we present next two projections of this vector field on two planes: the equatorial plane of the sphere in Fig. 3, and on one vertical plane passing through both poles.We present on the two rows of Fig. 4 these two projections for the axisymmetric case of this section, in absence of noise, as discussed next.On each row, in the first frame, the component of the vector field pushing the system state in the radial direction is presented, so that in red regions the state is driven outwards towards infinity, and in blue regions the state is driven towards the origin.In the second frame, the component of the vector field acting in the tangential direction is presented, so that red regions push the system counterclockwise and blue regions clockwise.In the third frame, the same information is presented in terms of the streamlines in the plane.Repellors in the plane are marked with empty circles and attractors in the plane are marked with filled circles.The same structure of this Figure is used later in Figs. 6, 9, 10, and the axes of the frames share the same radius limit at A = 3.2, so that they can be compared.In all these frames the azimuthal coordinate starts at three o'clock as in Fig. 1, and in the vertical cuts only the right half is presented, because in this plane it is always symmetric.In the equatorial plane of Fig. 4 we observe how the streamlines are almost parallel to radii passing through the origin, which is due to the azimuthal component of the field being significantly weaker than its radial counterpart.This will play a role when we perturb the symmetry in the next section.In the vertical plane, we observe that for this specific choice of flame saturation the spinning solution is an attractor of the system.This is not always the case, because also standing solutions can exist and be attractors, depending on how the flame saturates in the nonlinear regime [5].In these figures, higher values of the flame effective gain β push the solutions further away from the origin.Similary, higher values of the acoustic damping α push the system towards the origin.

The effect of a small loss of rotational symmetry
Due to small manufacturing scatter, small differences in the mass-flow rates of air and fuel are expected between the flames.In thermoacoustics, these differences affect the flame response.One way to account for this is to perturb to the effective flame gains β m , m = 0, . . ., M − 1 some additive white gaussian noise scalars with a small standard deviation, here chosen equal to 1/10th of the baseline value β sym .One realization of the resulting slightly asymmetric pattern is presented in Fig. 5 and is considered next.Of the pattern, only the mean value β sym and the 2n-azimuthal component have an effect on the system dynamics 2 , so that their sum is presented in the Fig. 5 with red markers as an equivalent pattern.The 2n = 2 locations where this equivalent pattern is maximum are presented in Fig. 5 at θ = θ (2n) ≈ 3π/4 and at θ = θ (2n) + π.The resulting vector field is presented in Fig. 6.The dynamics in 2 this can be observed in the governing equations (2) and discussed in [7] in the general nonlinear case and in [8] for the linear term Figure 4: Axisymmetric case for an annular combustor with M = 18 burners exhibiting an azimuthal mode of order n = 1.All the symbols are introduced in Fig. 3.In the top 3 frames, a cut on the equator plane of Fig. 3 is presented.In the bottom 3 frames a cut passing through the north and south poles of the sphere of Fig. 3 is presented.On each row the first frame presents the component of the vector field in the radial direction A, the second frame presents the component in the tangential direction, and the last frame presents the streamlines of the resulting 2D vector field.In the first two columns the contour levels where the components are zero are presented with black lines.These two sets of black lines are reported also in the third column, where their intersections are the fixed points of the system.We draw with a filled circle fixed points that are attractors in the plane, and with an empty circle points that are either saddles or repellors.the equatorial plane are quite different from the axisymmetric case of Fig. 4: only 4 fixed points survive, and of these 2 are saddles and 2 are attractors in the plane.The 2 attractors are azimuthally located approximately at the angle θ (2n) corresponding to to the 2n maxima of the linear coefficients of Fig. 5 just described.
This result matches qualitatively the experimental evidence presented in Fig. 7, and past evidence from the literature: Fig. 10 of Worth & Dawson [14], a similar pdf of the orientation angle presented by Ghirardo & Bothien [13] in their Fig. 6 after the addition of the dampers, and Fig. 6 of Bourgouin et al. [24].In all four cases there are two azimuthal locations (actually one location, and the same location rotated by π) where the antinode of the standing component of the pressure field is more likely to be found.These results apply regardless of whether the system is noisy or not, i.e regardless of the intensity σ of the background noise.Technically this is the case because from the equations this strong change in the topology of the deterministic vector field in the equatorial plane does not depend on σ.
Of the many vertical planes, we present in Fig. 6 the cut at the angle nθ 0 = θ (2n) .In this vertical plane the solutions that were perfectly spinning in the axisymmetric case of Fig. 4 are now pushed slightly off from the poles, but still spinning to a very good approximation, consistently with [8].We also find that the vector field is very close to the axisymmetric case, suggesting that existing theoretical results and the N 2n stability criteria for standing solutions [5] derived for perfectly symmetric configurations are robust against small asymmetries like this one.

The effect of the nature angle on the HRR response
Saurabh et al. [10,11] force a single flame with simultaneous axial and transverse excitation as in annular combustors.They show that the gain and the phase of the flame response to a spinning wave can significantly depend on whether the spinning wave rotates clockwise or anticlockwise in the annulus.They attribute this to a loss of reflection symmetry on a plane orthogonal to the direction of the spinning wave passing through the burner axis.Nygård et al. [12] force the annular combustor of Fig. 1 with loudspeakers so that the acoustic field is to a good approximation spinning clockwise or anticlockwise.They measure the integrated HRR and show that it differs between the two cases, as rediscussed here in Fig. 8.We can model this in the equations, by making the effective gains of the flames β m Figure 6: Effect of a small loss of rotational symmetry on an annular combustor with M = 18 burners exhibiting an azimuthal mode of order n = 1.All the symbols are introduced in Fig. 3, and the structure of the plot is the same of Fig. 4. The loss of symmetry considered here is modelled as some random scatter in the response of the flames, as presented in Fig. 5.In comparison to the axisymmetric case of Fig. 4, there are now just 4 fixed points in the xy plane, of which two are repellors and two are attractors in the plane.The vertical plane is very similar to the axisymmetric case.depend on the nature angle χ.We propose to write: to model the experimental evidence that when the system is spinning anticlockwise (2χ close to π/2) the HRR The resulting linearly stable system was then forced with the loudspeakers to push the system state to a spinning wave, first in the anticlockwise direction and then in the clockwise direction, approximately at the same level of acoustic amplitude.One period of the fluctuating HRR Q − Q of one burner representative of all others is then reconstructed for the two cases, normalized by the mean HRR Q.The timeseries of the two cases have been adjusted to start at zero because it was not possible to synchronize them with the occurrence of the acoustic pressure maximum at the same burner position.Of the two states, the flame response is larger when the system is spinning anticlockwise.More information can be found in [12].
Figure 9: Effect on the system dynamics of the dependence of the HRR on the nature angle, for γ = 0.14.Only the vertical plane is presented because the equatorial cut matches exactly the axisymmetric case of Fig. 4. The addition of noise is considered next in Fig. 10 response is larger than when the system is spinning clockwise (2χ close to −π/2).The coefficient γ models the strength of this effect.For this analysis the effect of ( 7) on the HRR is linear in the amplitude of oscillation A, because γ directly perturbs the gains β m , and we fix γ = 0.14 to match the experimental evidence that the HRR for the clockwise forced case is 23% smaller than the one for the anticlockwise forced case, as presented in Fig. 8.The resulting flow and discussion are presented in Fig. 9.In comparison to the baseline case of Fig. 4, the spinning solution in the northern hemisphere has moved to a larger amplitude, and the spinning solution in the southern hemisphere has moved to smaller amplitudes.To match the noisy experimental data, the background noise is now switched on in the model by Figure 10: Same as in Fig. 9, with the same limits on the colormaps, but now accounting also for the effect of the background noise σ, set here to 0.09.Simulation results of this case are presented in Fig. 11.tuning σ = 0.09, affecting the flow as in Fig. 10.Noise pushes the two attractors away from the poles, as proved by [7].The value of σ is chosen to match the statistics of the experiment.We run a simulation starting from an anticlockwise initial condition, and present in Fig. 11 a comparison with the experimental data.We observe how the two probability density functions qualitatively agree, have a similar width and tail to the lower values, but the simulation peaks at a location closer to the north pole at π/2.Regarding the effect of the noise intensity σ, we mention that in the numerical model it is possible to shift the position of the maximum while maintaining the solution in the northern hemisphere by increasing both the noise intensity σ and the coefficient γ, or both σ and the linear growth rate of the system ν 0 introduced in §3.We did not attempt a systematic optimization study in the space {ν 0 , σ, γ} to improve this result.
Other plausible reasons for the maximum pdf peak shift include a partial loss of symmetry as discussed in §4, which is not accounted for in this section but hinted at by the experimental result of Fig. 7, and the nonlinear saturation of the flame responses, which is at the moment unknown and modelled according to (6).
Figure 12: From an initial condition in the southern hemisphere close to the clockwise spinning solution, the system escapes to the north hemisphere in this specific simulation after about 7 500 periods, which is about 4.7 seconds for the experiment considered.
Another matter of interest is the fact that there exist two distinct attractors in Fig. 10.We find numerically that a point initialized close to the anticlockwise solution does not escape the northern hemisphere in the same time span measured in the experiment.Similarly, we find that a point initialized close to the clockwise solution escapes the southern hemisphere after a time that is small compared to the observation time of the experiment, as presented in Fig. 12.Both observations are typical of systems with a double well potential, where one well is deeper than the other.The double well is then called asymmetric, and the asymmetry is described in our case by the coefficient γ introduced in (7).The solution in the deeper well is the anticlockwise solution, which is more likely to be observed in the model and that is observed in the experiment, while in the shallower well we find the clockwise solution, which is not observed in the experiment.The two probabilities of escaping from one basin to the other are different.

Conclusions
We model an annular combustor with M = 18 burners, with a simple flame saturation model that neglects the HRR not in phase with the acoustic pressure.We show how a small loss of rotational symmetry pushes the pressure antinodes of the first order azimuthal instability towards two preferential azimuthal locations, one location and the same location rotated by π, as found in two different laboratory experiments and in one industrial scale engine.
We then model recent experimental evidence showing that the heat release response depends on the direction of the spinning wave, i.e. on the nature angle 2χ.For the annular experiment considered, the HRR is found to be larger for a wave rotating in anticlockwise direction than in clockwise direction.We model this effect and find that it pushes the anticlockwise solution to larger amplitudes and the clockwise solution to smaller amplitudes, while mantaining both solutions as attractors in the deterministic case, i.e in absence of background noise.When noise is considered, we present numerical evidence of how it is possible that the system escapes from the basin of attraction of the clockwise solution because of the noise, and that this can happen in a matter of seconds, and then easily go undetected in the experiment.We present numerical evidence that it is harder for the system to escape from the basin of attraction of the anticlockwise solution, which is the only one measured in the experiments.A comparison of the statistics of the nature angle 2χ between model and experiment shows a good qualitative agreement.

Figure 1 :
Figure 1: View of the annular rig from the top.The acoustic pressure field depends in this plane on the azimuthal angle θ and on the time t.The M = 18 burners, sketched with two concentric circles, are located at equispaced azimuthal angles θ m = 2πm/M, m = 0, . . ., M − 1.

Figure 2 :
Figure 2: Mean normalized chemiluminescence from the top.The M = 18 darker round regions correspond to the bluff-bodies of the burners.

Figure 3 :
Figure 3: Poincaré sphere representation of an azimuthal instability of order n, presented in(1).The radius A describes the amplitude of acoustic oscillation, the nature angle 2χ describes whether the system is spinning (at the poles) or standing (on the equator) and the angle nθ 0 describes the location of the pressure antinode of the standing component of the instability[13].

Figure 5 :
Figure 5: Burner to burner scatter can affect the response of the flames in an annular combustor.These slight changes can affect the effective gains of the flames, presented with black squares, which depart from the baseline symmetric value, presented with gray circles.

Figure 7 :
Figure 7: Probability density function of the orientation angle nθ 0 for the self excited case.The system lingers in the vicinity of two preferred orientations, one opposite of the other.

2 Figure 8 :
Figure8: Loudspeakers were connected to the walls of the combustion chamber by means of ducts designed to damp the thermoacoustic instability.The resulting linearly stable system was then forced with the loudspeakers to push the system state to a spinning wave, first in the anticlockwise direction and then in the clockwise direction, approximately at the same level of acoustic amplitude.One period of the fluctuating HRR Q − Q of one burner representative of all others is then reconstructed for the two cases, normalized by the mean HRR Q.The timeseries of the two cases have been adjusted to start at zero because it was not possible to synchronize them with the occurrence of the acoustic pressure maximum at the same burner position.Of the two states, the flame response is larger when the system is spinning anticlockwise.More information can be found in[12].

Figure 11 :
Figure 11: Comparison between simulations and experiments of the statistics of the nature angle 2χ.The experimental data refers to the self-excited case, ran at the same operating conditions of Fig. 8 but without loudpspeakers.Both histograms describe approximately 52 000 periods of acoustic oscillations.