Feedback control of unstable thermoacoustic modes in an annular Rijke tube

Control Engineering Practice 20 (2012) 770–782

Contents lists available at SciVerse ScienceDirect

Control Engineering Practice journal homepage: www.elsevier.com/locate/conengprac

Feedback control of unstable thermoacoustic modes in an annular Rijke tube Gregor Gelbert a,n, Jonas P. Moeck b, Christian O. Paschereit b, Rudibert King a a b

Chair of Measurement and Control, Berlin Institute of Technology, 10623 Berlin, Germany Institute of Fluid Dynamics and Engineering Acoustics, Berlin Institute of Technology, 10623 Berlin, Germany

a r t i c l e i n f o

a b s t r a c t

Article history: Received 9 September 2010 Accepted 26 March 2012 Available online 17 April 2012

Simulation and experimental results from an annular Rijke tube are presented. This system is a thermoacoustic surrogate system of an annular gas turbine combustor which, despite its simplicity, possesses the basic mechanisms to feature unstable azimuthal modes. A thermoacoustic network model is set up and used to derive low-order models for modal control of the system. The derived controllers are successfully applied in simulation and experiment. With the modal controllers, all unstable acoustic modes can be eliminated individually. A simultaneous use of all controllers results in a complete stabilization of the system. & 2012 Elsevier Ltd. All rights reserved.

Keywords: Combustion instabilities Azimuthal mode Annular combustion chamber Rijke tube Network model Model-based control

1. Introduction Combustion instabilities remain a serious problem in the development of modern gas turbines for aeroengines as well as for power generation applications. Lean premixed flames, introduced by the gas turbine industry to open the way for a significant reduction of nitric oxides emission, are particularly susceptible to this unsteady, generally undesirable phenomenon (Lieuwen & Yang, 2005). If a flame is perturbed by an acoustic wave, it responds with a fluctuation in heat release rate. This mechanism is exceedingly complex and may be related to a direct kinematic interaction, complex hydrodynamic processes (such as the excitation of flow instabilities), or the generation and propagation of fuel–air mixture ratio non-uniformities, among others (Lieuwen, 2003). The fluctuation in heat release, in turn, acts as a volume source and thus generates an acoustic wave. This acoustic wave is reflected from the combustion chamber boundaries and impinges again on the flame. If certain phase relationships prevail between the individual mechanisms and the damping is small, a mutual amplification occurs, and the perturbations grow until limited by nonlinear effects. Since these oscillations arise due to the interaction of unsteady heat release and the acoustic field in a resonator, this phenomenon is also more generally referred to as thermoacoustic instability. The necessary phase relation to promote an amplification of acoustic and heat release rate fluctuations is given by Rayleigh’s (1878) criterion, which states that the perturbations of pressure and heat release rate must have an

n

Corresponding author. Tel.: þ49 30 39978 9548; fax: þ49 30 39978 8108 E-mail address: [email protected] (G. Gelbert).

0967-0661/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.conengprac.2012.03.016

in-phase component to be destabilizing. Assessing this criterion in practical systems is often not possible because, in particular, the time delay associated with the generation of heat release rate perturbations through the acoustic waves is difficult to quantify due to the various transport mechanisms involved. This problem of predicting combustion instabilities and thus mitigate them in the design phase has motivated the development of various approaches to suppress pressure and heat release rate oscillations in unstable combustors. Passive control techniques, typically based on Helmholtz or quarter-wave resonators, aim to increase the damping in the system and thus inhibit the growth of perturbations (Noiray & Schuermans, 2012). These methods were shown to be effective and can be frequently found in full-scale engines (Bellucci, Flohr, Paschereit, & Magni, 2004). The major drawback of passive techniques is their limited flexibility since they are only effective in specific regions in the operating range (which they have been tuned to). Active control methods, in contrast, can be adapted to different operating conditions and thus do not suffer from this drawback. The simplest active control technique employs an open-loop forcing of the system (McManus, Poinsot, & Candel, 1993), typically with sinusoidal signals, to break the interaction between the acoustic field and the unsteady heat release rate in the flame. This type of control was successfully applied in various configurations. However, suitable control parameters (frequency and amplitude in the case of sinusoidal forcing) have to be determined in the experiment, and the fact that no information on the system state is utilized also limits the effectiveness. In closed-loop control applications to mitigate combustion instabilities, dynamic information on the system state is typically provided by a pressure sensor. The simplest feedback strategy is then to actuate

G. Gelbert et al. / Control Engineering Practice 20 (2012) 770–782

the acoustic field or the unsteady heat release rate in phase opposition so that the growth of disturbances is inhibited. Various closed-loop schemes have been applied in model simulations and test-rig experiments (Candel, 2002), including elementary empirical controllers based on proportional or phase-shifted feedback (McManus et al., 1993), model-based methods using reducedorder models of the combustor acoustics and the flame response (Campos-Delgado et al., 2003; Gelbert, Moeck, Bothien, King, & Paschereit, 2008; Morgans & Dowling, 2007), and adaptive algorithms (Banaszuk, Ariyur, & Krstic´, 2004; Moeck, Bothien, Paschereit, Gelbert, & King, 2007; Riley, Park, Dowling, Evesque, & Annaswamy, 2004). The most common actuators employed in active control of combustion instabilities are loudspeakers and (mostly solenoid) fuel valves. While the former is superior in terms of linearity and bandwidth, unsteady fuel injection appears to be the only possible choice at full engine scale, i.e., high power density and elevated pressure. A recent review on closed-loop control of combustion instabilities is given by Dowling and Morgans (2005). The vast majority of previous work related to control of combustion instabilities is concerned with single-burner configurations, in which only one flame is located in a combustor set-up that is essentially purely longitudinal. However, this is in strong contrast to the combustion chambers found in modern stationary gas turbines and in aeroengines. These are of annular type, with a large number of burners (up to over 50) being mounted circumferentially and firing all into one ring-shaped chamber. Apart from multiple sources of heat release, the geometry of an annular combustion chamber is, from an acoustic point of view, completely different compared to that in a single-burner set-up. In the annular case, the largest dimension is the circumference of the chamber so that the acoustic modes associated with the lowest resonance frequencies are azimuthal, i.e., they represent approximately harmonic waves in the angular direction. The azimuthal order of these modes, m, corresponds to the number of periods along the circumference. Combustion instabilities observed in gas turbine engines are thus often associated with this type of acoustic mode (Lieuwen & Yang, 2005). Even if there is no direct interaction between adjacent flames, all unsteady heat sources are coupled dynamically through the circumferential acoustic field. Therefore, modeling and control of combustion instabilities in annular configurations is a significantly more complex problem compared to the single-burner case. In a single-burner set-up, for instance, one sensor–actuator combination is generally sufficient for control because the locations of the pressure nodes associated with purely longitudinal modes are fixed. The azimuthal modes in an annular chamber, in contrast, are degenerate due to the rotational symmetry of the problem. Accordingly, the eigenspace associated with a certain azimuthal mode is two-dimensional so that the circumferential oscillation pattern is not fixed. Using then only a single sensor and/or a single actuator permits the unstable mode to ‘‘evade’’ the control action. This is analogous to the case where only one damping device is used in an annular chamber (Stow & Dowling, 2003). The degenerate eigenvalue splits due to the broken symmetry, and only one of the split modes is stabilized. Despite these difficulties, active control of combustion instabilities in annular configurations has been achieved in previous work. Most studies are purely based on simulations; experimental applications have been restricted to open-loop or simple empirical closed-loop controllers. Seume et al. (1998) applied a feedback control scheme based on phase-shifted pressure feedback using modulated pilot fuel to suppress unstable azimuthal modes in a full-scale annular combustion chamber with 24 burners. They observed second- and fourth-order circumferential modes. The application of feedback control resulted

771

