## Abstract

This paper brings a novel mathematical perspective in assessing the rise of the secondary dynamic modes to prominence during the suppression of thermoacoustic instability. This phenomenon is observed by many earlier investigators; however, without a complete analytical reasoning. We consider a Rijke tube with both a passive Helmholtz resonator and an active feedback control to suppress instabilities. The core dynamics is represented as a linear time-invariant multiple time-delay system of neutral type. Parametric stability of the resulting infinite-dimensional dynamics is investigated using a recent analytical tool: cluster treatment of characteristic roots paradigm. This tool reveals the stability outlook of such systems exhaustively and non-conservatively in the parameter space of the system. First, we examine the stability with and without the Helmholtz resonator. We then select an unstable operation for the resonator-mounted Rijke tube, impose a time-delayed integral feedback control over it and reveal the stabilizing controller parameters using the cluster treatment of characteristic roots methodology. When high control gains are inappropriately selected, the new analytical procedure declares how the secondary dynamic modes of the system exhibit instability although the initially unstable mode is now stabilized. All of these stability assessments are cross-validated using experimental results from a laboratory-scale Rijke tube set-up.

## 1. Introduction

Thermoacoustic instabilities are self-excited oscillations in an acoustic enclosure which are caused by a coupling between the unsteady heat release and the surrounding pressure waves. Regionally confined heat release drives the acoustic waves which propagate within the enclosure and partially reflect from the boundaries. As they arrive back at the heating zone, they interact with the unsteady heat release disturbances again after some acoustic time delays. If the two effects, the acoustic perturbations and the perturbations caused by the heat release, synchronize they amplify and cause thermoacoustic instabilities. This famously known phenomenon, first hypothesized by Rayleigh [1], has been the mainstay of many investigations over the centuries. Thermoacoustic instabilities are observed in complex combustion systems such as gas turbines, aero-engines, afterburners [2,3], as well as relatively simple laboratory-scale abstractions such as a Rijke tube [4–8].

Suppression of thermoacoustic instabilities is mainly achieved via two approaches: passive and active. The passive approach aims to add acoustic damping to the system in order to dissipate the energy, for instance, with the aid of Helmholtz resonators or acoustic liners and alike [9,10]. When resonators are mounted on the combustor, an effective suppression is achieved at some narrow frequency ranges, which can be adjusted only by physically altering the resonator design (i.e. tuning [11]). Active approach, on the other hand, includes an additional feedback control loop in the system, with the use of sensors (e.g. microphone) and actuators (e.g. loudspeaker, fuel valve). An important advantage of the active control techniques is the tunability of the frequency ranges in which the operation remains effective.

In general, the passive absorbers perform well in attenuating the amplitude of pressure oscillations. However, they are typically unable to eliminate them completely; in other words, these devices may not suffice to stabilize the unstable thermoacoustic dynamics. In that regard, combination of active and passive approaches may achieve better results. The passive absorber removes the bulk of the acoustic energy and the active control reduces the residual pressure fluctuations to practically zero. As the active controller treats the already attenuated acoustic pressures (after the passive absorption), the energy required for it would be relatively small, which makes the operation highly desirable. An interesting study on such a deployment is [12] which also uses a Helmholtz resonator for passive absorption. In addition, a feedback control is used to optimally manipulate the amplitude and the frequency of the oscillations of the resonator’s back-plate. This operation is akin to structures which are known as ‘noise cancellation’ devices applied at multiple modes concurrently.

In this paper, we use a Rijke tube as the experimental platform. For control purposes, we use a passive Helmholtz resonator and separately deploy a feedback control loop consisting of a microphone and a loudspeaker (sensor and actuator, respectively). This combined construct of the Rijke tube together with both the active and passive control elements leads to a complex system dynamics which will be studied here.

Under certain assumptions, it has been shown that the mathematical model of the thermoacoustic dynamics in the Rijke tube reduces to a linear time-invariant multiple time-delayed system of neutral type (LTI-NMTDS) [4,5,7]. The infinite-dimensional nature of the dynamics is preserved in this time-delayed model and the structure can now benefit from the existing mathematical arsenal on TDS (e.g. [13]) for the prediction and the suppression of thermoacoustic instability.

