Numerical prediction of combustion instability limit cycle oscillations for a combustor with a long flame

A coupled numerical approach is investigated for predicting combustion instability limit cycle characteristics when the combustor contains a long flame. The test case is the ORACLES combustor, with a turbulent premixed flame a metre long: it exhibits limit cycle oscillations at ∼ 50 Hz and normalised velocity amplitude ahead of the flame of ∼ 0.29. The approach obtains the flame response to acoustic excitation using Large Eddy Simulations (LES), and couples this with a low-order wave-based network representation for the acoustic waves within the combustor. The flame cannot be treated as acoustically compact; the spatial distribution of both its response and the subsequent effect on the acoustics must be accounted for. The long flame is uniformly segmented axially, each segment being much shorter than the flow wavelengths at play. A series of “local” flame describing functions, one for the heat release rate response within each segment to velocity forcing at a fixed reference location, are extracted from the LES. These use the Computational Fluid Dynamics toolbox, OpenFOAM, with an incompressible approximation for the flow-field and combustion modelled using the Partially Stirred Reactor model with a global onestep reaction mechanism. For coupling with the low-order acoustic network modelling, compact acoustic jump conditions are derived and applied across each flame segment, while between flame segments, wave propagation occurs. Limit cycle predictions from the proposed coupled method agree well with those predicted using the continuous 1-D linearised Euler equations, validating the flame segmentation implementation. Limit cycle predictions (frequency 51.6 Hz and amplitude 0.38) also agree well with experimental measurements, validating the low-order coupled method as a prediction tool for combustors with long flames. A sensitivity analysis shows that the predicted limit cycle amplitude decreases rapidly when acoustic losses at boundaries are accounted for, and increases if combustor heat losses downstream of the flame are accounted for. This motivates more accurate determination of combustor boundary and temperature behaviour for thermoacoustic predictions. © 2017 The Authors. Published by Elsevier Inc. on behalf of The Combustion Institute. This is an open access article under the CC BY license. ( http://creativecommons.org/licenses/by/4.0/ )


Introduction
Lean premixed combustion is preferred in modern aero-engines and industrial gas turbines due to its effect in reducing NO x emissions [1] . However, lean premixed systems are highly susceptible to combustion instabilities, which originate from the coupling between acoustic disturbances and unsteady heat release rate fluctuations in the combustor [2,3] . Combustion instabilities may lead to early ageing of the combustion chamber, or, in extreme cases, to severe structural damage [2,3] . The ability to predict, and if necessary suppress, these thermoacoustic instabilities during the early design stage of a combustor would offer enormous benefits, but this still constitutes a challenge due to the complex mechanisms and combustor geometries involved [2,4] .
Approaches for numerically simulating and analysing combustion instabilities generally fall within two categories. The first involves direct numerical calculation of the coupled acoustic oscillations and unsteady flame heat release within the combustor, via complete 3D compressible Computational Fluid Dynamics (CFD) simulations [5] . For example, investigations were conducted on self-excited azimuthal modes using LES in a full scale helicopter combustion chamber [6] .
The second, much less time-consuming, approach is to decouple the calculations of the perturbed flame response and the acoustic waves. The unsteady heat release rate from the flame resulting from acoustic excitation is characterised via a flame transfer function (FTF) for linear analysis [7,8] or a flame describing function (FDF) for weakly nonlinear analysis [9,10] . It can be obtained from experiments [11][12][13] , analytical models [7,14,15] or numerical simulations [16][17][18][19] . The generation, propagation and reflection/transmission of the acoustic waves are captured by either a low-order acoustic network model or a Helmholtz solver. Both exploit the fact that the acoustic wave behaviour remains linear in gas-turbine combustors, even during limit cycle oscillations [20] . Low-order acoustic network models represent the combustor geometry as a network of simple geometry modules, each sustaining acoustic waves superimposed on a mean flow. The acoustic wave behaviour is assumed low-dimensional typically just longitudinal and circumferential waves, and the mean flow variables are assumed constant within each geometry module. Neighbouring modules are connected by acoustic transfer matrices [21][22][23][24] . The alternative, Helmholtz solvers, assume zero mean flow velocity and numerically solve the wave equation in the frequency domain (the Helmholtz equation) for general geometry and temperature variations [25][26][27] . They must account for convective effects via corrections.
Both low-order acoustic network models and Helmholtz solvers typically make the assumption that the flame can be considered compact compared to the wavelength of the dominant acoustic waves [25] . This enables the coupling of the acoustic and flame models via acoustic jump conditions across the flame [8,22,25] . When the length of the flame approaches the order of dominant wavelengths, this assumption breaks down and the flame cannot be treated as a single compact perturbation source or sink [28][29][30] . Any prediction of thermoacoustic behaviour will then need to account for the spatial variation in the heat release rate response to acoustic disturbances along the flame, and in the generation of new acoustic waves by the flame.
The spatial distribution of the flame's response introduces the concept of a "local" and continuously varying flame model, characterising the axially varying flame response to normalised velocity perturbations at a fixed reference position. However, current analyses using the continuous concept are mainly based on analytical models with many assumptions [28,29] . For example, distributed FTFs or FDFs with spatially uniform gains and local time delays determined using a "particle flight time" method [31] were incorporated in Helmholtz solvers [32,33] . Practical implementation considers a discrete series of FTFs or FDFs, applied at different downstream locations within the flame [34,35] . When either measuring or simulating local FTFs or FDFs, their accuracy will invoke a trade-off; increasing the number of discrete flame elements captures the non-compactness more accurately, but also deteriorates the signal to noise ratio (SNR) as each flame zone contains less fluctuating heat release rate.
Simulating the response of the flame requires time-intensive, high fidelity numerical simulations which capture the unsteady flow and flame structures of turbulent combustion [25,36] . Recent work has used "low-Mach number" or "incompressible" large eddy simulations to achieve this [17][18][19]37,38] , based on the assumption that the flame is largely unaffected by compressibility effects and responds mainly to hydrodynamic fluctuations [39,40] . This offers the benefit of allowing a larger time step (now restricted by the inverse of the flow speed rather than the speed of sound). These previous simulations assume an acoustically compact flame, whose heat release rate is extracted by spatially integrating the entire reacting zone. The present work also employs "incompressible" 1 LES, but uses this as the basis for extracting local FDFs for a succession of axial flame segments. The obtained local FDFs are substituted into a low-order acoustic network model, where they appear as discrete compact flame elements. It is shown rigorously that compact acoustic jump conditions can be applied across each discrete flame segment. The mean flow variables (such as axial velocity, temperature, speed of sound) are permitted to vary across flame elements. This coupled approach, invoking a segmented flame, local FDFs and a series of compact acoustic jump conditions, is then used to predict combustion instability, including nonlinear features such frequency and amplitude of limit cycle oscillations.
The target case of the present study is a rectangular dump combustor, containing a perfectly premixed flame, named One Rig for Accurate Comparisons with Large-Eddy Simulations (ORACLES) and developed in a European research programme framework [42] . The combustion chamber has a length of 2 m, with the flame extending downstream approximately half the length of the chamber. The combustor is unstable; low frequency limit cycle oscillations at f ∼ 50 Hz occur, with a normalised velocity amplitude of A ∼ 0.29 upstream of the flame. Note that although the flame length L f is approximately half the combustor length, it is the flame length compared to the acoustic wavelength c b, avg / f in the combustion chamber that matters for acoustic compactness, with c b, avg being the average speed of sound in the combustion chamber. At the low frequency of the oscillations, this ratio is f L f / c avg ∼ 1 / 15 ( c b, avg ≈ 750 m s −1 (see Fig. 9 (b))), implying that flame non-compactness, while being higher than typical, is not extremely high. Furthermore, it should be noted that there are not many instability benchmark experiments for combustors with long flames. The present work will investigate the benefits of flame segmentation for limit cycle predictions at this acoustic compactness ratio.
The objectives of the present study are to: (1) adapt a low-order thermoacoustic network model to account for a segmented flame and series of local FDFs; (2) to determine the local FDFs using LES for the first time; (3) to validate the limit cycle characteristics predicted using the coupled approach, by comparing with those determined using a reference technique implementing the 1-D linearised Euler equations (LEEs); (4) to validate the coupled limit cycle predictions by comparing with available experimental results; (5) to investigate the sensitivity of the limit cycle predictions to the (unknown) level of acoustic damping at the combustor boundaries.
The remainder of the paper is organised as follows. The ORA-CLES combustor and the main experimental results are described in Section 2 . The low-order network representation of the combustor geometry, boundaries and segmented flame, along with the reference LEEs method are presented in Section 3 . The numerical details of the LES used to obtain the local FDFs are presented in Section 4 . Section 5 presents the mean flow and thermodynamic properties of the combustor as predicted by the LES, the local FDFs for the flame segments and the continuous FDF for the LEEs method. It then presents predictions of the limit cycle oscillations. The sensitivity of these predictions to the acoustic damping at the combustor boundaries is presented in Section 5 . Conclusions are drawn in the final section.

