Time domain simulations of nonlinear thermoacoustic behaviour in a simple combustor using a wave-based approach Journal of Sound and Vibration

Lean premixed combustion chambers are susceptible to combustion instabilities arising from the coupling between the heat release rate perturbations and the acoustic disturbances. These instabilities are not desirable and knowledge of this complex mechanism is necessary in order to prevent or at least suppress them. A low order model is developed comprising a linear acoustic network and an improved analytical form for a flame describing function (FDF). The latter includes both the saturation of the amplitude of heat release rate perturbations and the change of phase lag relative to oncoming acoustic velocity fluctuations when the instability grows into a limit cycle. A stability map is constructed by moving the flame along the Rijke tube based on the eigenvalues resolved from the network. The acoustic model is then converted into the time domain and combined with the flame describing function to determine the evolutions of the heat release rate disturbances and velocity perturbations within the tube. It is shown that this method can be used to capture some quite intricate nonlinear behaviour of combustion instabilities and the results in the time domain are consistent with those predicted in the frequency domain. & 2015 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). model capture and explain the most significant features of nonlinear behaviour within combustion chambers. A low order network model was developed linking the nonlinear flame model with linear acoustic wave descriptions of a simple combustion chamber. This reduced order model was then converted into the time domain to reproduce some previously observed nonlinear features of combustion instabilities. The instantaneous amplitude of flow velocity disturbances was calculated with a proposed method and this value then used to determine the instantaneous gain and time delay of the flame describing function. Results showed this time domain can reproduce an important range of features nonlinear behaviour single tone instability, “ triggering , “ mode switching ” and frequency shift during the limit cycle.