Among many alternatives, the feedback control law we select in this investigation is time-delayed integral (TDI) type, which has also been used in other investigations. Although the time delays are generally undesirable in controllers, it has been shown in various studies that when they are carefully selected, these delays could improve the closed-loop performance considerably [14,15]. A similar control logic, so-called the phase-shift control, has also been widely used earlier to combat thermoacoustic instabilities [5,16]. In all of these earlier studies; however, a targeted suppression effort was conducted on a particular unstable mode exclusively. Therefore, there is always a possibility of inviting the instabilities at secondary frequencies corresponding to other dynamic modes in the system which are ignored [5,17]. We prevent such a risk in this study, by holistically synthesizing the feedback controller to achieve the asymptotic stability of the closed-loop system dynamics including *all* the infinitely many modes. This is the key mathematical strength in this text.

The paper starts with the derivations of the LTI-NMTDS model of the Rijke tube dynamics containing a resonator and a feedback control loop. The stability of the combined dynamics is analysed using a mathematical paradigm called the cluster treatment of characteristic roots (CTCR) [18]. It reveals the conditions for stable and unstable operations in the domain of the system parameters, exhaustively and non-conservatively. First, we present these results for the uncontrolled Rijke tube with and without the resonator. Then, we choose an unstable operation for the resonator-mounted tube and attempt to stabilize the dynamics using a feedback controller. The controller parameters that can achieve this stabilization are exhaustively identified. These declarations yield the highlight theme of the paper by displaying how certain control settings could invite instability in higher order dynamic modes of the system, while the originally unstable mode is being stabilized. This fact was observed in many earlier studies without a justifying analytical framework, which is covered in this paper. Finally, we verify all the analytical findings with experimental results over a laboratory-scale Rijke tube.

## 2. Mathematical model and dynamic stability

In this section, we discuss the mathematical model of the thermoacoustic dynamics in a Rijke tube combining the Helmholtz resonator and the feedback control effects. The dynamics in the Rijke tube is governed by the first principles of conservation of mass, momentum and energy [16,19]. In the following derivations, we use a general notation of

When linearized under the listed assumptions, the pressure and velocity fluctuations, respectively, *f*(*x*,*t*) and *g*(*x*,*t*) represent the acoustic waves travelling downstream and upstream in the tube, respectively ([5] and as in [20]). They are denoted as *f*_{i}(*t*) and *g*_{i}(*t*) in figure 1*a*, with the subscript *i*=1,…,7 indicating the cross sections along the tube that the functions are evaluated at. *x* is the position along the tube, and *x*=0 corresponds to the location of the heater. The upstream and downstream tube ends are marked as *x*=−*x*_{u} and downstream *x*=*x*_{d}, where *x*_{u} and *x*_{d} denote the distances between ①–② and ③–⑦, respectively. The Helmholtz resonator is mounted at *x*=*x*_{h} and the microphone is located at *x*=*x*_{m} to record the downstream pressure fluctuations. The microphone senses the sound pressure signal which is fed through the controller and linear amplifier, and finally back into the system as a controlling pressure wave via a loudspeaker at the upstream end of the tube.

The corresponding block diagram to this system is given in figure 1*b*. In this representation, the system variables are the acoustic waves given in the Laplace domain as *F*_{i}(*s*) and *G*_{i}(*s*), *i*=1,…,7 that propagate upwards and downwards along the tube, respectively. The time delays

### (a) Thermoacoustic dynamics in Rijke tube