Experimental test case
The ORACLES combustor considered in the present study is schematically illustrated in Fig. 1 (a), with more information available in [42,43] . The test rig contains three main sections: (1) two straight settling chambers followed by smoothly converging units in which air and propane, supplied at the left end, are fully mixed; (2) two separate narrow rectangular channels of length l 02 = 3 . 270 m, used to obtain two fully developed turbulent flows. These channels are separated from one another by a 10 mm thick plate and both feed (3) a single combustion chamber of length l 3 = 2 m. The chamber consists of a sudden expansion of h = 29 . 9 mm in the positive and negative height directions. The test rig features rectangular cross-sections with width W = 0 . 1505 m.
The combustor is operated at atmospheric pressure with the mean temperature of the fresh mixture being 276 K and its equivalence ratio being φ = 0 . 75 . The flow in the two narrow rectangular channels is the same, with bulk velocity U bulk = 11 m s −1 . Two flames attached to the upstream end of the combustion chamber extend downstream with length ∼ 1 m. (Note that the flame length has not been measured experimentally, only predicted by LES in previous [34] and current work). The combustor exhibits combustion instability leading to limit cycle oscillations with frequency ∼ 50 Hz and normalised axial velocity amplitude ahead of the flame ( x = −5 h ) of ∼ 0.29 [34,42] . Neither the mean temperature in the combustion chamber nor the flame's response to oncoming acoustic excitation have been experimentally measured. In the present work, this information is determined by the LES and presented in Section 4 .

Model of the combustor
The acoustic wavelength at the combustion instability frequency of 50 Hz (the average speed of sound in the combustion chamber approximately equals 750 m s −1 , and the average acoustic wavelength is thus ∼ 15 m in the reacting flow) is much larger than the cross-sectional size of the combustor, allowing the acoustic waves to be treated as 1-D plane waves. The combustor geometry can then be simplified for acoustic treatment, with the upper and lower parts combined into a single stream. The length of the converging unit is much smaller than the dominant acoustic wavelength, allowing simplification to two connected straight ducts with heights H 1 and H 2 , and lengths 0 . 175 m , respectively [44,45] , where H ( x ) denotes the axial height profile of the convergent unit. The effective lengths of the straight upstream settling chamber and subsequent narrow channel thus become l 1 = l 011 + δl 1 = 0 . 34 + 0 . 245 = 0 . 585 m and l 2 = l 02 + δl 2 = 3 . 270 + 0 . 175 = 3 . 445 m , respectively, as sketched in Fig. 1 (b). The final simplified combustor thus only contains three rectangular elements; its dimensions are summarised in Table 1 .
The entrance of the settling chamber is treated acoustically as a closed end, leading to the assumption that the acoustic velocity perturbation is zero ( u 1 = 0 ), and the pressure reflection coefficient is R 1 = 1 . We note that for non-zero Mach numbers, a closed end does not give zero acoustic energy flux [46] , rather zero-mass-flow-rate-perturbation ( ˙ m = 0 ) [47] , leading being defined as the ratio of reflected to incident waves). The settling chamber mean Mach number has a very low value of M 1 = ū 1 / c 1 = 2 . 6 / 333 ≈ 0 . 008 , making the pressure reflection coefficient for zero-mass-flow-rate-perturbation R 1 = 0 . 984 . We note that although a closed upstream end is a good approximation to the constricted fuel-air entrance, the real boundary condition will be more complicated, and the porous plugs and grids shown in Fig. 1 (a), will lead to further acoustic losses [48][49][50] . Unfortunately, R 1 has not been experimentally characterised. In the first set of calculations, we thus use R 1 = 1 ; and in Section 5.5 , we investigate the sensitivity of the thermoacoustic predictions to R 1 . The downstream end of the combustion chamber is taken as the end-point of the acoustic network, with an "open" pressure reflection coefficient of R 2 = −1 . Again, for non-zero Mach numbers and isentropic flows, zero acoustic flux follows from the zero-enthalpy-perturbation boundary condition ( p / ρ + ū u = 0 ) [47] . However, as this pressure reflection coefficient R 2 is known to drop in absolute value with frequency [15,46,51] , and any entropy perturbations present will render the zero-enthalpy-perturbation condition invalid, the actual pressure reflection coefficient will be more complicated. In the preliminary calculation, we thus use R 2 = −1 and then go on to study the sensitivity to R 2 in Section 5.5 .

