Low order modelling of thermoacoustic instabilities in aero-propulsion engines

Aswathy Surendran Department of Engineering Physics and Computation, Technical University of Munich, 85747 Garching, Germany School of Chemical and Physical Sciences, Keele University, Sta ff ordshire ST5 5BG, United Kingdom

Charles Boakes Department of Mechanical Engineering, Imperial College London, London SW7 2AZ, United Kingdom

Dong Yang Department of Mechanics and Aerospace Engineering, Southern University of Science and Technology, Shenzhen, 518055, PR China

Aimee S. Morgans 1

Department of Mechanical Engineering, Imperial College London, London SW7 2AZ, United Kingdom

ABSTRACT In the present work, we investigate the thermoacoustic stability behaviour of a combustor-heat exchanger assembly pertinent to the Synergetic Air Breathing Rocket Engine (SABRE) of Reaction Engines Ltd. Since the thermoacoustic behaviour of heat exchangers is poorly understood, we have approximated it to a combination of heat transfer / heat sink and aeroacoustic scattering (dissipation) mechanisms. The ϵ - NTU method is used to characterise the low frequency, linearised unsteady heat transfer (heat exchanger transfer function) behaviour in compact heat exchangers. The aeroacoustic response of the tube row is evaluated from the linearised conservation equations, wherein the hydrodynamics of the bias flow in the tube row is described as being similar to that of an isentropic contraction followed by a sudden expansion. The stability predictions are carried out through low order modelling of the combustor and subsequent eigenvalue analysis. The flame dynamics is approximated using the n − τ law, with tau being a parameter. Stability predictions show that the thermoacoustic response of the heat exchanger was dominated by aeroacoustic dissipation. Only one mode was predicted to be unstable and can be stabilised by moving the heat exchanger upstream. Further work is required to better characterise the flame and the heat exchanger transfer functions.

1 a.morgans@imperial.ac.uk

a slaty. inter.noise 21-24 AUGUST SCOTTISH EVENT CAMPUS O ¥, ? GLASGOW

1. INTRODUCTION

Thermoacoustic instability occurs when there is positive feedback loop or interaction between the acoustic perturbations of the combustion system and the unsteady heat release rate fluctuations, leading to low frequency high amplitude pressure fluctuations that can be detrimental to the combustion system [1]. In the present work, we investigate the thermoacoustic behaviour of the pre-burner assembly in SABRE (Synergetic Air Breathing Rocket Engine), which is the world’s first pre-cooled, air-breathing rocket engine that is being developed by Reaction Engines Ltd. U.K. [2]. SABRE (see Figure 1) has the potential to facilitate single stage to orbit space vehicles and the next generation of supersonic passenger aircraft.

Pre-burner

Ramjet Rocket engine

Nacelle

Supersonic air intake cone

Nozzle

Compressor Turbine

SABRE

Pre-cooler

Figure 1: Schematic of SABRE (Reproduced from www.reactionengines.co.uk ).

There are two modes of operation for SABRE: (i) air-breathing mode (Mach 5.4) where liquid H 2 is burnt with air from the surrounding atmosphere and (ii) rocket mode where the fuel is burnt with the stored oxidiser. In air-breathing mode, the incoming air decelerates from supersonic to subsonic speeds after passing through a shock array at the supersonic inlet cone. This sudden deceleration generates significant amount of heat, which is then transferred to a helium cycle through the pre-cooler heat exchanger assembly. The heat thus removed is then used to drive the turbo-compressor unit. At low air-speeds, when the heat generated is not su ffi cient to drive the air compressor, a secondary combustor called pre-burner (Figure 2) is utilised. In addition to a heat source, the pre-burner also has a heat exchanger to supply heat to the helium cycle.

Round duct Matrix burner Heat exchanger

Choked outlet

H 2 flames

Inlet

Rectangular duct Transition duct

Figure 2: Schematic of a pre-burner combustor with heat exchanger.

Heat exchangers can act as both heat sinks and acoustic dampers (due to scattering) and therefore they can influence the thermoacoustic behaviour of combustion units that incorporate them. We

expect that the heat exchanger assembly inside the pre-burner could be used as a potential mitigator of thermoacoustic modes. The structure of the paper is as follows: the modelling aspects of the pre-burner is provided in Section 2. This section describes the modelling of individual acoustic components of the pre-burner and is followed by the section on stability predictions (Section 3). The stability predictions and associated discussion are provided in Section 4, with the conclusions given in Section 5.

2. MODELLING OF THE PRE-BURNER

The acoustic field (given in Equation 1) inside the pre-burner is assumed to consist forward and backward travelling planar pressure waves denoted by ˆ p + and ˆ p − [3] respectively.