All the cross-sectional variables in figure 1*b* are interconnected through some causal relations. Next, we present a series of such sectional relationships starting with the interaction of acoustic waves at the heating zone (*x*=0) and resonator zone (*x*=*x*_{h})
**0**_{2} is a vacuous matrix of size 2×2. The descriptors of **H**_{h} (transfer matrix at the heating zone) and **H**_{r} (transfer matrix at the Helmholtz resonator attachment) can be found in appendices A and B, respectively. The acoustic pressure variables are connected to one another through some time delays as indicated in figure 1*b*. These relations are captured, using the same set of variables as in equation (2.3), in the form of some transport delays
*R*_{u} and *R*_{d}. In addition, the feedback control signal *W*(*s*) (Laplace domain representation of *w*(*t*)) acts as a control input to the system. This function denotes the speaker-induced additive acoustic pressure fluctuation at the upstream end of the tube and it is incorporated into the system dynamics as follows:
**R** is written in terms of the output vector in (2.3) as
*W*(*s*), as it is a critical element in the overall system behaviour.

### (b) A distinguishing philosophy for the feedback control loop and relevant discussions

We consider a broadly used feedback control configuration as shown in figure 1*b*. The control command *W*(*s*) should be looked at as an ensemble of several dynamical components: downstream acoustic pressure sensor, the feedback control logic, linear power amplifier and the loudspeaker (actuator). Such a combined effect can be modelled as follows:
*K*(*s*) denotes the actual control logic which is user defined, whereas *L*(*s*) represents the combined transfer function of all the remaining components in the feedback loop, such as the microphone, linear amplifier and the loudspeaker. *P*_{6}(*s*) is the combination of *F*_{6}(*s*) and *G*_{6}(*s*) as per (2.1) at the cross section ⑥. In this investigation, we adopt a TDI control logic, which consists of a control gain *K* and a control delay *τ*_{5}:
*τ*_{5} in (2.8) is characteristically selected as a ‘phase shift operator’ when the dynamic suppression is aimed at a particular mode and the corresponding frequency. This logic comes with a controversial proposition. Such modal suppression attempts need to first identify the mode in which the tube dynamics will exhibit instability, and this knowledge is gained *a priori* to the control implementation. Once this troublesome frequency is identified, the relevant phase shift is determined by selecting an appropriate delay *τ*_{5}. In this study, we ignore the interlinking of the control delay with the targeted frequency of instability. This point constitutes the highlight theme of the paper, that the system at hand is infinite-dimensional, thus it has infinitely many modes, and any one of them can exhibit instability while the conventional controller is occupied with suppression of a specific frequency. To preclude such an occurrence the dynamic suppression (a.k.a. stabilization) has to be performed in a holistic fashion and to assure stability in all the infinitely many modes of the system, not only a few. As such the delay in (2.8) has to be conceived as a completely free and independent control parameter to be selected for the design of the control logic. Just to give the reader a preamble, what provides the means to this objective here is the mathematical tool for assessing the complete stability of LTI-NMTDS, so-called the CTCR paradigm.

One critical element in (2.7) is *L*(*s*) as it represents the dynamic properties of many components. To identify *L*(*s*) with some level of fidelity, we conduct a set of open-loop system identification tests using a data acquisition system (dSpace MicroAutoBox), linear amplifier, loudspeaker (New Wave Audio, IC-8S) and the microphone. The microphone is positioned facing the loudspeaker and the set-up is structured as shown in figure 2. Note that the Rijke tube-related elements are kept completely out of this exercise. A chirp signal is provided as the input with a unity gain, covering a frequency range of 10–1000 Hz for a duration of 100 s. The data acquisition takes place at a sampling rate of 10 kHz. A transfer function is conceived using empirical and iterative model matching trials to represent this set-up, such that the frequency response characteristics concur between those obtained from experimental data and from the model. The resulting transfer function for the given system is

The comparison of Bode plots using the experimental data and the model transfer function of (2.9) is shown in figure 3. The identified transfer function matches very well with the experiment within the range of 150–600 Hz, although there is an acceptable level of deviation at frequencies outside this range. A higher order, but more representative transfer function with better fidelity could still be obtained. This would, however, increase the computational overhead further. For the objectives of the paper it is found to be sufficient to adopt *L*(*s*) as in (2.9).

### (c) Closed-loop system dynamics