Modifying the low-order acoustic network model for axially segmented flames
The low-order acoustic network modelling tool OSCILOS (the open source combustion instability low-order simulator), is used. This is being developed in the authors' group [15,24] and has been validated against experiments [19] . The combustion system is represented as a network of simple connected acoustic modules, each corresponding to a certain component of the system [21,52] . Within each module, two planar acoustic waves propagate, one in either direction, at the speed of sound with respect to the mean flow. Wave strengths between modules are tracked using the linearised flow conservation equations, and the acoustic boundary conditions are applied at either end. For the current system, the settling chamber, the narrow channel and the downstream half of the combustion chamber ( x ∈ [1,2]

Wave propagation in the j th module
According to linear acoustic theory, all flow and thermodynamic variables (e.g. p ) can be represented as the sum of a mean value (e.g. p ) and a perturbation (e.g. p ). The latter is assumed small compared to the relevant mean value. For modules in which the mean thermodynamic properties and axial flow velocity can be considered uniform, the acoustic field is then the superposition of downstream and upstream propagating plane waves, which can be characterised by considering both the pressure fluctuation, p , and velocity fluctuation, u , at a given location. Assuming that all fluctuating quantities have a time dependence in the form a = ˆ a e i ωt , ω being complex angular frequency, the pressure ˆ p j , velocity ˆ u j , and entropy ˆ s j amplitudes in element j (with inlet at x + j−1 and outlet at x − j as sketched in Fig. 2 ) can be expressed as: where are three time delays. ρ is the density and c the speed of sound.

Wave strength jump conditions between elements with no flame
For an abrupt area increase, such as that which occurs between the narrow rectangular channel and the combustion chamber, the mass and energy fluxes are conserved across the interface. (Flames are considered downstream of the abrupt increase). The momentum flux increases by the axial force on the walls [22] , to account for the flow separation and hence acoustic energy losses that occur at such area increases. This leads to three jump conditions linking waves either sides of the interface x j : where j = A j+1 / A j . A j is the cross-sectional surface area of the module j . M = ū / c is the Mach number and γ the ratio of specific heats.
For an abrupt area decrease, such as that which occurs between the settling chamber and the narrow channel, the mass and energy flux are conserved across the interface, and the flow can also be assumed isentropic [22] . The second equation in the matrix equation (2) is then replaced by: Readers should refer to [24] for more detailed derivation of these equations and the relations of the mean properties.

Low-order network representation of a segmented flame
Flames which are sufficiently long cannot be treated as a single compact flame sheet like those in [19,27] . Any prediction of thermoacoustic behaviour needs to account for the variation in the heat release rate response to acoustic perturbations along the flame, and subsequently for the axial variation in generated flow perturbations. In this work, we segment the long flame axially into a discrete series of uniform-width flame segments, each with axial width short enough to be represented as a compact flame sheet. This allows us to implement a series of discrete compact flame models for the response of the unsteady heat release rate to incoming excitation, and to apply a series of compact acoustic jump conditions for the subsequent effect on the flow perturbations. Between neighbouring compact flame sheets, the flow is considered uniform and without heat release rate.
In order to rigorously derive the appropriate form of the acoustic jump relations, we start with the 1-D linearised mass, momentum and energy conservation equations for a duct of uniform cross-section, in which viscous terms and heat diffusion are neglected [53] : Here ˆ ˙ q is the cross-sectionally integrated heat release rate perturbation, U = ρc ˆ u and S = ρc 2 ˆ s /c p . It should be noted that the mean thermodynamic and flow properties (e.g. p , ū and T ) used above are cross-sectionally averaged values and depend only on axial location, x . (In this work they are extracted from LES results.) γ is considered only a function of mean temperature, and is calculated within OSCILOS [54] , assuming values for air for simplicity.
We now would like to show that the continuous long flame in Fig. 3 (a) can be uniformly axially segmented to a series of discrete compact flames (schematically shown in Fig. 3 (b)) for thermoacoustic analysis under certain conditions. For ducts with uniform or slightly changing mean properties, pressure and velocity perturbations can be considered to be the superposition of a downstream and upstream waves [55,56] , e.g., as described in Eq.
(1) . These two waves both propagate at around the speed of sound. The entropy waves propagate at the bulk mean flow velocity [52] . We now consider the zone [ x j,c , x j+1 ,c ] of width l in Fig. 3 (a). When l c / f and l ū / f, axial changes in acoustic waves and entropy waves remain small [52] . The axial integral of the linearised energy equation ( Eq. (6) ) along it can be simplified to: where ( a ) j+1 ,c j,c represents the axial averaged value of a within the We now consider the simplified configuration sketched in Fig. 3 (b), with a compact flame, denoted by j , located at the centre of the zone [ x j,c , x j+1 ,c ] . The distance between neighbouring compact flames is l , and the mean thermodynamic and flow properties are considered uniform between neighbouring flames, equalling axially averaged values within this zone: The mean pressure, velocity and temperature upstream of the flame j are denoted p j , ū j and T j , respectively, with values downstream denoted p j+1 , ū j+1 and T j+1 . There are thus three regions within [ x j,c , x j+1 ,c ] : a compact flame sheet separating two uniform flow regions, with jumps of mean thermodynamic properties across the compact flame. The axial integral of the linearised energy equation ( Eq. (6) ) along [ x j,c , x j+1 ,c ] can also be simplified to: When changes in mean thermodynamic and flow properties are small within [ x j,c , x j+1 ,c ] , the energy equation (9) can be considered consistent with Eq. (7) , with the same true for the mass ( Eq. (4) ) and momentum ( Eq. (5) ) equations. It should be noted that the condition l c / f and l ū / f can be used as the criterion of the maximum flame segment width or the minimum flame segment number, which is further discussed in Section 5 . Thus a discrete sequence of compact flames can be used to approximate the continuous long flame in the thermoacoustic analysis. The spatial intervals [ x + j−1 , x − j ] can be considered flame-free modules in the thermoacoustic analysis, within which the wave propagation transfer matrix ( Eq. (1) ) applies. As the mean temperature varies smoothly with axial location within these modules, the three time delays are better calculated using to jump conditions ( Eq. (2) ) to account for the heat release rate perturbation from the compact flame j , leading to: It should be noted that the flame segmentation method was proposed by the authors in [34,35] for thermoacoustic prediction of the ORACLES combustor using the low-order code LOTAN [57] , although full details on how the series of sub-flames was treated were not given. Low-order network modelling assuming zero-mean-flow was also used by Leandro and Polifke [58] , which dealt with the distributed heat release by segmenting the flame to a sequence of linear n − τ flame models; the criterion for flame segmentation and the treatment of the mean thermodynamic properties were not detailed. For low Mach number situations, the jump conditions ( Eq. (10) ) can be simplified to the version retaining Mach number terms up to the first order as shown in the Appendix A or those in [59,60] , or the zero-mean-flow jump conditions [25] ; a recent analytical investigation [61] proposed that the generation of entropy waves by a compact perfectly premixed flame can be neglected for vanishing flow Mach numbers. However, the current work proposes a general model suitable for non-compact flames or heaters, applicable not only to low-Mach number flow situations. The entropy waves can be important in certain configurations and complete jump conditions will then be necessary.