p ′ = ˆ p + e i ω ( t − x / ( c + ¯ u )) + ˆ p − e i ω ( t + x / ( c − ¯ u )) , (1)

The flow properties such as pressure, velocity, temperature and density are assumed to be formed of a mean part, denoted by an overbar (¯), and a small fluctuating part, denoted by a prime ( ′ ), i.e.,

p = ¯ p + p ′ ; u = ¯ u + u ′ ; T = ¯ T + T ′ ; ρ = ¯ ρ + ρ ′ . (2a-d)

The fluctuating quantities when expressed in terms of the acoustic pressure waves ˆ p ± 1 , 2 , are given by

p ′ 1 , 2 = ˆ p + 1 , 2 + ˆ p − 1 , 2 ; u ′ 1 , 2 = ˆ p + 1 , 2 − ˆ p − 1 , 2 ¯ ρ 1 , 2 c 1 , 2 ; ρ ′ 1 , 2 =

 ˆ p + 1 , 2 + ˆ p − 1 , 2 c 2 1 , 2

 − s ′ 1 , 2 ¯ ρ 1 , 2

c p 1 , 2 (3a-c)

where c is the speed of sound, c p is the specific heat at constant pressure, s ′ 1 is the incoming entropy wave in region 1 and s ′ 2 is the outgoing entropy wave in region 2. T ′ is expressed in terms of p ′ and ρ ′ . The entropy waves have the form s ′ = ˆ s e i ω ( t − x / ¯ u ) (4)

where t is the time and x is the location at which the quantity is computed.

The pre-burner shown in Figure 2 can be split into its individual acoustic components as shown in Figure 3. The pre-burner can be modelled as a network of acoustic components like inlet / outlet, straight and area varying ducts, the flame and the heat exchanger. Each of these components can be described through a transfer matrix [4] that relates the upstream and downstream acoustic fields. If the upstream quantities are denoted by subscript 1 and the downstream quantities are denoted by 2, the transfer matrix that relates the up- and downstream acoustic and entropy field can be written as

ˆ p + 2 ˆ p − 2 ˆ s 2

ˆ p + 1 ˆ p − 1 ˆ s 1











 (5)

T 11 T 12 T 13 T 21 T 22 T 23 T 31 T 32 T 33

When we have multiple components, the acoustic fields upstream (denoted by 1) and downstream (denoted by N ) can be described by a single transfer matrix [ T ] 1 → N which is obtained by multiplying the transfer matrices of individual components, as in Equation 6

ˆ p + N ˆ p − N ˆ s N

ˆ p + 1 ˆ p − 1 ˆ s 1

ˆ p + 1 ˆ p − 1 ˆ s 1

 (6)



 = [ T ] 1 → N



 = [ T ] N − 1 → N · · · [ T ] 2 → 3 [ T ] 1 → 2



The inlet and the outlet boundaries are described by reflection coe ffi cients that relate the upstream and downstream travelling waves. The reflection coe ffi cient at the inlet is given by Equation 7,

Heat exchanger

Choked outlet

Inlet Matrix burner

Straight duct

Transition duct

Straight duct

Area varying

Inlet Duct FTF

HTF Duct Outlet

Duct

Figure 3: Components of the pre-burner.

whereas the reflection coe ffi cient at the downstream end is given by Equation 8. Since entropy waves are convected downstream, the reflection at the downstream end will have contribution from entropy waves.

ˆ p + 1 = R 1 , p ˆ p − 1 (7) ˆ p − N = R N , p ˆ p + N + R N , s ˆ s N (8)

In the present study, we assume that the inlet is a rigid wall with R 1 = 1 and the outlet is acoustically compact and choked. Subsequently, we use the reflection coe ffi cients proposed by Marble and Candel [5].

R 2 , p = 1 − ( γ − 1) M / 2

1 + ( γ − 1) M / 2 (9)

R 2 , s = − M / 2 1 + ( γ − 1) M / 2 (10)

The transfer matrices for the straight and area varying ducts can be found by linearising the steady conservation equations for mass, momentum and energy equations, substituting the relations in Equation 3a-c and rearranging the terms to form transfer matrices. For example, the transfer matrix for the straight duct would be

ˆ p + 2 ˆ p − 2 ˆ s 2

e − i ω L / ( c + ¯ u ) 0 0

ˆ p + 1 ˆ p − 1 ˆ s 1











 (11)

0 e i ω L / ( c − ¯ u ) 0

0 0 e − i ω L / ¯ u

