Thermo-acoustic instability: Theory and experiments

June 28-30, 2015. Ann Arbor, MI, USA Proceedings Proceedings of of the the 12th 12th IFAC IFAC Workshop Workshop on on Time Time Delay Delay Systems Systems Proceedings of the 12th IFAC Workshop on Delay June 28-30, 2015. Ann Arbor, MI, USA Proceedings of the 12th IFAC Workshop on Time Time DelayatSystems Systems Available online www.sciencedirect.com June 28-30, 2015. Ann Arbor, MI, USA June 28-30, 2015. Ann Arbor, MI, USA June 28-30, 2015. Ann Arbor, MI, USA

ScienceDirect IFAC-PapersOnLine 48-12 (2015) 075–080 and experiments Thermo-acoustic instability: Theory Thermo-acoustic instability: Theory and experiments Thermo-acoustic instability: Theory and experiments Thermo-acoustic instability: Theory and experiments Thermo-acousticUmut instability: Theory and Zalluhoglu*, Nejat Olgac** experiments

Umut Umut Zalluhoglu*, Zalluhoglu*, Nejat Nejat Olgac** Olgac** Umut Umut Zalluhoglu*, Zalluhoglu*, Nejat Nejat Olgac** Olgac** *University of Connecticut, Storrs, CT 06269-3139 USA (e-mail: [email protected]). *University of Connecticut, Storrs, CT 06269-3139 USA (e-mail: [email protected]). ** University of Connecticut, Storrs, CT 06269-3139 USA (e-mail: [email protected]) *University of Connecticut, Storrs, CT 06269-3139 USA (e-mail: [email protected]). *University of Storrs, CT USA (e-mail: [email protected]). *University of Connecticut, Connecticut, Storrs, CT 06269-3139 06269-3139 USAUSA (e-mail: [email protected]). ** University of Connecticut, Storrs, CT 06269-3139 (e-mail: [email protected]) ** ** University University of of Connecticut, Connecticut, Storrs, Storrs, CT CT 06269-3139 06269-3139 USA USA (e-mail: (e-mail: [email protected]) [email protected]) ** University of Connecticut, Storrs, CT 06269-3139 USA (e-mail: [email protected]) Abstract: This article deals with thermo-acoustic instability (TAI) on a simple Rijke tube, an extensively Abstract: This articletest deals with thermo-acoustic instability (TAI)model on a simple Rijke tube, and an extensively studied combustion platform. We revisit theinstability mathematical of the Rijke dynamics show that Abstract: This deals with (TAI) on tube, an Abstract: This article articletest deals with thermo-acoustic thermo-acoustic instability (TAI)model on aaa simple simple Rijke tube, and an extensively extensively Abstract: This article deals with thermo-acoustic instability (TAI) on simple Rijke tube, an extensively studied combustion platform. We revisit the mathematical of the dynamics show that when linearized it falls into the class of linear time-invariant neutral time delay systems. In orderthat to studied combustion test platform. We revisit the mathematical model of dynamics and studied combustion test into platform. We of revisit thetime-invariant mathematicalneutral model time of the thedelay dynamics andInshow show that studied combustion test platform. We revisit the mathematical model of the dynamics and show that when linearized it falls the class linear systems. order to analyze its stability characteristics, we implement a recent mathematical paradigm called the cluster when linearized linearized it it falls falls into into the the class class of of linear linear time-invariant time-invariant neutral neutral time time delay delay systems. systems. In In order order to to when when linearized it falls into the (CTCR). class linear neutral timeparadigm delay systems. In to analyze itsofstability stability characteristics, weofimplement implement recent mathematical called theorder cluster treatment characteristic roots It can time-invariant reveal the stable operating regions incalled the space of the analyze its characteristics, we aaa recent mathematical paradigm the cluster analyze its stability characteristics, we implement recent mathematical paradigm called the cluster analyze its stability characteristics, we implement a recent mathematical paradigm called the cluster treatment of characteristic roots (CTCR). It can reveal the stable operating regions in the space of the operationalofparameters, suchroots as the geometric dimensions of stable the Rijke tube. The results of this prediction treatment characteristic (CTCR). It can reveal the operating regions in the space of the treatment ofparameters, characteristic (CTCR). It reveal stable operating regions in the space of treatment characteristic roots (CTCR). It can can reveal the the stable operating regions in of thethis space of the the operational suchroots ascompared the geometric dimensions offindings. the Rijke Rijke tube. The The results prediction mechanismofparameters, are very favorably with experimental operational such as the geometric dimensions of the tube. results of this prediction operational parameters, such the dimensions the operational parameters, such as ascompared the geometric geometric dimensions of offindings. the Rijke Rijke tube. tube. The The results results of of this this prediction prediction mechanism are very favorably with experimental mechanism are very favorably compared experimental findings. © 2015, IFAC Federation of with Automatic Hosting by Elsevier Ltd. All rights reserved. mechanism are very compared with experimental  Control)findings. mechanism are(International very favorably favorably compared with experimental findings.  linearized dynamics following the procedures described in 1. INTRODUCTION  linearized dynamics following described in Dowling (1997). This leads the to aprocedures linear time-invariant, linearized dynamics following the procedures described 1. INTRODUCTION linearized (1997). dynamicsThis following the procedures described in in 1. INTRODUCTION linearized dynamics following the procedures described in Dowling leads to a linear time-invariant, INTRODUCTION A critical problem 1. in the modern-day gas turbines is the multiple time delay system, which is infinite dimensional. 1. INTRODUCTION Dowling (1997). This leads to aa linear time-invariant, Dowling (1997). This leads to linear time-invariant, Dowling (1997). This leads to a linear time-invariant, A in turbines is time itdelay which is infinite dimensional. thermo-acoustic instability (TAI), whichgas is also known as multiple Furthermore, falls system, into a notoriously cumbersome special A critical critical problem problem in the the modern-day modern-day gas turbines is the the multiple system, which dimensional. A in modern-day gas turbines is the multiple time time itdelay delay system, which is is infinite infinite dimensional. A critical critical problem problem in the the modern-day gas turbines is old the multiple time delay system, which is infinite dimensional. thermo-acoustic instability (TAI), which is also known as Furthermore, falls into a notoriously cumbersome special combustion instability. Although TAI is a 1½ century sub-class which is known as neutral type. The stability thermo-acoustic instability instability (TAI), (TAI), which which is is also also known known as as Furthermore, it into cumbersome special thermo-acoustic Furthermore, it falls falls into aaa notoriously notoriously cumbersome special thermo-acoustic instability (TAI), which is 1½ alsocentury known as analysis Furthermore, it falls into notoriously cumbersome special combustion instability. Although TAI is a old sub-class which is known as neutral type. The stability problem, especially last two decades have seen numerous of this mathematically complex dynamics is a combustion instability. instability. Although Although TAI TAI is is aa 1½ 1½ century century old old sub-class which is known as neutral type. The stability combustion sub-class which is known as neutral type. The stability combustion instability. Although TAI is a 1½ century old sub-class which is known as neutral type. The stability problem, especially last two decades have seen numerous analysis of this mathematically complex dynamics is a scientific advances in its prediction and control (Dowling challenging task. In earlier investigations the Nyquist problem, especially especially last last two two decades decades have have seen seen numerous numerous analysis mathematically complex is problem, analysis of of this this mathematically complex dynamics dynamics is aaa problem, especially last two decades havecontrol seen numerous analysis of this mathematically complex dynamics is scientific advances in its prediction and (Dowling challenging task. In earlier investigations the Nyquist 1997; Annaswamy et al., 2000; Selimefendigil et al., 2013; Criterion is used to tackle this problem for a given system scientific advances advances in in its its prediction prediction and and control control (Dowling (Dowling challenging task. earlier investigations the scientific challenging task.to In Intackle earlier investigations the Nyquist Nyquist scientific in its prediction andemanates control (Dowling challenging task. In earlier investigations the 1997; et Selimefendigil et is used this problem for aa given system Olgac et advances al., 2015a). Its2000; complexity from2013; the Criterion (Dowling, 1997; Dowling and Morgans, 2005; OlgacNyquist et al., 1997; Annaswamy Annaswamy et al., al., 2000; Selimefendigil et al., al., 2013; Criterion is used to tackle this problem for given system 1997; Annaswamy et al., 2000; Selimefendigil et al., 2013; Criterion is used to tackle this problem for a given system 1997; Annaswamy et al., 2000; Selimefendigil et al., 2013; Criterion is used to tackle this problem for a given system Olgac et al., 2015a). Its complexity emanates from the (Dowling, 1997; Dowling and Morgans, 2005; Olgac et al., dynamic coupling between acoustic and thermal events in an 2015b). It is, however, not practical if one aims to establish Olgac et et al., al., 2015a). 2015a). Its Its complexity complexity emanates emanates from from the the (Dowling, 1997; Dowling and Morgans, 2005; Olgac et Olgac (Dowling, 1997; Dowling and Morgans, 2005; Olgac et al., al., Olgac etcoupling al., 2015a). Its acoustic complexity emanates fromin the (Dowling, 1997; Dowling and Morgans, 2005; Olgac et al., dynamic between and thermal events an 2015b). It is, however, not practical if one aims to establish enclosure (a.k.a. the combustor). A regionally-confined stability maps in the space of the system parameters. dynamic coupling coupling between between acoustic acoustic and and thermal thermal events events in in an an 2015b). It is, however, not practical if one aims to establish dynamic 2015b). It is, however, not practical if one aims to establish dynamic coupling between acoustic and thermal events in an 2015b). It is, however, not practical if one aims to establish enclosure (a.k.a. the combustor). regionally-confined maps in the space of the system parameters. unsteady release a series ofA waves within stability enclosureheat (a.k.a. thedrives combustor). Aacoustic regionally-confined stability maps in space system parameters. enclosure (a.k.a. the combustor). A regionally-confined stability mapsand in the the space of of the the system has parameters. The systems mathematics literature a strong evidence enclosure (a.k.a. thedrives combustor). Aacoustic regionally-confined stability maps in the space of the system parameters. unsteady heat release a series of waves within the combustor. These waves return and influence their own unsteady heat heat release release drives drives aa series series of of acoustic acoustic waves waves within within The systems and mathematics literature has aa strong unsteady of advances in the time delay systems (TDS) in the evidence last two unsteady heat release drives a series of acoustic waves within The systems and mathematics literature has evidence the waves influence their source (the heatThese release) afterreturn some and acoustic reflections and The systems and mathematics literature has aa strong strong evidence the combustor. combustor. These waves return and influence their own own The systems and mathematics literature has strong evidence of advances in the time delay systems (TDS) in the last two the combustor. These waves return and influence their own decades (Richard, 2003). Many numerical stability analysis the combustor. These waves return and influence their own of advances in the time delay systems (TDS) in the last two source (the heat after some acoustic and time delays 1878). mechanism of advances in the time delay systems (TDS) in the last two source (the (Rayleigh, heat release) release) afterSuch somea regenerative acoustic reflections reflections and decades of advances in the time delay systems (TDS) in the last two (Richard, 2003). Many numerical stability analysis source (the heat release) after some acoustic reflections and tools are developed (Engelborghs et al., 2002; Breda et al., source (the heat release) after some acoustic reflections and decades (Richard, 2003). Many numerical stability analysis time delays (Rayleigh, 1878). Such a regenerative mechanism may ultimately cause instability, which is reflected in the decades (Richard, 2003). Many numerical stability analysis time delays delays (Rayleigh, (Rayleigh, 1878). 1878). Such Such aa regenerative regenerative mechanism mechanism tools decades (Richard, 2003). Many numerical stability analysis are developed (Engelborghs et al., 2002; Breda et al., time 2005; Vyhlidal and Zitek, 2009), which approximate the time delays (Rayleigh, 1878). Such a regenerative mechanism tools are developed (Engelborghs et al., 2002; Breda et al., may ultimately instability, which reflected in form pressure cause fluctuations with growing tools developed (Engelborghs et al., 2002; Breda et al., may of ultimately cause instability, which is isamplitudes reflected within in the the characteristic tools are are developed (Engelborghs et al., 2002; Breda et al., 2005; Vyhlidal and Zitek, 2009), which approximate the may ultimately cause instability, which is reflected in the root distribution of TDS for given delay may ultimately cause instability, which is reflected in the 2005; Vyhlidal and Zitek, 2009), which approximate the form of pressure fluctuations with growing amplitudes within the combustion chamber. Very importantly, this regenerative 2005; Vyhlidal and Zitek, 2009), which approximate the form of of pressure pressure fluctuations fluctuations with with growing growing amplitudes amplitudes within within characteristic 2005; Vyhlidal and Zitek, 2009), which approximate the root distribution of TDS for given delay form compositions. Since they determine stability in given a pointwise form of pressure fluctuations withofgrowing amplitudes within root distribution of TDS for delay the Very this regenerative mechanism is anchamber. inherent feature the overall and characteristic characteristic root distribution of TDS for given delay the combustion combustion chamber. Very importantly, importantly, thisdynamics regenerative characteristic root distribution of stability TDS for given delay compositions. Since they determine in aa stability pointwise the combustion chamber. Very importantly, this regenerative fashion, a high computational demand arises when is the combustion chamber. Very importantly, this regenerative compositions. Since they determine stability in pointwise mechanism is an inherent feature of the overall dynamics and is not induced by an exogenous dynamic enforcement. As we compositions. Since they determine stability in a pointwise mechanism is is an an inherent inherent feature feature of of the the overall overall dynamics dynamics and and fashion, compositions. Since they determine stability in a pointwise a high computational demand arises when stability is mechanism of concern in a broader parametric space (such as delays). mechanism is an inherent feature of the overall dynamics and fashion, a high computational demand arises when stability is is induced by enforcement. As we will theexogenous unsteady dynamic heat release simply acts fashion, aa high computational demand arises when stability is is not notexplain inducedlater, by an an exogenous dynamic enforcement. Aslike we of fashion, high computational demand arises when stability is concern in a broader parametric space (such as delays). is not induced by an exogenous dynamic enforcement. As we Lyapunov–Krasovskii approaches are also developed for the is not induced by an exogenous dynamic enforcement. As we of concern in a broader parametric space (such as delays). will explain later, the unsteady heat release simply acts like an amplifier on this regenerative effect. of concern in a broader parametric space (such as delays). will explain explain later, later, the the unsteady unsteady heat heat release release simply simply acts acts like like stability of concern in a broader parametric space (such as delays). Lyapunov–Krasovskii approaches are also developed for the will determination of TDS (Gu and Niculescu, 2006). will explain on later, the unsteady heat release simply acts like Lyapunov–Krasovskii an effect. approaches are are also also developed developed for for the the Lyapunov–Krasovskii approaches an amplifier amplifier on this this regenerative regenerative effect. Lyapunov–Krasovskii approaches are also developed for the stability determination of TDS (Gu and Niculescu, 2006). an on effect. To better understand the TAI phenomenon, many researchers However, they provide only conservative results and are an amplifier amplifier on this this regenerative regenerative effect. stability determination of TDS (Gu and Niculescu, 2006). stability determination of TDS (Gu and Niculescu, 2006). stability determination TDSselection (Gu andofNiculescu, 2006). To better understand the TAI phenomenon, many researchers they provide conservative results and are focus on much simpler thermoacoustic devices, one of which However, strongly dependent onofonly the the LyapunovTo better understand the TAI phenomenon, many researchers However, To better understand the TAI many researchers However, they they provide provide only only conservative conservative results results and and are are To better understand thethermoacoustic TAI phenomenon, phenomenon, many researchers However, they provide only conservative results and are focus on much simpler devices, one of which strongly dependent on the selection of the Lyapunovis the Rijke tube (Rijke, 1859). It was first discovered in 1859 Krasovskii functionals. Alternatively a the recent Cluster focus on much simpler thermoacoustic devices, one of which strongly dependent on the selection of Lyapunovfocus on much simpler thermoacoustic devices, one of which strongly dependent on the selection of the Lyapunovfocus on occupies much simpler devices, one ofinwhich strongly dependent on the selection of Lyapunovis the Rijke tube (Rijke, 1859). It was first discovered 1859 Cluster Alternatively aa the recent and some thermoacoustic prestigious combustion research labs Krasovskii Treatment offunctionals. Characteristic Roots (CTCR) paradigm is is thestill Rijke tube (Rijke, 1859). It was first discovered in 1859 Krasovskii functionals. Alternatively Cluster is Rijke tube (Rijke, 1859). It first in Krasovskii functionals. Alternatively a recent recent Cluster is the the Rijke tube (Rijke, 1859). It was was first discovered discovered in 1859 1859 Krasovskii functionals. Alternatively a recent Cluster and still occupies some prestigious combustion research labs of Characteristic Roots (CTCR) paradigm is Treatment in its purest form (Matveev 2003; Dowling and Morgans, deployed here to remedy such Roots shortfalls (Olgacparadigm and Sipahi, and still occupies some prestigious combustion research labs Treatment of Characteristic (CTCR) is and occupies some prestigious combustion research labs Treatment of to Characteristic Roots (CTCR) paradigm is anditsstill still occupies some prestigious combustion research labsa Treatment of Characteristic Roots (CTCR) paradigm is in purest form (Matveev 2003; Dowling and Morgans, deployed here remedy such shortfalls (Olgac and Sipahi, 2005). The device simply consists of a cylindrical tube with 2005; Sipahi and Olgac, 2006). CTCR creates exhaustive and in its purest form (Matveev 2003; Dowling and deployed here to remedy such shortfalls (Olgac and in purest form (Matveev 2003; Dowling and Morgans, Morgans, deployed hereand to Olgac, remedy2006). such CTCR shortfalls (Olgac and Sipahi, Sipahi, in its itssource purest formsimply (Matveev 2003; Morgans, deployed here to remedy such shortfalls (Olgac and Sipahi, 2005). The device consists of cylindrical with aa non-conservative 2005; Sipahi Sipahi creates exhaustive and heat placed inside, as shown inaaDowling Fig. 1. andtube “stability maps” increates the domain of and the 2005). The device simply consists of cylindrical tube with 2005; and Olgac, 2006). CTCR exhaustive 2005). The simply consists of cylindrical tube 2005; Sipahi and and Olgac, Olgac, 2006). CTCRincreates creates exhaustive and 2005). The device device simply consists ofinaaFig. cylindrical tube with with aa system 2005; Sipahi 2006). CTCR exhaustive and heat source placed inside, as shown 1. non-conservative “stability maps” the domain of the parameters.“stability maps” in the domain of the heat source placed inside, as shown in Fig. 1. non-conservative heat placed inside, Fig. It commonly in the in community non-conservative non-conservative “stability maps” maps” in in the the domain domain of of the the heatissource source placed accepted inside, as as shown shown incombustion Fig. 1. 1. “stability system parameters. system parameters. It accepted community main contributions of this paper are the reformulation of (Dowling, 1997; Selimefendigil et al.,combustion 2013) that the thermo- The system parameters. It is is commonly commonly accepted in in the the combustion community system parameters. It is commonly accepted in the combustion community It is commonly accepted in the combustion community The main contributions of this paper systems are the reformulation reformulation of (Dowling, 1997; Selimefendigil et al., 2013) that the thermothe Rijke tube model from a control perspective and acoustic dynamics remains within linear realm until main contributions of paper are of (Dowling, 1997; 1997; Selimefendigil Selimefendigil et et al., al., 2013) 2013) that that the the thermothermo- The The main tube contributions of this this paper systems are the the reformulation reformulation of (Dowling, The main contributions of this paper are the of (Dowling, 1997; Selimefendigil et al., unstable, 2013) thathowever, the thermothe Rijke model from a control perspective and acoustic dynamics remains within linear realm until implementation of the CTCR paradigm to predict detailed instabilities set in. When it becomes the the Rijke tube model from a control systems perspective and acoustic dynamics dynamics remains remains within within linear linear realm realm until until the Rijke tube model from a control systems perspective and acoustic the tubethemodel fromstability control systems perspective and acoustic dynamics remains within linear realm until detailed implementation of the CTCR paradigm to predict TAI.Rijke First, unique maps of the system are instabilities set it becomes the pressure fluctuations grow (whichhowever, is a typical implementation ofa the the CTCR paradigm to predict predict instabilities set in. in. When When itexponentially becomes unstable, unstable, however, the detailed detailed implementation of CTCR paradigm to instabilities set in. When it becomes unstable, however, the implementation of the CTCR paradigm to predict detailed instabilities set in. When it becomes unstable, however, the TAI. First, the unique stability maps of the system are created. Then a series of experiments are performed to pressure fluctuations grow exponentially (which is a typical linear system property) to the levels which invite TAI. First, the unique stability maps of the system are pressure fluctuations fluctuations grow grow exponentially exponentially (which (which is is aa typical typical TAI. TAI. First, thea unique unique stability maps of ofarethe theperformed system are are pressure First, the stability maps system created. Then series of experiments to pressure fluctuations grow exponentially (which is a typical validate these maps. linear to levels created. Then Then aa series series of of experiments experiments are are performed performed to to nonlinearities they eventually into which a limit invite cycle created. linear system systemandproperty) property) to the thesettle levels which invite linear property) to the levels which invite created. Thenmaps. a series of experiments are performed to validate these linear system system property) to2000). thesettle levels which invite nonlinearities and they eventually into a limit cycle validate these maps. behaviour (Annaswamy et al., nonlinearities and and they they eventually settle settle into into a limit cycle cycle validate maps. of the following segments: Section 2 The textthese consists nonlinearities nonlinearities and they eteventually eventually settle into aa limit limit cycle validate these maps. behaviour (Annaswamy al., 2000). The text consists of following 2 behaviour (Annaswamy et al., 2000). contains the mathematical of the segments: Rijke tube Section dynamics. The text consists of the themodel following segments: Section behaviour (Annaswamy et al., 2000). In this paper, to detect the “onset” of instability, we take the The text consists of behaviour (Annaswamy et al., 2000). following segments: Section 222 The texttheconsists of the themodel following segments: Section contains mathematical of the Rijke tube dynamics. In this paper, to detect the “onset” of instability, we take the contains the mathematical model of the Rijke tube dynamics. In this paper, to detect the “onset” of instability, we take the contains the mathematical model of the Rijke tube dynamics. In to detect In this this paper, paper, to 2015 detect the the “onset” “onset” of of instability, instability, we we take take the the 75 contains the mathematical model of the Rijke tube dynamics. Copyright © IFAC 2405-8963 © 2015, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Copyright IFAC 2015 75 Copyright © IFAC responsibility 2015 75 Peer review© under of International Federation of Automatic Control. Copyright © 75 Copyright © IFAC IFAC 2015 2015 75 10.1016/j.ifacol.2015.09.356

IFAC TDS 2015 76 June 28-30, 2015. Ann Arbor, MI, USA

Umut Zalluhoglu et al. / IFAC-PapersOnLine 48-12 (2015) 075–080

waves f i (t ) and g i (t ) . The subscript i  1,2,3,4 denotes the cross-sections in Fig. 1, where these functions are evaluated. The time delays  u  xu / c ,  d  xd / c correspond to the travel times of acoustic waves moving at the speed of sound. In order to derive causal relationships between each variable,

Section 3 reviews the highlights of the CTCR paradigm. Section 4 is devoted to the experimental validation of theoretical findings. 2. THE MATHEMATICAL MODEL The mathematical model described in this section carries the underpinning physics of the Rijke tube under commonly applied assumptions (Dowling, 1997, Dowling and Morgans, 2005, Olgac et al., 2014):

we start with expressing

 f 2 (t )   f1 (t   u )        g 3 (t )   g 4 (t   d )  f1

e  u s

xd

g1

f3

f2

e  d s

Heat Release Zone

Ru

x Downstream

g 3 (t )T as the delayed

version of  f1 (t ) g 4 (t )T as

(i) the airflow is induced by natural buoyancy, therefore has low velocity and negligibly small Mach number; (ii) the heating zone is narrow compared to the tube length; (iii) the average (mean) flow quantities such as density  and the speed of sound c are assumed as constant along the tube; (iv) the acoustic wave propagation is taken as 1-D event;

④

 f 2 (t )

e  u s

g2

(3)

f4

Rd g3

e  d s

g4

f g Fig. 2. Block diagram representation of Rijke tube dynamics In Laplace domain, (3) transforms into

③ Heating Zone

 e  u s  F2 ( s)   F ( s)     T( s) 1 , T( s)    0  G4 ( s)   G3 ( s)  

0 Upstream

②

①

(4)

Next, we identify the dynamics representing the thermoacoustic coupling. The conservation of mass, momentum and energy equations, when used with (1) and (2) across the heating zone results in

-xu

Fig. 1. Schematic representation of Rijke tube dynamics

 1  f (t )   ~0  1  g (t )   1  1  1 1  2    Q(t )  (5) 1  2    1       f 3 (t )       1   1  g 3 (t )   Ac     1   1 

Differently in this paper we reformulate the model from a control systems perspective. The dynamics of the Rijke tube is governed by the first principles of conservation of mass, momentum and energy. When linearized under the listed assumptions, the pressure and velocity fluctuations in the tube obey the linear wave equation. They can be expressed as

~ p( x, t )  f (t  x / c )  g (t  x / c ) ~ u ( x, t )  [ f (t  x / c )  g (t  x / c )] /  c

0  e   d s 

where  is the heat capacity ratio, A is the cross-sectional ~ area of the tube and Q (t ) is the unsteady heat release. For brevity, derivation of (5) is suppressed here, but interested reader may refer to Dowling (1997) for details. ~ The unsteady heat release Q (t ) is causally influenced by the velocity fluctuations at the heating zone u~ (t )

(1) (2)

where overscript ~ denotes the non-steady (i.e., fluctuating) components of the quantity  . f ( x, t ) and g ( x, t ) represent the acoustic waves traveling upwards and downwards in the tube respectively. As shown in Fig. 1, x is the position along the tube, where x  0 corresponds to location of the heater. At the tube ends x   xu and x  xd , where xu and xd denote the distances between ①-② and ③-④ (the upstream and downstream sides from the heating zone, which are represented by subscripts u and d), respectively.

2

(Selimefendigil et al., 2011). Following Gelbert et al., (2012), we take it as a first order transfer function in Laplace domain ~ Q(s) / u~2 (s)  a /(bs  1) (6) where a and b act as gain and time constant, affected by mean heat release ( Q ), mean flow velocity ( u 2 ) and the heater geometry. Evaluating (2) at cross-section ② and taking its Laplace transform gives

The block diagram of the overall dynamics is given in Fig. 2. In this representation, the system variables are the acoustic

u~2 (s)  F2 (s)  G2 (s)/ c 76

(7)

IFAC TDS 2015 June 28-30, 2015. Ann Arbor, MI, USA

Umut Zalluhoglu et al. / IFAC-PapersOnLine 48-12 (2015) 075–080

77

Substituting (6) and (7) in Laplace transform of (5), one obtains

 u and  d . To simplify the further operations, it is meaningful to introduce a delay vector 1, 2   (2 u , 2 d )

ˆ ( s) G2 ( s)   Y ˆ ( s) F2 ( s)  X  F ( s)   G (s)   3   3  ˆ ˆ where X( s) and Y( s) are

and express the characteristic equation as

1  ˆ ( s)   1 a / A c 2 X     1 bs  1 

(8)



 Rd a1    e

1 1     ˆ a / A c 2 1  Y ( s)   1    1  bs  1   1  

 G2 ( s)   F ( s)  ˆ ( s) 1 Y ˆ ( s)    Z( s) 2 , Z( s)  X F ( s ) G ( s )  3   3  Next, we express the relationship between g 2 (t )



 Ru a  1 e

 1s



2

 2 Ac 

(15)

It is a quasi-polynomial involving two independent delays (1, 2 ) , which are related to ( xu , xd ) by   2 x / c . They vary with the tube length and heater location. The rest of the parameters in (15) are taken as constants. 3. STABILITY ANALYSIS METHODOLOGY

(9)

The Rijke tube dynamics falls into the class of linear time invariant multiple time delay systems with two delays, as discussed in Sec. 2. Its asymptotic stability is completely determined by the roots of its characteristic polynomial in (15). In general form, it can be written as

f 3 (t )T

f 4 (t )T similar to (3) and (4) as

 G ( s)   g1 (t )   g 2 (t   u )   G1 (s)           T( s) 2   f 4 (t )   f 3 (t   d )   F4 (s)   F3 ( s) 

 2 s

 Ru Rd 2 Ac 2   a  a e  ( 1  2 ) s  a(1   )  0

1   1    1 

In order to get the transfer matrix relating the two vectors, (8) can be rewritten as

and g1 (t )



CE( s,1, 2 )  2 Abc 2  Ru Rd e ( 1  2 ) s  1 s

(10)

CEs, 1, 2  

n

m

l

 Pijk e( j k )s si  0 1

2

(16)

i 0 j 0 k 0

The boundary conditions at cross-sections ① and ④ can be characterized by the reflection coefficients Ru and Rd . Both have values close to -1 due to the atmospheric pressure condition at the tube ends. This leads to the following expression

R  F1 ( s)   G ( s)     R 1 , R   u  G4 ( s)   F4 ( s)   0

0   Rd 

where n is the order of the characteristic polynomial, m and l are the highest order of commensuracy of delays  1 and  2 respectively. Notice that the highest order term in (16) involves delays, thus it is representative of a neutral timedelay system. Next, we give two theorems that determine the stability of this class of systems.

(11)

Theorem 1: The stability of a neutral time delay system requires the dynamics governed by its associated difference equation to be stable. Therefore as a necessary condition, all the zeros of

Now that all the transfer matrices are obtained, (4), (9), (10) and (11) can be combined to close the loop in Fig. 2 as

 F1 (s)   F (s)     RT(s)Z(s)T(s) 1   G4 (s)   G4 ( s) 

m

D ( s,  1 ,  2 ) 

(12)

1

2

(17)

have to lie in the left half of the complex plane. This is known as the strong stability condition (Hale and Verduyn Lunel, 2002). (13) Theorem 2: The necessary and sufficient condition for the asymptotic stability of (16) is that all its infinitely many zeros should lie in the open left half plane (Niculescu, 2001, Hale and Verduyn Lunel, 2002). The proof is again suppressed here.

In (13), M(s) represents the overall system matrix, determinant of which gives the characteristic equation of the system. Utilizing the Sylvester’s determinant theorem (Akritas et al., 1996) the characteristic equation can also be written as

CE(s,. u , d )  det[I  RT2 (s)Z(s)]

 Pnjk e( j  k )s  0 j 0 k 0

Then the return difference can be obtained as follows

 F ( s)  M( s) 1   0, M( s)  I  RT( s)Z( s)T( s)  G4 ( s) 

l

The conditions of Theorems 1 and 2 can be evaluated using the CTCR paradigm as explained in Olgac and Sipahi (2005) and Olgac et al., (2008). Apparently, if the strong stability condition in Theorem 1 is not fulfilled, there is no need to follow the CTCR procedure. If it is satisfied however, delay dependent stability posture of the system can be determined by CTCR. Details on the paradigm are described comprehensively in Sipahi and Olgac (2006) and Fazelinia et

(14)

Notice that the exponential terms in (14) originate from

T 2 ( s) , which contains twice the acoustic travel delays

77

IFAC TDS 2015 78 June 28-30, 2015. Ann Arbor, MI, USA

Umut Zalluhoglu et al. / IFAC-PapersOnLine 48-12 (2015) 075–080

al., (2007). In order to convey some highlights and definitions, a brief review is presented next.

the imaginary root set S remains invariant from kernel to offspring. The set  constitutes the complete distribution of

For a linear time-invariant time delay system characterized by (16) to switch the stability posture, the characteristic roots should cross the imaginary axis, say at  i for some delays

(1, 2 ) where the characteristic equation (16) has root sets containing at least one pair of imaginary roots, S .

τ  2  . We use

The root tendency (RT) for an imaginary characteristic root s  S  , along a delay axis (say  1 ) is defined by

τ,  notation to indicate their causal

correspondence. In order to assess the stability of the system exhaustively, all delay compositions for which (16) has an imaginary root should be determined. Let’s denote the complete set of such imaginary root frequencies with Ω , and the corresponding root set with S as



Ω   CE( s   i, τ)  0, τ  2  ,    S  s   i   Ω 



RT

   s   Re  sgn  s i    1  s i 

1

(22)

This property indicates the evolution direction of the imaginary root  i (to the left or the right-half of the complex plane) as we pass through a hypercurve by increasing one delay (  1 ) infinitesimally.

(18)

It is trivial to observe from (16) that an imaginary root will be repeated

Proposition 2. At s  S  as  1 (or  2 ) increases at a point on a kernel and its corresponding offspring, the root tendency

2 2   (1 j , 2k )  1  j , 2  k , j  0,1, ..., k  0,1, ... (19)    

, (or RT s 2 ), remains invariant so long as the other delay  2 (or  1 ) is kept fixed. This feature, in essence, declares the stabilizing (or destabilizing) transitions across the regional boundaries defined by  .

s   i corresponding to infinitely many times as

τ  2 

RT

These trajectories divide the delay domain into encapsulated regions in which the number of unstable roots, NU, is fixed. Consequently, any stability reversal can only occur at the boundaries of these regions. The system is declared stable when NU  0 . One needs to determine the boundaries of these infinitely many regions exhaustively for a complete stability declaration. The following two propositions bring a discipline to this chaotic looking picture.



s

Based on these propositions of CTCR, one can establish the stability outlook of the system by generating its stability map in (1, 2 ) domain. To determine the complete set of kernel and offspring hypercurves, we use the Spectral Delay Space (SDS) method (Fazelinia et al., 2007; Gao et al., 2014). A new coordinate set, (1, 2 )  (1, 2) is defined in this approach, which confines the kernel and offspring hypercurves in 2-dimensional blocks of 2 length. The building block (BB) contains the building hypercurves, which later transform-to kernel hypercurves in (1, 2 ) domain. The blocks that contain the reflection hypercurves are then simply obtained by stacking the BB. They transform to the offspring hypercurves in (1, 2 ) domain. Once the complete set  is obtained, RT is evaluated and the stable/unstable regions on the stability map is revealed.

Proposition 1. There are only a manageably small number of trajectories in (1, 2 ) space called the kernel hypercurves (Sipahi and Olgac, 2006) 2 2   0   τ τ,  , τ  2  ,   Ω, 0  1  ,0  2   (20)   

Notice that 0 represents the loci of the smallest  1 and  2 combinations corresponding to each imaginary root frequency   Ω . All the remaining hypercurves are created from this set, 0 , utilizing the point-wise translation property in (19) for j, k  0 . They are called the offspring

4. MAIN RESULTS In this section, the stability of the Rijke tube is assessed using the CTCR paradigm and these analytical results are compared with several experimental tests. Three cylindrical glass tubes of different lengths are heated with an electrical coil. Its power is supplied through a variac. The temperatures below and above the heater are measured via two thermocouples. Sound pressure is measured at the downstream end with a piezoelectric microphone. The data acquisition of the measured pressure took place at 8000 Hz sampling rate.

hypercurves, and denoted by  jk where j and k identify the jth and kth generation offspring in  1 and  2 , respectively. Consequently, the complete set of kernel and offspring hypercurves,  , becomes

  0  j 1  k 1 jk

1

(21)

The operational parameters that represent our experiment are measured and taken as Ru  Rd  0.93 ,   1.4 ,

Any point on the trajectories of kernel, 0 , imposes its

  Ω signature identically onto its offspring,  jk . Thus,

A  7.5 cm 2 ,   1.2 kg/m3 , c  340 m/s . Benefiting from Gelbert et al., (2012), the heat release dynamics are taken as 78

IFAC TDS 2015 June 28-30, 2015. Ann Arbor, MI, USA

Umut Zalluhoglu et al. / IFAC-PapersOnLine 48-12 (2015) 075–080

a  200 and b  0.002 . Substituting these parameters in the characteristic equation given in (15), we obtain the following normalized form: CE( s, 1 , 2 )  (1  0.865e ( 1  2 ) s ) s  178.78 e  1 s  178.78 e  2 s  598.71e  ( 1  2 ) s  692.23  0

The agreement between experimental TAI detection overlaid on the stability chart of CTCR is remarkable, especially for the intervals A1B1, A2B2 and A3B3. These intervals correspond to the case when the heater is located at the lower half of the tubes ( xu  xd ) . In accordance with Matveev, (2003) and Raun et al., (1993), it is not possible to destabilize the Rijke tube when the heater is located at the midpoint ( xu  xd ) or when it is close to tube ends. Some minor disagreements are observed when the heater is located at the upper half of the tube ( xu  x d ) , such as C2D2 (for 1.02 m tube) and C3D3 (for 1.22 m) intervals. At such locations, the instability is driven by higher modes of the dynamics, instead of the fundamental (first) mode (Raun et al., 1993). Such discrepancies are attributed to the lack of nonlinear damping features that were unaccounted for in the mathematical model. For instance, Matveev, (2003) reports that the reflection of the acoustic waves at the tube ends are in fact weaker for the secondary (higher) modes. In addition, the pressure oscillations at higher frequencies tend to have lower signal to noise ratio, which are more difficult to record. Our experimental equipment might have failed to capture such sound levels.

(23)

First, the strong stability condition in Theorem 1 should be checked. The coefficient of the highest order s term in (23) corresponds to the associated difference equation:

D(s,1, 2 )  1  0.865e ( 1  2 ) s  0

79

(24)

which implies | e ( 1  2 ) s |  1.156  1 , meaning that all the infinitely many roots of (24) have Re(s) < 0. Thus, the strong stability condition is satisfied. Next, to generate the stability map of the Rijke tube in ( 1 , 2 ) domain, the CTCR procedure described in Sec. 2 is followed. The SDS representation of (23) with the respective building (red) and reflection (blue) hypercurves is obtained in Fig. 3. Again, they carry the kernel and offspring signature respectively.

Fig. 3. SDS representation of the system. Building block (BB) is in black solid frame.

Fig. 4. The CTCR generated stability map of the Rijke tube, compared with experimental results. Next, the modal characteristics of the Rijke tube are studied during an unstable operation. In Fig. 4 it is declared that TAI occurs for the 1.22 m length tube when the heater is located at ( xu , xd )  (0.20 m,1.02 m) . At this operating point, the recorded sound pressure is shown in Fig. 5a, where the instability and nonlinear limit cycle can be observed. As the zoomed-in figure in Fig. 5b shows, at the onset of instability the pressure oscillations grow exponentially. This is linear system behaviour, thus it can be claimed that the linearized mathematical model described in Sec. 2 is sufficient to predict TAI. A fast Fourier transformation (FFT) is performed at the onset of instability and power corresponding spectral power density is obtained as in Fig. 5c. In Fig. 5d, the approximate characteristic root locations of (23) are shown, generated with the Quasi-polynomial Mapping based Rootfinder (QPmR) algorithm developed by Vyhlidal and Zitek (2009). The imaginary parts of the unstable roots match well with the detected frequencies in Fig. 5c, supporting the accuracy of the mathematical model and stability analysis.

The CTCR generated stability map is given in Fig. 4. To facilitate visualization it is displayed in ( xu , xd ) space rather than delay space ( 1 , 2 ) . The kernel and offspring hypercurves are coloured in red and blue respectively. The stable regions are shaded in grey. Several experimental tests are conducted using three different tube lengths ( xu + x d = 0.51 m, 1.02 m and 1.22 m), which are represented by black parallel lines. Any point on these lines depicts a certain heater location determined by xu and xd . To test the structure at any heater location, we energize the heating element, and wait until the temperature difference below and above the heater reach steady state. If the recorded pressure perturbations attain a noticeable level, the system is declared unstable. The unstable declared operating points are marked with green, while the stable ones are marked with black circles.

79

IFAC TDS 2015 80 June 28-30, 2015. Ann Arbor, MI, USA

Umut Zalluhoglu et al. / IFAC-PapersOnLine 48-12 (2015) 075–080

Gu, K., Niculescu, S. I. (2006). Stability Analysis of Timedelay Systems: A Lyapunov Approach. In Advanced Topics in Control Systems Theory, Springer London, 4, 139-170. Fazelinia, H., Sipahi, R., Olgac, N. (2007). Stability Robustness Analysis of Multiple TimeDelayed Systems Using “Building Block” Concept, IEEE Trans. Automat. Control, 52, 799-810. Engelborghs, K., Luzyanina, T., Roose, D.. 2002. Numerical bifurcation analysis of delay differential equations using DDE-BIFTOOL. ACM Trans. Math. Software, 28, 1-21. Hale, J.K., Verduyn Lunel, S.M., 2002. Strong stabilization of neutral functional differential equations, IMA J. Math. Control Info. 19, 5-23. Matveev, K.I., (2003). Thermoacoustic instabilities in the Rijke tube: experiments and modeling, PhD Thesis, California Institute of Technology, Pasadena, CA. Niculescu, S. I., (2001). Delay effects on stability: A robust control approach, Springer. Olgac, N., Sipahi, R. (2005). The cluster treatment of characteristic roots and the neutral type time-delayed systems, ASME J. Dyn. Sys. Measur. Control, 127, 88-97. Olgac, N., Vyhlídal, T., Sipahi, R. (2008). A new perspective in the stability assessment of neutral systems with multiple and cross-talking delays, SIAM J. on Control and Optim., 47, 1, 327-344. Olgac, N., Zalluhoglu, U., Kammer, A.S. (2014). Predicting thermoacoustic instability; a novel analytical approach and its experimental validation, J. Propul. Power, 30, 4, 1005-1015. Olgac, N., Zalluhoglu, U., Kammer, A.S. (2015a). A new perspective in designing delayed feedback control for thermo-acoustic instabilities (TAI), Combust. Sci. Tech., 187, 5, 697-720. Olgac, N., Cepeda-Gomez, R., Zalluhoglu, U., Kammer, A.S. (2015b). Parametric investigation of thermoacoustic instability (TAI) in a rijke tube: a time-delay perspective, Int. J. Spray Combust. Dyn., 7, 1, 39-68, 2015. Raun, R.L., Beckstead, M.W., Finlinson, J.C., Brooks, K.P., (1993). A review of Rijke tubes, Rijke burners and related devices, Prog. in Energy and Comb. Sci., 19, 313-364. Rayleigh, J.W.S. (1878). The explanation of certain acoustic phenomena, Nature 18, 319-321. Richard, J.P. (2003). Time-delay systems: an overview of some recent advances and open problems. Automatica, 39, 1667-1694. Rijke, P.L., (1859). Notice of a new method of causing a new vibration of the air contained in a tube open at both ends, Philos. Mag. Series 4, 17, 419-422. Selimefendigil, F., Sujith, R.I., Polifke, W. (2011). Identification of heat transfer dynamics for non-modal analysis of thermoacoustic stability, Appl. Math. Comput., 217, 5134–5150. Sipahi, R., Olgac, N. (2006). A unique methodology for the stability robustness of multiple time delay systems, Sys. Control Lett., 55, 819-825. Subramanian, P., Mariappan, S., Sujith, R.I., Wahi, P. (2013). Subcritical bifurcation and bistability in thermoacoustic systems, J. of Fluid Mech., 715, 210-238. Vyhlídal, T., Zítek, P., (2009). Mapping based algorithm for large-scale computation of quasi-polynomial zeros, IEEE Trans. Automat. Control, 54, 171-177.

Fig. 5. (a) Time trace of the recorded sound pressure; (b) zoomed-in view of the onset of instability; (c) power spectral density of the interval in (b); (d) QPmRgenerated dominant poles; all for 1.22 m length tube with the heater located at ( xu , xd )  (0.20 m,1.02 m) . 6. CONCLUSIONS In this paper, thermo-acoustic instability phenomenon in a Rijke tube is studied. First, it is shown that the mathematical model of the Rijke tube falls under the linear time-invariant neutral type multiple time delay systems class. Then utilizing the cluster treatment of characteristic roots (CTCR) paradigm, the thermo-acoustically stable and unstable operating conditions are determined. The results are then verified on an experimental Rijke tube set-up. The authors plan to demonstrate a part of this experimental study during the 12th IFAC Workshop on Time Delay Systems as well. REFERENCES Akritas, A. G.; Akritas, E. K.; Malaschonok, G. I. (1996). Various proofs of Sylvester's (determinant) identity. Math. Comput. in Simul., 42, 4: 585-593. Annaswamy, A. M., Fleifil, M., Rumsey, J. W., Prasanth, R., Hathout, J. P., Ghoniem, A. F. (2000). Thermoacoustic instability: model-based optimal control designs and experimental validation. IEEE Trans. Automat. Control, 8, 6, 905-918. Breda, D., Maset, S., Vermiglio, R., (2005). Pseudospectral differencing methods for characteristic roots of delay differential equations. SIAM J. Sci. Comput., 27, 482-495. Dowling, A.P. (1997). Nonlinear self-excited oscillations of a ducted flame, J Fluid Mech., 346, 271–290. Dowling, A.P. and Morgans, A.S. (2005). Feedback control of combustion oscillations, Annu. Rev. Fluid Mech., 37, 151182. Gao, Q., Zalluhoglu, U., Olgac, N. (2014). Investigation of local stability transitions in the spectral delay space and delay space, ASME J. Dyn. Sys. Measur. Control, 136, 051011-1. Gelbert, G., Moeck, J.P., Paschereit, C.O., King R., (2012). Feedback control of unstable thermoacoustic modes in an annular Rijke tube, Control Eng. Prac. 20, 770-782.

80