Compact local flame describing functions
To close the thermoacoustic system, it is necessary to link the heat release rate perturbation, ˙ Q j , of each compact flame segment to the acoustic velocity perturbations at a fixed reference position. The local flame describing function for flame j can be expressed as: where ˆ u ref /U bulk is the normalised velocity amplitude at the reference position, ˙ Q the mean total (not local) heat release rate, and Fig. 1 ) is taken as the reference position, due to the availability of experimental data here. It can be calculated using: It should be noted that when a fully compressible solver is used, the flame responds both to hydrodynamic and acoustic disturbances [4] . Acoustic waves (or velocity perturbations) propagate in the CFD domain, with the propagation speed (or the time delay) relying on the local bulk flow velocity and speed of sound. When the distance between the reference point and the sub-flames is not small compared to the wavelength, the wave amplitudes alter along the combustor in a way which depends on the boundary conditions [62] . The distance between the reference point and flame thus should be sufficiently small when using fully compressible solvers [63] . However, when an incompressible solver is used, acoustic terms are removed from the conservation equations. Density changes are uniquely associated with temperature changes, but not acoustic perturbations [4,41] . The acoustic wave feedback into the CFD domain from the boundaries is thus suppressed [38] , and the effect of boundary conditions on the flame response is negligible. The change in amplitude of the fresh mixture velocity perturbation is small during convection from the reference point to the sub-flames. An incompressible solver is used in the present study, and the flame model is thus believed to be robust to the boundary conditions.

Determining the thermoacoustic modes
The pressure and velocity perturbations at the upstream end of the thermoacoustic network shown in Fig. 1 (b) can be related via the pressure reflection coefficient, Substituting the local discrete flame describing functions ( Eq. (11) ) into the modified jump conditions ( Eq. (10) ), and combining the transfer matrices for wave propagation ( Eq. (1) ) and boundary conditions, the governing matrix equation for the system becomes: where α j = F j ˙ Q / (U bulk A j ) and α j = 0 when the interface j is not a flame. D is a square matrix of size (3 N + 1) × (3 N + 1) . The thermoacoustic modal frequencies of the system, ω = ω r − i σ, occur when det(D ) = 0 . These frequencies are determined numerically, yielding an angular frequency ω r = 2 π f and a growth rate σ . The stability of the system is defined by the sign of the growth rate, with a positive value corresponding to an unstable thermoacoustic mode.

Reference predictions using the linearised Euler equations (LEEs)
A first validation of the low-order coupled approach, in which segmented discrete flame models and acoustic jump conditions are embedded in a low-order thermoacoustic network model, is performed using the one-dimensional linearised Euler equations (LEEs, Eqs. (4) - (6) ). The settling chamber and narrow channel upstream of the flame are treated using wave-based models, with the LEEs used only in the combustor. The two computational domains are connected by matching the impedance, Z 3 = ˆ , and the entropy fluctuation, S 3 (x + 2 ) , immediately downstream of the abrupt area increase interface. The LEEs implement the unsteady heat release rate response of the long flame using a continuous local flame describing function. This describes the local heat release rate response to oncoming acoustic disturbances: (14) In practice, the heat release rate perturbation, ˆ ˙ q, is difficult to directly extract from numerical or experimental data, as it relates to such a small flame area that the signal to noise ratio becomes too small. In this work, it is reconstructed based on the segmented flame method; a full description of this is contained in Section 5.3 . In order to implement the LEEs, spatial discretisation of the three equations ( Eqs. (4) -(6) ) along the 2-meter-long combustion chamber is via a second order finite difference scheme on a uniform grid containing 20 0 0 points. The main thermoacoustic mode predicted using this method is used as the reference against which predictions using the coupled low-order segmented flame method are compared.