Now that the feedback-controlled pressure wave *W*(*s*) is expressed in terms of the acoustic system variables, we substitute (2.7) in (2.5) and get
**I**_{4} is 4×4 identity matrix and **M** represents the entire system matrix. Equation (2.12) presents a classical eigenvalue problem and **M** leads to a determinant which forms the characteristic equation of the system

### (d) Stability analysis of linear time-invariant multiple time-delayed system

In its general form, the characteristic equation for linear time-invariant multiple time-delay systems can be written as
*i* and *k*. *z*_{i} represents the number of delay compositions associated with *s*^{i} term. The time-delay vector is defined as *m* denotes the number of delays. If *s*^{n} multiplies some terms which also contain delays, that is, if there is a *k*∈(1,2,…,*z*_{n}) such that *α*_{nk}≠0 and

### Theorem 2.1 ([23])

*A necessary condition for asymptotic stability of LTI-NMTDS such as* (2.14) *is* *It declares that the discrete operator in* (*2.15*) *has its entire infinite spectra located in*

### Theorem 2.2 ([23])

*The necessary and sufficient condition for any LTI-TDS to be stable is that (i) the infinite spectra of its characteristic equation, such as in* (*2.14*), *lies in* *and (ii) χ*_{CE}(*s,*** τ**)

*is bounded away from*

Theorem 2.1 is an apparent precondition to theorem 2.2. It enforces the requirement (ii) in theorem 2.2 [24]. If the difference operator (2.15) is stable, we only need to determine the conditions for which all the characteristic roots of (2.14) to lie in

Considerably large scholarly work appeared especially in the last six decades on the time-delay systems primarily addressing their stability characteristics (see a survey [13]). Some numerical tools and approximate root finding methods are also developed as in [25]. These routines operate on systems for a fixed set of delays and approximate the dominant characteristic root locations of the infinite-dimensional LTI-MTDS in

In this study, we focus on the stability of (2.13) in the six-dimensional parametric space of *x*_{u}≤*x*≤*x*_{d}), or the feedback control formation. In such parametric spaces a point-wise stability search would simply be prohibitive due to computational limitations. Instead, we adopt a remedial procedure here, which is called the CTCR paradigm. The CTCR declares the stability boundaries exhaustively and non-conservatively in the space of such system and operating parameters. This capability generates the so-called *stability maps* in a desirably select parametric space, which reveal the stable and unstable operating conditions. The interested reader is referred to [18], on pp. 329–332, and references therein for further details of the paradigm.

## 3. Main results

In this section, stability properties of various systems represented by the characteristic equation (2.13) are investigated analytically using the CTCR paradigm. We approach the key objectives of this paper through some systematic stages.

(1) First, the stability of the uncontrolled Rijke tube without the resonator is studied. The heater is intentionally positioned in the tube such that the operation becomes unstable.

(2) Then, CTCR analysis of the tube with a Helmholtz resonator is performed. The resonator design is kept suboptimal, again on purpose and its location on the tube is intentionally selected such that the system still operates in an unstable mode even after including the resonator.

(3) On this unstable system, the TDI control logic is deployed and the overall system stability variations are studied for the controlled dynamics. The controller parameters that stabilize the system are identified exhaustively, again using the CTCR-generated stability maps.

(4) At the very end of this chain of tests, we demonstrate the key contribution: a given feedback control logic, if not selected pursuant to a holistic analytical methodology, may in fact excite the secondary modes of the system, while suppressing the unstable mode of the dynamics. This critical nuance was extensively observed in earlier investigations; however, a crisp analytical reasoning was never found until this study.

All these new revelations of CTCR are validated step by step using the results from the corresponding experimental tests which are conducted on a laboratory-scale Rijke tube. A picture of the set-up is provided in figure 4 and the critical parameters of the test platform are provided in table 1.

### (a) Uncontrolled and passively controlled Rijke tube