in a reduction of the pressure oscillation amplitude by 17 dB. Recent studies employed more sophisticated control schemes but are exclusively restricted to simulations. Schuermans (2003) and Schuermans, Bellucci, and Paschereit (2003) developed a threedimensional network modeling approach based on a modal expansion of the Greens function associated with the Helmholtz equation. This type of representation proved to be particularly useful for control studies since state-space models of the acoustic input–output properties are immediately obtained. A model for an annular combustion chamber was set up, and H1 loop-shaping was used to suppress unstable azimuthal modes in simulations. To reduce the size of the plant model and the number of inputs and outputs, a modal decomposition along the circumference was used. Similar modeling techniques are used in the present work; for the experimental application of the controller, a circumferential modal decomposition of the sensor and actuator signals is exploited as well. Based on a time domain network model, Schuermans, Paschereit, and Monkewitz (2006) performed simulations of thermoacoustic instabilities in a rotationally symmetric annular combustion chamber. They observed that, while growing, the unstable mode was of standing type (resulting from axisymmetric white noise excitation). However, the final limit cycle solution was always of rotating type (with no preference of either spinning direction). It was argued that the saturation nonlinearity in the flame destabilizes the standing-wave mode at finite amplitude and therefore promotes the traveling-wave solution. Krstic and Banaszuk (2006) developed an adaptive scheme to control instabilities in a dual-input dual-output system, which represented two acoustic modes in the model of an annular combustion chamber. Morgans and Stow (2007) and Illingworth and Morgans (2010) used a low-order acoustic network representation to develop model-based and adaptive control strategies for thermoacoustic instabilities in an annular model combustor. In the simulations, they were able to suppress decoupled azimuthal modes in a rotationally symmetric set-up but also coupled modes resulting from non-identical burners. To the authors’ knowledge, no model-based control scheme has been used so far to suppress thermoacoustically unstable azimuthal modes in an experiment. This has to be attributed to the high complexity of a full annular combustion chamber. Therefore in this work, a surrogate system is presented which is of significantly lower complexity than an annular combustion chamber but has all essential ingredients to exhibit thermoacoustic instabilities associated with circumferential acoustic modes. This surrogate system is an annular extension of the well-known Rijke tube (Heckl, 1990; Raun, Beckstead, Finlinson, & Brooks, 1993), the simplest device exhibiting thermoacoustic oscillations. In its most basic form, the Rijke tube consists only of a constantcross-section duct with a heating grid mounted inside. Due to the dependence of the convective heat transfer on the local flow velocity, acoustic waves cause a perturbation in the heat release rate, and thus thermoacoustic instability may occur. It is the simplicity of the Rijke tube that makes it attractive for fundamental investigations on thermoacoustic instabilities and their control. The surrogate system, which is described in detail in the next section, will be referred to as an annular Rijke tube, and this mock-up will be used to demonstrate suitable modeling and control methodologies.

2. Experimental set-up The annular Rijke tube has 12 tubes connected to an annular duct at their downstream end. Fig. 1 shows a schematic of the setup with the basic geometrical dimensions. The annular duct is

772

G. Gelbert et al. / Control Engineering Practice 20 (2012) 770–782

Fig. 2. Photograph of a heating grid (top view).

Fig. 1. Schematic of the annular Rijke tube. Top: top view with microphones (Mic) and speakers (Sp). The microphones 2, 4, 6, 8, 10 and 12 are mounted at the tubes under the speakers. Bottom: cut plane with the basic dimensions in mm.

605 mm in length and has a mean diameter of 720 mm with a ‘‘hub–tip ratio’’ of 0.8. The tubes have an inner diameter of 60 mm, and all parts are made of aluminum with a wall thickness of 10 mm. As in a conventional Rijke tube (Heckl, 1990; Raun et al., 1993), the sources of mean and unsteady heat release are electrically driven heating grids. In each of the 12 tubes, one heating grid is mounted (see Fig. 1). Obviously this type of heat source is much simpler than a premixed flame, however, both exhibit a qualitatively similar linear response to axial flow perturbations – lowpass characteristics with an associated phase-lag (see, e.g., Ducruix, Durox, & Candel, 2000; Lighthill, 1954). The axial extent of the grids is 15 mm. This is much smaller than the wavelengths associated with the relevant mode frequencies (below 1000 Hz, see Section 4). The grids can thus be considered as compact heat sources whose spatial details are not important for the generation of acoustic waves but only their integral heat release rate. A photograph of one of the heating grids is shown in Fig. 2. The grid has a meandering pattern and is made from flat wire. Three independent DC sources (EA-PS 9080-50) are used to power the heating grids, each with a maximum output of 1.5 kW, supplying four grids in parallel. It was taken special care that all (custommade) heating grids are identical so that the nominal system would be as symmetric as possible. The relative difference in the cold electrical resistance of the 12 heating grids was less than 0.5%.

A flame in an annular combustion chamber of a gas turbine is subject to different types of velocity disturbances. The effect of axial fluctuations on the unsteady heat release response has been well studied in single-burner configurations, and this effect is also present in an annular chamber. In contrast to the single-burner case, a flame in an annular configuration is also subject to transverse velocity disturbances associated with the circumferential acoustic modes (O’Connor & Lieuwen, 2011). The latter effect is not present in the annular Rijke tube considered here because the sources of unsteady heat release are located in the tubes, in which only axial waves are present. It is assumed that this effect is generally dominant, as was found in large-eddy simulations of an annular combustion chamber by Staffelbach, Gicquel, and Poinsot (2009). No external source is used to drive a mean flow in the annular set-up. As in the original Rijke tube, the mean flow is solely convection induced. This restricts the parameter space (power input and mean flow velocity cannot be varied independently) but allows for essentially noise-free measurements. Moreover, the acoustic boundary conditions are well defined. The mean flow induced by the steady heat input is very small, of the order of 1 m/s in the present case. The Mach number (ratio of mean flow velocity to speed of sound) is therefore of the order of 10  3. Accordingly, any effects resulting from a non-zero Mach number are neglected in the acoustic model of the system (Section 3). Resulting from the order of rotational symmetry of the configuration, all eigenvalues with azimuthal orders unequal to 0 or multiples of 6 are degenerate (Perrin & Charnley, 1973). The eigenspace associated with this type of eigenvalue is two-dimensional. A suitable basis in the case of an annular duct is given by fcos mf,sin mfg, which represents two standing modes orthogonal to each other. In the following, the corresponding modes are referred to as cosine mode and sine mode of azimuthal order m. The rotating-mode state, which may appear as limit-cycle solution (Schuermans et al., 2006), is constructed from a superposition of both modes, their amplitudes being shifted by p=2 in temporal phase. Forty-eight microphones in total can be flush-mounted to the walls of the annular duct and the tubes. However, in this study, only 12 microphones (G.R.A.S. 40BP), one mounted on each tube upstream of the heating grids, were used (see Fig. 1). This measurement set-up allows a decomposition of the pressure field up to a circumferential mode order of m¼6. For the sixth-order azimuthal mode, however, only the cosine part can be detected. With reference to Fig. 1, it is easy to see that sin 6fj ¼ 0, where fj ðj ¼ 1, . . . ,12Þ correspond to the angular locations of the microphones. Fortunately, the sine component of this mode is of no importance due to the discrete rotational symmetry of order 12. It has pressure nodes at the angular locations of the tubes and therefore does not couple to the heat release fluctuations. In fact,

G. Gelbert et al. / Control Engineering Practice 20 (2012) 770–782

finite element computations showed that the two m¼ 6 modes of first axial order are separated by several hundred Hz. To facilitate feedback control, six identical speakers are mounted to the annular duct at equidistant angles (starting at f ¼ 301) and equal axial locations, see Fig. 1. With this actuator set-up, not all azimuthal modes up to order 6 can be controlled independently, which follows directly from the Nyquist criterion. Forcing an m¼5 mode results, for instance, in an excitation of m¼1 modes as well. The sine part of the m¼6 mode is not controllable, but as stated above, this mode is of no relevance. The cosine mode of azimuthal order 6 is important, as modeling and experimental results in later sections will show. This mode is controllable, but the actuation also forces the axisymmetric mode (m¼0). As mentioned in the introduction, a modal control scheme will be employed to suppress thermoacoustic oscillations in the experimental set-up. In this approach, the controller acts only on a single circumferential mode and is completely oblivious of the dynamics of the modes of other circumferential orders. Due to the insufficient circumferential actuator resolution, the controller may thus destabilize modes through backfolding, in principle. However, by the use of 12 microphones, the feedback loop of lower-order aliases is not closed. For example, although actuating an m¼6 mode with six speakers simultaneously forces an axisymmetric mode (m¼0), no feedback for the latter is established because the 12 microphones ensure sufficient circumferential resolution to separate the two. This is illustrated in Fig. 3, which depicts the principal azimuthal pressure variation of the cosine m¼6 mode with actuator and sensor locations. It is evident that the actuation pattern for cosine mode m¼6 is identical to that of the axisymmetric mode m ¼0, and therefore, forcing of only one of the two is not possible. However, the two modes can be well distinguished with the available sensor resolution. When controlling cosine mode m¼ 6, the axisymmetric mode will be forced in the same way, but since the two modes can be identified with the modal decomposition of the sensor signals, the latter mode is not part of the feedback loop. This holds not true for higher-order aliases, but in this case, the resonance frequencies of the targeted mode and the aliases are sufficiently spaced apart so that unstable feedback is highly unlikely. Although loudspeakers are frequently used as actuators for control of combustion oscillations in test-rig experiments, they cannot be used in real applications at elevated pressure and high power density because they lack the necessary actuation authority and robustness against the harsh environment (Dowling & Morgans, 2005). More realistic actuation devices in combustion systems are based on a modulation of the fuel flow, thus exploiting the stored chemical energy. An analogous type of actuation in the Rijke tube would be to modulate the heating power of the

Fig. 3. Principal azimuthal pressure pattern of cosine mode m¼ 6 with locations of actuators and sensors.

773

grids. While this is in principle possible, the heating grids in the present study cannot be used for this task. Considering the wire as a lumped capacity, the characteristic time for heat transfer and temperature equilibration is found to be of the order of 1 s. Therefore, when modulating the supply power at frequencies in the hundreds of Hz, which are the typical instability frequencies in this system, the wire temperature will remain essentially constant so that no unsteady heat transfer to the fluid occurs. The controllers were implemented on a dSPACE DS1103 PPC control board that was programmed using Simulink and the RealTime Interface of dSPACE. Unsteady pressure data were acquired with a custom-made amplifier–low-pass filter combination and National Instruments hardware (cDAQ-9172 with four NI 9215 simultaneous-sampling modules). Data were acquired at a sampling frequency of 10 kHz, and the controllers were run at a sampling frequency of 8192 Hz.

3. System model Fig. 4 shows a block diagram of the model of the annular Rijke tube. The model is set up as a network of acoustic elements, and the response of each block is represented by a state space system. Accordingly, the model is a time-domain model. In comparison to the cut plane in Fig. 1, the system is rotated by 901 in the clockwise direction. The cold upstream part is on the left and the hot downstream part is on the right-hand side. Except for the actuator voltage vact , all inputs and outputs are 12-dimensional vectors, which contain the acoustic pressures p or velocities u in the 12 tubes at the corresponding axial position. Since six speakers are mounted to the system, the actuator voltage vact is a sixdimensional vector. In order to access the pressures at the microphone positions pmic , the tubes upstream of the heating grids are split into an impedance Zopen ends and an admittance Atubes . Since the communication between the tubes is only taking place via the downstream annulus, the state space systems and accordingly the transfer matrices of the tubes, heating grids, and loudspeakers all have block-diagonal structure (e.g. Atubes ¼ diag½Atube1 ,Atube2 , . . . , Atube12 ). Table 1 provides an overview of the elements of the network model. The column ‘‘Quantity’’ shows how many of these elements are included in the corresponding block. The column ‘‘States’’ gives the number of states of the block, which is calculated from the number of elements multiplied by the order of each element. For example, the element Z open end models the impedance of a tube from the open end to a microphone with a state space system of order 7. The block Zopen ends consists of 12 of these elements and thus has the order 12  7¼ 84. Interconnection of the subsystems to a single state space description of the complete system can be conveniently achieved with Matlab routines. From the sum of the order of the subsystems, the order of the complete model adds up to 976.

Fig. 4. Block diagram of the network model. Each block is represented by a state space system.

774

G. Gelbert et al. / Control Engineering Practice 20 (2012) 770–782

Table 1 Elements of the network model. Element

Description

Quantity

States

Z open end Atube Heating grid Tube downstream Zannulus Gspeaker

Impedance of a tube from the open upstream end to the microphones Admittance of a tube from the microphones to the heating grids Transfer function of a heating grid Tube from the heating grids to the annulus Impedance of the annulus with 18 inputs and 12 outputs Transfer function of the speakers

12 12 12 12 1 6

12  7¼ 84 12  8¼ 96 12  1¼ 12 12  16¼ 192 550 6  7 ¼42

Complete model

6 inputs (actuator voltages), 12 outputs (pressure at the microphones)

976

Fig. 6. Pressure distribution of three modes in an axial cut plane of the annular Rijke tube from a 3D FEM calculation. Red and blue correspond to high and low acoustic pressure, respectively. Left: cosine mode m¼ 5, middle: sine mode m¼ 5, right: cosine mode m¼ 6. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this article.) Fig. 5. Schematic diagram of an annular duct. Consideration as a SISO system with a velocity excitation uk on the area Ak and a pressure response pj.

" with Ae ¼

The elements of the network and their derivation are explained in more detail in the following.

3.1. Impedance of the annular duct The impedance of an arbitrary acoustic system can be calculated by means of a modal expansion of the Greens function associated with the Helmholtz equation. A detailed description of the approach can be found in the work of Schuermans (2003) and Schuermans et al. (2003). If the largest extension of the area Ak (see Fig. 5) on which the system is excited is small compared to the acoustic wavelength, the SISO frequency response from a velocity excitation at xk to a pressure measurement at xj is given by 1 X p^ j ce ðxj Þce ðxk Þ : ¼ iocAk Hjk ðioÞ ¼ u^ k L ð o2 ioae o2e Þ e¼0 e

ð1Þ

Here and for the rest of this section, the pressure is scaled with the characteristic impedance rc, the product of density and speed of sound. In Eq. (1), o is the angular frequency, i is the imaginary unit and p^ j and u^ k are the frequency domain representations of the acoustic pressure at xj and the acoustic velocity at xk . The righthand side of the equation is an infinite sum over all acoustic modes ce and eigenfrequencies oe of the system, which are ordered such R 2 that oe o oe þ 1 . Furthermore Le ¼ V 9ce 9 dV is a normalization factor and ae the modal damping coefficient. Depending on the frequency range of interest, the sum can be truncated at a sufficiently high value. For simple geometries, such as tubes and annular ducts, eigenfrequencies and modes can be calculated analytically. For systems of arbitrary complexity, they can be obtained with standard finite element method (FEM) software. Each modal component e of Eq. (1) is represented by a secondorder state space system of the form x_ e ðtÞ ¼ Ae xe ðtÞ þ be,k uk ðtÞ, pj,e ðtÞ ¼ cTj,e xe ðtÞ þ djk,e uk ðtÞ

0

oe

oe

ae

cTj,e ¼ ½0 cce ðxj Þ,

"

# ,

be,k ¼

0

Ak

Le

#

ce ðxk Þ ,

djk,e ¼ 0,

ð2Þ

which can be easily checked with the relation He ðioÞ ¼ cTe ðioIAe Þ1 be þde . The complete pressure at xj then results P from the sum over the acoustic modes, i.e., pj ðtÞ ¼ e pj,e ðtÞ ¼ P T e ðcj,e xe ðtÞ þdjk,e uk ðtÞÞ. Taking into account Eþ1 eigenvalues, K inputs and J outputs, the state space model of the MIMO system follows from superposition of all modes and all inputs 2 3 2 32 3 2 32 3 b0;1 . . . b0,K A0 x0 u1 x_ 0 6 ^ 7 6 7 6 7 6 7 6 & ^ 54 ^ 7 & 4 5¼4 54 ^ 5 þ 4 ^ 5, bE,1 . . . bE,K AE xE uK x_ E 2

3 2 T p1 c1;0 6 ^ 7 6 ^ ¼ 4 5 4 pJ cJ,0

 & 

3 2 x0 0 7 6 6 ^ 54 ^ 7 þ 5 4^ cJ,E xE 0

cT1,E

32

 & 

3 32 u1 0 76 ^ 7 ^ 54 5: uK 0

ð3Þ

The dynamic matrix A of the element has dimension 2(Eþ1)  2(Eþ1), B is a 2(Eþ1)  K and C a J  2(Eþ1) matrix. Since the radial extent of the annulus is small compared to its axial and circumferential length, radial modes, with a dependence of the velocity and pressure distribution in the radial direction, do only appear at very high frequencies. For this reason, twodimensional eigenfunctions are sufficient for modeling the annulus. Using two-dimensional eigenfunctions, the index e in Eqs. (1) and (2) is a two-dimensional multi-index e ¼ fm,lg, with m as the azimuthal and l the longitudinal order of the mode. In the model of the annulus, 13 azimuthal (m ¼0y12) and 11 longitudinal (l¼ 0y10) modes are taken into account. This results in a state space system of order 550. The axial shape of the eigenfunctions used in the annulus is determined by the downstream boundary condition. Since the annulus is open on its downstream side, p ¼0 was used as boundary condition. For illustration, Fig. 6 shows the pressure distributions of three modes of the annular Rijke tube in an axial plane. The results are from a three-dimensional FEM calculation of the Helmholtz equation for the complete system, and the plane was chosen to be

G. Gelbert et al. / Control Engineering Practice 20 (2012) 770–782

775

in the annulus. The three modes shown are the ones that were recognized to be thermoacoustically unstable in the experiment (see Section 6.2). The modes have the azimuthal order cosine m¼5, sine m¼5, and cosine m¼6. The axial pressure distribution (not visible in the figure) corresponds in each case to an acoustic half-wave, with pressure nodes at the tube inlets and the outlet of the annulus. It can be seen that the modes do not vary in the radial direction, which indicates that the two-dimensional ansatz for the eigenfunctions in the annulus is appropriate. Furthermore, it is evident that the extrema of cosine mode m¼6 coincide with the azimuthal position of the tubes.

RLS ðHeÞ ¼ expðHe2 =22i  0:6133HeÞ from their calculations. The first term in the exponential determines the magnitude of the reflection coefficient, which decreases with increasing frequency. The second term in the exponential acts as a time delay (with dead-time 2  0.6133 R/c) and thus determines the phase. This term is also referred to as ‘‘end correction’’ because it has the effect to shift the pressure node out of the tube and thus makes it appear longer by an additional length of 0.6133  R. The impedance corresponding to this reflection coefficient is calculated from a Taylor expansion of Z ¼ ð1 þRÞ=ð1RÞ assuming small He and gives Z LS ðHeÞ ¼ ðHeÞ2 =4 þ i0:6133He.

3.2. Tubes

3.3. Heating grids

Since the tubes have an inner diameter of 60 mm, only the plane mode can propagate within the frequency range of interest ( o 1 kHz). All non-planar modes are cut-off, i.e., they decay exponentially in the axial direction (Pierce, 1981). The first azimuthal mode starts to propagate in the tubes at frequencies above 3 kHz. Therefore, the acoustic field in the tubes can be treated as one-dimensional. In principle, the impedance of a tube can be obtained from its acoustic modes as shown in the last subsection. However, a disadvantage of the modal approach is that the calculated state space systems only represent the acoustic resonance frequencies accurately but not the acoustic antiresonances. Although this error is getting smaller with increasing number of modes taken into account in the element, it gives rise to inaccuracies in the frequency response of the complete network model. For the correct frequency response of an acoustic network, the correct frequencies of the element’s antiresonances are as important as the correct frequencies of the resonances. For example, if two elements with correct resonance frequencies but inaccurate antiresonances are linked, the resulting compound element will have incorrect resonances and antiresonances. For this reason, a different approach was chosen for modeling the tubes. A plane acoustic field can be represented by a superposition of a downstream traveling wave f and an upstream traveling wave g. Neglecting mean flow and acoustic damping, these waves travel with the speed of sound c and without attenuation. If L is the length of the tube, the time the waves need to pass through the tube is given by L/c. Thus, the following relations between the waves on the left (index l) and right (index r) end of the tube are found:

In each tube, the relations between the up- and downstream acoustic variables (pus ,uus and pds ,uds ) are determined by the dynamic response of the heat source in that tube. In a general linear framework, this relation can be expressed by a 2  2 transfer matrix in frequency domain " # " # p^ ds p^ us ¼ TðioÞ : ð6Þ u^ ds u^ us

f^ r ¼ eioL=c f^ l

FðioÞ ¼

ð4Þ

and ioL=c

g^ l ¼ e

g^ r :

ð5Þ

To generate state space models of the tubes, Pade´ approximants are used to substitute the time delays. The longer the tube, the higher the order of the chosen approximation (see Table 1). The relationship of these waves with the acoustic pressure and the velocity used in the network model is given by the equations p ¼ f þ g and u ¼ f g. Special care was given to the boundary condition of the upstream tubes. The idealized impedance (Z¼p/u) of an open end is zero (because of p¼ 0 at the boundary), and the reflection coefficient is minus one. A reflection coefficient of minus one corresponds to full reflection in the tube and no radiation to the outside, which is an approximation and becomes increasingly inaccurate for higher frequencies. This is why a more realistic boundary condition was used instead. Levine and Schwinger (1948) calculated the frequency-dependent reflection coefficient for an unflanged open pipe end. For low Helmholtz numbers He ¼ oR=c, where R is the radius of the tube, follows

Since 12 identical heating grids are used, the response characteristics in each tube are assumed to be identical, too. For the type of compact heat source considered here, and for a vanishing mean flow Mach number, it is reasonable to assume that the pressure loss is negligible. Furthermore, the heat source induces a jump in acoustic velocity via (Lieuwen, 2003) Aðu^ ds u^ us Þ ¼

g1 ^ q, gp0

ð7Þ

where q^ and g denote the heat release rate and the ratio of specific heats, respectively; A is the cross-sectional area, and p0 represents the mean pressure. The heat release transfer function is introduced as the normalized response of heat release rate fluctuations to normalized perturbations in the upstream velocity FðioÞ ¼

q^ u us : u^ us q

ð8Þ

As proposed by Lighthill (1954) for the unsteady heat transfer from a circular wire, the heat release transfer function is modeled as a first order low-pass 1=2 , 1 þ ioct

ð9Þ

with phase parameter ct . The steady-state gain 1=2 corresponds to the linearized square-root dependence of the heat transfer coefficient on the approach flow velocity. For small wire Strouhal numbers odw =u us (dw is the wire diameter), Lighthill suggests ct ¼ 0:2dw =u us . From direct numerical simulations of the unsteady heat transfer from a circular wire, Kopitz (2008) found a larger constant of proportionality between ct and dw =u us of 0.4. The latter value is used in the following. Also, since the heating grids in the present set-up are made from flat wire, the corresponding hydraulic diameter (0.49 mm) for dw is used. With a measured mean flow velocity of 0.67 m/s, this yields ct ¼ 0:29 ms. In any case, the value of the phase parameter must be considered an estimate, subject to some inaccuracy. In addition to the uncertainty in the constant of proportionality, the mean velocity is not constant over the cross section, which is another source of uncertainty. Although the dynamical response of a heating grid is qualitatively similar to that of a premixed flame, there are some significant quantitative differences. The time delay of a premixed flame is typically of the order of a few milliseconds. This is an order of

776

G. Gelbert et al. / Control Engineering Practice 20 (2012) 770–782

magnitude larger than the time delay associated with the phase lag of the heating grid response (ct ). Also, the temperature increase across a lean premixed flame is of the order of 1500 K, whereas the heating grid used here raises the temperature only by about 80 K. Since the jump in acoustic velocity across the heat source is proportional to the relative temperature increase [see Eq. (10) below], the energy addition to the acoustic field through the unsteady response of the heating grid is comparatively small. Due to small damping in the system, linearly unstable thermoacoustic modes are present nonetheless; however, the growth rates are small, of the order of 0.1 s  1. To connect the upstream and downstream parts of the model, the downstream pressure is required as input to the heating grid transfer matrix (see Fig. 4). This is achieved by simply interchanging the up- and downstream pressure. Furthermore, the mean heat release in Eq. (8) is given by q ¼ r us u us Acp ðT ds T us Þ, where T us and T ds are the temperatures up- and downstream of the heater element, respectively, and cp is the specific heat at constant pressure. After some elementary manipulation, the transfer matrix of the heat source can be written as 3" " # 2 1 # x 0 p^ us p^ ds   4 5 ¼ : ð10Þ 0 1 þ TT ds 1 F u^ ds u^ us

Fig. 7. Bode diagram of the transfer function of excitation at speaker 1 to pressure measurement at microphone 2 ðT ¼ pmic ð2Þ=vact ð1ÞÞ. Top: magnitude, bottom: phase (wrapped).

us

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi In Eq. (10), x ¼ ðrcÞus =ðrcÞds ¼ T ds =T us is the ratio of characteristic impedances up- and downstream of the heater element. The temperatures are set according to measurements from the experiment, i.e., T us ¼ 296 K and T ds ¼ 374 K are used. Downstream of the heating grids, there is a weak decrease in temperature due to heat loss to the walls. Since this variation is not taken into account in the model, an average value is used for the downstream temperature. Corresponding to these temperatures, the speed of sound in the up- and downstream elements of the network are set to cus ¼ 345 m=s and cds ¼ 388 m=s. With the final model of the annular Rijke tube, both the hot system, with powered heating grids, and the cold system, with inactive heating grids, can be represented. Switching from hot to cold just requires setting T ds ¼ T us . In this case, the transfer matrix of the heating grid (10) becomes the identity matrix, and the speed of sound in the network’s downstream elements is identical to that in the upstream elements. As in the experimental set-up, the hot model is an unstable system, whereas the cold model is a stable system. Hence, by turning on the heating grids in the model, some eigenvalues move from the stable half-plane to the unstable half-plane. 3.4. Transfer function of the speakers The linear frequency responses of the annular Rijke tube from each speaker to each microphone were measured. Doing these measurement in a straightforward way, i.e., simply exciting the input and measuring the output, is only possible in stable systems. For this reason, measurements were done on the cold Rijke tube (heating grids turned off). For the excitation, a 32 times repeated sine sweep was used that linearly ran through frequencies from 160 to 1500 Hz within 1 s. Thus, the length of each dataset is 32 s. The measurements showed that the transfer functions from a speaker to a microphone all had a similar frequency response for actuator–sensor pairs with the same relative positions. For example, the transfer function from speaker 1 to microphone 1 is identical to that from speaker 1 to microphone 3 or from speaker 5 to microphone 9 or 11 (see Fig. 1), and so on. This was expected as the system is symmetric and six nominally identical speakers are used. Fig. 7 shows in black the measured frequency response from speaker 1 to microphone 2, which is mounted in the tube under the speaker (see Fig. 1).

In the model of the annular Rijke tube, the transfer function of a speaker relates the applied input voltage to the wall normal acoustic velocity generated in the hole the speaker is mounted in. This wall normal velocity is, in turn, an additional input to the impedance of the downstream annulus. Comparing the model’s frequency response without speaker transfer function (using the identity matrix for GSpeakers , see Fig. 4) with the measured frequency response, an additional high-pass characteristic and an additional constant phase decay were observed in the measurement. (Because of the propagation time of the acoustic waves from the speaker to the microphone, the model also shows a constant phase decay without the transfer function of the speaker.) This is reasonable since every speaker is a high pass and the amplifiers of the speakers may have small time delays. By looking at the differences between model and measurement, a second-order high pass with a corner frequency of fc ¼200 Hz, a gain of K ¼ 400 m=sV, and a time delay of T0 ¼0.6 ms was identified as speaker transfer function  2 io GSpeaker ðioÞ ¼ K eioT 0 : ð11Þ io þ2pf c With the identification of the speaker transfer function, the network model of the annular Rijke tube is completed. The frequency response of the complete model from speaker 1 to microphone 2 is shown along with the measurement in Fig. 7. For the frequency regime up to 600 Hz, the model reproduces the measured frequency response well, although the resonance of the model at 407 Hz is hardly visible in the measurement. This resonance is associated with a mode in the annulus with the order m¼1 in the azimuthal and l ¼1 in the longitudinal direction. Also, in the higher frequency regime, some differences between model and measurement are apparent. In the model of the annulus, a heuristic approach for the modal damping coefficient of mode ml

aml ¼ aoml with a ¼ 0:02

ð12Þ

is used. The damping of a mode thus simply scales with its eigenfrequency. By adjusting the individual modal damping coefficients according to the measurement, a better match between measurement and model could be achieved, but this was not the primary goal of this work. In the next section it will be shown that the modal frequency responses of the relevant modes are captured with good accuracy by the model.

G. Gelbert et al. / Control Engineering Practice 20 (2012) 770–782

4. Modal frequency response and control A more distinct frequency response compared to the excitation with one speaker is obtained by modal excitation of the system. For this purpose, cosine and sine modal base vectors for the actuation are defined as 2 3 cosðmfact,1 Þ 6 cosðmf 7 1 6 act,2 Þ 7 6 7 eact ð13Þ m,c ¼ 7 ^ lm 6 4 5 cosðmfact,6 Þ and 2 eact m,s ¼

sinðmfact,1 Þ

3

6 sinðmf 7 1 6 act,2 Þ 7 6 7 6 7 ^ lm 4 5 sinðmfact,6 Þ

ð14Þ

with

lm ¼

( pffiffiffi 6, pffiffiffiffiffiffiffiffiffi 6=2

m ¼ 0; 3,6, . . . , otherwise:

In these equations, m is the azimuthal order, fact,j the azimuthal angle of speaker j (see Fig. 1), and lm is a normalization coefficient. The normalization coefficients are introduced in order act act act to define an orthonormal base, i.e., eact m,c  en,c ¼ em,s  en,s ¼ dmn , in which dij is the Kronecker delta. In the same way, 12-dimensional base vectors (em,c and em,s ) for the modal decomposition of the 12 pressure readings are defined. The cosine and sine modal coefficients of the pressure are then calculated as p m,c9s ðtÞ ¼ eTm,c9s  pmic ðtÞ,

ð15Þ

where pmic ðtÞ is the vector with the pressure readings of the 12 microphones. With these definitions, the modal transfer functions can be measured. The system is excited in the shape of an azimuthal mode, i.e.  vact ðtÞ ¼ eact m,c9s v m,c9s ðtÞ

ð16Þ

is used as input, and the modal coefficient of the pressure is calculated using Eq. (15). The modal transfer function relates the modal input v m,c9s with the modal pressure coefficient p m,c9s , i.e., T m,c9s ðioÞ ¼ p^ m,c9s =v^ m,c9s . For the excitation, the same sine sweep signal as described in the previous section is used.

777

Fig. 8 shows in black the frequency response of the azimuthal modes m¼5 and cosine mode m ¼6 of the cold system (heating grids turned off). As stated in Section 2, sine and cosine mode m¼5 form a degenerate pair and thus have identical responses. The first resonance peak of modes m¼5 is at 522 Hz and the one of cosine mode m¼ 6 is at 525 Hz. In both cases, the phase drops by 1801 at the resonances, which indicates that the modes are stable. On the other hand, in the hot annular Rijke tube (heating grids turned on), the modes sine and cosine m¼ 5 and cosine m¼6 are unstable. Measurements in the limit cycle showed a double resonance peak at 558 and 562 Hz, resulting from these modes. Since the eigenfrequencies of the acoustic modes linearly scale with the speed of sound which, in turn, scales with the squareroot of the temperature, the resonance frequencies of the hot system are noticeably higher as the ones of the cold system. In order to design a model-based controller that can be used to stabilize the hot system in the experiment, a low-order model describing the instabilities is necessary. The complete model of the annular Rijke tube with its 976 states is too large for direct controller synthesis. Furthermore, since only some acoustic modes become unstable, not the complete system needs to be controlled. As the system is axially symmetric, the azimuthal modes are decoupled and separate controllers for each unstable azimuthal mode may be set up (modal control). In this way, the order of the model for controller synthesis can be considerably reduced. The models of the modes m¼5 and cosine m¼6 were obtained from a modal projection of the complete model. The input of the complete model is replaced with the modal input according to Eq. (16), and the new output, the modal pressure coefficient, follows from a modal projection according to Eq. (15). If the state space matrices of the complete model with input vector vact ðtÞ and output vector pmic ðtÞ (see Fig. 4) are given by A, B, C and D, the modal model for the azimuthal cosine mode m with input v m,c ðtÞ and output p m,c ðtÞ reads as _ ¼ AxðtÞ þ Beact  xðtÞ m,c v m,c ðtÞ,  p m,c ðtÞ ¼ eTm,c CxðtÞ þ eTm,c Deact m,c v m,c ðtÞ:

ð17Þ

The modal models for the azimuthal sine modes are obtained in the same way, using the sine base vectors for the projection. For all modes with azimuthal orders unequal to 0 or multiples of 6 (degenerate modes, see Section 2), it makes no difference whether sine or cosine base vectors are used for the projection; both

Fig. 8. Bode diagram of the modal transfer functions T m ðioÞ ¼ p^ m,c9s =v^ m of the cold system. Left: sine/cosine mode m¼ 5, right: cosine mode m¼ 6; top: magnitude, bottom: phase (wrapped).

778

G. Gelbert et al. / Control Engineering Practice 20 (2012) 770–782

projections result in the same model. In the following, mode m¼5 therefore refers to the sine and cosine component of the mode. Initially, the system matrix of this model is the same as that of the complete model, and accordingly, also the order is the same (976). However, since the azimuthal modes are decoupled in the complete model, most of the states in Eq. (17) become redundant and can be removed by calculation of a minimal realization. The frequency response of the modal models of the cold system for mode m ¼5 and cosine m ¼6 is shown in Fig. 8 along with the measurement. There is good agreement with the measurement, although the damping is slightly too large around the second resonance and too small in the higher frequency regime. As discussed in Section 3, this could be adjusted by manually adapting the damping coefficients of the model. However, this was not done as model accuracy was deemed sufficient. In the hot case, the calculation of a minimal realization resulted in a model of order 81 for mode m ¼5 and a model of order 59 for mode m¼6. The frequency response of these models is shown in black in Fig. 9. The first resonance of mode m¼5 is at 558 Hz and the one of mode m¼6 is at 561 Hz. This is in very good agreement with resonances that were measured in the annular Rijke tube in the limit cycle (558 and 562 Hz). In contrast to the cold case (Fig. 8), the phase now increases by 1801 at these resonances. This indicates that these resonances are unstable or, in other words, that a pair of unstable eigenvalues belongs to each of it. For controller synthesis, the gain K in the transfer function of the speakers [see Eq. (11)] was set to unity. This scaling has the advantage that all variables of the model have a similar magnitude. For this reason, the frequency response in Fig. 9 is shifted by 52 dB compared to the one with the original gain (K ¼ 400 m=sV). The constant decay in the phases of both modes follows from the time delay of 1.21 ms that is inherent to the system. This time delay has two contributions – one from the time delay in the transfer function of the speakers [0.6 ms, see Eq. (11)] and the other from the propagation time of the acoustic waves from the speakers to the microphones (45 mm at 345 m/s plus 186 mm at 388 m/s ¼0.61 ms, see Fig. 1). Before the modal models of the hot system were used for controller synthesis, their order was further reduced using balanced truncation (Laub, Heath, Paige, & Ward, 1987). This resulted in models of order 12 for both modes. The frequency response of the reduced models is also shown in Fig. 9,

demonstrating good agreement with the original models. From both reduced models, a linear-quadratic regulator with weighting of the output was calculated (lqry command in Matlab) and a Kalman state estimator. The controller synthesis was done in continuous time domain, and since the two models have a quite similar frequency response, the same weights can be used for both controllers. The input and output weights of the regulator were chosen to be unity, i.e., Rlqry ¼ Q lqry ¼ 1. This simple choice was possible due to the aforementioned scaling of the input. In the Kalman estimator, the spectral density of the input noise was set to Q kalm ¼ 10 000 and the one of the measurement noises was set to Rkalm ¼ 1. This choice follows from the assumption that the measurement is accurate and greater uncertainty in the model. Fig. 10 shows the frequency responses of the controllers, the closed loops and, for comparison, again the plants. The transfer function of the controller was calculated by combining the state regulator and the observer in one block. If this block is denoted by CðioÞ and the transfer function of the plant PðioÞ, the complementary sensitivity, which describes how measurement noise is processed in the closed loop, is given by Scomp ðioÞ ¼

CðioÞPðioÞ : 1 þCðioÞPðioÞ

ð18Þ

For both modes, it can be observed that the phase of the closed loop now drops again by 1801 at the first resonance, which indicates that the controller stabilizes its mode. At the resonances, the magnitude of the controllers starts to decline. This shows that the controllers essentially affect only the unstable resonances. The schematic in Fig. 11 shows how the modal controllers are implemented. Input to each controller is the modal coefficient of the pressure and output is the modal actuation variable. Since the modal model of mode m¼5 is the same for the sine and the cosine component, also the same controller can be used to control the sine and cosine part of this mode. The modal pressure coefficients are obtained from an online modal decomposition according to Eq. (15), and the input vector for the speakers, vact ðtÞ, follows from a modal recombination according to Eq. (16). The recombination results from a superposition over the unstable modes, i.e., the modal actuation variables multiplied with their act   base vector are summed vact ðtÞ ¼ eact m5,c v m5,c ðtÞ þ em5,s v m5,s ðtÞ þ  eact v ðtÞ. m6,c m6,c

Fig. 9. Bode diagram of the modal transfer functions of the hot system. Models from modal projection (black) and reduced models of order 12. Left: sine/cosine mode m¼ 5, right: cosine mode m¼6; top: magnitude, bottom: phase.

G. Gelbert et al. / Control Engineering Practice 20 (2012) 770–782

779

Fig. 10. Bode diagrams of the reduced models of the hot system (open loop), the closed loops, and the controllers. Left: sine/cosine mode m¼ 5, right: cosine mode m¼6; top: magnitude, bottom: phase.

Fig. 11. Schematic representation of the modal control approach

5. Single-mode model and modal model from projection In the present work, the modal models of the unstable modes were derived from a modal projection of the complete model. This approach was described in the preceding section [see Eq. (17)]. Schuermans (2003) and Schuermans et al. (2003) also used modal control to stabilize an unstable annular thermoacoustic system, but the manner with which the modal models are derived is different. In this section, the main differences between the modal model used by Schuermans et al., termed single-mode model in the following, and the modal model used in the present work, termed modal model from projection, are discussed. Thereby, also interesting basic characteristics of the complete network model, termed complete model in the following, are described. To obtain a modal model for a single azimuthal mode m, Schuermans et al. set up a network model that relates the modal amplitudes of the signals rather than the pressure and velocity signals themselves. The basic structure of the model remains identical to that of the complete model (see Fig. 4) but the annulus is represented with a state space system that only contains the azimuthal mode m of interest [see Eq. (3)]. The annulus is then connected with only one tube and one heat source. As a result, not only the order of the annulus is reduced, but also all states from the remaining K 1 tubes and heat sources (in the present case K ¼12) can be omitted. Using such a singleazimuthal-mode approach to obtain modal models for the unstable modes of the experiment (sine and cosine mode m¼5 and cosine mode m¼6), it was realized that the eigenvalues of the derived modal models did not accurately match the corresponding eigenvalues of the complete model. While investigating the source of this discrepancy, the following observations with respect to the complete model were made. If, for example, the state space model of the annulus contains azimuthal modes only up to the order M¼3 (i.e., m¼0, 1, 2, 3), also the complete model (annulus with tubes) contains azimuthal modes only up to this order. This is not surprising since the tubes only communicate among each other via the annulus. The modes

of the complete model represent acoustic modes of the complete system (and not acoustic modes of the subsystems). If another mode is added in the annulus (M¼4), additional eigenvalues appear in the complete model that are associated with the azimuthal modes m¼4 of the complete system. Initially, the addition of a higher mode does not affect the eigenvalues of the lower modes (m¼0, 1, 2, 3). In other words, the complete model with M ¼4 contains exactly the same eigenvalues as the complete model with M¼3 and some additional eigenvalues; the complete model with M¼ 5 contains all eigenvalues of the complete model with M¼4; and the complete model with M ¼6 contains all eigenvalues of the complete model with M ¼5. This behavior changes if the order of modes taken into account in the annulus is chosen to be larger than half the number of attached tubes, i.e., M 46 in the present case. As a result of the order of discrete rotational symmetry of the system, the higher modes also affect the eigenvalues of the lower modes. For instance, mode m¼7 couples with mode m ¼5 and modifies the associated eigenvalues, mode m¼ 8 couples with mode m¼4, mode m¼9 couples with mode m¼3 and so on. It is observed that the coupling and the associated effect on the eigenvalues is weaker for lower modes (m¼0, 1, 2) compared to higher modes (m ¼3, 4, 5). Generalizing this observation, the following can be said for a network model of an annular combustion system with K tubes connected to the annulus. If the model is intended to quantitatively predict the behavior of azimuthal modes with an order close to half the number of tubes (i.e., m ¼ K=21, K=22,K=23), azimuthal modes with an order larger than half the number of tubes (m ¼ K=2 þ1,K=2 þ 2,K=2 þ3) should be included in the model of the annulus. According to its construction, the single-mode model cannot take into account any coupling with other azimuthal modes. For this reason, such a modal model correctly reproduces the eigenvalues of the complete model only if there is no coupling between the azimuthal modes in the complete model, or in other words, if the order of azimuthal modes taken into account in the annulus of the complete model is lower than M oK=2. Summarizing the previously described observations, the single-mode model is only suitable to describe modes of low azimuthal order (m oK=4, in the present case m¼0, 1, 2). As Schuermans et al. (2003) consider a system with K ¼24 tubes connected to the annulus and the azimuthal modes m¼ 1 and m¼2 are unstable in their simulation, this restriction causes no problems in their work. However, due to the modal coupling that is introduced by the finite number of tubes, the single-mode

780

G. Gelbert et al. / Control Engineering Practice 20 (2012) 770–782

model is inaccurate for modes of the order K=4 r m rK=2, and it fails for higher modes (m 4 K=2). In Fig. 12, the low-frequency eigenvalues of the system with powered heating grids are shown. The chosen presentation differs from the standard used in control engineering. Compared with the standard, the figure is rotated by 901 in the clockwise direction, and negative frequencies are not shown. However, all eigenvalues appear as complex conjugate pairs. The squares are the eigenvalues of the complete model with 13 azimuthal and 11 longitudinal modes taken into account in the annulus (i.e., modes up to order M¼ 12 in the azimuthal and L¼10 in the axial direction). Each of the eigenvalues is associated with an acoustic mode of the complete system. The eigenvalues that are numbered in the figure all belong to modes that have the shape of a half-wave in the axial direction (as described at the end of Section 3.1), and the number indicates the azimuthal order of the corresponding mode (from m¼ 0 to m¼6). The three eigenvalues that are not numbered in the upper right corner of the lower frame all belong to the next group of modes, which have the shape of a full wave in the axial direction. The upper frame in the figure shows in detail the eigenvalue of mode m¼4 (at 547 Hz) and the unstable eigenvalues (negative damping coefficient) of modes m¼5 (at 558 Hz) and m¼ 6 (at 561 Hz). The eigenvalues of modesm¼1,2,3,4 and 5 all appear twice – a manifestation of the degenerate eigenspace resulting from rotational symmetry. The eigenvalue of mode m¼ 6 is distinct and corresponds to the cosine part only. The circles represent the corresponding eigenvalues of the single-mode model. For modes m¼0 (at 162 Hz), m¼1 (at 234 Hz), and m¼2 (at 370 Hz), the divergence from the complete model is small, but for modes m¼3, 4, 5, and 6, the agreement is unsatisfactory; mode frequencies and damping are overpredicted. Thus, in this example, the single-mode model fails to predict that modes m ¼5 and m¼6 become unstable when the heating grids are powered. The lowest eigenvalues of the m¼5 and m ¼6 modal models from projection are shown as crosses in Fig. 12. Since these modal models are extracted from the complete model [see Eq. (17)], their eigenvalues always exactly match the corresponding

Fig. 12. Low-frequency eigenvalues of the models with powered heating grids. Squares: complete model with 13 azimuthal and 11 longitudinal modes taken into account the annulus, circles: eigenvalues of the single-mode models, crosses: eigenvalues of the modal models of mode m ¼5 and m¼ 6 from a projection of the complete system.

eigenvalues of the complete model, regardless of the number of modes taken into account in the annulus. The only drawback of the method is the relatively high computational cost that is associated with the calculation of the minimal realization.

6. Results 6.1. Simulation results In this section, results from a simulation of the complete network model (Section 3) connected to the modal controllers described in Section 4 are shown. One intention of the simulation was to test whether the modal controllers, which were designed on the basis of the reduced-order models, were able to stabilize the original, complete network model. For the simulation, a Simulink model was set up, and two additional elements were included. To account for the nonlinear saturation of the real system, saturation blocks were introduced at the input to the heat release transfer functions. Input of these transfer functions are the upstream velocities [see Eq. (8)], and the reason for introducing the saturation blocks at this location is the following. The pressure fluctuations associated with combustion instabilities typically do not exceed a few percent of the mean pressure. Linearity of the acoustic field can therefore be assumed even on the limit cycle (Dowling, 1997). In contrast, the unsteady heat release response is mainly a hydrodynamic effect and therefore starts to saturate at velocity fluctuation levels comparable to the mean velocity. This is true for flames as well as heating grids (Balachandran, Ayoola, Kaminski, Dowling, & Mastorakos, 2005; Heckl, 1990). In order to approximately match limit-cycle amplitudes observed in the experiment, the saturation level was set to 0.25 m/s. Furthermore, source terms were introduced by superimposing the downstream velocity of the heating grids with uncorrelated white noise of small amplitude. These signals ensure that the instabilities start growing at the beginning of the simulation, and they also simulate the measurement noise. In the last section, it was shown that sine and cosine mode m¼5 and cosine mode m¼ 6 become unstable when the heating grids are turned on. Fig. 13 shows the acoustic pressure at microphone 3 and the modal amplitudes of the three unstable modes over time. At the beginning of the simulation, no controller is active, and all three unstable modes start growing. After 2.8 s, the amplitudes of cosine and sine mode m¼5 start decreasing

Fig. 13. Results from a simulation of the complete network model. Top: pressure  at microphone 3, below: modal amplitudes pðtÞ of the unstable modes cosine m¼5, sine m ¼5, and cosine m¼ 6. All pressures are given in kPa. The vertical dashes indicate when the controllers are activated.

G. Gelbert et al. / Control Engineering Practice 20 (2012) 770–782

while cosine mode m¼6 still increases until it reaches the limit cycle bounded by the saturation in the heat source. The emergence of a limit-cycle oscillation can be explained qualitatively in the following way. All effects responsible for a reduction of the oscillation energy (radiation from the tubes, damping in the annular duct) are linear in the acoustic variables and thus in the system state. On the other hand, the oscillation energy is increased by the heat release–acoustic interaction, represented through the transfer function of the heating grid. Since limiters were added at the input of the heat release transfer function, there is now an amplitude dependence in this mechanism that, from a describing function point of view, essentially decreases the gain for larger oscillation levels (Dowling, 1997). In the case of linear instability, the oscillation energy input is larger than the damping at small amplitudes. However, since the input increases sublinearly at larger amplitudes due to the saturation, it will coincide with the damping at some finite oscillation level. This represents an equilibrium between oscillation energy gain and loss and therefore defines the amplitude of the limit-cycle oscillation. The mechanism leads to a coupling between the unstable modes, such that the mode that reaches the limit cycle prevents the growth of the other modes. Since mode m¼ 6 is now dominant, the controller for this mode is turned on first, at t¼4 s. The control action immediately stabilizes the mode. With the modal amplitude of mode m¼ 6 also the pressure at microphone 3 decreases. Since the system is now no longer limited by the saturation, sine and cosine mode m¼5 start growing again. In the following, a spinning mode m¼5 with cosine and sine part shifted by p=2 in time develops. At 9.0 s the controller of cosine mode m¼5 is turned on additionally and stabilizes this mode. Now only sine mode m ¼5, which is a standing mode, is left, and a slight increase of its amplitude is observed. At 11 s, also the controller for the sine component is turned on and mode m¼5 is suppressed. The system is completely stabilized and thus also the pressure at microphone 3 reduces to the noise level. 6.2. Experimental results Fig. 14 shows results from measurements at the annular Rijke tube with powered heating grids. At the beginning, no controller is active, and the system is in the limit cycle. In contrast to the simulation, not only cosine mode m¼ 6 is dominating the system, but also sine and cosine mode m ¼5 have significant contributions. The reason for this difference might be a weak coupling between the modes, such that mode m¼ 6 excites the two m¼5 modes and vice versa. In the model, such a reciprocal excitation does not exist. On the other hand, the mechanism of reciprocal suppression, observed in the simulation, also exists. Accordingly, in the limit cycle state, the modal amplitudes are also coupled by the nonlinear saturation of the system. The two opposed coupling mechanisms lead to a permanent ‘‘struggle’’ between the modes. This can be observed in the magnification of the first 4 s on the right-hand side of Fig. 14. For example, the amplitude of sine mode m ¼5 has two maxima in the first 0.5 s (indicated by vertical dashes in the figure). A decrease in the amplitude of mode m¼6 can be observed when mode m ¼5 is large and vice versa. A similar phenomenon is clearly visible between 2 and 2.6 s. Whenever cosine mode m¼6 has a maximum in amplitude, cosine mode m ¼5 has a minimum. At t¼ 11.2 s, the controller of cosine mode m ¼5 is activated and stabilizes this mode immediately. Concurrently, the amplitude of sine mode m¼ 5 starts increasing and reaches larger maximum values as before. Also, the pattern of the pressure measured at microphone 3 changes since now two standing waves are dominating the system. At 21.4 s, the controller of sine

781

Fig. 14. Experimental results. Top: pressure at microphone 3, below: modal  amplitudes pðtÞ of cosine mode m¼ 5, sine mode m¼ 5 and cosine mode m¼ 6. All pressures are given in kPa. Left: full plot, right: magnification of the first 4 s. Vertical dashes correspond to specific instances in time referred to in the discussion in the text.

mode m ¼5 is activated in addition and stabilizes this mode. At the same time, cosine mode m ¼6 grows and reaches larger amplitudes as before. Finally, at 30.8 s, also the controller of cosine mode m¼6 is turned on, thereby stabilizing the complete system. The pressure oscillations at microphone 3 completely disappear. As in the simulation, the sequence that is used to switch on the controllers is of no importance. As soon as a controller is activated, the corresponding mode is stabilized, and when the controller is switched off, the mode starts to grow again until it reaches the limit cycle. A complete stabilization of the system is only achieved with all three controllers running in parallel. The spectra of the acoustic pressure for the case without and with control are shown in Fig. 15. With the three controllers, the peak amplitude is reduced by more than 60 dB.

7. Discussion and conclusion Model-based control was successfully applied to suppress thermoacoustic instabilities in an annular Rijke tube. For timedomain simulations and controller synthesis, a thermoacoustic network model was set up. Because of its high order, the model could not be directly used to build a controller. To reduce the order for controller synthesis, a modal control approach was used. Here, multiple controllers run in parallel, each one is controlling one unstable azimuthal mode.

782

G. Gelbert et al. / Control Engineering Practice 20 (2012) 770–782

Fig. 15. Spectra of the acoustic pressure at microphone 3 without and with control. Experimental results.

The modal models of the unstable modes were obtained from a modal projection of the complete model. From the modal models of mode m ¼5 and cosine m¼ 6 of the hot system, reduced-order models were calculated. These reduced models were used to design a linear-quadratic regulator for each unstable mode. The modal control approach made it possible to control all three unstable modes individually, both in the simulation and in the experiment. Running all three controllers in parallel, a complete stabilization of the system was achieved, reducing the peak oscillation amplitude by more than 60 dB in the experiment. Besides the use of the model for controller synthesis and timedomain simulations, it can be further employed for a detailed analysis of the dynamic properties of the system. As shown in Fig. 12, each low frequency eigenvalue of the model can be related to an acoustic mode of the complete system. Calculating the eigenvalues from the state space system with 976 states requires only a few seconds on a standard PC. Thus, the model is applicable for studying the effects of parameter variations. For example, the effects of geometrical asymmetries as well as asymmetries in the heat release response of the heating grids can be easily studied with the model. The frequency responses of the modal models of the cold system (heating grids turned off) were compared with the frequency responses from the measurements and showed good agreement. This good agreement is an indicator for the high quality of the model. Although a real gas turbine is in many aspects more complex than the generic configuration considered in the present work, the same modeling methodology can be applied. During the design stage of a new combustor, such models can be used to predict the thermoacoustic behavior of the engine. Time consuming and costly fixes on prototypes suffering from thermoacoustic problems could be reduced or even avoided in this way (Banaszuk, Mehta, & Hagen, 2007).

Acknowledgments Financial support from the German Research Foundation through the collaborative research center Control of complex turbulent shear flows (SFB 557) is gratefully acknowledged. The authors also wish to thank the anonymous reviewers for their thorough review, which helped to significantly improve an earlier draft. References Balachandran, R., Ayoola, B. O., Kaminski, C. F., Dowling, A. P., & Mastorakos, E. (2005). Experimental investigation of the nonlinear response of turbulent premixed flames to imposed inlet velocity oscillations. Combustion and Flame, 143(1–2), 37–55.

Banaszuk, A., Ariyur, K. B., & Krstic´, M. (2004). An adaptive algorithm for control of combustion instability. Automatica, 40, 1965–1972. Banaszuk, A., Mehta, P., & Hagen, G. (2007). The role of control in design: From fixing problems to the design of dynamics. Control Engineering Practice, 15(10), 1292–1305. Bellucci, V., Flohr, P., Paschereit, C. O., & Magni, F. (2004). On the use of Helmholtz resonators for damping acoustic pulsations in industrial gas turbines. Journal of Engineering for Gas Turbines and Power, 126, 271–275. Campos-Delgado, D. U., Schuermans, B. B. H., Zhou, K., Paschereit, C. O., Gallestey, E., & Poncet, A. (2003). Thermoacoustic instabilities: Modeling and control. IEEE Transactions on Control Systems Technology, 11(4), 429–447. Candel, S. (2002). Combustion dynamics and control: Progress and challenges. Proceedings of the Combustion Institute, 29(1), 1–28. Dowling, A. P. (1997). Nonlinear self-excited oscillations of a ducted flame. Journal of Fluid Mechanics, 346, 271–290. Dowling, A. P., & Morgans, A. S. (2005). Feedback control of combustion oscillations. Annual Review of Fluid Mechanics, 37(2), 151–182. Ducruix, S., Durox, D., & Candel, S. (2000). Theoretical and experimental determinations of the transfer function of a laminar premixed flame. Proceedings of the Combustion Institute, 28(1), 765–773. Gelbert, G., Moeck, J. P., Bothien, M. R., King, R., & Paschereit, C. O. (2008). Model predictive control of thermoacoustic instabilities in a swirl-stabilized combustor. AIAA Paper 2008-1055. Heckl, M. A. (1990). Non-linear acoustic effects in the Rijke tube. Acustica, 72, 63–71. Illingworth, S. J., & Morgans, A. S. (2010). Adaptive feedback control of combustion instability in annular combustors. Combustion Science and Technology, 182(2), 143–164. Kopitz, J. (2008). Kombinierte anwendung von str¨ omungssimulation, netzwerkmodellierung und regelungstechnik zur vorhersage thermoakustischer instabilit¨ aten. ¨ Munchen. ¨ Ph.D. thesis, Technische Universitat Krstic, M., & Banaszuk, A. (2006). Multivariable adaptive control of instabilities arising in jet engines. Control Engineering Practice, 14(7), 833–842. Laub, A., Heath, M., Paige, C., & Ward, R. (1987). Computation of system balancing transformations and other applications of simultaneous diagonalization algorithms. IEEE Transactions on Automatic Control, 32(2), 115–122. Levine, H., & Schwinger, J. (1948). On the radiation of sound from an unflanged circular pipe. Physical Review, 73(4), 383–406. Lieuwen, T. (2003). Modeling premixed combustion–acoustic wave interactions: A review. Journal of Propulsion and Power, 19(5), 765–781. Lieuwen, T. C., & Yang, V. (Eds.). (2005). Combustion instabilities in gas turbine engines. Progress in astronautics and aeronautics (Vol. 210). AIAA, Inc. Lighthill, M. J. (1954). The response of laminar skin friction and heat transfer to fluctuations in the stream velocity. Proceedings of the Royal Society of London A, 224, 1–23. McManus, K. R., Poinsot, T., & Candel, S. (1993). A review of active control of combustion instabilities. Progress in Energy and Combustion Science, 19, 1–29. Moeck, J. P., Bothien, M. R., Paschereit, C. O., Gelbert, G., & King, R. (2007). Twoparameter extremum seeking for control of thermoacoustic instabilities and characterization of linear growth. AIAA Paper 2007-1416. Morgans, A. S., & Dowling, A. P. (2007). Model-based control of combustion instabilities. Journal of Sound and Vibration, 299, 261–282. Morgans, A. S., & Stow, S. R. (2007). Model-based control of combustion instabilities in annular combustors. Combustion and Flame, 150, 380–399. Noiray, N., & Schuermans, B. (2012). Theoretical and experimental investigations on damper performance for suppression of thermoacoustic oscillations. Journal of Sound and Vibration, 331(12), 2753–2763. O’Connor, J., & Lieuwen, T. (2011). Disturbance field characteristics of a transversely excited burner. Combustion Science and Technology, 183(5), 427–443. Perrin, R., & Charnley, T. (1973). Group theory and the bell. Journal of Sound and Vibration, 31(4), 411–418. Pierce, A. D. (1981). Acoustics: An introduction to its physical principles and applications. New York: McGraw-Hill. Raun, R. L., Beckstead, M. W., Finlinson, J. C., & Brooks, K. P. (1993). A review of Rijke tubes, Rijke burners and related devices. Progress in Energy and Combustion Science, 19(4), 313–364. Rayleigh, J. W. S. (1878). The explanation of certain acoustical phenomena. Nature, 18, 319–321. Riley, A. J., Park, S., Dowling, A. P., Evesque, S., & Annaswamy, A. M. (2004). Advanced closed-loop control of an atmospheric gaseous lean-premixed combustor. Journal of Engineering for Gas Turbines and Power, 126, 708–716. Schuermans, B. (2003). Modeling and control of thermoacoustic instabilities. Ph.D. thesis, E´cole Polytechnique Fe´de´rale de Lausanne, Switzerland. Schuermans, B., Bellucci, V., & Paschereit, C. O. (2003). Thermoacoustic modeling and control of multi-burner combustion systems. ASME Paper 2003-GT-38688. Schuermans, B., Paschereit, C. O., & Monkewitz, P. (2006). Non-linear combustion instabilities in annular gas-turbine combustors. AIAA Paper 2006-0549. Seume, J. R., Vortmeyer, N., Krause, W., Hermann, J., Hantschk, C.-C., Zangl, P., et al. (1998). Application of active combustion instability control to a heavy duty gas turbine. Journal of Engineering for Gas Turbines and Power, 120, 721–726. Staffelbach, G., Gicquel, L. Y. M., & Poinsot, T. (2009). Large eddy simulation of selfexcited azimuthal modes in annular combustors. Proceedings of the Combustion Institute, 32(2), 2909–2916. Stow, S. R., & Dowling, A. P. (2003). Modelling of circumferential modal coupling due to Helmholtz resonators. ASME Paper GT2003-38168.