Introduction
Combustion instability has been a major issue in the design stage of the industrial gas turbine for decades. It arises due to the coupling between the unstable combustion process and acoustic disturbances within the combustion chamber. The mechanism can be briefly described as: acoustic noise with a broad frequency bandwidth is produced during the combustion process [1]. These sound waves propagate inside the combustion chamber, interact with the boundaries and return back to the combustion zone with a time delay that depends on the size of the combustion chamber, disturbances of speed of sound and impedances at the boundaries of the combustion chamber. These pressure oscillations generate in turn perturbations of the flow field by modifying the local flowrate, reactant compositions or thermodynamic properties in the flame region producing heat release rate disturbances [2,3]. When these disturbances are synchronised, they amplify leading to an increase of acoustic energy in the system and a resonance is generally observed at specific tones. These self-sustained instabilities are more likely to happen in lean premixed combustion systems, which are preferred to minimise the emission of NO x for modern industrial gas turbine design [4].
Combustion instabilities are generally not desirable because they may lead to an early ageing of the combustion chamber or in extreme cases to severe structural damage [5]. The stability of a combustion chamber is determined by the balance between the energy gain from the heat released from unsteady combustion and the dissipation due to the viscous thermal damping [6,5], radiation from the boundaries [7] and various relaxation processes in flows with particles or droplets [2], which are considered to be proportional to the acoustic disturbance level [8][9][10]. Flame perturbations arise in different ways and originate mainly from the convection of hydrodynamic perturbations [2] or disturbances in the fuel and air injection supplies [11][12][13]. These mechanisms are susceptible to acoustic disturbances, which may lead to flame wrinkles that are convected along the flame front and modify the flame surface area due to the unbalance between the local flame speed and the flow velocity [14,15,5]. Improvement in the knowledge of the response of flames to acoustic disturbances is thus one of the main approaches in the strategy for reducing or eliminating combustion instabilities. For weak perturbations, the flame transfer function (FTF) has been used for decades to describe the linear response of heat release rate perturbations to disturbances in the velocity u 0 or the equivalence ratio ϕ 0 . The link can be mathematically expressed in the Laplace-domain: whereT u ðsÞ andT ϕ ðsÞ denote the responses of the flame to velocity perturbationsûðsÞ and incoming mixture inhomogeneitiesφðsÞ respectively. In practical combustion systems, the complex functionsT u andT ϕ vary with the modulation frequency and can be described in terms of a gain and phase lag. The gain, specifying the amplification effect of the flame on the incident acoustic flow perturbations, features the shape of low-pass filter and is cut-off at a certain frequency related to the shape of the steady flame and other properties [16][17][18][19]. The phase lag, also represented as the time lag of heat release perturbations with respect to oncoming flow perturbations, depends on certain quantities such as local thermodynamic conditions, turbulence properties [20], modulation levels [21,18,22] and thermal boundary conditions [13].
In particular for premixed flames, this time lag indicates the time necessary for an incoming or an impinging vortex convecting from the injector tip to the flame front, and is considered constant at low frequencies for small perturbations [17,18]. The transfer function is then combined with the acoustic response of the combustion chamber to analyse thermoacoustic stability and predict the eigenfrequencies and corresponding growth rates of potential instabilities [23,24]. In practice, most flames feature a nonlinear response to flow disturbances [25,16,18,21,26]. Saturation of the heat release rate amplitude or a phase lag change relative to the acoustic pressure (and hence a change in the Rayleigh source term 〈p 0 _ q 0 〉 driving combustion instabilities [27,28]) with the modulation level have been experimentally identified in many cases [29,30,21,22]. For laminar conical confined flames, the phase lag is almost constant for larger modulations revealing that the time delay is inversely proportional to the undulation frequency [21]. The nonlinear behaviour of four kinds of premixed flames with different geometries are summarised in [18] and it is found that the flame stretch becomes more progressive as the perturbation level increases, leading to a reduction of heat release rate disturbances compared to that anticipated by the linear theory. Kashinath et al. [31] attribute the nonlinearity to the unbalance between the convection speed of disturbances and the mean flow velocity and [32] finds that saturation is enhanced when the former is smaller than the latter. Perturbations at the flame holding point are convected along the flame front and produce additional wrinkles [33] and the saturation of the movement of holding point is identified as another reason for the nonlinearity of a 'V' shape premixed flame in a duct [34,35]. This is also found in swirling flames, in which the centre of flame exhibits a saturation process for high forcing levels [22]. Furthermore, the fluctuation in the flame angle may also make contributions to the saturation of flame wrinkle [36,26] and produces nonlinearities in the heat release rate perturbations. A limit cycle is thus established due to the balance between the saturation of heat release rate perturbations and acoustic losses [10].
It is clear that the linear flame transfer function is not enough to anticipate all of the characteristics of combustion instabilities; in particular it cannot predict the saturation amplitude or the time to limit cycle. The concept of a "flame describing function" (FDF) has been proposed, treating the response of heat release rate as a linear regime for small perturbation levels and a weakly nonlinear regime for larger disturbance amplitudes [8,30]. "Describing function" methods have already been used in nonlinear control community for decades due to their ease of use [37,38]. They assume that the nonlinearity can be described as a transfer function, which relies on the ratio of the magnitude of the fundamental term in the Fourier series of the output wave to the sinusoidal input wave. This limits their application to systems in which the magnitudes of harmonics are relatively weak compared to that of fundamental frequency component. This concept is used for example by Dowling [8,34] to account for the fact that the instantaneous heat release cannot go negative and is incorporated into a linear thermoacoustic network model for plane acoustic waves to predict the amplitude of the oscillations and the time evolution of the disturbances in further studies [39,40]. The experimentally determined FDF for the premixed flames of a matrix burner was recently combined with an acoustic network to predict the stability map, in which the eigenfrequency is varied by changing the length of combustion chamber [30,41,42]. This method was then extended to predict the evolution of the eigenfrequencies and corresponding growth rates with the increase of modulation level for turbulent swirling flames [43,44]. This method can explain phenomena such as "hysteresis", "triggering", "mode switching", frequency shift during the growth of instabilities and two mode limit cycles [30,43,42,45] limited by the experimental cost and the transient growth under certain conditions may not be observed due to the limit of instrument's response time.
This difficulty can be overcome numerically. A combined analytical/numerical approach is employed in this work to systematically analyse different features of nonlinear behaviour of combustion instabilities, based on a proposed analytical form for a flame model accounting for both the nonlinearities of the gain and the time delay. This nonlinear flame model is applied to a Rijke tube where the perturbations can be treated as longitudinal waves to simplify the complex nonlinear mechanism. The rest of the paper is organised as follows. The acoustic network, the nonlinear flame model and the improved boundary conditions treating the tube ends as dampers depending on frequencies are presented in Section 2.

Acoustic wave equations
Analysis is carried out on a simple model combustion chamber which is schematically described in Fig. 1. The configuration consists of a cylindrical Rijke tube with both ends open to the atmosphere. For a horizontal tube, the mean flow inside the tube is from external blowers [46,47], while for a vertical Rijke tube, the flow comes from natural convection arising from mean combustion products [48]. Denoting the distance along the tube axis by the vector x, the inlet and outlet of the tube are at x ¼ x 0 ¼ 0 and x ¼ x 2 ¼ l respectively. A premixed gaseous flame located at x ¼ x 1 is used as the heat source. The fully mixed fresh gases are ignited and rapidly turn to burned gases. In the analysis, the following assumptions are implemented: (1) The envisaged frequencies are assumed sufficiently small to consider the combustion zone to be "compact" compared to the acoustic wavelength and to only take into account the longitudinal waves. (2) The fluids before and after the combustion zone are assumed to be perfect gases. (3) The dissipation of acoustic waves throughout the tube is negligible and acoustic damping only happens at the ends of the tube. (4) There is no noise produced by the entropy waves formed during the unstable combustion processthese waves are assumed to leave the tube without interaction with the flow at the end of tube. (5) The flame is assumed to be always stabilized at the burner outlet. Flame intrinsic instabilities [49,50] and irregular response to strong disturbances [51,52], which may cause that the final state of the unstable system is not a limit cycle [53,45], have not been accounted for in this work. Now consider weak disturbances induced upstream of the combustion zone. Acoustic waves propagate in both directions. The upstream and downstream regions are indicated by the subscripts 1 and 2. Taking the flow to be composed of a steady uniform mean flow (denoted as ) and small perturbations (denoted as 0 ), it is thus possible to write out the pressure and velocity in the upstream and downstream regions [54]: u and c represent the longitudinal velocity and speed of sound respectively. A þ k and A À k (k ¼ 1; 2) denote the amplitudes of the inward and outward propagating acoustic waves respectively. Their link can be described by the reflection coefficient at the boundary. The pressure reflection coefficients at the inlet and outlet are characterised by R 1 and R 2 respectively. The entropy wave may also generate inward travelling acoustic waves at the outlet and make contributions to the inward acoustic wave. For example, this effect for chocked exits has been considered in [55]. However, these acoustic waves can be neglected in this work as the interaction between the entropy wave and weak subsonic flow will be very small [56]. In order to relate the upstream and downstream acoustic waves, the mass, momentum and energy conservation equations across the flame zone (x ¼ x 1 ) are combined with the perfect gas equation to obtain two equations which are independent of density ρ 2 which contains the entropy component (see for example [8]):

Flame model
In this work, only the transfer functionT u ðsÞ is accounted for. For industrial combustion chambers, the impedance of the fuel feed line is generally much larger than that of the air line. When the pressure disturbances interact with the injection unit (fuelþair), only the air line will respond significantly to these disturbances. There is thus a relation between the mixture composition and velocity disturbances which can be expressed as ϕ 0 =ϕ ¼ Àu 0 =u [55], or ϕ 0 =ϕ ¼ Àk 1 u 0 =u=ð1 þk 1 u 0 =uÞ proposed in [9], where k 1 is a coefficient of the mixture strength model. Theoretical models have been devised to analytically describe the transfer functionT u ðsÞ. For example, the G-equation model obtains the flame shape by using a kinematic model for the flame position response to oncoming velocity perturbations, and has been studied for various flame shapes: conical, V-shape, M-shape and matrix flames [14,59,34,15,60,35,13,31], with their detailed derivations collected in Lieuwen's book [61]. In this description, perturbations in the flame position are generated by velocity disturbances normal to the mean flame front position. Flow perturbations generated at the burner outlet are transported by the mean flow and act as a convective wave along the axial direction [14,60]. The flame surface area and hence the heat release rate perturbations can be found from the flame shape. A simpler model, which empirically captures the flame response shape for weak velocity disturbances (u 0 1 % 0) predicted by the G-equation for conical flames [60], is the famous n-τ model [62] filtered by a first-order low pass filter [59,8]: where ω c ¼ 2πf c indicates the angular cut-off frequency of the first-order filter, τ f represents the time delay and a f denotes the gain and is equal to 1 in this work. For the laminar conical premixed flame, f c % 3s L =2πr, where s L denotes the flame speed and r is the radius of burner nozzle [15]. This value increases if the flame is replaced by V-shape, M-shape or matrix flames [18,60].
The time delay τ f of the laminar conical premixed flame is equal to H=3u 1 for lower frequency band and increases with frequency, while the V-shape flame features a constant time delay of 2H=3u 1 [18], where H specifies the height of the flame. For a turbulent swirled flame, the profile of flame transfer function cannot be described by a simple n-τ model filtered by a low pass filter. The axial and swirled velocity induced vorticity fluctuations must be accounted for at high frequencies in the flame model, this is for example the one used by [63]. However, in our work, we only consider the laminar conical premixed flame and the proposed flame model is enough to capture the flame response. Time delay dispersion is not accounted for in this laminar flame model although it is an important mechanism for turbulent flames and creates uncertainties in the control of combustion instabilities [64].
A nonlinear flame describing function depends on s and velocity ratioû 1 =u 1 and it is assumed here can be decoupled as where the superscript^indicates the signal amplitude. The nonlinear function Lðû 1 =u 1 Þ describes the saturation of heat release rate with velocity perturbationsû 1 =u 1 , and Lðû 1 =u 1 Þ ¼Gðû 1 =u 1 ; 0Þ. A new mathematical link is proposed in this paper: α and β are two coefficients which determine the shape of the nonlinear model. For weak velocity disturbances, the link between b _ q = _ q andû 1 =u 1 is linear. For example, the linear region for the first case isû 1 =u 1 A ½0 0:4, with corresponding proportional coefficient L % 1 in the right figure. One may denote the upper limit of the linear region by the velocity disturbance ratio ðû 1 =u 1 Þ s . When the level of velocity disturbances exceeds ðû 1 =u 1 Þ s , L decreases and the heat release rate perturbations begin to saturate. It can be found that the saturation limit ðû 1 =u 1 Þ s is mainly determined by the coefficient α; ðû 1 =u 1 Þ s decreases as α increases towards unity. The smoothness of the saturation corner is controlled by the coefficient β. There are also some simpler models, such as that proposed by [65], which is the third-order Taylor expansion of a hyperbolic tangent function. However, this model becomes inaccurate at very larger velocity ratios. In following time domain simulations and our recent contributions to the feedback control of combustion instabilities [66], the response of the flame to very large velocity disturbances will also be accounted for. The proposed model in this paper guarantees saturation for larger velocity perturbation ratios.
One may also introduce a simple nonlinear model of the time delay, using the mathematical description: where τ 0 f means the time delay whenû 1 =u 1 ¼ 0 and τ f N is a time delay to describe the change of τ f as L changes. The final flame describing function thus can be expressed as Assuming the mean heat release rate per unit surface area to be _ q ¼ ρ 1 u 1 C p ðT 2 À T 1 Þ [54] and C p to be constant for the temperature band envisaged in this paper, it is thus possible to substitute the flame model into Eq. (3) The eigenvalues s k ¼ σ k þ j2πf k ðk ¼ 1; 2; …Þ can be obtained by solving the equation detðZðsÞÞ ¼ 0, where detðZðsÞÞ denotes the determinant of the coefficient matrix of the above equation.
In the following analysis, the parameters used for the analysis and time domain simulations are listed in Table 1 unless otherwise stated. The natural frequency of the cold Rijke tube with perfect acoustic reflections at the boundaries is given by f ⋆ k ¼ kc 1 =2l ¼ 173:6k Hz ðk ¼ 1; 2; …Þ. These values are used to normalise the frequency in later work.