We start with the uncontrolled plain Rijke tube, without a Helmholtz resonator. The characteristic equation of the dynamics in the absence of the feedback loop and the resonator can be obtained by simply substituting *K*(*s*)=0 and *e*^{−2(τ1+τ2+τ3+τ4)s}|=1.16>1 yields infinitely many zeros all of which lie in *χ*_{D}(*s*,** τ**)<0, is satisfied. Next, the asymptotic stability of the complete system given in equation (3.1) is tested as per theorem 2.2, using the CTCR paradigm. Its stability map is created in the selected parametric space of tube length (i.e.

*a*. This display is unique to the CTCR paradigm as explained earlier.

In figure 5*a*, the shaded region represents stable operating conditions. That is, when Rijke tube length and heater location composition (*L*,*x*_{u}) is selected outside this shaded area, the pressure oscillations are expected to grow in amplitude exponentially (a linear system characteristic). When such an unstable test is conducted the growing amplitudes eventually settle into a limit cycle (a nonlinear system characteristic) [7,16]. The colour-coded curve that separates the stable and unstable regions represents the loci of (*L*,*x*_{u}) compositions, for which (3.1) has a pair of characteristic roots on the imaginary axis. The system is marginally stable along this curve, and the colour coding represents the imaginary root locations (resonant frequencies) in Hertz. Furthermore, the number of unstable roots (NU) of the characteristic polynomial (3.1) is noted in each region. When *NU*=0, all the characteristic roots are in

As given in table 1, the tube length in our experiment is *L*=0.508 m. In the following tests, we intentionally fix the heater location at *x*_{u}=0.102 m (point A in figure 5*a*), which makes the thermoacoustic dynamics in the system unstable. On this unstable Rijke tube, we investigate the effects of a common passive control approach by mounting a Helmholtz resonator. The effectiveness of the resonator in dynamic suppression is directly correlated to the proximity of the frequency at which the resonator offers maximum damping and the frequency of unstable pressure oscillations [27]. We select the geometric dimensions of the resonator (see appendix B and table 1) such that its resonant frequency (254.7 Hz) is considerably different from the frequency of unstable pressure oscillations at point A (349.1 Hz). This indicates a sub-optimum resonator design for passive control of the unstable Rijke tube dynamics. Taking *K*(*s*)=0 (i.e. no feedback control) and forming the **H**_{r} matrix as in Appendix B, characteristic equation (2.13) exhibits a new form
*D*_{2}(*s*,** τ**)=

*D*

_{1}(

*s*,

*τ*). Therefore, the necessary condition of theorem 2.1,

*χ*

_{D}(

*s*,

**)<0, is automatically satisfied. Fixing the tube length at**

*τ**L*=0.508 m, the stability map for (3.3) is generated using the unique procedure of CTCR once more, but this time in the domain of

*b*. Similar to figure 5

*a*, the shaded region is stable (corresponding to

*NU*=0) and the resonant frequencies are colour-coded on the curves separating the stable and unstable regions in figure 5

*b*. While the heater is at

*x*

_{u}=0.102 m (point A in figure 5

*a*), we attach the Helmholtz resonator at

*x*

_{r}=0.330 m, which corresponds to point B in figure 5

*b*. The new operating configuration purposely selected once again within the unstable region, but with some different characteristics as discussed next.

To compare the two different unstable operating conditions at points A and B, we deploy the Quasi-Polynomial mapping-based Rootfinder (QPmR) algorithm [25]. QPmR is a numerical tool that approximates the dominant characteristic roots of the system characteristic equation (3.3). Note again that QPmR is a ‘point-wise’ stability analysis tool in contrast with ‘global’ CTCR methodology that produces the stability maps in a broad range of system parameters. A brief comparison of the point-wise numerical tools and the CTCR paradigm was offered earlier in §2d. Using QPmR on equations (3.1) and (3.3) for points A and B, respectively, we obtain the dominant root distributions as presented in figure 6*a*. Note that there is only a pair of characteristic roots in *a*. To verify this, the pressure oscillations are recorded experimentally (at 10 kHz sampling rate) as the system operates at these points. In figure 6*b*, these time traces of pressure oscillations for unstable operating points A and B are displayed (without and with the resonator, respectively). The fundamental frequencies of these two pressure readings are extracted at the onset of the respective exponential growth phases, performing sectional fast Fourier transformations (FFT). The corresponding power spectral densities are shown in figure 6*c*, with the experimentally measured centre frequencies at *f*^{A}_{ex}=349.1Hz, *f*^{B}_{ex}=175.8Hz. These frequencies are also calculated using the dominant roots of the conceived models as shown in figure 6*a*, *f*^{A}_{mod}=340.8Hz, *f*^{B}_{mod}=175.3Hz. The disagreements between the analytically obtained values and those from the experiments are very small, which gives a strong encouragement to the procedure we present here.

As an important observation, the amplitudes of pressure oscillations after the limit cycle behaviour sets in figure 6*b*, are much smaller for operating point B when compared with those of point A. One can claim that the addition of the Helmholtz resonator to the Rijke tube as a passive control element reduces the residual oscillatory energy, which is exhibited by lower limit cycle amplitudes. It is obvious that, the resonator is absorbing some energy, but without being able to stabilize the system as the CTCR stability maps declared in figure 5*b*. Therefore, the proposed analytical methodology used here provides a definitive new insight to the stability posture of the complex dynamics.

### (b) Feedback-controlled Rijke tube

Finally, we take the unstable operating conditions of B, and introduce a feedback loop using the control law (2.8) on it. The objective is to stabilize the unstable dynamics of the resonator-mounted Rijke tube, without targeting any modal frequency. Once again, the present approach decouples the treatment from some particular modes of the system. Therefore, it is a holistic procedure of stabilization in the truly infinite-dimensional sense.

When we consider the controller gain and feedback delay (*K*,*τ*_{5}) as the only free parameters and rewrite (2.13) using the model in (2.9), it exhibits a sixth order LTI-NMTDS
*K*,*τ*_{5}) form the domain where the stability maps will be explored using the CTCR paradigm.

