Feedback control of combustion instabilities from within limit cycle oscillations using H∞ loop-shaping and the ν-gap metric
===============================================================================================================================
* Jingxuan Li
* Aimee S. Morgans
## Abstract
Combustion instabilities arise owing to a two-way coupling between acoustic waves and unsteady heat release. Oscillation amplitudes successively grow, until nonlinear effects cause saturation into limit cycle oscillations. Feedback control, in which an actuator modifies some combustor input in response to a sensor measurement, can suppress combustion instabilities. Linear feedback controllers are typically designed, using linear combustor models. However, when activated from within limit cycle, the linear model is invalid, and such controllers are not guaranteed to stabilize. This work develops a feedback control strategy guaranteed to stabilize from within limit cycle oscillations. A low-order model of a simple combustor, exhibiting the essential features of more complex systems, is presented. Linear plane acoustic wave modelling is combined with a weakly nonlinear describing function for the flame. The latter is determined numerically using a level set approach. Its implication is that the open-loop transfer function (OLTF) needed for controller design varies with oscillation level. The difference between the mean and the rest of the OLTFs is characterized using the *ν*-gap metric, providing the minimum required ‘robustness margin’ for an H∞ loop-shaping controller. Such controllers are designed and achieve stability both for linear fluctuations and from within limit cycle oscillations.
* combustion instabilities
* flame describing function
* *ℋ*∞ loop-shaping
* *ν*-gap metric
* level set approach
* low-order controller
## 1. Introduction
It is desirable to operate modern industrial gas turbines and aeroengines under lean premixed combustion conditions in order to reduce NO*x* emissions. However, lean premixed systems are highly susceptible to combustion instabilities arising from the coupling between the heat release rate perturbations and the acoustic disturbances within the combustion chamber [1]. These instabilities are nearly always undesirable as they lead to large oscillation amplitudes and damaging structural vibrations of the combustion chamber [2].
The stability of a combustion chamber is determined by the balance between the energy gained from the flame/acoustic interactions and various dissipation processes [1,3,4]. Active feedback control can be used to interrupt the coupling between the acoustic waves and unsteady heat release and prevent or suppress instability. The design of most types of controller requires prior knowledge of how the combustor responds to actuation—the ‘open-loop transfer function’ (OLTF) between monitored thermodynamic properties within the combustion chamber (typically one or more pressure measurements) and the actuator signal(s) [4].
Although linear models may be accurate at low perturbation levels, the dynamics of real unstable combustors become dominated by nonlinear mechanisms once amplitudes grow sufficiently. In many combustors, including gas turbine and aeroengine combustors, it is notable that even during large amplitude oscillations the behaviour of the acoustic waves remains linear, with nonlinearity arising almost entirely from the way in which the flame’s unsteady heat release responds to oncoming flow disturbances [5–7]. (Acoustic waves in rocket engines are, however, typically nonlinear [8].) This flame nonlinearity is ‘weak’ such that a ‘describing function’ provides a good model of the nonlinear flame response. This assumes that the flow forcing amplitude level is a quasi-steady parameter, and that the dominant frequency of the unsteady heat release rate response matches the incoming flow forcing frequency, but with a gain and phase shift that depend on the forcing level as well as the frequency. This gives rise to a ‘family’ of frequency responses that depend on forcing level. It has been shown that this weak nonlinearity can cause saturation into limit cycle either via saturation of the heat release rate amplitude or via a change in phase between the heat release rate and acoustic pressure as the modulation level increases [6,9]. This may cause the dominant unstable frequency to change, which, in turn, alters the OLTF and complicates the controller design. Although some recent work has identified nonlinear states other than limit cycle oscillations as resulting from combustion instability [10–12], the focus of this work will be limit cycle oscillations, which are by far the most common reported nonlinear state.
The design of linear feedback controllers is typically based on the OLTF corresponding to small, linear disturbances. Such feedback controllers are then activated from within nonlinear limit cycle oscillations, during which the linear OLTF will not be valid. It is typically observed that such controllers successfully stabilize oscillations despite this paradox [13–15]. However, their effectiveness is not guaranteed, and if they fail then there has, until now, been no framework for systematically improving their design.
In this work, the concept of a flame describing function (FDF) is combined with the *ν*-gap metric [16] in order to design controllers which are guaranteed to achieve stabilization across all oscillation levels, from small linear fluctuations to large limit cycle oscillations. In modern robust control systems theory, the concept of a *ν*-gap between open-loop plants represents the most useful measure of distance between systems when one is concerned specifically with applying feedback to those systems. The *ν*-gap metric fits very naturally into H∞ loop-shaping controller design, providing a bound on the minimum required ‘robustness stability margin’ for an H∞ loop-shaping controller [16,17]. Modern linear robust control, of which H∞ loop-shaping is one methodology, has been successfully used in diverse fields, such as control of combustion instabilities [18–20], control of developing flow [21], voltage converter design [22] and vehicle oscillation damper optimization [23], to achieve control even in nonlinear systems with limit cycles. H∞ loop-shaping combines classical loop-shaping, which can obtain tradeoffs between good performance and robust stability, with modern H∞ optimization in order to guarantee closed-loop stability and a level of robust stability at all frequencies [16,17]. It can be applied to plants with high dimensions and is computationally inexpensive [19], at least for plants of low and moderate order. The resulting controllers generally have orders comparable to those of the plants.
It is typical to approximate the OLTF as a high-order rational transfer function (RTF) for the purposes of controller design, meaning that high-order controllers will be generated. These are less reliable to implement than low-order controllers. As a final step, this work proposes a low-order controller design method, which combines H∞ loop-shaping, the *ν*-gap metric and a low-order fitting procedure.
The remainder of the paper is organized as follows. Section 2 presents the ducted laminar flame combustor model used as a test case in the paper, including the coupling of linear modelling for the acoustic waves with the nonlinear flame response. Section 3 demonstrates the determination of the nonlinear FDF from simulations using a level set approach (LSA), ensuring that our flame model is grounded in flow physics. Section 4 analyses the resulting nonlinear thermoacoustic instabilities for different flame positions. Section 5 presents the design of linear controllers which are guaranteed to be stabilizing from within limit cycle oscillations, using H∞ loop-shaping and the *ν*-gap metric. This includes the design of standard high-order controllers and a method for obtaining low-order controllers. The viscothermal damping of acoustic waves in ducts is presented in appendix, and conclusions are drawn in the final section.
## 2. Combustor model
A geometrically simple model combustor is considered for this study. The geometrical simplicity permits us to embed the advanced modelling features that capture the important phenomena at play in real combustors. This includes linear acoustic wave behaviour with realistic boundary and viscothermal losses, and complex nonlinear flame kinematics for sufficiently high amplitude perturbations.
The model considered is a ducted laminar flame combustor, often known as a Rijke tube. It is shown schematically in figure 1, comprising a cylindrical Rijke tube with both ends open to the atmosphere. Denoting the distance along the duct axis by the vector *x*, the tube inlet and outlet are at *x*=*x*=0 and *x*=*x*2=*l*, respectively. The inner diameter of the tube is *D*, and a laminar premixed methane–air conical flame is located at *x*=*x*1=*x**f*. The Bunsen-type burner has an axisymmetric geometry with outlet diameter *d*. A loudspeaker is located at the bottom of the tube as an actuator for control, with a pressure sensor located a distance *x*=*x**r* downstream of the combustion zone. The modelling analysis makes the following assumptions: (i) the combustion zone is considered ‘compact’ compared with the acoustic wavelength, and we limit attention to low frequencies where we need to account only for longitudinal waves [24]. (ii) The fluids before and after the flame are assumed to be perfect gases [24]. (iii) There is no noise produced by the entropy waves formed during the unstable combustion process—these waves are assumed to leave the tube without interaction with the flow at the end of tube [25]. (iv) The flame remains stabilized at the burner outlet. Flame intrinsic instabilities [26] and irregular response to strong disturbances [27] are not accounted for.
[Figure 1.](http://rspa.royalsocietypublishing.org/content/472/2191/20150821/F1)
Figure 1.
Schematic view of the Rijke tube and the feedback control configuration.
The flow can be taken to be composed of a steady uniform mean flow (denoted ¯ ) and small perturbations (denoted ′). Harmonic time variations are considered for which all fluctuating variables have the form a′=a^ eiωt, where *ω* is complex angular frequency. The regions upstream and downstream of the flame are indicated by the subscripts 1 and 2, respectively. Then, the pressure, *p*, and longitudinal velocity, *u*, in the two regions can be expressed in terms of the amplitudes An+ and An− of the upward and downward propagating acoustic waves, respectively, their wavenumbers, k±=k¯±+Δk±, where k¯±=ω/(c¯±u¯) and Δ*k*± are wavenumber corrections for the viscothermal damping effects (see appendix), and the speed of sound, *c*, and density, *ρ* [28]. The pressure, *p*, and longitudinal velocity, *u*, in the two regions are expressed as pn(x,t)=p¯n+p^n(x) eiωt=p¯n+(An+eikn+(x−xn−1)+An−e−ikn−(x−xn−1)) eiωt(n=1,2)andun(x,t)=u¯n+u^n(x) eiωt=u¯n+1ρ¯nc¯n(An+eikn+(x−xn−1)−An−e−ikn−(x−xn−1)) eiωt.}2.1The acoustic wave strengths, An+ and An−, are related by pressure reflection coefficients1 at the tube boundaries—the inlet and outlet pressure reflection coefficients are characterized by *R*1 and *R*2, respectively.
For an open duct, radiation of sound from the duct ends can be accounted for using frequency-dependent reflection coefficients that act as low-pass filters [29]. The tube also appears acoustically longer than its physical length, and this end correction can be accounted for. The reflection coefficient for an unflanged pipe end can be approximately expressed in the Laplace domain as [29,30] R=−(1+2τd2s2)e−0.61τdsfor 2πfτd≪1.2.2where τd=D/c¯. It was confirmed in [31] that the gain of this reflection coefficient decreases with both duct diameter and frequency. The reflection coefficients at the inlet and outlet of the Rijke tube are denoted by *R*1 and *R*2, respectively.
In order to relate the acoustic wave strengths across the flame, the linearized mass, momentum and energy conservation equations across the flame zone (*x*=*x*f) are combined with the perfect gas equation. Denoting the jump across the flame as [⋅]12, the amplitude of the flame’s heat release rate perturbation as Q˙^ and the cross-sectional area of the duct as *S*, this yields two equations which are independent of the downstream density, *ρ*2, which contains an entropy component [32] S2S1[p^(xf)]12+ρ¯1u¯1[u^(xf)]12+[u¯]12(ρ¯1u^1(xf)+p^1(xf)c¯12u¯1)=02.3and γγ−1[S(p¯u^(xf)+p^(xf)u¯)]12+ρ¯1u¯1S1[u¯u^(xf)]12=Q˙^.2.4Note that upstream of the flame, *S*1=*π*(*D*2−*d*2)/4. *γ* is the ratio of specific heats of the gases, considered constant throughout the duct owing to the small temperature jump across the flame. Substituting the linearized pressure and velocity expressions from equation (2.1) together with the boundary conditions into the above two conservation equations, the governing equations linking the acoustic wave strengths before and after the flame are obtained in terms of *s*, the Laplace transform variable [33] (where sn=s+Λns, *n*=1,2, the second term accounting for viscothermal damping with the coefficient *Λ**n* defined in equation (A 2)), time delays τn±=(xn−xn−1)/(c¯n±u¯n) and Ξ=(S1/S2)(ρ¯2c¯2)/(ρ¯1c¯1). The governing equation is expressed as [1−R100M1e−τ1+s1M2eτ1−s1−M3−M4M5e−τ1+s1−M6eτ1−s1−1ΞM71ΞM800−R2e−τ2+s2eτ2−s2]⏟M0[A1+A1−A2+A2−]⏟A=[00−γ−1c¯1S1Q˙^0].2.5The matrix coefficients, Mn, are given in appendix. In order to close equation (2.5) that contains five unknowns (An± and Q˙^), a flame model relating Q˙^ to the acoustic fluctuations is needed. In this work, we use an FDF, whose form is derived numerically in §3.
Although the combustor configuration contains only two acoustic modules (the two regions upstream and downstream of the flame), the system is complex owing to the frequency-dependent pressure reflection coefficients, viscothermal damping terms and the presence of flame nonlinearity. A numerical solver is used to solve the nonlinear eigenproblem [34]. The solver is part of *OSCILOS*, an open source combustion instability low-order network simulation tool. It represents complicated combustor geometries as a network of simple connected modules, modelling the acoustic wave behaviour analytically using linear wave-based methods, and able to incorporate a wide range of linear and nonlinear flame models. It can predict thermoacoustic modal frequencies, growth rates, mode shapes and the time evolution of disturbances, and has been validated against experiments [35].
In this work, a relatively low cost numerical method, compared with large eddy simulations [35,36], is used to determine the FDF of the laminar conical premixed flame. This numerical method is described in §3.
## 3. Determining the flame describing function
In the model combustor, a laminar conical premixed flame within the duct generates a mean and fluctuating heat release rate. The heat release rate is taken to be directly proportional to the flame’s surface area, which means that an accurate flame model for insertion into equation (2.5) can be found by capturing the kinematic response of the flame to oncoming acoustic disturbances.
The instantaneous location of a moving flame front subjected to oncoming flow perturbations can be captured using the LSA proposed by Markstein [37]. This kinematic model is also known as the *G*-equation model, and treats the flame as an extremely thin interface (corresponding to the iscontour *G*=0), separating fresh reacting flow (denoted by *G*<0) from the burned gases (*G*>0). The flame front is assumed to propagate normal to itself with the velocity **U**⋅**n**−*s*L, where **U** is the fresh reacting flow velocity vector and in this work is two-dimensional, **U**=(*u**x*,*u**r*), owing to the axisymmetric configuration. *s*L is the flame burning speed and **n**=**∇***G*/|**∇***G*| is the normal vector. The transport equation for *G* is [38] ∂G∂t+U⋅∇G=sL|∇G|.3.1
The mutual interactions of the flame front kinematics and the hydrodynamic flow field are not accounted for [39], meaning the fresh mixture flow field above the burner outlet and inside the flame can be treated using a simple model which is independent of the flame evolution. It has been shown that when the *G*-equation model is combined with a simple incompressible velocity perturbation model in the fresh mixture, complex flame front evolution phenomena, including the ‘pinch-off’ that occurs for strong flow perturbations, can be captured [40,41]. This model was recently validated using direct numerical simulations [42], and is expressed as ux=u¯1+u1′=u¯1+u^1cos(kx−ωt)andur=ur′=kr2u^1sin(kx−ωt).3.2The composition of the fresh reacting flow is considered constant, and the dependence of the flame burning speed on the flame front curvature is taken to be sL=sL0(1−Lκ), where *s*L0 is the burning speed of a laminar planar flame, L is the Markstein length, which depends on mixture equivalence ratio *ϕ* and thermal conditions [43] (and hence is constant in this work), whereas *κ*=**∇**⋅**n** is the local flame curvature [37,43].
The *G*-equation (equation (3.1)) combined with the incompressible flow velocity model (equation (3.2)) and the flame curvature effect is solved numerically using a fifth-order weighted essentially non-oscillatory (WENO) scheme [44] for spatial discretization and a third-order total variation diminishing (TVD) Runge–Kutta scheme for time integration [45]. The spatial and temporal resolutions are fixed in all calculations: Δ*x*=Δ*r*=5×10−3*d* and Δt=5×10−4d/u¯1, respectively, where *d* is the diameter of the burner outlet. To significantly reduce the computational cost, a local LSA is employed, only accounting for the grids around the flame front [46]. The boundary condition for the centreline of the flame is symmetry—flame axisymmetry is exploited to simulate only half of the flame. The flame base is considered always attached to the burner lip, with *G* at the flame base reinitialized to 0 at every simulation time step. No boundary condition is needed for other boundaries, because the local LSA is implemented.
For premixed flames, the heat release rate can be calculated using Q˙(t)=2πρ¯1sL0ω˙T∬A(1−Lκ)|∇G|δ(G) dr dx,3.3where ρ¯1 denotes the mean density of the fresh reacting flow, ω˙T indicates the heat release rate per volume from combustion and A indicates the whole space of the LSA simulation. *δ*(*G*) is the delta function, representing integration over the flame front, which is performed using a high-order scheme [47].
The weakly nonlinear flame model to be extracted from these simulations is in the form of an ‘FDF’, expressed as F(u^1/u¯1,s)=Q˙^1(u^1/u¯1,s)/Q˙¯u^1(s)/u¯1,3.4where Q˙^1(u^1/u¯1,s) is the Laplace transform of the fundamental term of heat release rate perturbations. Note that this can be thought of as an ‘input amplitude dependent’ transfer function or frequency response.
To validate the LSA, a flame for which experimental data are available is firstly simulated. A laminar conical methane–air premixed flame with burner outlet diameter *d*=20 mm, equivalence ratio *ϕ*=1.0, mean velocity u¯1=1.5 m s−1, laminar flame speed *s*L0=0.368 m s−1 and Markstein length L=1 mm [43] is considered. Figure 2 compares the simulated flame front location from the LSA with a four-colour Schlieren image from experiments, when the perturbation frequency and level are *f*p=100 Hz and u^1/u¯1=0.1, respectively. The flame front is seen to respond strongly even for a relatively weak perturbation level. The prediction matches the experimental result with the flame cusps accurately captured. The evolution of the flame front throughout an entire forcing cycle for modulation at *f*p=50 Hz and u^1/u¯1=0.15 is shown in figure 3*a*. At t/T=12, where T=1/fp, pinch-off occurs and a flame pocket is produced. The frequency spectrum of the heat release rate perturbations is shown for 50 Hz forcing in figure 3*b*. Although the nonlinear flame response causes harmonics at high frequencies, the fact that the dominant flame response frequency exactly matches the forcing frequency validates the ‘weakly’ nonlinear assumption of the FDF.
[Figure 2.](http://rspa.royalsocietypublishing.org/content/472/2191/20150821/F2)
Figure 2.
Comparisonbetween locations of the flame front deduced using LSA (*a*) and Schlieren techniques [48,49] (*b*) of a laminar conical premixed flame when the flame is modulated at *f*p=100 Hz for a perturbation level u^1/u¯1=0.1. t/T=120, where T=1/fp. In the Schlieren image, the yellow part represents the flame front, and the red part indicates the hot plume surrounding the flame. (Online version in colour.)
[Figure 3.](http://rspa.royalsocietypublishing.org/content/472/2191/20150821/F3)
Figure 3.
(*a*) Flame front at four instants in a forcing cycle. t/T=14,12,34 and 1. *f*p= 50 Hz. u^1/u¯1=0.15. (*b*) Spectra of u1′/u¯1 (represented by the dashed line and corresponding to the left vertical axis) and Q˙′/Q˙¯ (represented by the continuous line and corresponding to the right vertical axis).
The FDF, F, calculated using the LSA is compared with experimental results for a similar laminar conical premixed flame [50], with almost the same burner outlet diameter. The simulation and experimentally measured FDFs are shown in figure 4. The FDF has the shape of low-pass filter at constant u^1/u¯1, with nonlinearity evident at higher frequencies—the gain of FDF decreases with increasing u^1/u¯1, but with the response linear at low frequencies even for a large u^1/u¯1=0.214. The phase dependence on forcing amplitude begins at a higher frequency than the gain dependence. The gain and phase lag of the FDF both generally decrease with modulation level, which is consistent with the experimental results. It can clearly be concluded that the LSA captures both the flame front motion and the main features of the FDF of laminar conical premixed flames.
[Figure 4.](http://rspa.royalsocietypublishing.org/content/472/2191/20150821/F4)
Figure 4.
FDF of a laminar conical premixed methane–air flame deduced from the LSA (*a*) and experiments [50] (*b*). Top graph: gain. Bottom graph: phase. Mean flow velocity is u¯1=2.12 m s−1 and equivalence ratio is *ϕ*=1.08. Burner diameter is *d*=22 mm. (Online version in colour.)
Simulation is then performed for the combustor test-case flame. For a laminar conical premixed flame, the cut-off frequency of the FDF is inversely proportional to the diameter of the burner outlet [6,51,52]. To guarantee that a sufficiently large heat release feeds the unstable thermoacoustic mode, the diameter of the burner outlet is reduced to *d*=6 mm (meaning that the cut-off frequency is approximately three times that of the FDF shown in figure 4), with the other flow parameters unchanged. Simulations are performed for modulation frequencies ranging from 10 to 500 Hz in intervals of 10 Hz, and with normalized velocity perturbation levels ranging from 0.025 to 0.4 in intervals of 0.025. The FDF is shown in figure 5. The shape remains similar, with the response generally falling off with both frequency and modulation level. Although |F| does not strictly decrease with increasing u^1/u¯1 above 300 Hz, the accuracy of the flow velocity perturbation model (equation (3.2)) breaks down at high frequencies. With the dominant unstable frequency in the Rijke tube being 210 Hz, the FDF results can be considered accurate below 300 Hz. The FDF provided by the LSA simulations (denoted FLSA) is discrete in both frequency and forcing amplitude and is limited to the frequency range [10, 500] Hz. The FDF can be considered a set of flame transfer functions (FTFs)2 for different u^1/u¯1, each of which depends only on *s*. A fitting procedure is performed on each FTF using the Matlab command *fitfrd*3 . The resulting state-space models for different u^1/u¯1, each of order 14, are combined to give a ‘fitted FDF’, denoted F. The maximum error for each fitting procedure is smaller than 0.005, for *f*∈[10,500] Hz and u^1/u¯1∈[0.025,0.4].
[Figure 5.](http://rspa.royalsocietypublishing.org/content/472/2191/20150821/F5)
Figure 5.
The gain (*a*) and phase (*b*) of the FDF determined by the LSA for a laminar conical premixed methane–air flame. u¯1=1.5 m s−1. *d*=6 mm. *ϕ*=1.0. (Online version in colour.)
## 4. Simulation and analysis of the unstable combustor
By inserting the fitted FDF F into the acoustic wave equation (2.5), the thermoacoustic stability of the combustor is characterized by a matrix equation involving the acoustic matrix, **M******, from equation (2.5), a matrix, **M****1**, containing the FDF and the acoustic wave strength vector, A, from equation (2.5) (M0+M1)A=0,where M1=QF(u^1/u¯1,s)=[00000000e−τ1+s1−eτ1−s1000000]4.1where Q=T¯2/T¯1−1 and T¯ denotes the mean flow temperature. The thermoacoustic modes are given by the values of *s*=i*ω*=*σ*+i*ω*i at constant u^1/u¯1 for which det(M0+M1)=0, where det denotes the matrix determinant, *ω*i=2*πf* is real angular frequency and *σ* is the growth rate. The stability of the mode is defined by the sign of the growth rate, with a positive value corresponding to an unstable mode. This calculation is performed within *OSCILOS*. The combustor parameters considered for the analysis are listed in table 1.
View this table:
[Table 1.](http://rspa.royalsocietypublishing.org/content/472/2191/20150821/T1)
Table 1.
Parameters used in the analysis. They are fixed, unless otherwise stated.
It is found that only the mode with the eigenfrequency *f*≈210 Hz has the potential for instability (this corresponds to the cold flow fundamental frequency of 173 Hz for perfect acoustic reflections from the boundaries). This is consistent with both the FDF gain fall-off and the acoustic losses increasing with frequency, having a stabilizing effect on higher-order modes. A growth rate map for this unstable mode as a function of both flame position within the duct and flame velocity perturbation level is shown in figure 6, with only positive (unstable) growth rates shown. The most dangerous unstable position is a quarter of the way along the duct, which is generally the situation in a Rijke tube [53]. For each flame location, limit cycle oscillations will be established at the velocity fluctuation level corresponding to the uppermost contour, as this is where the growth rate has dropped to zero. The flame location, *x*f=0.2 m, is chosen for feedback controller design, corresponding to the limit cycle being established at a relatively large u^1/u¯1. Figure 7 shows the trajectories of the growth rate and eigenfrequency with increasing u^1/u¯1 for this flame location. For weak perturbations, the growth rate is strongly positive, implying a rapidly growing exponential envelope, e.g. for u^1/u¯1=0.025, the growth rate is 46 s−1, implying exponentially growth with an envelope of exp(46t). This growth rate decreases with u^1/u¯1, becoming zero when u^1/u¯1=0.18, hence this is the velocity fluctuation amplitude at which we expect limit cycle oscillations to occur.
[Figure 6.](http://rspa.royalsocietypublishing.org/content/472/2191/20150821/F6)
Figure 6.
Contour plot of the positive growth rate *σ* of the unstable mode with *x*f and u^1/u¯1. (Online version in colour.)
[Figure 7.](http://rspa.royalsocietypublishing.org/content/472/2191/20150821/F7)
Figure 7.
Trajectories of eigenvalues with increasing u^1/u¯1. *x*f=0.2 m. The mode with frequency ∼100 Hz is introduced by the time delay of the FDF and is always stable. (Online version in colour.)
## 5. Feedback control from within limit cycle oscillations
### (a) Open-loop transfer function
Actuators for feedback control of combustion instabilities generally fall within two categories. The first seeks to modify the fuel–air composition, normally via a fuel valve, whereas the second uses acoustic actuators to modify the pressure waves inside the combustor directly [1,18]. This work assumes the latter—these are generally better suited to light duty combustors. A loudspeaker at the bottom end of the duct is assumed, injecting a pressure signal *p*L. A microphone mounted *x**r* = 0.5 m downstream of the flame measures the pressure at that location. To simplify the analysis, we do not account for the instrument transfer functions of the loudspeaker and microphone [13,54]: we consider that the electrical driving signal to the loudspeaker *I*in equals *p*L, and the microphone electrical output signal *I*out equals the measured pressure signal *p**r*. Furthermore, we assume that the fuel–air ratio is unaffected by the acoustic waves, maintaining the applicability of the FDF calculated in §3. For model-based feedback controller design (figure 8), it is necessary to have a model for the OLTF from *p*L to *p**r*: P(s)=pr(s)pL(s)=B(M0+M1)−1C,5.1where B=[0,0,e−τr+s2,eτr−s2]andC=[1,0,0,0]T5.2and the time delays corresponding to the microphone position are τr±=(xr−xf)/(c¯2±u¯2).
[Figure 8.](http://rspa.royalsocietypublishing.org/content/472/2191/20150821/F8)
Figure 8.
The sketch of the negative feedback control system.
Owing to the flame nonlinearity, P(s) is nonlinear and depends on u^1/u¯1. Following weakly nonlinear theory, it can be cast into a set of linear OLTFs, P(s)∈P(s), with a different *P*(*s*) for each value of u^1/u¯1. The OLFTs for different u^1/u¯1 when *x*f=0.2 m are shown in figure 9.
[Figure 9.](http://rspa.royalsocietypublishing.org/content/472/2191/20150821/F9)
Figure 9.
Evolutions of the gain (*a*) and phase (*b*) of the OLTF *P* with *f* for different u^1/u¯1. (Online version in colour.)
Each modal peak is associated with a phase change of *π*, implying a second-order transfer function contribution. The direction of the phase change provides stability information—a phase increase in *π* corresponds to an unstable conjugate pair of poles, whereas a phase decrease in *π* indicates a stable conjugate pair. From figure 9, normalized velocity perturbation amplitudes, u^1/u¯1, of 0–0.15 correspond to unstable OLTFs, whereas amplitudes of 0.2 and 0.25 correspond to stable OLTFs. This is consistent with the previous finding that limit cycle oscillations occur at u^1/u¯1=0.18.
It is thus clear that a controller design based on just one of these linear OLTFs cannot guarantee the entire set of plants. In order to guarantee stability from within limit cycle oscillations, a robust controller must
* — reduce the size of perturbations when the system is in limit cycle, with the effect of reducing u^1/u¯1. The set of plants in this regime appear open loop stable and
* — result in closed loop stability for all plants. This requires a sufficiently large stability margin.
### (b) H∞ loop-shaping and the *ν*-gap metric
H∞ loop-shaping is used to design a robust stabilizing controller. H∞ loop-shaping combines classical loop-shaping design with guaranteed robustness at all frequencies. It involves designing a pre- or post-compensator to ‘shape’ the singular values (i.e. the gain for single-input–single-output (SISO) systems) of the open-loop system, so that they have a high gain where disturbance rejection is important, and a low gain where robustness to multiplicative plant uncertainty is required (see references [17,55] for an explanation of different types of uncertainty from a control systems perspective). A stabilizing controller is then synthesized using Matlab (e.g. using the *ncfsyn* command), which minimizes the H∞ norm of the closed-loop transfer function matrix (see equation (5.5))—this norm also provides a measure of the system robustness. If it is suitably small, then a robust stabilizing controller has been found. The gain of the shaped transfer function will then not be significantly affected by the controller, but the phase will be altered so as to achieve stability [17]. If the norm is not suitably small then another iteration on the choice of the pre- or post-compensator must be performed. The underpinning mathematics can be found in [17].
The shape of the open-loop plant *P* is first modified using pre- and post-compensators *W*1 and *W*2 (figure 8). For an SISO system, modifying either one of these has the same effect. The H∞ controller K∞ is then synthesized in Matlab, with the final feedback controller *K*, being K=W1K∞W2.
The *ν*-gap metric provides a useful measure of the distance in the state space between two open-loop plants with respect to how they behave when connected in a unity feedback control loop [16]. For SISO systems, it can be computed from the frequency response: the *ν*-gap, *δ**ν*, between plants *P*1 and *P*2 is given by δν(P1,P2)=supωψ(P1(iω),P2(iω)),provided the winding number condition is satisfied4,5.3where4 ψ(P1(iω),P2(iω))=|P1(iω)−P2(iω)|(1+|P1(iω)|2)1/2(1+|P2(iω)|2)1/2.5.4
The *ν*-gap metric fits very naturally into H∞ loop-shaping control design, which returns both a feedback controller design, and the corresponding ‘stability margin’, b∞, for that controller/plant pairing. The stability margin b∞ for a shaped plant *P**W*=*W*1*P*1*W*2 and the resulted controller K∞ is represented as [17] b∞=∥[IPW](I−K∞PW)−1[IK∞]∥∞−1,5.5where ∥⋅∥∞ is the infinity norm. If the controller is applied to a different plant (e.g. *P*2), separated from the initial by *ν*-gap *δ**ν*, as long as b∞>δν, then the H∞ loop-shaping controller is also guaranteed to stabilize this second system.
### (c) H∞ loop-shaping controller design
H∞ loop-shaping controller design will need to be designed based on a nominal linear plant. The differences between this nominal plant and all of the possible OLTFs will then be characterized using the *ν*-gap metric. The controller design can be schematically explained by the flowchart shown in figure 10. For combustion instabilities, nonlinearities are introduced by the FDF, which can be treated as a discrete set of FTFs for different velocity perturbation levels u^1/u¯1. From figure 6, limit cycle oscillations are generally established for u^1/u¯1≤0.2. We thus calculate the mean FTF, F0, by averaging over the FTFs for velocity perturbation levels u^1/u¯1≤0.25. The mean FTF is fitted to a state space form in order to construct the mean OLTF (indicated by *P*) for feedback controller design.
[Figure 10.](http://rspa.royalsocietypublishing.org/content/472/2191/20150821/F10)
Figure 10.
Flowchart of the H∞ loop-shaping controller design. Note that the post-compensator *W*2=*I*. (Online version in colour.)
The presence of time delays in the wave-based models (e.g. see equations (4.1), and (5.2)) means the system is infinite dimensional. In order to implement H∞ synthesis, the mean OLTF needs to be approximated by a RTF which will be finite dimensional. The exponential terms could be replaced by Padé approximants [56], but the presence of complicated damping mechanisms means that this is not straightforward. In this work, a frequency-domain fitting procedure using the Matlab command *fitfrd* is applied to obtain a RTF approximation to the mean OLTF. As the FDF is only valid over the [10, 500] Hz frequency range, we fit the OLTF over a restricted frequency range [0, 1000] Hz, and also add a weighted low-pass filter to the RTF to reduce its gain at high frequencies. The resulting RTF, represented by *P**F*, typically has an order of 16. A comparison of the original infinite-dimensional and RTF-fitted plants is shown in figure 11*a*.
[Figure 11.](http://rspa.royalsocietypublishing.org/content/472/2191/20150821/F11)
Figure 11.
Comparisons between the original plant and fitted plant. (*a*) High-order fitting of *P*. (*b*) Low-order fitting of *P**W*=*P**W*1. W1=WL2 (see equation (5.7)). (Online version in colour.)
As *P**F* is SISO, we let the post-compensator *W*2=*I* and design the pre-compensator *W*1. A first-order low-pass filter with a cut-off frequency 300 Hz is chosen as *W*1, to drop the gain at high frequencies to provide robustness to multiplicative plant uncertainty. The *ν*-gaps between the plant *P**F* and the plants *P* for different u^1/u¯1 can be identified from the frequency variations of *ψ* (see equation (5.4)), shown in figure 12*a*. Large bumps occur at each mode, with the deviation reaching a maximum at the unstable frequency—210 Hz. The maximum value of the *ν*-gap across all different u^1/u¯1 levels is max(δν)=0.18.
[Figure 12.](http://rspa.royalsocietypublishing.org/content/472/2191/20150821/F12)
Figure 12.
Plotof *ψ* with *f* for different u^1/u¯1. (*a*) High-order fitting. *ψ*(*P**F*(i*ω*),*P*(i*ω*)). (*b*) Low-order fitting. *ψ*(*P**WF*(i*ω*),*P**W*(i*ω*)). The dashed line represents the stability margin of the controller: b∞. (Online version in colour.)
A H∞ loop-shaping synthesis is then performed, using the Matlab command *ncfsyn* in the Robust control toolbox. The designed controller has an order of 16, and its stability margin is b∞=0.41, which is significantly larger than max(δν) above. This suggests that the controller should guarantee closed loop stability of all possible plants.
The performance of the negative feedback controller can be assessed by calculating the poles of the closed-loop transfer function P(s)/(1+K(s)P(s)), where P(s) is the original OLTF (equation (5.1)). Figure 13*a* shows the poles locations for different u^1/u¯1 in the presence of the feedback controller. The growth rates of all modes are now stable, confirming closed loop stability for all u^1/u¯1.
[Figure 13.](http://rspa.royalsocietypublishing.org/content/472/2191/20150821/F13)
Figure 13.
Trajectories of eigenvalues with increasing u^1/u¯1 for the feedback control system. (*a*) High-order controller deign. (*b*) Low-order controller design. (Online version in colour.)
Of course, for sufficiently large values of u^1/u¯1, the combustor exhibits limit cycle oscillations and the linear OLTF appears stable even without feedback control. For these conditions, the sensitivity transfer function provides an essential measure of the controller performance. S(s)=11+K(s)P(s).5.6When |*S*(i2*πf*)|<1, the feedback system attenuates disturbances compared with the open loop system, which corresponds to the amplitude of limit cycle oscillations and hence u^1/u¯1 being reduced. Thus, assuming that the limit cycle frequency matches or is close to the frequency of the unstable mode, *f*⋆, the requirement now becomes |*S*(i2*πf*⋆)|<15 . Figure 14*a* shows the plots of |*S*| versus frequency for different u^1/u¯1. The unstable frequency is from 205 to 215 Hz in the amplitude range 0≤u^1/u¯1≤0.2. The fact that |*S*(i2*πf*⋆)|<1 at these frequencies means that the current controller will reduce the limit cycle oscillation amplitude. For sufficient attenuation, the plant will appear unstable (see the yellow and orange colour lines in figure 9) and the feedback controller will then stabilize the system.
[Figure 14.](http://rspa.royalsocietypublishing.org/content/472/2191/20150821/F14)
Figure 14.
Evolution of |S| with frequency for different u^1/u¯1. (*a*) High-order controller design. (*b*) Low-order controller design. (Online version in colour.)
### (d) Low-order H∞ loop-shaping controller design
It has been shown that H∞ loop-shaping combined with the *ν*-gap metric can be used to design a robust controller that guarantees stabilization from within the limit cycle. The designed controller typically has an order comparable to that of the shaped plant [55]. High-order controllers are difficult to implement and exhibit lower reliability than low-order controllers. Although controller order reduction methods, such as balanced truncation methods [57], can be used on the plant or the controller, this work proposes an alternative approach to directly reduce the order of plant, provided that the maximum *ν*-gap (max(δν)) is smaller than the stability margin b∞.
The low-order method swaps the order of the shaping and fitting procedures. The plant *P* is firstly shaped with *W*1 using an ‘aggressive’ low-pass filter, which does not significantly alter the gain at low frequencies near the unstable mode, but reduces the gain at frequencies above the unstable mode, to provide robustness to fitting errors at these frequencies. A low-order fitting procedure is then implemented on the shaped plant *P**W*, with the fit over low frequencies near the unstable mode prioritized. A fourth-order low-pass filter is chosen as the pre-compensator *W*1 to shape the original plant, given by W1=WL2, where WL(s)=ωL2s2+2ξLωLs+ωL2,5.7where *ω*L=2*π* × 200 s−1 and ξL=2/2. *W*1 has a gain of 0.5 at 200 Hz, and so does not significantly alter the gain at the unstable frequency. The gain rapidly falls off to 0.025 at 500 Hz.
The controller designed by H∞ loop-shaping synthesis generally has a dimension comparable to that of the shaped plant, and so a fitted plant with low dimension is desirable. The shaped system contains only one unstable mode and exhibits fast gain fall-off at higher frequencies. The fitting procedure is performed on the shaped plant *P**W*=*P**W*1—this allows the dimension of the fitted RTF to be significantly reduced and yields a RTF with an order of 4. Figure 11*b* compares the original, *P**W*, and fitted, *P**WF*, plants. The fitted RTF captures the shape of original plant at low frequencies and has the shape of low-pass filter at high frequencies. Figure 12*b* shows the evolution of *ψ*(*P**WF*(i*ω*),*P**W*(i*ω*)) with frequency. Compared with the previous approach (figure 12*a*), there are larger deviations at high frequencies ( *f*≥250 Hz) where the low-order fitting is less accurate. However, owing to the shape of *W*1, relatively small values of *ψ*(*P**WF*(i*ω*),*P**W*(i*ω*)) are guaranteed for all frequencies and the maximum of *ν*-gap for all u^1/u¯1 is max(δν)=0.34.
A controller with an order of 3 is then obtained from the H∞ loop-shaping synthesis. The final controller K=W1K∞ has an order of 7 in the denominator and 3 in the numerator of its RTF, which is straightforward to implement. The stability margin is b∞=0.52, which is larger than max(δν).
The locations of the closed-loop poles for different u^1/u¯1 are shown in figure 13*b*. The unstable mode is stabilized for all u^1/u¯1, and the modes at high frequencies are not disturbed. However, the growth rate reduction for the unstable mode is less than for the high-order controller (figure 13*a*). This is due to the larger *ν*-gap caused by the low-order fitting. The gain of the sensitivity function *S*(*s*) with frequency is shown in figure 14*b* and is seen to be smaller than unity near the unstable mode frequency. This, combined with the fact that b∞>max(δν), means that the low-order controller is also guaranteed to be stabilizing from within limit cycle oscillations.
The robustness of the low-order controller design to changes in operating condition is investigated by moving the flame position along the combustor. Figure 15 shows that the growth rate of the dominant mode becomes negative for all unstable configurations. Thus, the controller designed stabilizes all possible flame locations within this combustor.
[Figure 15.](http://rspa.royalsocietypublishing.org/content/472/2191/20150821/F15)
Figure 15.
Plotof the growth rate of the unstable mode (*f*≈210 Hz) with and without control for different flame positions *x*f and different velocity perturbation levels u^1/u¯1. The dot markers represent the average growth rates. (Online version in colour.)
The controller implementation in time domain simulations for the design condition of *x*f=0.2 m is shown in figure 16 (more information on the time domain simulation method is found in [31,34,58]). The normalized velocity perturbations before the flame, u1′/u¯1, and the actuation signal from the loudspeaker *p*L are shown both before feedback control is activated and after it has been switched on. Prior to control, small disturbances are seen to grow rapidly, with a limit cycle established approximately beyond *t*=0.4 s. When the controller is switched on at *t*=0.5 s, the velocity perturbations decrease rapidly to the weak background white noise. The loudspeaker signal, *p*L, is initially large in order to attain control and decreases to a much smaller value in order to maintain control.
[Figure 16.](http://rspa.royalsocietypublishing.org/content/472/2191/20150821/F16)
Figure 16.
Evolution of *p*L and u1′/u¯1 with time. The controller is switched on when *t*=0.5 s. *x*f=0.2 m.
## 6. Conclusion
This article has presented a method of designing linear feedback controllers that are guaranteed to stabilize combustion instability, whenever they are activated. The guarantees are subject to some very reasonable assumptions on the nature of the flame nonlinearity that would be satisfied very widely in practice. A single feedback controller design then can guarantee stability whether control is activated within the regime of small linear disturbances and exponential growth or within the regime of limit cycle oscillations in which the amplitude has saturated.
The method combines the FDF approach for weakly nonlinear flame responses with H∞ loop-shaping for controller design and the *ν*-gap metric for characterizing nonlinearity as deviations from the ‘average’ linear model. The FDF gives rise to a ‘family’ of linear OLTFs corresponding to different flame perturbation levels. The maximum *ν*-gap between these and the average or nominal plant transfer function sets the minimum robustness stability margin required of the controller. An H∞ loop-shaping controller is then designed which has a stability margin exceeding this, and which also attenuates fluctuations within limit cycle.
The control strategy has been demonstrated on a representative model of a nonlinear combustion system comprising an unsteady laminar flame inside a Rijke tube. Linear plane waves within the duct are assumed to suffer damping owing to end radiation and viscothermal wall damping mechanisms. The nonlinear kinematic response of the flame to flow disturbances is modelled using an LSA, which combines the *G*-equation kinematic flame model with a two-dimensional flow velocity perturbation model and a flame front curvature-dependent burning speed. Simulations implementing this LSA have been validated by comparison with experimental data, and used to determine the FDF of the laminar conical premixed flame.
The calculated FDF of the target flame has been successfully incorporated into a recently developed code, *OSCILOS*, to calculate the eigenvalue evolutions of the nonlinear system with increasing flame velocity perturbation levels. The corresponding ‘family’ of OLTF between the monitored pressure perturbations downstream of the flame and the external pressure perturbations from an actuator at the entrance of the Rijke tube, has been established for the purpose of control.
A linear mean OLTF has been selected as the target plant for the controller design by averaging the FDF over a large range of velocity perturbation levels. The *ν*-gap metric has been applied to quantify the differences between the set of OLTFs for different perturbation levels and the selected plant. This metric fits naturally into the H∞ loop-shaping framework to provide the minimum required stability margin of the designed robust controller. An initial H∞ loop-shaping controller has been designed and shown to give effective robust control performance. Subsequently, a low-order H∞ loop-shaping controller based on an aggressively low-pass filtered weighted plant has been designed. Both controllers are guaranteed to stabilize the combustor for all possible nonlinear flame responses, as long as the *ν*-gap is smaller than the controller stability margin of the controller designed by H∞ loop-shaping synthesis. The resulting controller has been applied to simulations of the unsteady combustion system in the time domain and has been indeed found to suppress the combustion instability from within the limit cycle.
## Data accessibility
This paper has no data.
## Author' contributions
J.L. and A.S.M. conceived the mathematical models, interpreted the computational results, and wrote the paper. J.L. performed most of the simulations in consultation with A.M. All authors gave final approval for publication.
## Competing interests
We declare we have no competing interests.
## Funding
This work is supported by the European Research Council (ERC) Starting Grant ACOULOMODE (2013–2018).
## Acknowledgements
The authors thank the anonymous reviewers for their valuable comments and suggestions to improve the quality of the paper.
## Appendix A
**(a) Viscothermal damping of acoustic waves in a duct**
The propagation of sound waves in a duct has been found to be dissipated owing to viscothermal losses at the duct walls [30,59]. The forward and backward wavenumbers, k¯±=−is/c¯(1±M¯), are modified to k±=k¯±+Δk±, where the correction approximately equals Δk±=−is1/2c¯(1±M¯)ΛA 1with the coefficient *Λ* expressed as Λ=2ν1/2D(1+γ−1Pr1/2)A 2where *ν* and *P**r* are the kinematic viscosity and Prandtl number of gases, respectively. *D* is the diameter of the duct. *γ* is the ratio of specific heats. In the calculation, *P**r*=0.71 is considered constant, both for air and combustion products of hydrocarbon flames. *ν* can be calculated by ν=μ/ρ¯, with the dynamic viscosity *μ*air of air calculated by the Sutherland’s Law [60], with further corrections expressed as μair={μref(TTref)3/2Tref+CST+CSfor T∈[100,1000] Ka1T+a0for T∈[1000,3000] K,where *μ*ref=1.716×10−5 kg m−1 s−1, *T*ref=273.15 K, *C**S*=110.4 K, *a*1=2.653×10−8 and *a*=1.573×10−5. It was found that the dynamic viscosity of the hydrocarbon–air combustion products *μ*prod differed little from that of air, with a minor correction depending on the equivalence ratio *ϕ* [61], and is written as μprod=μair1+0.027ϕ,for T∈[500,4000] K, ϕ∈[0,4].A 3
The exponential component in propagating waves can be changed to exp(±ik±Δx)=exp(∓τ±(s+Λs1/2)),A 4where τ±=Δx/(c¯±u¯) and Δ*x* is the path length of the travelling acoustic wave.
**(b) The expressions of Mn**
M1=1+S1S2(M¯1(2+M¯1)−c¯2c¯1M¯2(1+M¯1))M3=1+M¯2M2=1−S1S2(M¯1(2−M¯1)−c¯2c¯1M¯2(1−M¯1))M4=1−M¯2M5=1+γM¯1+(γ−1)M¯12M7=1+γM¯2+(γ−1)M¯22M6=1−γM¯1+(γ−1)M¯12M8=1−γM¯2+(γ−1)M¯22,where M=u/c¯ is the Mach number.
## Footnotes
* 1 The pressure reflection coefficient is defined as the ratio of reflected to incident acoustic waves at the boundary [24].
* 2 For small enough flow disturbances, the flame’s unsteady heat release rate generally responses linearly and can be described by a flame transfer function F(s)=Q˙^/Q˙¯/(u^1/u¯1)[40].
* 3 *fitfrd* returns a RTF of the best fit of a set of frequency-response data.
* 4 The winding number condition is: wno(1+P∼2(iω)P1(iω))+η(P1(iω))−η(P2(iω))=0, where wno is the winding number, P∼ is the complex conjugate transpose of *P* and *η*(*P*) represents the number of open right-hand-plane poles of *P*[16].
* 5 It should be noted that the condition |*S*(i2*πf*⋆)|<1 is the sensitivity at the specific frequency, *f*⋆, that is relevant to reducing amplitude. It cannot be used to determine whether the plant over all frequencies is stabilized by the controller.
* Received November 30, 2015.
* Accepted June 10, 2016.
* © 2016 The Authors.
[http://creativecommons.org/licenses/by/4.0/](http://creativecommons.org/licenses/by/4.0/)
Published by the Royal Society under the terms of the Creative Commons Attribution License [http://creativecommons.org/licenses/by/4.0/](http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, provided the original author and source are credited.
## References
1. Candel S. 2002 Combustion dynamics and control: progress and challenges. Proc. Combust. Inst. 29, 1–28. ([doi:10.1016/S1540-7489(02)80007-4](http://dx.doi.org/10.1016/S1540-7489(02)80007-4))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1016/S1540-7489(02)80007-4&link_type=DOI)
2. 1. Lieuwen T,
2. V. Yang
(eds). 2005 Combustion instabilities in gas turbines, operational experience, fundamental mechanisms, and modeling. Progress in Astronautics and Aeronautics, vol. 210. Reston,VA: AIAA, Inc.
3. Chu BT. 1964 On the energy transfer to small disturbances in fluid flow (part I). Acta Mech. 1, 215–234. ([doi:10.1007/BF01387235](http://dx.doi.org/10.1007/BF01387235))
4. Dowling AP, Morgans AS. 2005 Feedback control of combustion oscillations. Annu. Rev. Fluid Mech. 37, 151–182. ([doi:10.1146/annurev.fluid.36.050802.122038](http://dx.doi.org/10.1146/annurev.fluid.36.050802.122038))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1146/annurev.fluid.36.050802.122038&link_type=DOI)
5. Balachandran R, Ayoola BO, Kaminski C, Dowling AP, Mastorakos E. 2005 Experimental investigation of the nonlinear response of turbulent premixed flames to imposed inlet velocity oscillations. Combust. Flame 143, 37–55. ([doi:10.1016/j.combustflame.2005.04.009](http://dx.doi.org/10.1016/j.combustflame.2005.04.009))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1016/j.combustflame.2005.04.009&link_type=DOI)
6. Noiray N, Durox D, Schuller T, Candel S. 2008 A unified framework for nonlinear combustion instability analysis based on the flame describing function. J. Fluid Mech. 615, 139–167. ([doi:10.1017/S0022112008003613](http://dx.doi.org/10.1017/S0022112008003613))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1017/S0022112008003613&link_type=DOI)
7. Schimek S, Moeck JP, Paschereit CO. 2011 An experimental investigation of the nonlinear response of an atmospheric swirl-stabilized premixed flame. J. Eng. Gas Turbines Power 133, 101502. ([doi:10.1115/1.4002946](http://dx.doi.org/10.1115/1.4002946))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1115/1.4002946&link_type=DOI)
8. 1. Harrje DT,
2. FH. Reardon
(eds). 1972 **Liquid propellant rocket combustion instability**. Washington, DC: NASA.
9. Poinsot T, Veynante D, Bourienne F, Candel S, Esposito E, Surget J. 1989 Initiation and suppression of combustion instabilities by active control. Proc. Combust. Inst. 22, 1363–1370. ([doi:10.1016/S0082-0784(89)80147-X](http://dx.doi.org/10.1016/S0082-0784(89)80147-X))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1016/S0082-0784(89)80147-X&link_type=DOI)
10. Kabiraj L, Sujith RI. 2012 Nonlinear self-excited thermoacoustic oscillations: intermittency and flame blowout. J. Fluid Mech. 713, 376–397. ([doi:10.1017/jfm.2012.463](http://dx.doi.org/10.1017/jfm.2012.463))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1017/jfm.2012.463&link_type=DOI)
11. Boudy F, Durox D, Schuller T, Candel S. 2013 Analysis of limit cycles sustained by two modes in the flame describing function framework. CR Mech. 341, 181–190. ([doi:10.1016/j.crme.2012.10.014](http://dx.doi.org/10.1016/j.crme.2012.10.014))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1016/j.crme.2012.10.014&link_type=DOI)
12. Waugh IC, Kashinath K, Juniper MP. 2014 Matrix-free continuation of limit cycles and their bifurcations for a ducted premixed flame. J. Fluid Mech. 759, 1–27. ([doi:10.1017/jfm.2014.549](http://dx.doi.org/10.1017/jfm.2014.549))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1017/jfm.2014.549&link_type=DOI)
13. Morgans AS, Dowling AP. 2007 Model-based control of combustion instabilities. J. Sound Vib. 299, 261–282. ([doi:10.1016/j.jsv.2006.07.014](http://dx.doi.org/10.1016/j.jsv.2006.07.014))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1016/j.jsv.2006.07.014&link_type=DOI)
14. Morgans AS, Stow SR. 2007 Model-based control of combustion instabilities in annular combustors. Combust. Flame 150, 380–399. ([doi:10.1016/j.combustflame.2007.06.002](http://dx.doi.org/10.1016/j.combustflame.2007.06.002))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1016/j.combustflame.2007.06.002&link_type=DOI)
[Web of Science](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=000248989800008&link_type=ISI)
15. Illingworth SJ, Morgans AS. 2010 Adaptive feedback control of combustion instability in annular combustors. Combust. Sci. Technol. 182, 143–164. ([doi:10.1080/00102200903258773](http://dx.doi.org/10.1080/00102200903258773))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1080/00102200903258773&link_type=DOI)
16. Vinnicombe G. 2001 **Uncertainty and feedback: H∞ loop-shaping and the *ν*-gap metric**. London, UK: Imperial College Press.
17. McFarlane D, Glover K. 1992 A loop shaping design procedure using H∞ synthesis. IEEE Trans. Autom. Control 37, 759–769. ([doi:10.1109/9.256330](http://dx.doi.org/10.1109/9.256330))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1109/9.256330&link_type=DOI)
[Web of Science](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=A1992HW45900006&link_type=ISI)
18. Campos-Delgado DU, Schuermans BBH, Zhou K, Paschereit CO, Gallestey EA, Poncet A. 2003 Thermoacoustic instabilities: modeling and control. IEEE Trans. Control Syst. Technol. 11, 1–15. ([doi:10.1109/TCST.2002.808096](http://dx.doi.org/10.1109/TCST.2002.808096))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1109/TCST.2002.808096&link_type=DOI)
19. Chu Y-C, Glover K, Dowling AP. 2003 Control of combustion oscillations via H∞ loop-shaping, *μ*-analysis and integral quadratic constraints. Automatica 39, 219–231. ([doi:10.1016/S0005-1098(02)00226-1](http://dx.doi.org/10.1016/S0005-1098(02)00226-1))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1016/S0005-1098(02)00226-1&link_type=DOI)
20. Yuan X-C, Glover K. 2013 Model-based control of thermoacoustic instabilities in partially premixed lean combustion – a design case study. Int. J. Control 86, 2052–2066. ([doi:10.1080/00207179.2013.830337](http://dx.doi.org/10.1080/00207179.2013.830337))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1080/00207179.2013.830337&link_type=DOI)
21. Lauga E, Bewley T. 2004 Performance of a linear robust control strategy on a nonlinear model of spatially developing flows. J. Fluid Mech. 512, 343–374. ([doi:10.1017/s0022112004009711](http://dx.doi.org/10.1017/s0022112004009711))
22. Weiss G, Zhong Q-C, Green T, Liang J. 2004 H∞ repetitive control of DC-AC converters in microgrids. IEEE Trans. Power Electron. 19, 219–230. ([doi:10.1109/TPEL.2003.820561](http://dx.doi.org/10.1109/TPEL.2003.820561))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1109/TPEL.2003.820561&link_type=DOI)
23. Lefebvre D, Chevrel P, Richard S. 2003 An H∞ based control design methodology dedicated to the active control of vehicle longitudinal oscillations. IEEE Trans. Control Syst. Technol. 11, 948–956. ([doi:10.1109/TCST.2003.815552](http://dx.doi.org/10.1109/TCST.2003.815552))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1109/TCST.2003.815552&link_type=DOI)
24. Poinsot T, Veynante D. 2005 **Theoretical and numerical combustion**. Philadelphia, PA: R.T. Edwards Inc.
25. Marble FE, Candel SM. 1977 Acoustic disturbance from gas non-uniformities convected through a nozzle. J. Sound Vib. 55, 225–243. ([doi:10.1016/0022-460X(77)90596-X](http://dx.doi.org/10.1016/0022-460X(77)90596-X))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1016/0022-460X(77)90596-X&link_type=DOI)
[Web of Science](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=A1977ED87100005&link_type=ISI)
26. Lieuwen T. 2003 Modeling premixed combustion-acoustic wave interactions: a review. J. Propul. Power 19, 765–781. ([doi:10.2514/2.6193](http://dx.doi.org/10.2514/2.6193))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.2514/2.6193&link_type=DOI)
27. Bourehla A, Baillot F. 1998 Appearance and stability of a laminar conical premixed flame subjected to an acoustic perturbation. Combust. Flame 114, 303–318. ([doi:10.1016/S0010-2180(97)00323-4](http://dx.doi.org/10.1016/S0010-2180(97)00323-4))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1016/S0010-2180(97)00323-4&link_type=DOI)
28. Dowling AP. 1995 The calculation of thermoacoustic oscillations. J. Sound Vib. 180, 557–581. ([doi:10.1006/jsvi.1995.0100](http://dx.doi.org/10.1006/jsvi.1995.0100))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1006/jsvi.1995.0100&link_type=DOI)
29. Levine H, Schwinger J. 1948 On the radiation of sound from an unflanged circular pipe. Phys. Rev. 73, 383–406. ([doi:10.1103/PhysRev.73.383](http://dx.doi.org/10.1103/PhysRev.73.383))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1103/PhysRev.73.383&link_type=DOI)
[Web of Science](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=A1948UB29500017&link_type=ISI)
30. Peters M, Hirschberg A, Reijnen A, Wijnands A. 1993 Damping and reflection coefficient measurements for an open pipe at low Mach and low Helmholtz numbers. J. Fluid Mech. 256, 499–534. ([doi:10.1017/S0022112093002861](http://dx.doi.org/10.1017/S0022112093002861))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1017/S0022112093002861&link_type=DOI)
31. Li J, Morgans AS. 2015 Time domain simulations of nonlinear thermoacoustic behaviour in a simple combustor using a wave-based approach. J. Sound Vib. 346, 345–360. ([doi:10.1016/j.jsv.2015.01.032](http://dx.doi.org/10.1016/j.jsv.2015.01.032))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1016/j.jsv.2015.01.032&link_type=DOI)
32. Dowling AP. 1997 Nonlinear self-excited oscillations of a ducted flame. J. Fluid Mech. 346, 271–290. ([doi:10.1017/S0022112097006484](http://dx.doi.org/10.1017/S0022112097006484))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1017/S0022112097006484&link_type=DOI)
33. Morgans AS, Annaswamy AM. 2008 Adaptive control of combustion instabilities for combustion systems with right-half plane zeros. Combust. Sci. Technol. 180, 1549–1571. ([doi:10.1080/00102200802125719](http://dx.doi.org/10.1080/00102200802125719))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1080/00102200802125719&link_type=DOI)
34. Li J, Yang D, Luzzato C, Morgans AS. 2014 OSCILOS: the open source combustion instability low order simulator. Technical Report, Imperial College London. See [http://www.oscilos.com](http://www.oscilos.com).
35. Han X, Li J, Morgans AS. 2015 Prediction of combustion instability limit cycle oscillations by combining flame describing function simulations with a thermoacoustic network model. Combust. Flame 162, 3632–3647. ([doi:10.1016/j.combustflame.2015.06.020](http://dx.doi.org/10.1016/j.combustflame.2015.06.020))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1016/j.combustflame.2015.06.020&link_type=DOI)
36. Han X, Morgans AS. 2015 Simulation of the flame describing function of a turbulent premixed flame using an open-source LES solver. Combust. Flame 162, 1778–1792. ([doi:10.1016/j.combustflame.2014.11.039](http://dx.doi.org/10.1016/j.combustflame.2014.11.039))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1016/j.combustflame.2014.11.039&link_type=DOI)
37. Markstein G. 1964 *Nonsteady flame propagation*. Oxford, UK: Pergamon Press.
38. Dowling AP. 1999 A kinematic model of a ducted flame. J. Fluid Mech. 394, 51–72. ([doi:10.1017/S0022112099005686](http://dx.doi.org/10.1017/S0022112099005686))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1017/S0022112099005686&link_type=DOI)
39. Luzzato CM, Assier RC, Morgans AS, Wu X. 2013 Modelling thermo-acoustic instabilities of an anchored laminar flame in a simple lean premixed combustor: including hydrodynamic effects. In *Proc. of the 19th AIAA/CEAS Aeroacoustics Conf., Berlin, Germany, 27–29 May*, pp. 1–15.
40. Schuller T, Ducruix S, Durox D, Candel S. 2002 Modeling tools for the prediction of premixed flame transfer functions. Proc. Combust. Inst. 29, 107–113. ([doi:10.1016/S1540-7489(02)80018-9](http://dx.doi.org/10.1016/S1540-7489(02)80018-9))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1016/S1540-7489(02)80018-9&link_type=DOI)
41. Schuller T. 2003 Mécanismes de Couplage dans les Interactions Acoustique-Combustion. PhD thesis, École Centrale Paris.
42. Kashinath K, Hemchandra S, Juniper MP. 2013 Nonlinear thermoacoustics of ducted premixed flames: the influence of perturbation convection speed. Combust. Flame 160, 2856–2865. ([doi:10.1016/j.combustflame.2013.06.019](http://dx.doi.org/10.1016/j.combustflame.2013.06.019))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1016/j.combustflame.2013.06.019&link_type=DOI)
43. Gu X, Haq M, Lawes M, Woolley R. 2000 Laminar burning velocity and Markstein lengths of methane–air mixtures. Combust. Flame 121, 41–58. ([doi:10.1016/S0010-2180(99)00142-X](http://dx.doi.org/10.1016/S0010-2180(99)00142-X))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1016/S0010-2180(99)00142-X&link_type=DOI)
[Web of Science](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=000085534000003&link_type=ISI)
44. Osher S, Fedkiw R. 2003 **Level set methods and dynamic implicit surfaces**. New York, NY: Springer.
45. Gottlieb S, Shu C-W. 1998 Total variation diminishing Runge–Kutta schemes. Math. Comput. 67, 73–85. ([doi:10.1090/S0025-5718-98-00913-2](http://dx.doi.org/10.1090/S0025-5718-98-00913-2))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1090/S0025-5718-98-00913-2&link_type=DOI)
[Web of Science](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=000071947700004&link_type=ISI)
46. Peng D, Merriman B, Osher S, Zhao H, Kang M. 1999 A PDE-based fast local level set method. J. Comput. Phys. 155, 410–438. ([doi:10.1006/jcph.1999.6345](http://dx.doi.org/10.1006/jcph.1999.6345))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1006/jcph.1999.6345&link_type=DOI)
[Web of Science](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=000083590300009&link_type=ISI)
47. Wen X. 2009 High order numerical methods to two dimensional delta function integrals in level set methods. J. Comput. Phys. 228, 4273–4290. ([doi:10.1016/j.jcp.2009.03.004](http://dx.doi.org/10.1016/j.jcp.2009.03.004))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1016/j.jcp.2009.03.004&link_type=DOI)
48. Li J, Richecoeur F, Schuller T. 2013 Reconstruction of heat release rate disturbances based on transmission of ultrasounds: experiments and modeling for perturbed flames. Combust. Flame 160, 1779–1788. ([doi:10.1016/j.combustflame.2013.03.014](http://dx.doi.org/10.1016/j.combustflame.2013.03.014))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1016/j.combustflame.2013.03.014&link_type=DOI)
49. Li J, Durox D, Richecoeur F, Schuller T. 2015 Analysis of chemiluminescence, density and heat release rate fluctuations in acoustically perturbed laminar premixed flames. Combust. Flame 162, 3934–3945. ([doi:10.1016/j.combustflame.2015.07.031](http://dx.doi.org/10.1016/j.combustflame.2015.07.031))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1016/j.combustflame.2015.07.031&link_type=DOI)
50. Durox D, Schuller T, Noiray N, Candel S. 2009 Experimental analysis of nonlinear flame transfer functions for different flame geometries. Proc. Combust. Inst. 32, 1391–1398. ([doi:10.1016/j.proci.2008.06.204](http://dx.doi.org/10.1016/j.proci.2008.06.204))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1016/j.proci.2008.06.204&link_type=DOI)
51. Schuller T, Durox D, Candel S. 2003 A unified model for the prediction of laminar flame transfer functions: comparisons between conical and V-flame dynamics. Combust. Flame 134, 21–34. ([doi:10.1016/S0010-2180(03)00042-7](http://dx.doi.org/10.1016/S0010-2180(03)00042-7))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1016/S0010-2180(03)00042-7&link_type=DOI)
52. Preetham X, Santosh H, Lieuwen T. 2008 Dynamics of laminar premixed flames forced by harmonic velocity disturbances. J. Propul. Power 24, 1390–1402. ([doi:10.2514/1.35432](http://dx.doi.org/10.2514/1.35432))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.2514/1.35432&link_type=DOI)
53. Raun RL, Beckstead MW, Finlinson JC, Brooks KP. 1993 A review of Rijke tubes, Rijke burners and related devices. Prog. Energy Combust. Sci. 19, 313–364. ([doi:10.1016/0360-1285(93)90007-2](http://dx.doi.org/10.1016/0360-1285(93)90007-2))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1016/0360-1285(93)90007-2&link_type=DOI)
54. Evesque S, Dowling AP, Annaswamy AM. 2003 Self-tuning regulators for combustion oscillations. Proc. R. Soc. Lond. A 459, 1709–1749. ([doi:10.1098/rspa.2002.1085](http://dx.doi.org/10.1098/rspa.2002.1085))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1098/rspa.2002.1085&link_type=DOI)
55. Zhou K, Doyle JC, Glover K. 1996 **Robust and optimal control**. Upper Saddle Rive, NJ: Prentice Hall.
56. Brezinski C. 1996 Extrapolation algorithms and Padé approximations: a historical survey. Appl. Numer. Math. 20, 299–318. ([doi:10.1016/0168-9274(95)00110-7](http://dx.doi.org/10.1016/0168-9274(95)00110-7))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1016/0168-9274(95)00110-7&link_type=DOI)
57. Gugercin S, Antoulas AC. 2004 A survey of model reduction by balanced truncation and some new results. Int. J. Control 77, 748–766. ([doi:10.1080/00207170410001713448](http://dx.doi.org/10.1080/00207170410001713448))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1080/00207170410001713448&link_type=DOI)
[Web of Science](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=000222500300006&link_type=ISI)
58. Li J, Morgans AS. 2014 Model based control of nonlinear combustion instabilities, In *The 21st Int. Congress on Sound and Vibration, Beijing, China, 13–17 July*, pp. 1–8.
59. Kirchhoff G. 1868 Uber den einfluss der wärmteleitung in einem gas auf die schallbewegung. Ann. Phys. 134, 177–193. ([doi:10.1002/andp.18682100602](http://dx.doi.org/10.1002/andp.18682100602))
60. Sutherland W. 1893 The viscosity of gases and molecular force. Philos. Mag. Ser. 5 36, 507–531. ([doi:10.1080/14786449308620508](http://dx.doi.org/10.1080/14786449308620508))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1080/14786449308620508&link_type=DOI)
61. Mansouri SH, Heywood JB. 1980 Correlations for the viscosity and Prandtl number of hydrocarbon-air combustion products. Combust. Sci. Technol. 23, 251–256. ([doi:10.1080/00102208008952416](http://dx.doi.org/10.1080/00102208008952416))
[CrossRef](http://rspa.royalsocietypublishing.org/lookup/external-ref?access_num=10.1080/00102208008952416&link_type=DOI)