Boundary conditions
A plane acoustic wave propagating along the tube will be at least partially reflected at the end. The reflection and transmission at the end can be characterised by the pressure reflection coefficient R. Herein, we only account for the situation of an open tube although extension to other types of tube boundaries would be straightforward. Sound waves corresponding to the dominant mode can be completely reflected if the wave length is much larger than the diameter of the tube, which means that R ¼ À1. Moreover, radiation of sound from the end must be accounted for once the diameter is increased. Furthermore, it has previously been shown experimentally that the magnitude of the reflection coefficient may exceed unity when the Mach number is considerable [69,58]. Despite the likely importance of the boundary behaviour to thermoacoustic stability, accounting for the end radiation in many previous Rijke tube studies has been very simplistic [57,48]. In this section, we seek to investigate its effect more carefully. It is shown in [7] that the reflection coefficient acts as a low pass filter and its amplitude may attain a maximum value of unity at low frequency. The tube also appears to be acoustically longer than its physical length and the end correction must be accounted for in calculating the eigenfrequency. The reflection coefficient for an unflanged pipe end can be approximately expressed in the Laplace domain as [7,69] for 2πf τ r {1. τ r ¼ r=c and r denotes the radius of the tube. An end correction of 0:61r is included in the expression. For a flanged pipe end, the reflection coefficient is changed to [69] RðsÞ ¼ Àð1 þτ 2 r s 2 Þe À 1:64τr s It is clear that both the radius and flow temperature affect the pressure reflection coefficient. The magnitude of the reflection coefficient acts as a damping factor and decreases with radius and frequency. This is further validated by the examples shown in Fig. 4. The eigenfrequencies decrease due to the increase in the end correction. The growth rate of the higher mode is more sensitive to the change of radius. This means that thermoacoustic instabilities are less likely to occur at higher frequencies in the Rijke tube. Note that the acoustic losses at the flanged end are larger than those at the unflanged end.

Time domain simulations of nonlinear combustion instabilities
A commonly used approach for analysing nonlinear thermoacoustic instabilities in the time domain is the Galerkin discretization, which involves expanding the pressure perturbation as a Galerkin series [70][71][72]. The partial differential thermoacoustic wave equation is then reduced to a set of ordinary differential equations with coefficients dependent on the nonlinear acoustic loss or heat sources [73][74][75]. However, this approach is based on a continuous spatial approximation [31] and is limited to simple combustor configurations, although recent efforts have accounted for the temperature jump across Table 1 Parameters used in the analysis and time domain simulations. They are fixed unless otherwise stated. the flame [76,77]. Other discretization methods such as the modal expansion technique, have been used to produce a statespace network for the time domain simulations [64]. The wave-based approach is the one used in the present work. It has the benefits of extending easily to combustion chambers with varying cross sectional area and annular geometries and can account for a wide range of acoustic boundary conditions [23,40]. It can account for discontinuous acoustic velocity and temperature distributions using jump equations [54,78] (or see Eq. (2b)). Time domain simulations are based on the inverse Laplace transform of Eq. (8) so that the network can be mathematically expressed as The exponential terms are transformed to time delays. The reflected acoustic waves (A þ 1 ðtÞ and A À 2 ðtÞ) at the ends of tube are calculated by the convolution of the inverse Laplace transform of the reflection coefficients and the incident acoustic waves (A À 1 ðt À τ À 1 Þ and A þ 2 ðt Àτ þ 2 Þ). The tube ends are considered as unflanged and Eq. (9) is used as the pressure reflection coefficient. However, this expression only works for lower frequencies. The magnitude of jRj ¼ j1 À 2π 2 τ 2 r f 2 j decreases to zero at the cut-off frequency f R;c ¼ ffiffiffi 2 p c=2πr and then increases to infinity with frequency. Such inappropriate boundary conditions at higher frequencies may cause high frequency unstable modes leading to errors in the time domain simulations. To overcome this problem, a second-order low pass filter with a time delay is used to fit the frequency response of Eq. (9) in the frequency band ½0; 10f The magnitudes of these new reflection coefficients decrease to zero at high frequencies and can be safely used in time domain simulations as this paper considers only cases in which the first two acoustic modes can be unstable. Their inverse Laplace transforms are used in Eqs. (11b) and (11c). However, it should be noted that Eq. (9) is still used in frequency domain calculations since the cut-off frequency f R;c c 10f rate σ k and angular frequency ω k ¼ 2πf k respectively: where B k is the amplitude of kth mode and θ k denotes the phase. The value ð P n k ¼ 1 B 2 k Þ 1=2 is used to measure the amplitudê u 1 of the combined signal u 0 1 ðtÞ in this work. The eigenfrequency of the kth mode can be approximately expressed as f k ¼ kf 1 in the Rijke tube [79]. Denoting the period of the first mode as τ ¼ 1=f 1 , one can obtain the RMS value along the time interval ½t À τ=2; t þτ=2: Assuming that σ k τ-0, the above equation can be simplified as An external acoustic signal p L (t) from a loudspeaker located at the inlet of the Rijke tube is used to provide an additional disturbance. This signal is used to produce a very weak disturbance ( À 100 dB) at the beginning of the simulation to stimulate any instabilities of the system. It is also used to investigate the "triggering" of combustion instabilities.
An open loop transfer function (OLTF), which is generally used for control purposes [48,66], is utilised here to validate the time domain simulation code. The OLTF is from the sound signal p L produced by the loudspeaker to the signal p r measured by a reference microphone located at x ¼ x r ðx r 4 x 1 Þ. The profile of the OLTF changes with the amplitude of the inlet velocity perturbationsû 1 =u 1 due to nonlinearity. Note that the assumption of "weak nonlinearity" has been made: this assumes that the output harmonic content contains the same frequency as the input harmonic content, but with a gain and phase shift that is amplitude, as well as frequency dependent. An analytical OLTF can be easily obtained from Eq. (8): where the matrix W is expðÀτ 1 sÞ X 12 À Y 12R2 ðsÞ expðÀτ 2 sÞ X 21 À Y 21R1 ðsÞ expðÀτ 1 sÞÀΩGðû 1 =u 1 ; sÞðR 1 ðsÞ expðÀτ 1 sÞÀ1Þ X 22 À Y 22R2 ðsÞ expðÀτ 2 sÞ " # The OLTF from time domain simulations can be calculated using where fft denotes the fast Fourier transform. Validation is carried out with a linearly stable case: nonlinearity does not need to be accounted for which means that L ¼ 1. A chirp signal with a frequency swept linearly from 10 Hz (0:06f ⋆ 1 ) to 2000 Hz (11:5f ⋆ 1 ) is generated by the loudspeaker. The corresponding pressure p r (t) is plotted in Fig. 6(a). The calculated OLTFs together with that deduced from Eq. (8) are shown in Fig. 6(b). The time domain simulation results match perfectly with those from the analytical model, showing that our time domain network accurately reproduces the behaviour of frequency domain calculations.
The coefficients of the flame model are listed in Table 1. The cut-off frequency of f c ¼50 Hz means that only the first two modes have the potential for instability. The time delay τ 0 f ¼ 4 ms and τ N f ¼ À3 ms which means that the time delay  decreases with the velocity perturbation level. This is consistent with that of the laminar conical premixed flame [21]. Following the experimental method used by [30,42,80,53], the flame is also moved along the Rijke tube by a step of 0.01 m to build a stability map for the first two modes, as shown in Fig. 7. Only positive values designating the unstable zone are presented. The higher modes are always stable. The velocity ratio ðû 1 =u 1 Þ LC k ðk ¼ 1; 2Þ when the limit cycle is established can be identified by the outer boundary of the contours. The contours, featuring the shape of "Arnold tongues" [81], are vertical and "straight" when τ N f ¼ 0 (as illustrated in Fig. 7(a)), but tend to bend and intrude the "territory" of other modes for increased values of τ N f (the effect of time delay nonlinearity is illustrated in Fig. 7(b)). x 1 is carefully chosen to select cases with different stability characteristics to illustrate different phenomena in combustion instabilities. They are summarised in Table 2. The calculation is stepped forward with a time step of Δt ¼ 1 Â 10 À 5 s. The calculation time period is from 0 to 5 s.

Simulation of instabilities exhibiting a single tone (cases 1 and 2)
We consider first cases in which combustion instabilities exhibit a pure tone. Case 1, for which the flame position is x 1 =l ¼ 0:10, is the first example. From the trajectories of growth rate and normalised eigenfrequency, f k =f ⋆ k ; ðk ¼ 1; 2Þ, witĥ u 1 =u 1 , plotted in Fig. 8(a), it can be seen that the growth rate of the second mode is always negative indicating that it is stable. However, part of the trajectory of the first mode has a positive growth rate. The growth rate equals 12 s À 1 when u 1 =u 1 ¼ 0, meaning that the first mode is linearly unstable and weak disturbances increase exponentially with an envelope of expð12tÞ. The growth rate decreases withû 1 =u 1 and becomes zero whenû 1 =u 1 ¼ 0:14, hence this is the velocity fluctuation amplitude at which we expect limit cycle oscillations to occur. Saturation has occurred despite very weak change in flame gain/response, L (see dashed line case in Fig. 2). The implied limit cycle is validated by calculations in the time domain. Fig. 9(a) shows the time evolution of the signals u 0 1 =u 1 and its amplitudeû 1 =u 1 .The value calculated from Eq. (14) matches perfectly with the amplitude of the oscillating signal. Asû 1 =u 1 increases, L decreases and the growth rate drops until it equals to zero which leads to the saturation of oscillations. The oscillating frequency can be resolved by short-term Fourier transform shown in the bottom graph of Fig. 9(a). It is also possible to substitute the time domain resolvedû 1 =u 1 into Eq. (8) to solve the first eigenfrequency in the frequency domain. The two frequencies match perfectly again validate the time domain code. The heat release rate perturbations not shown here are found to saturate before the velocity disturbances. This is a result of the form of the nonlinear model, in which the growth rate at which heat release rate perturbations saturate may still be positive. It is interesting to plot the normalised heat release rate perturbations _ q 0 = _ q as the function of u 0 1 =u 1 as shown in Fig. 10(a). The trajectories expand and approach the limit cycle as time approaches infinity. When we change the flame position, for example x 1 =l ¼ 0:67 for case 2, the characteristics of the combustion instabilities also change. From Fig. 8(b), we observe that the normalised velocity perturbation ðû 1 =u 1 Þ LC 1 when the limit cycle is established is increased to 0.55. Asû 1 =u 1 increases, the growth rate also increases and reaches the maximum value when   Fig. 9(b)) and the limit cycle trajectories (shown in Fig. 10(b)). It is seen from Fig. 10(b) that the phase lag between _ q 0 = _ q and u 0 1 =u 1 also changes during the establishment of limit cycle.

Simulation of combustion instabilities exhibiting mode switching (case 3)
The phenomenon of "mode switching" is investigated in this section. The flame position x 1 is carefully selected so that the combustion instability consists of two unstable modes. The trajectories of the growth rate of the first two modes are shown in Fig. 11 for case 3. Whenû 1 =u 1 o0:12, the second mode is dominant, but asû 1 =u 1 increases, the growth rate of the second mode decreases sharply, becoming zero whenû 1 =u 1 ¼ 0:14. The first mode then becomes the dominant mode, saturating whenû 1 =u 1 ¼ 0:41. These predictions are consistent with time domain simulations as shown in Figs. 12 and 13 (a). From t¼ 0 s to t¼2.2 s, the two signals establish a preliminarily limit cycle corresponding the second mode. With increasingû 1 =u 1 , the first mode grows and the instability slowly switches to a new limit cycle with a lower frequency and a larger amplitude.

Simulation of combustion instabilities exhibiting mode triggering (cases 4 and 5)
We now consider the nonlinearly unstable system, for example case 4 shown in Fig. 14(a). The system is linearly stable with respect to weak amplitude disturbances (û 1 =u 1 o 0:11) but becomes unstable when the magnitude of disturbances exceeds a threshold value (ðû 1 =u 1 Þ TG 1 ¼ 0:11). This phenomenon is known as "triggering" and combustion instabilities may be excited by low amplitude disturbances, such as the background noise within the combustion chamber [5,82]. Recently, [73][74][75]83] showed that the non-normal transient growth induced by small finite disturbances can also cause "triggering". However, non-normality is out of scope of this paper and we only focus on the nonlinearity of combustion instabilities.  It is now interesting to express the velocity disturbance as u 0 1 ¼ u 0 1;se þu 0 1;ex ¼ u 0 1;se þ p L =ðρ 1 c 1 Þ, where u 0 1;se indicates the self-excited velocity disturbance and p L represents the external acoustic disturbances from the loudspeaker. As such, when p L does not exceed the threshold ρ 1 u 1 c 1 ðû 1 =u 1 Þ TG 1 ¼ 156 Pa, the flame is linearly stable. This is for example illustrated with the case shown in Fig. 15. The signal u 0 1 =u 1 is always zero until the presence of an external disturbance at t¼ 0.95 s. The external disturbance p L is a Gaussian white noise in the time interval [0.95 1.05] s. The amplitude of p L is 71 Pa which is lower than the threshold 156 Pa, and thus all the modes are damped leading to the disturbances being attenuated. In the time interval [1.45 1.55] s, a second disturbance with a larger amplitude of 284 Pa is used to trigger instability. The external disturbance is strong enough to excite the instability and a limit cycle of the first mode is established.  "Triggering" also happens between different modes [42]. This is for example illustrated with case 5 shown in Fig. 14(b). The first mode is unstable and saturates whenû 1 =u 1 ¼ 0:10. The second mode is stable at the beginning. Its growth rate increases withû 1 =u 1 , turning positive whenû 1 =u 1 4 0:12. The growth rate of the second mode reaches its maximum value whenû 1 =u 1 ¼ 0:17, then decreases and becomes zero whenû 1 =u 1 ¼ 0:20. Without external excitations, self-sustained combustion instabilities may be generated and a limit cycle is established for the first mode. This is also validated by the time domain simulation shown in Figs. 16 and 13(b). Once an external disturbance p L with a sufficiently large amplitude is added, the existing balance will be broken and another limit cycle will be established with a higher eigenfrequency.

Conclusions
This paper has presented an improved analytical form for a flame describing function (FDF). It is based on an extension of the n-τ model with a first-order low-pass filter, modified to account for the nonlinear saturation of the heat release rate perturbation amplitude and the variation of the time delay with respect to velocity disturbances. The aim of this model is to capture and explain the most significant features of nonlinear behaviour within combustion chambers. A low order network model was developed linking the nonlinear flame model with linear acoustic wave descriptions of a simple combustion chamber. This reduced order model was then converted into the time domain to reproduce some previously observed nonlinear features of combustion instabilities. The instantaneous amplitude of flow velocity disturbances was calculated with a proposed method and this value was then used to determine the instantaneous gain and time delay of the flame describing function. Results showed that this time domain network can reproduce an important range of features of nonlinear behaviour such as single tone instability, "triggering", "mode switching" and frequency shift during the limit cycle. This nonlinear low order model provides a new platform for the analyses of combustion instabilities (the code is called OSCILOS and is available freely) and motivates further developments in more complex configurations.