The difference operator of (3.4) is again identical to (3.2), so the necessary condition in theorem 2.1, *χ*_{D}(*s*,** τ**)<0 holds. We fix the earlier mentioned parameters

*L*=0.508 m,

*x*

_{u}=0.102 m and

*x*

_{r}=0.330 m (

*x*

_{h}=0.076 m) to identify the operating conditions at point B in figure 5

*b*, but this time with feedback control. The microphone, which is used as the sensor to measure the pressure fluctuations, is placed at

*x*

_{m}=0.356 m. The stabilizing controller designs in the domain of (

*K*,

*τ*

_{5}) are exhaustively declared via CTCR as shown in figure 7

*a*along with the marginal stability frequencies at the stability boundaries (indicated by colour coding) and the regional NU as marked. We take a closer look at it next.

The origin of figure 7*a* corresponds to the uncontrolled dynamics at point B in figure 5*b*, which is unstable at a modal frequency of *f*^{B}_{mod}=175.3Hz. When the control as suggested in (2.8) is deployed using the parameters (*K*,*τ*_{5}) from within the shaded region (where NU=0), the system is stabilized. But a striking feature in this figure appears for high-gain control applications such as *K*>5000. If the controller parameters are not carefully selected, it is clear that, the system can be destabilized at some higher dynamic modes of the system, although the initial unstable dynamic mode remains stabilized. This feature has been observed by other researchers evidenced by the considerable literature [8] and it forms the key contribution of this study. Same discovery is mentioned without justification in [5] as: ‘Simple controllers involving just a fixed time delay of phase-shift introduce a new oscillation mode, which becomes unstable as gain is increased and gives rise to a new peak in the pressure spectrum’. The excitement in the present findings is that the CTCR paradigm provides a holistic methodology to identify the destabilizing design selections as opposed to those that impart stability for high-gain time-delayed controller so that undesired instabilities at other frequencies (i.e. higher modes) will not occur. This document provides a novel tool for the design of a stabilizing controller, without knowing the particular frequency (or the mode) at which the system tends to exhibit the instability.