Numerical method for flame simulation
Previous LES [64,65] generally extended downstream over a short zone covering only part of the long flame. This enabled highly resolved prediction of the flow field in the region where experimental data are available. The flame in the current case has a length of around 1 m [34] , and simulating such a short computational zone will not capture the response of the whole heat release rate. In this work, LES are performed over a longer region starting at x = −0 . 15 m and ending at x = 1 m, as sketched in Fig. 4 (a). The computational domain thus extends upstream into the two rectangular channels, and downstream a meter into the combustion chamber. It includes the reference position for velocity fluctuation measurements, x ref = −5 h = −0 . 1495 m, and nearly the whole heat release rate zone (at some moments, a very small part of the flame may leave the computational domain. However, the time-averaged heat release rate is found to be well within the computational domain). Based on previous research [64][65][66][67] , the computational domain is then meshed using structured hexahedral cells, as shown in Fig. 4 (b) and (c). The total number of mesh cells is 5.46 million. In the inlet region, the grid points are 180 × 100 × 70 (streamwise x , vertical y and lateral z , respectively), while in the combustor region, the grid points are 300 × 200 × 70. The first node of the mesh is always located 0.1 mm away from the solid walls. The maximum cell skewness is 0.157, sufficient to achieve a good mesh  quality. The histogram of the closest wall-normal distance of the nodes normalised by the viscous lengthscale, y + , is shown in Fig. 5 (a). As an artificial unsteady velocity is imposed at the inlet, the velocity gradient at wall is high there, leading to large value of y + . The majority (87%) of the near-wall cells are located within the viscous sublayer ( y + ≤ 5 ), and the averaged value of y + equals 2.6, again indicating that the mesh is sufficient fine for the present LES.
As the numerical methods used in the present LES are nearly the same as those in our previous paper [19] , we only briefly recap the numerical setup in the present work. The response of the turbulent flame to upstream forcing of the velocity perturbations is simulated using incompressible reacting flow LES, based on the open-source CFD toolbox, OpenFOAM [68] . The equation of state used in the present LES is expressed as ρ = p 0 / ( R g T ) , where the pressure p 0 remains constant based on the zero-Mach number assumption implemented in the asymptotic expansion of the Navier-Stokes equations [41] , and R g denotes the gas constant. The changes in density ρ are thus uniquely associated with changes in temperature T , and acoustic waves are totally removed from the flow simulation; this allows a larger time step not restricted to the inverse of the speed of sound and accelerates the computing speed [4] . Specifically, a modified version of the reacting-FOAM solver is applied, which has been used successfully in previous studies [19,[69][70][71] . The governing equations are the Favrefiltered Navier-Stokes equations of mass and momentum, plus the species mass fraction and energy conservation equations [4] . The code employs a collocated Finite Volume Method (FVM) in which the discretisation is based on Gauss' theorem together with a semiimplicit time-integration scheme. The Preconditioned-Conjugate-Gradient (PCG) algorithm is used for solving the pressure, and the Preconditioned-BiConjugate-Gradient (PBiCG) algorithm is applied for velocity, enthalpy and species. The algorithm for solving mass and momentum conservation equations is based on the PIM-PLE method (merged PISO-SIMPLE algorithms) suitable for transient simulations [68] . For temporal advancement, a second-order Crank-Nicolson scheme is applied for time derivatives, and convection terms are discretised using a second order central difference scheme. In order to close the conservation equations, Sub-Grid Scale (SGS) turbulence modelling is required, in this case using the one equation SGS kinetic energy transport model [72] . The Partially-Stirred Reactor (PaSR) combustion model [19,69,71,73] is used to model turbulence-combustion interactions; this has been validated for premixed and partially-premixed flames in previous investigations [19,69,71] . For the mixing time scale τ m in the PaSR  model, a newly developed model [74,75] is applied. Based on the time scale of the SGS model and the Kolmogorov time scale, the mixing constant is set to 3.0 to describe the correct features of the present long flame. A global one-step propane/air mechanism involving 7 intermediate species [76] is used. This was designed to model lean premixed propane/air flames and has been successfully applied for previous LES studies of combustion instability [76] . A relatively large computational time step is allowed, which is fixed as 0.01 ms in the present LES, resulting in the maximum Courant number being around 0.5.
Zero-gradient conditions are applied to the outflow velocity ( du / dx = 0 ) and scalars at combustion chamber exit (see Fig. 4 (a)). Except for the two inlets, the outlet, and the two lateral planes ( z direction), all the other boundaries are treated as solid adiabatic non-slip walls [64][65][66][67] , as it was stated in [42,43] that combustor walls are thermally insulated with a refractory material, although it was shown in [77,78] that heat losses (or changes in the wall temperature) have a strong impact on the flame response. The boundaries in the lateral ( z ) direction (i.e. two lateral planes) are set to periodic [64,66] .
The unsteady inflow velocities are prescribed at the two channel inlets. For the base flow without acoustic forcing, the timeaveraged bulk velocity at the reference location (or at the inlet which is very close to the reference location) is U bulk = 11 m s −1 , giving a Reynolds number of R e ≈ 25, 0 0 0 [42] . In order to determine the flame describing function, a harmonic forcing is superimposed on the time-averaged bulk flow and artificial turbulence fluctuations [79] at the inlet of the computational zone. This leads to a stream-wise velocity profile which can be expressed as: u (y, t ) = U (y ) ( 1 + A sin (2 π f t ) + w t (t) ) (15) where U ( y ) is a mean vertical velocity profile (fitted using a seventh order polynomial to experimental data at the reference location x = −5 h of the mid y − z plane), f is the forcing frequency and A the normalised forcing amplitude. A and f are varied independently in order to obtain the flame describing functions. w t ( t ) corresponds to the artificial turbulence fluctuations with a level of 5% based on the experimental result [42] .
The subgrid activity parameter, s , is also calculated, defined as [80] : where ε t and ε μ represent the turbulence and molecular dissipation with the Reynolds number respectively. s satisfies 0 ≤s < 1, with s = 0 corresponding to DNS and s = 1 to LES at an infinite Reynolds number. The value of s thus quantifies the importance of the SGS model. In the present LES, the distribution of s for a reactive operating condition ( A = 0 . 29 and f = 50 Hz) is shown in Fig. 5 (b). The value of s remains small within the upstream channels but becomes larger in the downstream regions, indicating that the mesh cells are well-refined near the inlet but become relatively coarse downstream, where the SGS model becomes more important in this highly turbulent region.
The predicted vertical profiles of the mean and fluctuating velocities are compared to experimental data [42] at multiple axial locations (see Figs. 6 and 7 ). The simulation predictions match the experimental measurements quite well, although slight discrepancies exist at downstream locations. The LES codes therefore capture the reacting flowfield sufficiently accurately to be used as a basis for calculating the flame describing functions, which are presented in the next section.

Mean thermodynamic properties in the combustion chamber
In the thermoacoustic analysis, the mean thermodynamic and flow properties in the combustion chamber are extracted from LES results. Figure 8 shows the spatial distribution of the time-averaged temperature and axial velocity in the upper half of the plane z = 0 for upstream velocity forcing at f = 50 Hz and A = 0 . 29 .
The flame forms a triangle shape, consistent with experimental results [42] . Its effective length is around 1 m. The spatial information is then averaged in y to yield axial distributions in mean temperature, T , and mean axial velocity, ū , shown in Fig. 9 (a). In the approximate flame region, x ∈ [0, 1] m, the mean temperature increases from 1011 K to 1854 K, while the mean axial flow velocity increases from 6.0 m s −1 to 39.4 m s −1 . In the first set of thermoacoustic predictions, it is assumed that the combustion chamber does not exhibit heat loss downstream of the flame, such that T = 1854 K and ū = 39 . 4 m s −1 in the region x ∈ [1, 2] m.
The mean axial temperature, T , is then substituted into OSCI-LOS to calculate the axial variation of the ratio of specific heats, γ , and speed of sound, c , in the flame region, shown in Fig. 9 (b).
The thermodynamic properties of air are used to represent those of the fresh mixture and burned gases. It should be noted that γ changes with mean temperature, decreasing along the flame zone from 1.335 to 1.301, while the speed of sound increases from 622.5 m s −1 to 832.1 m s −1 . Note that for the sequence of compact flame elements in the thermoacoustic network model, the mean flow and thermodynamic properties between flame elements are calculated by axial averaging using Eq. (8) .

Local flame describing functions of flame segments
In order to extract the local flame models for the axial series of flame segments, the computational domain in the combustion chamber ( x ∈ [0, 1] m) is equally segmented. As outlined in Section 3. The heat release rate, ˙ Q j , within each segment, j , is spatially integrated and the time evolution is recorded. The first flame segment is numbered 3, the numbers 1 and 2 representing the first two interfaces separating the settling chamber, the narrow channel and the combustion chamber (see Fig. 1 (b)). Because self-excited oscillations are known to occur at ∼ 50 Hz, forced simulations are performed at three frequencies -40 Hz, 50 Hz and 75 Hz. For each forcing frequency, five forcing amplitudes are performed:   Fig. 11 confirms that the main response frequency matches the forcing frequency -this was true for all cases investigated and validates the weakly nonlinear assumption of the flame describing function [10,12] .
In order to determine the flame describing functions (FDFs) of the series of discrete compact flames, it is necessary to extract the amplitude and phase lag of the heat release rate response to the pure sine velocity forcing signal, u ref /U bulk , for different velocity forcing amplitudes. In this work, the DFT at the forcing frequency,   e.g. the main peak value shown in Fig. 11 , is used as the amplitude of the heat release rate perturbation at the forcing frequency.
However, it is difficult to obtain an accurate phase lag directly from the cross-correlation between ˙ Q j and u ref /U bulk , especially in segments with lower heat release rate where the SNR is very small. A post-processing method, used previously in ultrasonic signal processing [81,82] , was applied: (1) the raw signal ˙ Q j is first filtered using a zero-phase-shift band-pass filter to eliminate the low frequency noise due to the bulk flow motion and the high frequency noise due to flame nonlinearity and small-scale turbulence. For instance, a digital Butterworth filter with a lower cut-off frequency of 30 Hz and an upper cut-off frequency of 70 Hz, featuring an attenuation of less than 1 dB in the pass-band and at least 10 dB in the stop-band is used for the 50 Hz forcing cases. The resulting signal, ˙ Q j, flt (represented by the thick line in Fig. 10 (a)) contains only the forcing frequency. (2) A Hilbert transform is then used to obtain the envelope of the signal ˙ Q j, evlp (represented by the dashed line in Fig. 10 (a)). (3) This envelope enables renormalisation of the signal, ˙ Q j, flt , yielding a constant amplitude signal, ˙ Q j, nml = ˙ Q j, flt / ˙ Q j, evlp , as shown in Fig. 10 (b). (4) A cross-correlation is then conducted between the normalised signal, ˙ Q j, nml , and the pure sine forcing signal, u ref /U bulk , to obtain an accurate phase lag at the forcing frequency.
The amplitude and phase lag of the heat release rate perturbations are then substituted into Eq. (11) to construct the flame describing function, which is a complex function with a gain and a phase: This process is repeated across flame segments, forcing frequencies and forcing levels to yield the overall FDF shown in Fig. 12 . The gain generally decreases with increasing forcing level, A , due to nonlinear saturation. It decreases rapidly between A = 0.1 and 0.2, but less rapidly as A is further increased. This is consistent with other numerical predictions [64] . The phase lag decreases with forcing level, which is consistent with the response of a triangle laminar flame [13] . The phase lag of the flame response is generally proportional to the phase, ϕ c = −2 π f x/U bulk , of the bulk flow convection [14] . This is also confirmed by the plots shown in Fig. 12 . Because combustion instability occurs at around 50 Hz, the thermoacoustic analysis is restricted to the frequency range [25,100] Hz. The gain and phase of the FDF across this frequency range is obtained by linear interpolation using:

Reconstruction of 1-D axially continuous flame describing function for the long flame
The reference LEEs method requires a 1-D axially continuous form of the FDF. This section presents a method to reconstruct this from the discrete FDFs already obtained. When the width of each flame segment is sufficiently small, the change in flame response over one segment will also be small; the Taylor expansions for the gain and phase of the continuous FDF will be dominated by terms of zeroth and first order. The gain and phase can thus be approximated as linear functions of x , allowing the continuous flame describing function, F c , to be obtained by linear interpolation of the discrete FDFs, F j , for the series of compact flames. We assume that F c in the zone [ x j − l 2 , x j + l 2 ] takes the form: where β 1 and β 2 are two linear coefficients for the gain and phase, respectively. Integrating F c along the zone [ x j − l 2 , x j + l 2 ] leads to: When β 2 l π , the above equation can be simplified to: As mentioned in the previous section, the phase lag of the discrete flame describing functions approximately corresponds to the bulk flow convection, which has wavelength ∼0.22 m at 50 Hz. The phase lag when N s = 40 is then β 2 l ∼ 2 π × 0 . 025 / 0 . 22 = 0 . 23 π π . The condition for Eq. (22) is satisfied and N 0 and ϕ 0 can be calculated using: In this work, the gain and phase of F c are interpolated based on N j /l and ϕ j respectively. Figure 13 shows the comparison between the integrated axially continuous FDF, F * j , and the discrete series of FDFs for the series of flame segments, F j , for f p = 50 Hz, A = 0 . 1 and different N s . The axial variations are nearly identical when N s = 40 , again indicating that this flame segment number is an appropriate choice. When coarser flame segments are used, e.g., N s = 10 , differences between F * j and F j cannot be neglected; this validating this method. The continuous flame describing function at other frequencies in the range [25,100] Hz are then linearly interpolated using Eqs. (18) and (19) .

Prediction of thermoacoustic modes and limit cycle characteristics
The discrete flame describing functions, F j , for the series of compact flame segments along with the mean properties determined in Section 5.1 are then substituted into the governing equation ( Eq. (13) ) to determine the thermoacoustic modes of the system. Because the pressure reflection coefficients at the inlet and outlet are not experimentally measured, it is assumed in the preliminary predictions that the pressure reflection coefficients at the inlet and outlet are R 1 = 1 and R 2 = −1 , respectively; acoustic damping at the boundaries can be imposed by reducing the absolute values of the two coefficients and sensitivities of the thermoacoustic predictions to the values of R 1 and R 2 are further conducted in Section 5.5 .
Predictions using the reference 1-D LEEs method implement the same boundary conditions, use the continuous mean thermodynamic properties determined in Section 5.1 and the reconstructed continuous flame describing function F c determined in Section 5.3 . These are all substituted into the LEEs network modelling (presented in Section 3.3 ). Figure 14 shows the predicted modal frequency and corresponding growth rate for different normalised velocity amplitudes, A . The modal trajectories with increasing A for both the low-order coupled methods and the reference LEEs method are almost the same. For the lowest amplitude of A = 0 . 1 , the growth rate of the main mode is clearly positive, indicating that the system is unstable at the corresponding frequency of f = 51 . 8 Hz. As A increases, the growth rate quickly decreases, with a limit cycle is established when the growth rate is zero. This occurs between amplitudes of 0.3 and 0.4, although it is noted that the growth rate is flat and close to zero between amplitudes of 0.3 and 0.4. Linear interpolation can then be used to predict the amplitude  at which zero growth rate occurs, this being the predicted limit cycle amplitude. For the present results, A LC = 0.38 with the corresponding frequency being f LC = 51.6 Hz. These are very close to those predicted using the reference LEEs method, validating the flame segmentation approach of the low-order coupled method. The predicted limit cycle frequency f LC matches the experimental result ∼ 50 Hz; however, the predicted limit cycle amplitude A LC is overestimated compared to the experimental value ∼ 0.29, most likely because nearly no acoustic damping is accounted for using the pressure reflection coefficients R 1 = 1 and R 2 = −1 in the preliminary calculation.
The pressure and velocity modeshapes predicted by the two methods are compared in Fig. 15 . Predictions using the loworder coupled method match those predicted using the reference method. There are two large velocity perturbation jumps at the two interfaces (  It is interesting to calculate the axial distribution of the Rayleigh source term, R c , for the series of compact flames. This is defined by [83][84][85] : where ϕ p and ϕ q are the phases of the pressure perturbation, ˆ p j , and the heat release rate perturbation, ˙ Q j , at x − j respectively. This criterion provides a measure of the correlation between pressure and heat release rate perturbations, is seen to be the driver of combustion instabilities when mean flow is not accounted for [85] . Figure 16 shows the axial distribution of Rayleigh criterion R c in the flame region when A = 0 . 1 . Unlike the uniform behaviour of short flames with length much smaller than the wavelength of dominant acoustic waves, the long flame in the present study features Rayleigh sources in some locations (e.g., x ∈ [0.43, 0.67] m) and sinks in others (e.g., x ∈ [0.31, 0.43] m). It is potentially because of these opposing effects that the overall growth rates may tend to be smaller for long flames than for short flames [19,86] .
As outlined in Section 3.2.3 , in order for the approach of the low-order coupled method to be justified, the width of the flame segment must be much smaller than the flow perturbation wavelengths. However, it is interesting to examine the performance of this method as the width of the flame segments increases. This can be achieved by combining neighbouring flame segments e.g. by pairing flame segments, the segment number N s reduces from 40 to 20 and the flame segment width increases from 0.025 m to 0.05 m. Figure 13 shows the discrete flame describing function F j for N s = 40 , 20 and 10, f = 50 Hz and A = 0 . 1 . The gains of F j are divided by their corresponding flame segment width l (or multiplied by the corresponding N s ) for comparisons. As l is increased, the SNR of the heat release rate is increased, enabling better capture of the FDF gains and phases at the forcing frequency. However, the gains are lower than those for N s = 40 (also for the rest of four forcing levels), because the heat release rates in the neighbouring segments are not necessarily in phase.
The predicted modal frequency and corresponding growth rate across five forcing levels and five flame segment numbers ( N s = 40 , 20, 10, 5 and 1) are collected in Table 2 . As the number of flame segments decreases, the frequency nearly does not change, the modal growth rate decreases across all forcing levels, and the predicted limit cycle amplitude, A LC , decreases. When N s ≤ 5, the axial variation of the flame describing functions is totally different Table 2 Predicted modal frequency and corresponding growth rate for different forcing levels and limit cycle frequency f LC and normalised velocity perturbation amplitude A LC using different methods. R 1 = 1 and R 2 = −1 . Constant temperature in the zone [1,2] m. indicates predictions by assuming that entropy perturbations generated from one flame segment are fully attenuated when they reach the following flame segment. from that when N s = 40 and the combustor is in fact predicted to be linearly stable rather than unstable. This confirms the need to adopt a segmented flame approach for coupled low-order prediction of this combustor. Although the limit cycle amplitude A LC predicted by using a smaller N s (e.g., N s = 20 or 10) is closer to the experimental result A LC ∼ 0.29, it should be noted that all these predictions assume that R 1 = 1 and R 2 = −1 , and so this does not mean that the flame segmentation approach using a smaller N s provides a better prediction of the limit cycle. Note also that the growth rate is flat and close to zero between the amplitudes of 0.3 and 0.4 (e.g., that shown in Fig 14 ). When more realistic acoustic losses are accounted for at the two boundaries, the predicted limit cycle amplitudes using the coupled method with N s = 40 and the reference LEEs better match the experimental result. It should be noted that limit cycle predictions (frequency ∼ 50 Hz and amplitude 0.37) were also conducted in [34] . However, LES were performed only at one forcing frequency (50 Hz) and the 1.1-m-long combustion chamber was equally segmented to 14 slices with no criterion for the flame segmentation. The discrete FDF was constructed by assuming that the gain and phase in each segment remain constant for all frequencies.
Predictions of thermoacoustic modes and limit cycle characteristics are also conducted by assuming that entropy perturbations generated from each flame segment are fully attenuated potentially due to the dissipation and shear dispersion mechanisms [87,88] when they reach the following segment in the coupled low-order segmented flame method. Small differences are found between predictions considering and neglecting oncoming entropy perturbations for each flame segment (see Table 2 ), and can be neglected because the flow Mach number is small in the present study. The effect of entropy generations by the current flame on the thermoacoustic predictions is thus small.

Sensitivity analysis of the effect of reflection coefficients and mean temperature distribution
In the previous section, the thermoacoustic predictions assume that there are nearly no acoustic losses at the two boundaries ( R 1 = 1 and R 2 = −1 ). However, as shown in Fig. 1 (a), porous plugs and grids are installed at the beginning of the settling chamber to organise the flow, and are likely to lead to acoustic losses. An analysis of the sensitivity of the predicted limit cycle characteristics to losses at the boundaries is performed by successively varying the pressure reflection coefficient magnitudes (while keeping their phase fixed). As R 1 is reduced from 1.05 to 0.8 for R 2 fixed at −1 (note that the absolute value of reflection coefficient may exceed unity for a finite Mach number flow [47,89] ), the growth rates across all five forcing amplitudes uniformly decreases by 8 rad s −1 . The effect on the predicted limit cycle amplitude is quite dramatic, as shown in Fig. 17 , A LC decreasing from 0.42 to 0.28. This is associated with the flat growth rate, close to zero, between the forcing amplitudes 0.3 and 0.4; it should be remembered that this is visible in Fig. 14 . The change in modal frequency is very small. A similar analysis for the outlet boundary condition, reducing | R 2 | from 1.05 to 0.8, is also performed. The growth rates across the five forcing levels again decrease, with the predicted limit cycle amplitude decreasing from 0.39 to 0.29.
Finally, the preliminary calculation also assumed that the mean temperature in the downstream half of the combustion chamber was constant. However, heat losses occur in practice and the mean temperature is likely to drop with axial distance beyond the flame zone [90][91][92] . The effect of a linearly reducing axial mean temperature profile in downstream half of the combustor, x ∈ [1, 2] m, is therefore now modelled, in order to evaluate its effect on the predicted thermoacoustics. This can be expressed as: where T * = 1854 K is the mean temperature at x = 1 m and (T TD − 1) T * is the gradient of the linear distribution. When T TD = 1 , the mean temperature does not change; when T TD = 0 . 5 , the mean temperature changes from T * at x = 1 m to T * / 2 at x = 2 m. A sensitivity analysis of the predicted limit cycle to the coefficient T TD is performed, in which the effect of the downstream temperature distribution region, x ∈ [1, 2] m, is calculated using the 1-D LEEs method, while the upstream part of the combustor is calculated using the low-order coupled method. Fig. 17 , shows that as T TD decreases, the growth rates across all five forcing levels increases and an increased limit cycle amplitude is predicted e.g. A LC increases from 0.38 when T TD = 1 to 0.41 when T TD = 0 . 5 . Energy transfer occurs during communication between the acoustic waves and the temperature distributed zone, with the mean temperature distribution acting as an acoustic energy source; this is consistent with findings in [53,56] . To summarise, the predicted limit cycle amplitude decreases when acoustic losses at the end boundaries are accounted for, and increases as the mean temperature fall-off with axial distance beyond the flame is accounted for. These factors strongly affect the predicted limit cycle amplitude, motivating better experimental characterisation of them for combustion rigs in which thermoacoustic oscillations are relevant.

Conclusions
Combustion instabilities are problematic for modern industrial gas turbines and aero-engines. The ability to accurately predict the limit cycle behaviour resulting from instabilities at the early design stage is highly desirable, but this constitutes a significant challenge due to complex interactions between turbulent flames, hydrodynamic flow features, acoustic disturbances and loss mechanisms, all in complex geometries. A promising approach for thermoacoustic prediction is to couple low-order acoustic network models with flame models derived from high fidelity simulations. However, such coupled low-order methods have traditionally treated the entire flame as a compact sheet, an approach which is not valid for combustors with sufficiently long flames. This work has presented a coupled low-order method for the thermoacoustic analysis of combustors containing long flames. The method axially segments the flame, and treats each segment as a compact flame sheet. This method has been shown to be rigorously justified, provided the flame segment is sufficiently narrow. This enables the coupling of the acoustic field with the long flame via a series of local flame describing functions (FDF), one for each flame segment, and a series of acoustic jump conditions across each flame segment. The mean thermodynamic and flow properties between the neighbouring compact flame sheets can be considered uniform, enabling an analytical description of wave propagation between compact flame elements. The thermoacoustic modal frequencies and growth rates can then be calculated using a dispersion relation for the entire network; for nonlinear flame elements, this will depend nonlinearly on the amplitude of velocity perturbations at a reference location ahead of the flame.
The target case considered was a rectangular dump combustor, containing a perfectly premixed flame -the well-known ORACLES combustor. The combustion chamber has a length of 2 m, with the flame extending down nearly half of this length. The combustor is unstable at the investigated operating condition.
The nonlinear response of the heat release rate to velocity excitation was characterised by an axial sequence of flame describing functions, these being obtained from large eddy simulations (LES) of the flame region. The LES was performed using the CFD toolbox OpenFOAM, employing an incompressible flowfield assumption for computational efficiency, and modelling the turbulent combustion using the Partial Stirred Reactor (PaSR) model with a global one-step reaction mechanism. The computational zone in the combustion chamber was axially segmented into 40 equal width segments, the heat release rate in each being spatially integrated. Forced simulations enabled the sequence of FDFs to be extracted, as a function of frequency, forcing amplitude and segment number.
These, along with the axial variation of mean flow and thermodynamic properties, were coupled within a low-order acoustic network model. The frequency and amplitude of the limit cycle oscillations resulting from instability were then predicted. They firstly agreed very well with those from a reference method based on the linearised Euler equations, validating the flame segmentation implementation. They also compared well to the experimentally measured limit cycle characteristics. The predicted limit cycle frequency and amplitude were f LC = 51.6 Hz and A LC = 0.38, compared to experimental measurements of f LC ∼ 50 Hz and A LC ∼ 0.29.
The levels of boundary acoustic damping in the experiments were not known, and the above predictions almost neglected them. A sensitivity analysis was therefore performed which showed that the predicted limit cycle amplitude rapidly decreased when boundary acoustic damping was accounted for. The predicted limit cycle frequency was almost unchanged. It was also shown that mean temperature fall-off downstream of the flame zone within the combustor led to an increase in the predicted limit cycle amplitude, indicating that this led to energy transfer to the acoustics.
Overall, this work has presented a justified, viable and efficient method for predicting the thermoacoustic behaviour of combustors with long flames. It has also presented findings which motivate accurate determination, by experimental characterisation, improved simulations or modelling, of boundary acoustic damping and combustor temperatures for combustors in which thermoacoustic behaviour is likely to be important.
where β j = c j+1 / c j and = ( γ j+1 -1)/( γ j -1). The first two simplified jump conditions agree with those in [60] when the ratio of heat capacity is assumed constant. For the abrupt area decrease where no heat addition presents, the mean temperature change ratio, T j+1 / T j − 1 , has the order of magnitude of M 2 j and is assumed to be zero for low Mach number situations [24] . The jump conditions can be simplified to: