Unified approach to low-order modeling of can-annular combustors

Tiemo Pedergnana 1 and Nicolas Noiray 2

CAPS Laboratory, Department of Mechanical and Process Engineering, ETH Zürich Sonneggstrasse 3, 8092 Zürich, Switzerland

ABSTRACT We present a unified, physics-based low-order modeling approach to can-annular combustors, combining a linear stability analysis and a nonlinear Langevin approach to unravel the highly complex, multi-physical dynamics in this type of system. In the linear part, Howe’s Rayleigh conductivity model is combined with the projected Helmholtz equation and a Bloch wave ansatz to arrive at a single equation for the frequency spectrum of an N-can combustor. Starting from first principles, we illustrate and give a physical explanation for the coupling-induced amplification and suppression of thermoacoustic instabilities. Adding another layer of complexity, we then take into account nonlinearities in the flame response and the aeroacoustic coupling. The nonlinear dynamics of a symmetric model combustor are explored by exploiting the gradient structure of the averaged slow-flow dynamics. We obtain exact analytical expressions for the steady-state statistics, highlighting the connection between di ff erent emergent patterns and the resistive coupling between the cans. We leverage our analysis to explain the intermittent energy transfer between Bloch modes observed in real-world gas turbines.

1. INTRODUCTION

Thermoacoustic instabilities (TIs) arise from the constructive interaction between the turbulent flame and the sound field in the (rigidly enclosed) combustor [1]. If not su ffi ciently damped, the resulting instability can lead to high-amplitude limit cycles. In industrial machines, such as stationary gas turbines, the resulting pulsations induce high-cycle fatigue cracks in the metal parts surrounding the enclosure, which cause down-time, fees, and potentially catastrophic system failures. In the 21st century, the modeling of TIs has drawn renewed interest due to emission limits on power generation systems and the resulting increased demand for lean-premixed combustion. Lean premixed flames are highly sensitive to flow perturbations, which increases critically their response to acoustic forcing. Also, the lack of dilution holes in a premixed combustor leads a significant reduction in acoustic damping. These two factors make combustion instability problems in lean-premixed gas turbines similar to those encountered in rocket engines. In both applications, these instabilities remain a key problem in the design process [2,3]. There is ample literature on TIs in silo-type (one single can) or annular combustors (see, for instance, [4–7] and [8–10], respectively). In contrast, this work presents physics-based low-order

1 ptiemo@ethz.ch

2 noirayn@ethz.ch

| dite! al nllaieey, inter. noise 21-24 AUGUST SCOTTISH EVENT CAMPUS O ? e GLASGOW

modeling techniques for can-annular gas turbines, complementing the theoretical, numerical and experimental studies performed at Siemens [11–14], General Electric [15–17] and Ansaldo Energia Switzerland [18, 19]. While the largest can-annular gas turbines produce more than 750 MW in combined cycles at more than 62% e ffi ciency, and can be supplied with sustainably produced H 2 (hydrogen) blended with natural gas, the silo-type and annular designs become uneconomical at this scale due to several engineering trade-o ff s. This report is based on a physics-based linear stability analysis [20] and its subsequent extension, a nonlinear deterministic and stochastic time-domain analysis [21]. The linear stability analysis was preceded and motivated by a theoretical and experimental study on the aeroacoustic response of a whistling side branch, a simple example of a self-sustained system [22]. A sketch of the modeled system is given in Figure 1 (a), showing the control volume C j (not to scale), the flames and the acoustic pressure at the turbine outlet of the j th can, η j . As depicted in Figure 1 (b), the full system can be reduced to a ring of coupled oscillators with amplitudes A j and phases ϕ j . The coupling arises from the acoustic-hydrodynamic interaction in the aperture, shown in Figure 1 (c). At high bulk flow speed U tot , the turbulent wake is idealized as two vortex sheets separated by the can spacing d . This work considers the weakly coupled regime, where both the can spacing d and the aperture width W are small compared to the can length L .