To demonstrate this crucial point, we have selected three control settings: (*K*,*τ*_{5})=(8000,1.8 ms), (*K*,*τ*_{5})=(8000,3.1 ms), (*K*,*τ*_{5})=(6000,2.6 ms), which are represented by points C, D and E in figure 7*a*. The NU of (3.4) declared by the CTCR for these cases are *NU*=4, *NU*=2 and *NU*=0, respectively. Therefore, we expect to observe two unstable modes at point C, one unstable mode at point D and stable operation at point E.

The QPmR-generated characteristic root distribution of equation (3.4) is plotted in figure 7*b*, for (*K*,*τ*_{5}) compositions at these points. The initial unstable system with a characteristic root pair in *f*^{B}_{mod}=175.3Hz in figure 6*a*) is now stabilized by selecting (*K*,*τ*_{5}) at point E. The remarkable feature appears at this junction. Although this unstable root pair can be brought to *K*,*τ*_{5}) at points C or D, the other characteristic roots (i.e. those of the higher order modes) of the system migrate to *a*,*b* shows that the system has two pairs of characteristic roots in *a priori* selection of targeted modal suppression forms the main highlight of the present investigation.

Figures 8, 9 and 10 display the experimental results corresponding to those points C, D and E in figure 7*a*. In panel (a) of these three figures, the microphone-detected pressure fluctuations are shown. The control is off at the beginning, so the pressure fluctuations start to grow exponentially and then they are set into a limit cycle in all three figures. Until the control is turned on, we observe the same scenario as for point B in figure 6*b*. The instant the control triggers is marked with a thin red line. From that moment on the pressure fluctuations start decaying in all three cases; only to start an exponential growth after some time in figures 8*a* and 9*a* (for operating points C, D). Although the initial unstable mode is stabilized for all three (*K*,*τ*_{5}) selections, when the control parameters are selected at points C, D, secondary dynamic modes of the system are destabilized. This phenomenon exhibits itself as the pressure fluctuations growing in amplitude at relatively higher frequencies. On the other hand, when the controller parameters are set to point E, the system stays stable following explicit exponential decay behaviour (figure 10*a*).

In figures 8*b*, 9*b* and 10*b*, zoomed versions of the pressure fluctuations just after the control is turned on are shown. Linear system behaviour of exponential growth in the first two, and exponential decay in the last are very clearly observed. In figures 8*c*, 9*c* and 10*c*, we display the control signals which are produced at the DAC (D/A conversion) port of the real-time control card. In figures 8*d*, 9*d* and 10*d*, the corresponding sound pressure levels of the pressure fluctuations are displayed immediately before and after the controller is triggered. For all three cases, sound pressure levels show a peak at around 178.2 Hz before control is turned on, which are measured after the first limit cycle sets. This is expectedly very close to *f*^{B}_{ex}=175.8Hz, which was measured at the onset of instability (during the exponential growth without the control). After control is triggered, figure 8*d* shows two peaks, which are at *f*^{C2}_{ex}=902.1Hz. When compared with the imaginary parts of the corresponding characteristic roots for point C in figure 7*b*, *L*(*s*) and the experimental data for frequencies over 600 Hz, are likely to be accountable for this discrepancy. Needless to say, the much earlier assumptions which are made during the creation of the mathematical model may also have an impact on such disagreement.

A similar comparison is made for point D, between figures 9*d* and 7*b*. In this case turning the control on destabilizes one higher frequency mode at *f*^{D}_{ex}=754.4Hz in figure 9*d*. The corresponding frequency extracted from the model in figure 7*b* is *f*^{D}_{mod}=733.4Hz. There is a much better match in this case between the experimentally measured and analytically obtained frequencies.

And finally at point E, the control action is expected to bring a complete stability, and it does not excite any other dynamic modes of the system. After the control is triggered and the pressure oscillations decay in amplitude, no exponential growth is detected (figure 10*a*). Accordingly, there is no noticeable peak in figure 10*d* after the control is turned on because the system is stable at all natural modes. The QPmR tableau in figure 7*b* shows no characteristic roots in *NU*=0 in figure 7*a*. Therefore, this controller configuration renders a stable operation indeed, as predicted by the CTCR-generated stability map in figure 7*a*.

## 4. Conclusion

This paper proposes an exhaustive stability analysis technique to the study of thermoacoustic instability in a Rijke tube combined with active and passive sound pressure suppression elements. It is shown that the mathematical model of the overall system falls into the linear time-invariant multiple time-delay systems class of neutral type. This representation preserves the well-known infinite-dimensional nature of the problem without focusing on a particular modal behaviour. Therefore, the problem is treated in a holistic fashion. CTCR paradigm is used to generate the needed stability maps of the system, which reveal stable and unstable operating conditions in the space of the system parameters such as geometric dimensions of the Rijke tube, passive resonator location or feedback controller parameters. It is shown that careful design of the controller suppresses the unwanted thermoacoustic oscillations in the Rijke tube. Destabilization of some secondary dynamic modes (or higher order modes) of the system is observed when very high controller gain is deployed. Although this phenomenon has been known for quite some time from earlier experimental investigations, a complete analytical justification of it in the broad parametric space was missing. This point constitutes the major contribution of this work, that CTCR paradigm provides an analytical framework to detect these interesting features over the space of the selected operating parameters. The analytical results are favourably supported with experimental findings conducted on a simple Rijke tube set-up.

## Data accessibility

The logic behind the Cluster Treatment of Characteristic Roots paradigm that generates the stability maps can be accessed through [18]. The Quasi-Polynomial mapping-based Rootfinder (QPmR) algorithm and its detailed implementation is available in [25].

## Authors' contributions

U.Z. derived the mathematical model, implemented the theory, carried out the experiments and helped draft the manuscript; N.O. conceptualized the connection between thermoacoustic instability and time-delay systems, coordinated the entire study, developed the theory and drafted the manuscript.

## Competing interests

Authors have no competing interests.

## Funding

This research is partially supported by the National Science Foundation grant CMMI-1462301 and UCONN (University of Connecticut) research excellence program.

## Appendix A. Heating zone transfer function

The transfer matrix **H**_{h} represents the effect of heat injection on the incoming and outgoing acoustic pressure fluctuations at *x*=0 as
*Z*_{2}=1/*b*. Here, *a* and *b* are the gain and time constant of the first-order flame transfer function between the velocity and heat release fluctuations which is broadly adopted in the Rijke tube literature [6–8]. *γ* is the heat capacity ratio and *A* represents the cross-sectional area of the Rijke tube. Interested readers may find further details on the derivation of this heating zone dynamics within the cited references.

(A2) is sufficient to capture the dynamic relation at the heating zone when an electrical coil is used as the heat source [7,8]. If a flame is used; however, this may introduce some complexities in the dynamics (e.g. the effect of hydrodynamic region [28]), increasing the order of the transfer matrix in (A2).

## Appendix B. Resonator zone transfer function

The transfer matrix **H**_{r} represents the effect of the Helmholtz resonator on the incoming and outgoing acoustic waves at location *x*=*x*_{h} as
**H**_{r} is expressed as follows:
*C*_{r}=*C*/*M*, *S*_{n} and *l*_{eff}=*l*_{n}+0.85 *d*_{n} (*d*_{n} is the neck diameter) are the cross-sectional area and the effective length of the resonator’s neck. *V* is the volume of the resonator’s cavity. Interested readers are directed to the two references above for further details on the derivation of the Helmholtz resonator dynamics.

One must take note that nonlinearities become dominant if the airflow has a non-negligible Mach number as clearly pointed out in [30]. In such cases, the simple linear mechanical vibration absorber analogy and (B2) may not hold anymore and nonlinear analysis may be required as in [27]. In this study; however, the focus is on the airflow created only by natural buoyancy, and with negligible Mach number as clarified at the beginning of §2. This feature is necessary in adopting the linearized dynamics given in (B2).

- Received March 11, 2016.
- Accepted June 3, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.