2.1. Heat Sources and Sinks The heat source i.e., the flame and the heat sink i.e., the heat transfer contribution from the heat exchanger are assumed to be acoustically compact. In order to obtain the transfer matrices for the flame and the heat sink, we linearise the conservation equations shown in Equations 12-14

ρ 1 u 1 = ρ 2 u 2 (12)

p 1 + ρ 1 u 2 1 = p 2 + ρ 2 u 2 2 (13) ρ 1 u 1 H 1 A + ˙ Q = ρ 2 u 2 H 2 A (14)

where H is the stagnation enthalpy, A is the cross-sectional area and ˙ Q is the unsteady heat release rate / heat transfer rate fluctuations depending on whether it is a flame or a heat exchanger. Typically, the fluctuations in ˙ Q are related to the fluctuations in the acoustic velocity. This relation is defined

through a transfer function (Equation 15) in the frequency domain. For the flames, this is called a flame transfer function (FTF) [6].

T F = ˆ˙ Q ′ / ¯˙ Q ˆ u ′ re f / ¯ u re f (15)

where ˆ [ ] ′ quantities are the frequency domain equivalents of the perturbed quantities.

In practical situations, perturbations in the unsteady release rate may be caused by perturbations in flame speed, mixture density, or flame area, or by acoustic interactions fuel mixing or injection, amongst other behaviours. A mechanistic analysis of the pre-burner flame is outside the scope of the present work and hence we assume an approximated FTF as per available literature. Often, it is observed that the phase of the FTF is linearly related to the frequency and that the gain exhibits a low pass filter behaviour [7–9]. A simplest heat release rate model that captures this observation is the n − τ model, where the heat release rate fluctuations are related to the upstream acoustic velocity fluctuations with a time-lag of τ i.e.,

˙ Q ′ ( t ) = n u ′ re f ( t − τ ) (16)

Therefore, the FTF for the pre-burner flame can be written as

FTF = A 1 1 + i ω

ω c e − i ωτ (17)

Since no data is available for the pre-burner flame, we assume f c = ω c / 2 π = 400Hz, τ ϵ [0.5, 5]ms and A = 1.

2.2. Heat Exchanger model Heat exchanger act both as heat sinks and acoustic scatterers. Due to lack of data and literature for modelling heat exchangers as thermoacoustic components, we approximate the heat exchanger to the geometry shown in Figure 4. Here, we assume that all of the heat is transferred across the heat sink that is positioned prior to the solid cylindrical part of the heat exchanger. The heat sink is described through a transfer function (HTF) following the description in Section 2.1 and the geometry is described as a combination of an isentropic contraction, followed by a straight duct and an abrupt expansion. All these components can be expressed in terms of their individual transfer matrices and the detailed derivation can be found in [10].

Isentropic

contraction

Abrupt expansion Heat source/sink

Figure 4: Approximated geometry for the heat exchanger.

The HTF for the heat sink is derived from the commonly used ϵ − NTU method [11] as shown by Boakes [10]. The heat transfer rate fluctuations for the heat sink, thus obtained is given by

! 1 − ¯ ϵ

Q ′

¯ Q = ρ ′

¯ ρ + u ′

h ′

¯ h − ρ ′

¯ ρ − u ′

¯ u + ln 1 1 − ¯ ϵ

! (18)

¯ ϵ

¯ u

where

C min ; h ′

¯ u + ρ ′

¯ h = 0 . 6 u ′

! ; (19a-c)

ϵ = 1 − e − NTU ; NTU = hA ht

¯ ρ

h is the heat transfer coe ffi cient, A ht is the surface area over which heat transfer occurs and C is the minimum of the hot ( h ) and cold ( c ) heat capacity rates given by C h , c =  ˙ mc p 

h , c .

Heat transfer is a slow process compared to combustion and therefore it is reasonable to approximate this as a quasi-steady process in the low frequency ranges. Substituting Equation 19a-c into Equation 18 and after neglecting density perturbations, we get the low frequency limit of HTF as

! 1 − ¯ ϵ

HTF f → 0 = 1 − 0 . 4 ln 1 1 − ¯ ϵ

¯ ϵ (20)

3. STABILITY ANALYSIS

Stability analysis of the pre-burner assembly is conducted using OSCILOS [12] which is an open source tool implemented in MATLAB for studying combustion instability in longitudinal and annular combustor geometries. This tool enables users to quickly build network models of combustion systems with simplified treatment of geometric features, flames and boundary conditions. Linear eigenvalues analysis is done in the frequency domain and the results obtained are in terms of the eigenfrequencies and mode shapes of acoustic modes of system.

Equations 6-8 are written in a matrix form in terms of the unknown variables: ˆ p + 1 , ˆ p − 1 , ˆ p + N , ˆ p − N , ˆ s N (ˆ s 1 = 0), as shown in Equation 21. The eigenvalues, or eigenfrequencies in the present case, are computed by letting the determinant of the square matrix to zero. The eigenfrequency is a complex quantity ω = ω R − i δ whose real part gives the natural frequency of the acoustic mode and the imaginary part gives the growth rate. For the sign convention used in this study, the acoustic mode is said to be stable if δ is negative and it is unstable if δ is positive. 







ˆ p + 1 ˆ p − 1 ˆ p + N ˆ p − N ˆ s N

T 11 T 12 − 1 0 0

T 21 T 22 0 − 1 0

= 0 (21)

T 31 T 32 0 0 − 1

1 − R 1 , p 0 0 0

0 0 − R N , p 1 − R N , s

4. RESULTS AND DISCUSSION

The schematic of the pre-burner assembly under consideration is shown in Figure 2. It consists of an inlet, a section of duct terminated by the matrix flame, a duct of transitioning cross-section terminated by the heat exchanger (two rows of tubes separated by a small gap), and a section of duct terminated by a choked outlet. The full description of the dimensions and the values of the various flow and heat exchanger properties used in the present analysis can be found in Boakes [10]. Only plane propagating duct modes below the cut-o ff frequency ( f c / o = c / 2 l ) are considered. f c / o corresponding to the duct between the heat exchanger and the choked outlet is evaluated to be 1100Hz when c = 770m / s and duct length l = 350mm.

The eigenfrequencies and growth rates of the first nine acoustic modes of the pre-burner are shown in Figure 5 for τ in the range [0.5,5]ms. The modes were identified by the number of nodes and anti- nodes in the mode shapes of the velocity and pressure perturbations. Not all modes were observed

for every value of τ . It can be observed that the modal growth rates are very sensitive to τ , but only Mode 8 is unstable for certain values. Mode 2 and 5 can be considered as borderline cases.

mode 9 stable unstable

cut-off

mode 8

mode 6

mode 7

mode 5

mode 4

mode 2

mode 3

mode 1

Figure 5: Frequencies and growth rates of the acoustic modes of the pre-burner at various τ in the range [0.5,5]ms.

The instability of Mode 8 can be explained by quantifying the sources and sinks of acoustic energy. This is done by taking a cyclic average of the acoustic energy generated by a component, averaged over one time period [13] as given by Equation 22.

S = I

S ¯ I . n dS (22)

where I is the intensity or flux of acoustic energy. The cycle averaged generation rates, normalised by the magnitude of the generation rate of the flame, are presented in Table 1 for Modes 2, 5 and 8. The sources / sinks considered are the flame, the heat transfer and aeroacoustic interaction at the heat exchanger, and the choked outlet. Acoustic sources are positive and acoustic sinks or dampers are negative. For τ = 1.5ms, the net acoustic energy values indicate that Modes 2 and 5 are stable while Mode 8 is unstable.

Table 1: Table of normalised cycle averaged generation for modes 2, 5 and 8 at τ = 1.5ms.

S flame / | S flame | S hx , thermal / | S flame | S hx , aero / | S flame | S outlet / | S flame | S total / | S flame |

Mode 2 1.0 7.4 × 10 − 4 -9.9 × 10 − 1 -2.0 × 10 − 1 -1.8 × 10 − 1

Mode 5 -1.0 1.5 × 10 − 4 -5.4 × 10 − 1 -9.5 × 10 − 2 -1.6

Mode 8 1.0 -7.1 × 10 − 7 -3.0 × 10 − 2 -1.4 × 10 − 1 8.3 × 10 − 1

It is apparent that the thermal generation or dissipation of the heat exchanger is negligible compared to its aeroacoustic dissipation. Moreover, the aeroacoustic damping in Mode 8 is an order

of magnitude smaller than those for Modes 2 and 5. This behaviour can be explained by the velocity mode shapes presented in Figure 6. Acoustic dissipation is driven by the velocity perturbations and for Modes 2 and 5 the heat exchanger is at a location where the velocity perturbations has large amplitudes, whereas for Mode 8, this is a position close to a node. Hence, there is very small dissipation. To confirm this, we move the heat exchanger 100mm upstream and the new mode shapes are shown in Figure 7. Table 2 presents the frequencies and growth rates with and without the proposed modification. The e ff ect on the frequencies is small, but the growth rates are di ff erent. Notably, Mode 8 is stabilised, and Mode 2 is closer to the stability margin.

Table 2: Table of Mode 2, 5 and 8 frequencies and growth rates for τ = 1.5ms.

Original Modified

f R [Hz] δ [rad / s] f R [Hz] δ [rad / s]

Mode 2 359 -50.7 357 -5.68

Mode 5 633 -301 632 -401

Mode 8 1020 108 1000 -19.3

flame HX3 choked outlet

Figure 6: Mode shapes for modes 2, 5 and 8, for τ = 1.5ms.

Figure 7: Mode shapes for modes 2, 5 and 8, for τ = 1.5ms. Here the heat exchanger is moved 100mm upstream.

5. CONCLUSIONS

The stability analysis of the pre-burner assembly in SABRE was carried out using the open source tool OSCILOS. To achieve this, the pre-burner was split into its individual acoustic components that were described in terms of transfer matrices which relate the upstream and downstream acoustic and entropy fields. There were uncertainties involved in the modelling of the transfer functions for the flame and the heat exchanger. This was resolved by describing the flame dynamics through the n − τ law and the heat transfer part of the heat exchanger through the ϵ − NTU method. The aeroacoustic part of the heat exchanger was resolved by deriving transfer functions from the linearised conservation equations. Eigenfrequencies and growth rates for nine modes of the pre-burner were evaluated and it was found that for certain τ values, Mode 8 was unstable. Furthermore, it was observed that the aeroacoustic damping contribution of the heat exchanger dominated over the heat transfer contribution. From the mode shapes of Mode 8, one could infer that the acoustic damping from the heat exchanger had negligible influence since the heat exchanger was positioned close to a node. To stabilise Mode 8, the heat exchanger was moved away from the node and the proposed modification could stabilise Mode 8. However, this change pushed Mode 2 closer to the unstable regime.

ACKNOWLEDGEMENTS

We gratefully acknowledge the financial support from the European Research Council (ERC) Consolidator Grant AFIRMATIVE (2018–2023, Grant Number 772080). Dr. Surendran is a EuroTechPostdoc fellow at Technische Universität München, co-funded by the European Commission under its framework programme Horizon 2020. Grant Agreement number 754462.

REFERENCES

[1] T. C. Lieuwen and V. Yang. Combustion Instabilities in Gas Turbine Engines: Operational Experience, Fundamental Mechanisms, and Modeling , volume 210 of Progress in Astronautics and Aeronautics . American Institute of Aeronautics and Astronautics, Inc., 2005. [2] Reaction Engines Ltd. https://reactionengines.co.uk . Accessed: 2022-05-03. [3] A. P. Dowling and S. R. Stow. Acoustic Analysis of Gas Turbine Combustors. Journal of Propulsion and Power , 19(5):751–764, 2003. [4] N. A. Cumpsty and F. E. Marble. The interaction of entropy fluctuations withturbine blade rows: a mechanism of turbojet engine noise. Proceedings of the Royal Society A , 357(1690):323–344, 1977. [5] F. E. Marble and S. M. Candel. Acoustic disturbance from gas non-uniformitiesconvected through a nozzle. Journal of Sound and Vibration , 55(2):225–243, 1977. [6] L. Crocco and S. Cheng. Theory of combustion instability in liquid propellant rocket motors. Technical Report 8, AGARDO-GRAPH by Butterworths Scientific Publications, 1956. [7] D. Kim, J. G. Lee, B. D. Quay, D. A. Santavicca, K. Kim, and S. Srinivasan. E ff ect of flame structure on the flame transfer function in a premixed gas turbine combustor. Journal of Engineering for Gas Turbines and Power , 132(2):021502, 2010. [8] F. Boudy. Nonlinear dynamics and control analysis of combustion instabilitiesbased on the Flame Describing Function (FDF) . PhD thesis, EcoleCentrale Paris, 2012. [9] A. Chaparro, E. Landry, and B. M. Cetegen. Transfer function characteristics ofblu ff -body stabilised conical v-shaped premixed turbulent propane-air flames. Combustion and Flame , 145(1):290–299, 2006.

[10] C. Boakes. Predicting Instabilities in a Rocket Engine Preburner. Master’s thesis, Imperial College London, 2019. [11] F. P. Incropera, D. P. Dewitt, T. L. Bergman, and A. S. Lavine. Foundations of Heat Transfer . John Wiley and Sons, 6th edition, 2013. [12] J. Li, D. Yang, C. Luzzato, and A. S. Morgans. Open source combustion instability low order simulator (OSCILOS). Technical report, Imperial College London, 2017. [13] C.L. Morfey. Sound transmission and generation in ducts with flow. Journal of Sound and Vibration , 14(1):37–55, 1971.