2. LINEAR STIBILITY

The linear stability analysis starts from the Helmholtz equation in the frequency domain [23]:

∂ 2 ˆ p ( j ) ( s , x )

! 2 ˆ p ( j ) ( s , x ) = − s γ j − 1

∂ x 2 − s

c 2 j ˆ Q ( j ) ( s , x ) in V j , (1)

c j

ˆ p ( j ) ( s , x ) ˆ u ( j ) ( s , x ) · n = Z ( j ) ( s , x ) on σ j . (2)

In Equation 1 and Equation 2 , ˆ p ( j ) and ˆ u ( j ) denote the acoustic pressure and velocity in the j th can, s = ν + i ω is the Laplace variable, where ω and ν are the angular frequency and growth rate of thermoacoustic oscillations at a frequency f = ω/ 2 π , respectively, i is the imaginary unit, γ j and c j are the specific heat ratio and the ambient speed of sound in the j th can, respectively, n is the outward facing normal vector to the boundary σ j , Z ( j ) is the acoustic impedance on σ j and ˆ Q ( j ) is the unsteady heat release rate fluctuations per unit volume in the flame region. By symmetry, we set γ j ≡ γ and c j ≡ c ∀ j . We assume that the TI is dominated by a single natural eigenmode ψ k with corresponding modal amplitude ˆ η ( j ) , k and eigenfrequency ω k . Expanding the acoustic pressure in terms of ψ k yields the projected Helmholtz equation:

ˆ p ( j ) ( s , x ) = ˆ η ( j ) , k ( s ) ψ k ( x ) . (3)

ˆ u ( j ) ( s , x ) = − ˆ η ( j ) , k ( s ) ∇ ψ k ( x )

s ρ , (4)

where ˆ η ( j ) , k is the dominant modal amplitude defined by

ˆ η ( j ) , k = s ρ c 2

γ − 1

1 V j Λ j

Z

V j ˆ Q ( j ) ( s , x ) ψ k ( x ) dV

s 2 + ω 2 k

ρ c 2

| ψ k ( x ) | 2

Z ( j ) ( s , x ) dS ! , j = 1 , ..., N , (5)

Z

− ˆ η ( j ) , k

σ j

where Λ j = 1 / V j R

C j | ψ k | 2 dV is the mode normalization factor and V j = Vol( C j ) is the volume of C j . By symmetry, V j ≡ V , Λ j ≡ Λ and Z ( j ) ≡ Z ∀ j . For simplicity, we drop the subscript k and assume a

first row of

W η j

turbine vanes

A a

can outlet

C j can combustion

x L

d

chamber . . . . . .

q j

flame

0

(a)

cross-section A

burner outlet

A j − 1 µ

ϕ j ( t ) A j ( t )

λ

U tot

0

resistive coupling

mean

reactive coupling

axial velocity turbulent

wake

vortex sheet

(b) (c)

d

Figure 1: (a) Sketch of the modeled can-annular combustor, showing the control volume C j , the acoustic pressure at the can outlet η j , the aperture width W , cross-section area A a and depth d , respectively, the can length L and cross section area A , respectively, as well as the coherent heat release fluctuations of the flame q j , which characterize the flame response to acoustic perturbations. (b) Ring of coupled oscillators with amplitudes A j and phases ϕ j with linear resistive coupling λ and the reactive coupling µ (neglected in this study). (c) Acoustic-hydrodynamic interaction in the aperture between the cans, the mean bulk flow speed of the combustion products U tot and the turbulent wake bounded by two vortex sheets. This figure was adopted from [21] with minor changes

pressure antinode at the turbine inlet and a proportional flame response with zero time delay. For the coupling, we use Howe’s model of the aeroacoustic response of a rectangular aperture [24]. Starting from the unsteady Bernoulli equation, the following equation can be derived: Z 1

! Z 1

h

− 1 ζ ′ ( µ ) G ( ξ, µ ) d µ + ( α + + α − ξ ) e iSt c ξ = 1 , (6)

− 1 ζ ′ ( µ ) { ln | ξ − µ | + L + ( ξ, µ ) } d µ − π St 2 c

W

where ξ is the normalized streamwise variable, ζ ′ is the normalized vortex sheet displacement, St c = ω c W / 2 U is the Strouhal number based on the complex frequency ω c = ω + i ν , ν is the growth rate, µ is an integration variable corresponding to ξ , α ± are constants of integration,

G ( ξ, µ ) = − H ( ξ − µ )( ξ − µ ) e iSt c ( ξ − µ ) , (7)

where H ( · ) is the Heaviside function and

(2 B / W ) 2 + ( ξ − µ ) 2 } + p

1 + ( W / 2 B ) 2 ( ξ − µ ) 2 − ( W / 2 B ) | ξ − µ | . (8)

L + ( ξ, µ ) = − ln { 2 B / W + p

From the solution of Equation 6 , the Rayleigh conductivity K R can be computed using

Z 1

K R ( ω c ) = − π B

− 1 ζ ′ ( µ, ω c ) d µ. (9)

2

From Equation 9 , the acoustic impedance Z = ρ cZ s is readily obtained:

Z s = [ ˆ p ]

ρ c ˆ u a = − sA a

cK R . (10)

Assuming a perfectly symmetric system, we drastically simplify our analysis by making a Bloch-wave ansatz, expressing the modal amplitudes of neighboring cans as follows:

ˆ η ( j + 1) = ˆ η ( j ) e − i 2 π b

N and ˆ η ( j − 1) = ˆ η ( j ) e i 2 π b

N , (11)

where b is the Bloch wavenumber. The associated discrete azimuthal waves are called Bloch modes. This ansatz reduces the problem of determining the linear stability of the thermoacoustic system to a single complex equation:

s 2 − 2 ν 0 s + ω 2 k − b 0 K R ( s ) sin 2 ( π b / N ) = 0 , (12)

where b 0 = 16 R c 2 / V > 0. Using the above results, a linear stability analysis is performed in [20]. In that study, it is observed that various geometric and acoustic parameters influence the stability of the system, leading to coupling-induced modal damping or instability, depending on the parameter regime. A physical explanation of the overall phenomenon of coupling-induced changes to the frequency spectrum is given in Figure 3 : the instantaneous pressure di ff erences at the coupling interfaces drive the aeroacoustic coupling which can damp or enhance TIs. In Figure 3 (a) and (b), b = 4, and the instantaneous pressure di ff erence is greater than in Figure 3 (c) and (d) ( b = 1), which leads to an instability for the higher-order Bloch mode, but not for the lower-order one. This interpretation is key to understanding that the can-to-can coupling has no e ff ect when all the cans are in phase, and that its e ff ect is the strongest for the push-pull mode, which features a phase di ff erence of π between neighboring cans.

3. NONLINEAR AND STOCHASTIC DYNAMICS

Motivated by the rich linear dynamics induced from the aeroacoustic coupling, the above model is extended to the time domain. The system is represented by an amplitude-phase system of 12 coupled oscillators at a fixed frequency ω 0 , which considers only resistive (damper-like) coupling. Expressing the acoustic pressure η j as η j = A j cos ( ϕ j + ω 0 t ), ˙ η j = − ω 0 A j sin ( ϕ j + ω 0 t ), the nonlinear and stochastic dynamics are given by the following set of Langevin equations:

˙ A j = − ∂

∂ A j V ( A m , ϕ m ) + ζ j , (13)

∂ ∂ϕ j V ( A m , ϕ m ) + χ j , (14)

A j ˙ ϕ j = − 1

A j

where j , m = 1 , . . . , N and the potential V is defined as follows:

N X

 − ν A 2 l 2 + κ A 4 l 32 − λ A l

2  A l + 1 cos( ∆ l + 1 ) − A l  − ϑ A l

h ( A 3 l + 1 + A 2 l A l + 1 ) cos( ∆ l + 1 )

V ( A m , ϕ m ) =

8

l = 1

2 ( A 2 l + 2 A 2 l + 1 ) i − Γ 4 ω 2 0 ln( A l )  , (15)

− A l A 2 l + 1 2 cos(2 ∆ l + 1 ) − A l

Bloch mode with b = 4

norm. acoustic pressure

1

ˆ η j / ˆ η j , max ˆ η j / ˆ η j , max ξ ξ

0

(a)

− 1

norm. vortex sheet displacement Re  ζ ′ ( t ) 

1

0

(b)

− 1

Bloch mode with b = 1

1

0

(c)

− 1

7

1

0

(d)

− 1

1 2 3 4 5 6 7 8 9 10 11 12 can number j

Figure 2: Acoustic-hydrodynamic interaction corresponding to a coupling-induced instability for the second set of parameters in [20]. Vanishing can spacing d = 0 was assumed. Shown in (a) and (b) are the normalized acoustic pressure distribution at a given time instant and the real part of the normalized vortex sheet displacement Re  ζ ′ e ( − i ω + ν ) t  at 4 equally spaced points in time during an acoustic cycle with period T = 2 π/ω , respectively, for a Bloch mode with b = 4. Shown in (c) and (d) is the same for a Bloch mode with b = 1. The normalized frequency spectrum ( ω/ω k , ν/ω k ) is (0 . 93 , 9 . 3 × 10 − 3 ) for b = 4 (coupling-induced instability) and (0 . 995 , − 1 . 06 × 10 − 2 ) for b = 1 (stable thermoacoustic mode)

where λ is the linear and κ the nonlinear resistive coupling, ∆ j = ϕ j − ϕ j − 1 is the phase di ff erence and ζ j and χ j are 2 N uncorrelated white noise sources of intensity Γ / 2 ω 2 0 . Equations 13 and 14 , which are derived in [21] from stochastic and deterministic averaging of a set of coupled Van der Pol equations, imply that the trajectories ( A j , ϕ j ), j = 1 , . . . , N , are stochastically perturbed and attracted towards lower values of the potential V defined in Equation 15 . The phase di ff erences of the fixed points, or Bloch modes, of the Langevin system are found to be

∆ ∗ b = 2 π b

N , (16)

where b is the Bloch wavenumber. It is shown in [21] that dissipative coupling generally favors low-order, and amplifying coupling high-order Bloch modes. Figure 3 shows an example of a nonlinear deterministic time trace for the parameters listed in [21] with G = κ Γ / 8 | ν | 2 ω 2 0 = 0. In this example, the resistive coupling is dissipative ( λ = ν/ 4) and the system converges to a low-order Bloch mode with ∆ ∗ − 1 = − 2 π/ N . For t →∞ , in the steady state, the probability density function of the system can be deduced from the stationary Fokker–Planck equation:

∂ x P ∞ ( x ) + Γ 4 ω 2 0

 ∂ V ( x )

∂ P ∞ ( x )

 T , (17)

0 = ∂

∂ x

∂ x

where the newly introduced variables x = ( U 1 , V 1 , . . . , U N , V N ) are related to ( η j , ˙ η j ) by η j = A j cos φ j

Bloch mode with b = − 1 π

2 A ∗ 0

2 A ∗ 0

phase di ff . ∆ j

amplitude A j

ac. press. η j

A ∗ 0

0

0

0 4995 4995 4995 Cycles Cycles Cycles 5000 5000 5000 − 2 A ∗ 0

− π

1 . 2 A ∗ 0

2 A ∗ 0

2 π

η j

A j

ϕ j

− 1 . 2 A ∗ 0

0

0

Figure 3: Example of a nonlinear deterministic time trace computed from Equations 13 and 14 for the parameters listed in [21] with G = 0. In this example, the resistive coupling is dissipative ( λ = ν/ 4) and the system converges to a low-order Bloch mode with ∆ ∗ − 1 = − 2 π/ N

and ˙ η j = − ω 0 A j sin φ j , where

A j = q

U 2 j + V 2 j (18)

φ j = arctan( V j , U j ) + ω 0 t . (19)

The solution to Equation 17 is

P ∞ ( η m , ˙ η m ) = N exp ( − 4 ω 2 0 V ( A m ( η m , ˙ η m ) , ϕ m ( η m , ˙ η m )) / Γ ) . (20)

Equation 20 implies that regions of low values of the potential V correspond to a high probability density P ∞ in the steady state. Figure 4 shows an example of a nonlinear stochastic time trace for the parameters listed in [21] with G = 2 . 88 × 10 − 1 . In this example, the resistive coupling is amplifying ( λ = − ν/ 4) and the system is attracted to the push-pull mode with ∆ j = ± π . In contrast to the deterministic system shown in Figure 3 , in the steady state, the system does not converge exactly to any Bloch mode. Instead, we observe intermittent energy transfer between di ff erent Bloch modes, reproducing results obtained from real-world gas turbines [18]. A more detailed discussion of the stochastic dynamics is given in [21].

4. CONCLUSIONS

This work presents a unified, physics-based approach to low-order modeling of can-annular combustors, covering linear, nonlinear and stochastic dynamics. Although not entirely general, the presented models can be generalized to an arbitrary degree, enabling realistic representations of experiments or simulations. The analysis above demonstrates that

– The aeracoustic coupling between the cans is driven by the instantaneous pressure di ff erence across the coupling apertures between neighboring cans

stochastic push-pull mode π

2 A ∗ 0

2 A ∗ 0

phase di ff . ∆ j

amplitude A j

ac. press. η j

A ∗ 0

0

0

0 4995 4995 4995 Cycles Cycles Cycles 5000 5000 5000 − 2 A ∗ 0

− π

1 . 2 A ∗ 0

2 A ∗ 0

2 π

η j

A j

ϕ j

− 1 . 2 A ∗ 0

0

0

Figure 4: Example of a nonlinear stochastic time trace computed from Equations 13 and 14 for the parameters listed in [21] with G = 2 . 88 × 10 − 1 . In this example, the resistive coupling is amplifying ( λ = − ν/ 4) and the system is attracted to the push-pull mode with ∆ ∗ 6 = ± π

– In the time domain, the nonlinear dynamics are well approximated by a Langevin system. The associated potential gives an accurate description of the dynamics in the deterministic limit

– Under the addition of noise, intermittent energy transfer between Bloch modes is observed. The Langevin potential predicts which state is most favored in the steady state

As an extension of the present framework, quantitative models may be derived in analogous fashion by using empirically calibrated versions of the coupling and the coherent flame response. The results presented here add to the informed design of novel mitigation measures against TIs in can-annular gas turbines.

REFERENCES

[1] J. W. S. Rayleigh. The explanation of certain acoustical phenomena. Nature , 18(455):319–321, 1878. [2] J.J. Keller. Thermoacoustic oscillations in combustion chambers of gas turbines. AIAA J. , 33(12):2280–2287, 1995. [3] T. Poinsot. Prediction and control of combustion instabilities in real engines. Proc. Combust. Inst." , 36(1):1–28, 2017. [4] L. Crocco. Aspects of combustion stability in liquid propellant rocket motors part I: fundamentals. low frequency instability with monopropellants. J. Am. Rocket Soc. , 21(6):163– 178, 1951. [5] J.J. Keller, W. Egli, and J. Hellat. Thermally induced low-frequency oscillations. Z. Angew. Math. Phys. , 36(2):250–274, 1985. [6] T. Schuller, T. Poinsot, and S. Candel. Dynamics and control of premixed combustion systems based on flame transfer and describing functions. J. Fluid Mech. , 894, 2020.

[7] W. Polifke. Modeling and analysis of premixed flame dynamics by means of distributed time delays. Progress in Energy and Combustion Science , 79, 2020. [8] N. Noiray and B. Schuermans. On the dynamic nature of azimuthal thermoacoustic modes in annular gas turbine combustion chambers. Proc. R. Soc. A , 469(2151), 2013. [9] G. Ghirardo and M.P. Juniper. Azimuthal instabilities in annular combustors: Standing and spinning modes. Proc. R. Soc. A , 469(2157), 2013. [10] Abel Faure-Beaulieu, Thomas Indlekofer, James R. Dawson, and Nicolas Noiray. Imperfect symmetry of real annular combustors: beating thermoacoustic modes and heteroclinic orbits. J. Fluid Mech. , 925:R1, 2021. [11] S. Bethke, W. Krebs, P. Flohr, and B. Prade. Thermoacoustic properties of can annular combustors. In 8th AIAA / CEAS Aeroacoustics Conference and Exhibit , pages 1–11, 2002. [12] Werner Krebs, Sven Bethke, Joachim Lepers, Patrick Flohr, Bernd Prade, Cli ff Johnson, and Stan Sattinger. Thermoacoustic Design Tools And Passive Control: Siemens Power Generation Approaches , pages 89–112. AIAA, 2005. [13] P. Kaufmann, W. Krebs, R. Valdes, and U. Wever. 3D thermoacoustic properties of single can and multi can combustor configurations. In Proceedings of the ASME Turbo Expo , volume 3A, pages 527–538, 2008. [14] F. Farisco, L. Panek, and J.B.W. Kok. Thermo-acoustic cross-talk between cans in a can-annular combustor. Int. J. Spray Combust. , 9(4):452–469, 2017. [15] K. Venkatesan, A. Cross, C. Yoon, F. Han, and S. Bethke. Heavy duty gas turbine combustion dynamics study using a two-can combustion system. In Proceedings of the ASME Turbo Expo , volume 4A, pages 1–9, 2019. [16] K. Moon, H. Jegal, C. Yoon, and K.T. Kim. Cross-talk-interaction-induced combustion instabilities in a can-annular lean-premixed combustor configuration. Combust. Flame , 220:178–188, 2020. [17] K. Moon, C. Yoon, and K.T. Kim. Influence of rotational asymmetry on thermoacoustic instabilities in a can-annular lean-premixed combustor. Combust. Flame , 223:295–306, 2021. [18] G. Ghirardo, C. Di Giovine, J.P. Moeck, and M.R. Bothien. Thermoacoustics of can-annular combustors. J. Eng. Gas Turbine Power , 141(1), 09 2019. [19] G. Ghirardo, J. P. Moeck, and M. R. Bothien. E ff ect of Noise and Nonlinearities on Thermoacoustics of Can-Annular Combustors. J. Eng. Gas Turbine Power , 142(4), 01 2020. [20] Pedergnana T. and Noiray N. Coupling-induced instability in a ring of thermoacoustic oscillators. Proc. R. Soc. A , 478(2259), 2022. [21] Pedergnana T. and Noiray N. Steady-state statistics, emergent patterns and intermittent energy transfer in a ring of oscillators. Nonlinear Dyn. , 2022. [22] T. Pedergnana, C. Bourquard, A. Faure-Beaulieu, and N. Noiray. Modeling the nonlinear aeroacoustic response of a harmonically forced side branch aperture under turbulent grazing flow. Phys. Rev. Fluids , 6(2), 2021. [23] George C. Maling. Simplified analysis of the Rijke phenomenon. J. Acoust. Soc. Am. , 35(7):1058–1060, 1963. [24] M. Howe. Influence of wall thickness on Rayleigh conductivity and flow-induced aperture tones. J. Fluid Struct. , 11(4):351–366, 1997.