Time domain modelling and stability analysis of complex thermoacoustic systems [original]

657

Time domain modelling and stability analysis of

complex thermoacoustic systems†

M R Bothien∗,J P Moeck, A Lacarelle,and C O Paschereit

Institute of Fluid Dynamics and Engineering Acoustics (Hermann-Föttinger-Institute), Technical University Berlin,

Berlin, Germany

The manuscript was received on 31 October 2006 and was accepted after revision for publication on 11 April 2007.

DOI: 10.1243/09576509JPE384

Abstract: A methodology allowing for a modular setup of complex acoustic systems is developed.

The transfer behaviour of the individual subsystems is formulated in time domain. Subsystem

descriptions can be obtained by analytical considerations, numerical methods, or experimental

data. Once the complex subsystems have been characterized experimentally, changes in sys-

tem geometry can be implemented easily by exchanging or adding subsystems. To validate the

modelling approach, experiments are conducted in an acoustic test rig with a combustor-type

geometry. Results are compared to predictions from the model, demonstrating accuracy in fre-

quency and time domain. Application to thermoacoustic instabilities arising in lean-premixed

combustion is given.The influence of a modified fuel distribution on an unstable operating point

of a lean-premixed combustor is studied and validated with experimental data. Additionally, a

study on the parameters governing the flame transfer function is performed to generate a stability

map of a model combustor. An advantage of the state-space approach is that stability of a ther-

moacoustic system can be determined by simply solving a matrix eigenvalue problem. This is in

strong contrast to the traditional approach, where the complete model is formulated in frequency

domain with infinite-dimensional transfer functions. The time domain approach is based on the

methodology presented by Schuermans et al. [1]. In contrast to their work, however, subsystems

are not obtained from modal expansions but are characterized by using system identification

techniques. Additionally, accuracy of the time domain model is verified by experiments.

Keywords: stability analysis, state-space system, thermoacoustics, time domain simulation

1 INTRODUCTION

As emission levels for modern gas turbines become

more and more restrictive, the gas turbine industry

stays abreast of changes by using lean-premixed com-

bustion. However, the leaner a combustion system

operates the more it is prone to suddenly occur-

ring large pressure oscillations. These so-called ther-

moacoustic instabilities arise from the interaction of

∗Corresponding author: Institute of Fluid Dynamics and Engineer-

ing Acoustics (Hermann-Föttinger-Institute), Technical University

Berlin, Müller-Breslau-Straße 8, 10623 Berlin, Germany. email:

mirko[email protected]

†Extended version of a paper originally presented at the 2006 Con-

ference on Modelling Fluid Flow (CMFF’06),Budapest University of

Technology and Economics, 6–9 Sepetember 2006.

unsteady heat release and the acoustic field in the

engine. If the two mechanisms constructively inter-

fere, high amplitude pressure pulsations occur, which

have a detrimental effect on the combustion pro-

cess [2]. Among others, Dowling and Hubbard [3]

expressed this fact in terms of an acoustic energy

equation, which states that perturbations grow if their

net energy gain from the combustion–acoustic inter-

action is greater than their losses across the system

boundaries and due to dissipation.

Thermoacoustic instabilities are a major issue in gas

turbines – both in stationary and aerospace applica-

tions but also in boilers and furnaces [2,4,5]. Strong

pressure oscillations cause higher emissions of pollu-

tants (NOxand CO), decrease the engine performance

and significantly increase the noise level. If the fluctu-

ations are high enough, they can even cause structural

damage [2,6,7].

JPE384 ©IMechE 2007 Proc. IMechE Vol. 221 Part A: J. Power and Energy

658 M R Bothien,J P Moeck,A Lacarelle,andCOPaschereit

To predict unreliable and unsafe operation, it is

necessary to describe how the system’s acoustic field

couples with the heat release. A model of the acous-

tic properties, supplemented by information on the

flame response, allows for prediction and control

of combustion instabilities. Thus, maintenance rates

can be extended and the operating process can be

improved. Furthermore, a model can also be used to

formulate design constraints to new burner develop-

ments.

As full-scale engine tests are extremely costly and

time-consuming, detailed investigations of thermoa-

coustic phenomena are commonly conducted in sin-

gle burner test rigs. Due to the high complexity of

geometry and interaction mechanisms in combustion

systems,a comprehensivedescriptionoftheiracoustic

field is cumbersome [1]. Although detailed numer-

ical simulations (taking into account reaction, flow,

and acoustics) were shown to be able to accurately

calculate the thermoacoustic interaction mechanisms

in single burner configurations [8], time expense,

and computing power requirements are very high.

Schuermans et al. [9] decreased the computational

effort by coupling low-order acoustic models, repre-

senting complex impedance boundary conditions, to a

computational fluid dynamics (CFD) solver. Especially

if stability has to be assessed for a high number of

parameter combinations (as shown in section 5), the

use of a CFD code is not practicable. Moreover, as a

CFD simulation is only able to observe the dominant

mode of the instability, it is usually not possible to

assess more linearly unstable modes [10]. Additionally,

marginally stable modes cannot be detected, which

is in contrast to the network approach described in

the following. Sattelmayer and Polifke [11,12] stated

that it is extremely difficult to derive suitable mea-

sures to improve the system stability using results

of large eddy simulation (LES) calculations, although

these are able to give detailed information on the

periodic reacting flow. Therefore, methods based on

analytical modelling will remain the most appropri-

ate stability assessment tool in the relatively near

future.

One further step to reduce computational effort is

to use low-order models to describe the whole com-

bustion system. The common approach is to divide

the total system into several subsystems, which are

described by transfer functions accounting for plane

wave propagation. Connection of these subsystems

results in a low-order representation of the total

system’s acoustic properties. The network approach

has the main advantages that information about

the system’s acoustics and insight into the phys-

ical processes can be obtained very quickly [13],

which is of major importance in early design phases.

Furthermore, changes in geometry or setup, as

for example the performance of different burners,

can easily be implemented in the model by sim-

ply exchanging a subsystem. Due to this modular

setup, it is possible to incorporate an experimen-

tally determined flame transfer function in the model.

The flame–acoustic interaction is crucial for the sta-

bility of thermoacoustic systems and it is far from

trivial to represent the linear flame response accu-

rately in a CFD computation. A drawback of this

approach is the limitation to plane wave acoustics.

To overcome this problem, multi-dimensional mod-

els have to be considered. Evesque and Polifke [13]

set up a two-dimensional low-order network, which

was validated against a three-dimensional finite ele-

ment method (FEM) solver showing good agree-

ment.

Traditionally, stability analysis of thermoacoustic

systems is accomplished by setting up a low-order

network model, which is then investigated in the fre-

quency domain. This is either done by solving the

system’s dispersion relation [14,15], or by graphical

methods (Nyquist plots) as described by Sattelmayer

and Polifke [11,12]. Both the graphical approach and

finding the roots of the dispersion relation are gen-

erally not trivial but rather time-consuming tasks,

especially if stability maps have to be generated [16].

In contrast to the traditional frequency domain

approach, another procedure is followed here, by

modelling the acoustic properties of the system in the

time domain. Therefore, all subsystems are character-

ized in state-space form. For the resulting state-space

network system, only a matrix eigenvalue problem has

to be solved to perform a stability analysis, which is

numerically straightforward.Schuermansetal. [7] also

described the system in state-space form, however,

they used modal expansion techniques to describe the

individual subsystems.

Once a model for the thermoacoustic proper-

ties of the combustion system is at hand, control

strategies can be implemented to suppress the feed-

back cycle. These strategies can either be passive by

implementing changes in the combustor geometry or

introduction of damping devices (e.g. Helmholtz or

λ/4-resonators) [17,18], or active by means of open-

or closed-loop control [3,19].

The remaining parts are structured as follows: sec-

tions 2 and 3 deal with the characterization of the sub-

systems in state-space form. The modelling approach

is validated against experimental data in section 4.

Application to thermoacoustic stability analysis is

then given in section 5. Here, results from frequency

domain and state-space calculations are also com-

pared. In section 6, a simulation model is used to

predict the stabilizing influence of a modified fuel

distribution in an experimental lean-premixed com-

bustor.

Proc. IMechE Vol. 221 Part A: J. Power and Energy JPE384 ©IMechE 2007

Modelling and analysis of complex thermoacoustic systems 659

2 STATE-SPACE REPRESENTATION OF ACOUSTIC

SUBSYSTEMS

As in the frequency domain transfer matrix approach

[20], it is assumed that only plane waves propagate

between two subsystems. This assumption is valid for

elements that are coupled via ducts, for which the fre-

quencies considered are below the cut-on frequency of

the first non-planar mode. Neglecting entropy waves,

this results in two degrees of freedom for all coupling

planes. The variables accounting for that can be either

acoustic pressure pand velocity vor the Riemann

invariants of the plane wave field, denoted here with f

and gfor downstream and upstream travelling waves,

respectively. The primitive variables are related to

the Riemann invariants by p=f+gand v=f−g.

Here (as in the following), the acoustic pressure has

been scaled with the characteristic impedance ρc, the

product of density and speed of sound.

Subsequently, acoustic systems are characterized

using the scattering form, i.e. the f-wave upstream

and the g-wave downstream are treated as inputs

and, accordingly, the f-wave downstream and the

g-wave upstream are obtained as outputs.The generic

acoustic element in state-space form can then be

written as

x=Ax+Bfu

gd

fd

gu=Cx+Dfu

gd

(1)

where subscripts u and d denote upstream and down-

stream positions, respectively. A,B,C, and Dare

time-invariant N×N,N×2, 2 ×N, and 2 ×2 matri-

ces, respectively, and xis the N-dimensional state

vector. The dimension of the state depends on the

acoustic subsystem to be characterized and has to

be chosen sufficiently large so that the system’s fre-

quency response is represented with acceptable accu-

racy in the frequency range considered. An inlet or exit

boundary condition can be written similar to equation

(1) as a single input single output (SISO) system.

3 SUBSYSTEM CHARACTERIZATION

The modular setup allows for different ways of describ-

ing the subsystems, i.e. by analytical considerations,

experimental data, or numerical methods (e.g. using

FEM [21]). The subsystems’ transfer functions relate

the incoming and outgoing acoustic waves. Usually,

acoustic subsystem characterization based on FEM

computations or experiments results in frequency

domain data. In the frequency domain the rela-

tion of the Riemann invariants is represented by the

scattering matrix S

fd

gu=S11 S12

S21 S22fu

gd(2)

where the elements Sij are complex-valued functions

of frequency.

Once a subsystem’s acoustic transfer function is

given in the frequency-dependent form, equation (2),

it can easily be transformed into the corresponding

state-space form (equation (1)). The elements of the

matrices of this state-space realization depend on how

the system’s state-space variables are defined. The

same input-output behaviour can be obtained from

different definitions of the system’s state variables. If

these variables are determined from the poles of a

partial fraction expansion of the elements Sij, which

is mostly the case in this article, the system matrix A

(equation (1)) assumes the so-called Jordan canonical

form. The diagonal elements then comprise the poles

of the partial fraction expansion. Detailed informa-

tion on how this and other state-space formulations

can be obtained from a transfer function in frequency

domain can be found in reference [22]. In the follow-

ing sections 4 and 5, more information on the specific

form ofthetransfer functionofeachindividual subsys-

tem is given. Here, either the scattering matrix, transfer

function, or state-space form is chosen, depending on

which formulation is more descriptive.

3.1 Analytical description

The description of elements with simple geome-

tries (e.g. ducts) is straightforward. Neglecting sound

absorption at the walls and in the fluid, a duct only

represents a time-delay (e−iωτ ) in the system, i.e. an

acoustic wave is transmitted with a delay, depending

on the elements length and the speed of sound.

It is not possible to exactly describe a system with

time-delay with a rational transfer function.Therefore,

it has to be approximated, e.g. with a Padé approxi-

mation of the exponential function. The order of the

Padé approximation, which is necessary for a suf-

ficient accuracy, depends on the desired frequency

range (100 Hz–1000 Hz in the case considered here)

and the time-delay. A duct with a length of 0.5 m at

ambient conditions requires at least an order of 8 in

the frequency band mentioned above.

3.2 Experimental description

Transfer functions of more complex subsystems can

be determined by means of experiments. For this pur-

pose, the element is regarded as a black box, for which

the four unknown elements of the scattering matrix

(equation (2)) have to be determined. This procedure

JPE384 ©IMechE 2007 Proc. IMechE Vol. 221 Part A: J. Power and Energy

660 M R Bothien,J P Moeck,A Lacarelle,andCOPaschereit

has been applied previously [23–25] and thus only a

short description is presented here.

The system given by equation (2) consists of two

equations. Thus, two linearly independent states are

necessary to solve for the four unknowns Sij. The

test states are created by successive acoustic exci-

tation on both sides of the element. To obtain the

Riemann invariants from sound pressure measure-

ments at least two microphones at two different axial

locations in a duct are necessary. In the presence of

measurement uncertainty, better results are obtained

by acquiring the pressure at more (n) positions. Hav-

ing data of more than two microphones, the problem

is over-determined. This approach is called the multi-

microphone method (MMM) [24]. The up- and down-

stream travelling waves are fitted to the measured

quantities by means of the least squares method

f

g=Z+⎡

⎢

⎣

⎤

⎥

⎦(3)

where Zis given by

Z=⎡

⎢

⎣

e−iωx1

c+u0eiωx1

c−u0

e−iωxn

c+u0eiωxn

c−u0

⎤

⎥

⎦

(4)

Here, xiis the distance of microphone ito the

reference plane, i.e. the location where fand gare

determined. Z+denotes the pseudo-inverse of Z

Z+=(ZTZ)−1ZT(5)

where the superscripts −1 and T denote the inverse

and the conjugate transpose of a matrix, respectively.

4 EXPERIMENTAL VALIDATION OF THE

STATE-SPACE APPROACH

4.1 Experimental setup

The network approach based on state-space models

described in the preceding section was validated in an

acoustic test rig. Figure 1 shows a schematic overview

of the experimental facility.

The model for the test rig is split into subsystems,

whose scattering matrices were determined exper-

imentally (burner, up-, and downstream end) and

analytically (ducts). Interconnection of the subsys-

tems yields a model for the test rig, which can be used

to predict its acoustic field.

The up- and downstream ends are equipped with

circumferentially mounted loudspeakers, which allow

for acoustic excitation of the rig. With these loud-

speakers, it is possible to adjust the linearly indepen-

dent test states that are necessary to experimentally

obtain a system’s scattering matrix. For the determi-

nation of the acoustic field, the test rig’s duct sections

are instrumented with 1/4 condenser microphones.

There are three microphones upstream of the burner

element and four downstream, thus allowing for the

application of the MMM.

4.2 Setup of the simulation model

4.2.1 Loudspeaker elements

The speaker elements’ transfer functions were deter-

mined experimentally. The elements have two inputs

and one output. In the case of the upstream element,

the g-wave and the loudspeaker voltage u1are the

inputs while the f-wave represents the output. In the

frequency domain this can be written as

f=Rg +Gu1(6)

where Ris the reflection coefficient and Gthe transfer

function of the loudspeaker’s controlling input volt-

age to the outgoing f-wave. Equation (6) implies the

assumption that the contribution of the reflection

coefficient and the loudspeaker transfer function can

be linearly superimposed. Results shown in section 4.3

confirm the validity of this assumption.

Similar to the approach described previously, two

excitation states have to be generated. This is achieved

using loudspeaker excitation at each side separately.

Downstream excitation (index A) yields

fA=RgA(7)

from equation (6) as u1is zero. This allows for direct

calculation of the reflection coefficient R. Excitation

Fig. 1 Schematic overview of the acoustic test rig splitted up into its subsystems

Proc. IMechE Vol. 221 Part A: J. Power and Energy JPE384 ©IMechE 2007

Modelling and analysis of complex thermoacoustic systems 661

Fig. 2 Loudspeaker transfer function (G) and reflection

coefficient (R) for measured values (black), model

with 30 states (grey), and model with 15 states

(grey dashed). Top: absolute value, middle: phase

of G, bottom: phase of R

with the upstream loudspeaker (index B) and solving

for Ggives

G=fB−RgB

(8)

where Ris calculated from equation (7).

With complex values for Rand Gcalculated at

discrete frequencies, frequency domain system iden-

tification algorithms [26,27] were used to obtain

state-space models. As stated before, the frequency-

dependent transfer functions are approximated by

partial fraction expansions. The number of poles used

for this expansion determines the number of states

of the state-space system. Figure 2 depicts gain and

phase of Rand G. Results for measured data as

well as for models identified with 30 states and 15

states are shown.

As can be seen in Fig. 2, measured and computed

curves for gain and phase of the reflection coefficient R

modelled with 30 states match almost perfectly. In the

frequency regime above 700 Hz, the model obtained

with 15 states shows minor deviations. In the case of

the loudspeaker’s transfer function G, excellent results

are obtained with 30 states. With 15 states, the results

are fair except for some frequency bands, in which G

has several local extremal values. Regarding the phase,

both models agree well with the measured values. The

larger slope of R, i.e. the larger time-delay, is due to

the position of the loudspeaker. The incident wave

travels to the upstream element’s end and back to

the reference plane, whereas the emitted wave of the

loudspeaker has to travel a much shorter distance.

Up to this point, only two SISO state-space systems

exist. These have to be concatenated to a double input

single output state-space system, which completely

describes the upstream end

x=AR0

0AGx+BR0

0BGg

u1

f=[CRCG]x+[DRDG]g

u1(9)

Indices Rand Gdenote the matrices of the reflec-

tion coefficient’s and the loudspeaker’s state-space

systems, respectively.

4.2.2 Artificial output element

The grey highlighted element in Fig. 1, is not part of

the real system but only an artificial component of the

model system. Through this, it is possible to access

the signal of any microphone by placing the artificial

output element (AOE) at the position of the desired

microphone. The element is modelled in state-space

as a static gain without additional dynamics

⎡

⎣

⎤

⎦=⎡

⎣

⎤

⎦fu

gd(10)

Besides the Dmatrix, all other matrices of the state-

space system (A,B, and C) are empty. The output

element has a physical length of zero. The up- and

downstream travelling waves are neither delayed nor

changed in amplitude and therefore, the output ele-

ment does not influence the acoustic behaviour of the

total system.

4.2.3 Interconnected complete system

The outputs of each subsystem have to be con-

nected with the inputs of its adjacent subsystems.

Linking together the subsystems’ state-space repre-

sentations results in one single state-space model for

the complete system. The interconnection is exem-

plarily shown for two adjacent subsystems of the

acoustic test rig.

The transfer behaviour of a subsystem defines the

relation betweenincidentandoutgoingwaves.Bycon-

necting the two subsystems‘duct’ and‘burner’ (Fig. 3),

one system is obtained with the Riemann invariants f1

and g3as inputs and g1and f3as outputs. Thus, the

two systems of equations describing the individual

subsystems are now combined to one, representing

the transfer behaviour of the connected system. This

combination is a matrix operation referred to as the

JPE384 ©IMechE 2007 Proc. IMechE Vol. 221 Part A: J. Power and Energy

662 M R Bothien,J P Moeck,A Lacarelle,andCOPaschereit

Fig. 3 Interconnection of two subsystems by means of

the Redheffer Star Product

Redheffer Star Product [7,22]. Applied to all subsys-

temsofFig.1,this methodyieldsone singlestate-space

expression of the total system. Interconnection of the

individual state-space systems can conveniently be

done with built-in MATLAB routines.

Interconnecting all subsystems resulted in a sys-

tem with an order of 200. To reduce the order of the

complete system, model reduction techniques were

used [28], allowing the elimination of all states with

negligible influence on the input–output behaviour.

Through this, the model’s order was reduced to 60

states exhibiting only a slight deviation in frequency

response compared to the higher order system.

4.3 Comparison of model and experiments

Top and middle frames in Fig. 4 depict the transfer

function of the input voltage of the upstream loud-

speaker u1to the second microphone downstream of

the burner for the reduced order system. This is the

position, where the AOE is located in the simulation

model.

Note that the grey curves in Fig. 4 show gain and

phase of the transfer function obtained from the inter-

connection of the subsystems’ state-space represen-

tations. The individual state-space formulations of the

seven subsystems (up- and downstream end, burner

element, ducts, and output element) were linked to

one state-space system.

Only for frequencies below 150 Hz, around 470 Hz,

and above 900 Hz, a slight deviation in magnitude is

observed. Regarding the phase, the results are even

better, showing anobservabledeviationonly at300 Hz.

It is important to note here, that the modelled fre-

quency response shown in Fig. 4 was not simply iden-

tified from the experimental data but was predicted

based on the individual subsystem models.

The interconnected state-space model of the

complete system – reduced to 60 states – was used

to simulate the pressure fluctuations at the defined

microphone position. For this purpose, the upstream

loudspeaker was driven with a signal, which consisted

of three sines (all the same in amplitude) with arbi-

trary frequencies set to 203, 440, and 580 Hz. White

Fig. 4 Transfer function of input voltage u1(upstream

loudspeaker) to sound pressure at position of out-

put element; measured (black), interconnected

model (grey). Top: magnitude, middle: phase,

bottom: normalized sound pressure amplitude

for measured (black) and predicted (grey) data

and their deviation

noise, band-pass filtered between 100 and 600 Hz,

superimposed these sines.

Time responses of measured and predicted pres-

sure fluctuations are also shown in Fig. 4 (bottom).

Both signals show nearly perfect correspondence in

phase and in amplitude. For clarity, their difference is

also plotted (same scale). A comparison of measured

and predicted data obtained with downstream excita-

tion or excitation on both sides exhibited equally good

results.

5 APPLICATION TO THERMOACOUSTIC

STABILITY ANALYSIS

5.1 Stability of thermoacoustic systems

To assess the linear stability of a thermoacoustic

system, the most commonly used procedure is the

following [1,11,12].

1. The system is divided into different elements (duct,

burner, etc.).

2. Transfer matrices for the individual elements are set

up in frequency domain.

3. All transfer matrices are assembled in a system

matrix accounting for proper matching conditions.

4. The system’s dispersion relation det(system

matrix) =0 is solved to obtain the system eigen-

values.

Proc. IMechE Vol. 221 Part A: J. Power and Energy JPE384 ©IMechE 2007

Modelling and analysis of complex thermoacoustic systems 663

5. Eigenvalues with negative imaginary parts indi-

cate instability (time dependence ∼eiωt), the cor-

responding real parts the frequency of the growing

oscillation.

Concerning this procedure, some remarks should be

made. In step 2, the transfer matrices are generally of

a non-rational form, since the acoustic wave propaga-

tion is governed by a partial differential equation. Also,

solving det (system matrix)=0 will only be possible,

if all transfer matrices are available in analytical form.

Except for academic examples this will not be the case.

At least some elements will have to be determined

from CFD or FEM computations or experimentally.

However, even if closed-form expressions would be

available for all transfer matrices involved, solving for

the system’s eigenvalues can be quite cumbersome.

The dispersion relation is in general a complex-valued

transcendental equation, rapidly increasing in length

with network size.

Thestabilityanalysis based onastate-space descrip-

tion follows a procedure similar to the one given above

for the traditional frequency domain approach. Con-

necting all subsystems results in a state-space system

that represents the dynamics of the complete system.

To assess stability, the temporal evolution of free oscil-

lations has to be considered, which is governed by the

dynamicsmatrix A.Similarto therootsofthetranscen-

dental dispersion relation in the traditional approach,

the eigenvalues of the Amatrix, s(say), govern the sta-

bility of the system. Re(s)>0 indicates an unstable

mode with frequency |Im(s)|/2π. This procedure was

proposed by Schuermans et al. [1].

From a numerical point of view, determining system

stability from the state-space approach is straightfor-

ward, since calculating the eigenvalues from a given A

matrix can be performed with standard linear algebra

routines. In contrast to that, solving the transcen-

dental dispersion relations resulting from frequency

domain representations of thermoacoustic systems

is usually not a trivial task, even for networks of

medium complexity. These dispersion relations gen-

erally have an infinite number of solutions with only

those in the low-frequency regime being of physi-

cal relevance. Evidently, the state-space approach can

only provide a finite number of eigenvalues equal

to the dimension of A. This is due to the fact that

the infinite-dimensional system has been reduced

to a finite-dimensional state-space. In practice, how-

ever, the analysis will be limited to the low-frequency

regime. One reason for this is that thermoacoustic

instabilities most commonly appear in the form of

low-frequency oscillations. In addition to that, the net-

work approach based on plane wave matching is only

valid for frequencies below the cut-on frequency of the

first non-planar mode.

Fig. 5 Model system for stability analysis

5.2 A thermoacoustic model system

To assess feasibility of the state-space approach with

regard to linear stability analysis, the stability of

a thermoacoustic model system is investigated. In

order to consider the accuracy of the state-space

approach, the results have to be compared to a cor-

responding frequency domain calculation and there-

fore, only elements with fully analytical description

are chosen.

The model system is shown in Fig. 5 and comprises

the following subsystems:

(a) constant inlet impedance;

(b) uniform duct;

(d) distributed time-lag model for the flame with

variable mean time-delay and variable variance;

(e) uniform duct;

(f) Levine–Schwinger outlet impedance.

For the inlet impedance, a value of Zin =9 is chosen,

corresponding to a partially reflecting boundary with

a reflection coefficient of 0.8. The duct connecting the

inlet and the burner has a length of Lc, a temperature of

Tc, and carries a mean flow with velocity U. All param-

eter values describing the model system can be found

in Table 1.

The burner is represented by the L–ζmodel (cf. [23])

for a compact element with the transfer matrix

T=1−ikLeff −ζM

0α(11)

Leff is the effective length of burner and area expan-

sion, ζthe loss coefficient, and kdenotes the wave

number of the plane wave field. The change in cross-

sectional area is given by α. It can be noted here, that

the upper right transfer function in equation (11) is

improper [22] and, therefore, cannot be given in state-

space form. However, the scattering form of the L–ζ

model has a fully proper representation.

For the flame, a time-lag model with a Gaussian dis-

tribution of time-delays [29] was chosen. In this case,

the transfer function coupling the velocities across the

flame reads

T22 =1+Th

−1e−iωτ e−1/2σ2ω2(12)

JPE384 ©IMechE 2007 Proc. IMechE Vol. 221 Part A: J. Power and Energy

664 M R Bothien,J P Moeck,A Lacarelle,andCOPaschereit

Table 1 Parameter values describing the model system

Parameter LcTcUL

eff ζα LhThr

Value 0.95 m 523 K 10 m/s 0.13 m 5 0.3025 1.3 m 1573 K 0.1 m

where τand σdenote mean and variance of the time-

delay distribution.

The duct coupling the flame with the outlet bound-

ary condition is assigned a length of Lhand a tem-

perature of Th. The outlet is modelled as a long wave

Levine–Schwinger exit condition [30], the (specific)

impedance being given by

Zout =1

4(kr)2+i0.6133kr (13)

where ris the outlet duct radius. The reflection coeffi-

cient is related to the impedance through the bilinear

transform R=(Z−1)/(Z+1). Thus, the condition

relating fand gat the boundary is a proper transfer

function and can be converted to a state-space real-

ization. The ducts are assumed to support undamped

acoustic wave propagation. As stated in section 3.1,

a state-space model of the ducts is determined using

Padé approximations.

To obtain stability characteristics, the procedure

outlined in section 5.1 is followed. The roots of the

dispersion relation are determined with an itera-

tive solver working on a uniformly distributed ran-

dom field of initial guesses in the complex ω-plane.

The search was restricted to a rectangular domain

with |Re(ω)|/2π<1000 Hz and |Im(ω)|<100 rad/s. As

mentioned before, practical relevance justifies the

restriction to the low-frequency regime. Also, highly

damped modes (those with large Im (ω)) are not of

interest and highly unstable ones are usually consid-

ered unphysical.

Figure 6 shows the complex eigenfrequencies of

the model system in the specified interval for two

sets of flame model parameters (σvariable, τ=7 ms),

Fig. 6 Eigenvalues from frequency domain and

state-space computation; ◦,: frequency

domain σ◦=2ms,σ=2.5 ms, +,×: state-space

σ+=2ms,σ×=2.5 ms

computed from frequency domain and state-space

approach (see figure caption). It can be clearly seen

that the frequency domain and the state-space cal-

culations give identical results for both parameter

combinations.

Furthermore, it can be observed how a slight

increase in the variance of the time-delay distribution

can lead to a stabilization of the system. For σ=2ms

a low-frequency mode close to 100 Hz has a negative

damping rate and is therefore linearly unstable. For a

25 per cent higher variance, the unstable mode and the

least stable mode are shifted to higher damping rates

resulting in a stable system (,×). Note also, how only

the lowest frequency modes are affected by a change

in the time-delay variance.

Havingshownthat thestabilityanalysisbased onthe

state-space model accurately reproduces the results

from the traditional frequency domain approach, this

method can be used to investigate the influence of the

flame model parameters in more detail.

Figure 7 illustrates the dependence of system sta-

bility on the flame model parameters as obtained

from the state-space description. Note that the sta-

bility problem has been solved 15 000 times in (σ,τ)-

parameter space in less than an hour. This would be

impossible to accomplish by using graphical meth-

ods based on Bode or Nyquist diagrams. Generating

a stability map based on solutions of the frequency

domain dispersion relation would be much more

time-consuming as well. The greyscale coding indi-

cates the growth rate of the least stable mode. The

thick black line in the plot denotes the stability border.

As can be seen, the highest growth rate can be found

for a fully concentrated (σ=0) time-delay of about

Fig. 7 Growth rate in 1/s of the least stable mode as a

function of time-delay mean τand variance σ

(equation (12))

Modelling and analysis of complex thermoacoustic systems 665

Fig. 8 Frequency in Hz of the least stable mode as a

function of time-delay mean τand variance σ

(equation (12))

2.5 ms. For this case, two more local maxima in growth

rate can be identified at τ=8.5 ms and τ=14.5 ms.

An increase in time-delay variance has a purely stabi-

lizing effect. For a fixed mean time-delay, the growth

rate of the most unstable mode decreases monotoni-

cally with variance. In fact, for σ>2.25 ms all modes

are stable. The stabilizing effect of the time-delay

variance is in accordance with Polifke et al. [31].

A contour plot of the frequency of the least stable

mode as a function of time-delay mean and variance is

shown in Fig. 8.The stability border is added again as a

thick black line. In the unstable region, low-frequency

modes in the range of 75–175 Hz are dominant. For

certain values of τand σisocontours of the domi-

nant mode’s frequency meet, indicating that there are

at least two unstable eigenvalues. Note that lines of

constant mode frequency have an essentially vertical

trend. This corresponds to the fact that the frequency

of the least stable mode does not change considerably

with time-delay variance [31].

With respect to computing time, the stability anal-

ysis based on state-space models was recognized to

be clearly more efficient than the frequency domain

approach. Additionally, using the latter for paramet-

ric studies is not straightforward, since for all sets of

parameters all roots of the dispersion relation in the

domain considered must be located. This can only

be ensured if a high number of initial guesses for

the iterative solver is chosen. However, this results

in a considerably higher computing time. Although

graphical methods are a useful tool to investigate ther-

moacoustic stability for fixed parameter sets, they are

not suitable for generating stability maps.

6 APPLICATION TO AN ATMOSPHERIC

COMBUSTOR TEST RIG

In this section, the state-space methodology is

applied to investigate the stability of an atmospheric

Fig. 9 Setup of combustion test rig with instrumen-

tation

combustor test rig with a swirl stabilized burner.

A schematic of the test rig is shown in Fig. 9.

The lean-premixed burner features aerodynamic

flame stabilization through a vortex breakdown

induced by strong swirl. Pressure field and OH∗-

chemiluminescence, an indicator of heat release [32],

can be monitored. To investigate its thermoacoustic

properties the test rig is equipped with loudspeak-

ers and microphones mounted up- and downstream

of the flame. A more detailed description of the test

rig can be found in references [33] and [34]. In order

to generate thermoacoustic instabilities, a resonance

tube was attached to the combustion chamber. This

resulted in a high amplitude quarter wave mode insta-

bility at ∼100 Hz for certain preheated, lean-premixed

operating conditions.

Stability analysis based on the state-space approach

was applied to investigate the impact of a fuel distribu-

tion modification on thermoacoustic system stability.

The fuel distribution was modified as such, that the

majority of the fuel (natural gas) was directed towards

the outer flame region. This resulted in an enriched

outer combustion zone, whereas the inner region was

made leaner. For more details see reference [35].

The impact of the fuel distribution modification on

the flame response was taken into account by means

of the flame transfer function. Flame responses for

both, the baseline and the modified case were mea-

sured at the operating conditions where stability was

investigated. For the modified distribution, the phase

response was quite similar to the standard case, how-

ever, in the frequency range of 80–150 Hz, the ampli-

tude response was considerably lower [16,35].

In addition to the experimentally determined flame

responses, a parametric model for the upstream

reflection coefficient was identified and included in

the model. This was considered necessary since the

upstream acoustic boundary condition was rather

complex due to the preheater and air supply geome-

try. The complete model was composed of 11 elements

and had a state dimension of 120.

System eigenfrequencies for the baseline and the

modified case are shown in Fig. 10. The baseline

666 M R Bothien,J P Moeck,A Lacarelle,andCOPaschereit

Fig. 10 Complex eigenfrequencies of system model; :

baseline, ◦: modified

case has an eigenvalue with negative imaginary part

close to 100 Hz, indicating an unstable mode. For the

enriched outer combustion zone, the initially unsta-

ble eigenvalue, representing the quarter wave mode,

moves across the real axis to the upper half plane. This

implies that the system is now linearly stable.

The stabilizing effect of the enriched outer com-

bustion zone could be verified in the experiment.

Frequency spectra of the unsteady combustor pres-

sure for the standard and the modified case are shown

in Fig. 11.

The spectrum for the standard case exhibits a dom-

inant peak at 100 Hz corresponding to the unstable

quarter wave mode. In contrast to that, the spectrum

for the modified case shows a significant reduction of

the spectral peak amplitude (more than 20 dB). It is

difficult to judge if the system is actually stable since

the pressure amplitudes at the quarter wave frequency

were still high. Despite that, high pulsation ampli-

tudes can also occur in a stable lightly damped system,

which is driven by noise [36]. Compared to the stan-

dard case, the stabilizing effect of the modified fuel

distribution, as predicted by the model, can be clearly

seen in the experimental data. It should be also men-

tioned that enriching the outer combustion zone did

not have a stabilizing effect for all operating conditions

that were tested [35].

Fig. 11 Measured amplitude spectra of the unsteady

combustor pressure for the baseline and the

modified case

7 SUMMARY

A modular approach for modelling complex thermoa-

coustic systems based on state-space models was

presented. Frequency and time domain results for

a purely acoustic model system were compared to

experimental dataandproved themethodtobe of high

accuracy.

Additionally, application to stability analysis of ther-

moacoustic systems was investigated. Results were

compared to a standard frequency domain analy-

sis of the same system. Exactly the same results

were obtained with a fraction of computational effort

compared to that of the standard procedure.

The impact of a modified fuel distribution on ther-

moacoustic stability of an experimental combustor

was investigated using the state-space methodol-

ogy. Based on the model, a stabilizing effect was

predicted, which was confirmed by experimental

data.

REFERENCES

1 Schuermans, B., Bellucci, V., and Paschereit, C. O.

Thermoacoustic modeling and control of multi burner

combustion systems. In Proceedings of the ASME Turbo

Expo, Atlanta, USA, 2003, paper 2003-GT-38688.

2 Lieuwen,T.C.andYang,V.(Eds) Combustioninstabilities

in gas turbine engines. In Progress in astronautics and

aeronautics, vol. 210, 2005 (AIAA, Inc., Reston).

3 Dowling, A. P. and Hubbard, S. Instability in lean pre-

mixed combustors. Proc. Instn Mech. Engrs, Part A:

J. Power and Energy, 2000, 214, 317–332.

4 Dowling, A. P. and Morgans, A. S. Feedback control of

combustion oscillations. Annu. Rev. Fluid Mech., 2005,

37, 151–182.

5 Seume, J. R., Vortmeyer, N., Krause, W., Hermann,

J., Hantschk, C.-C., Zangl, P., Gleis, S., Vortmeyer,

D., and Orthmann, A. Application of active combus-

tion instability control to a heavy duty gas turbine.

Trans. ASME, J. Eng. Gas Turbines Power, 1998, 19(6),

549–566.

6 Sewell, J., Sobieski, P., and Beers, C. Application of

continuous combustion dynamics monitoring on large

industrial gas turbines. In Proceedings of the ASME

Turbo Expo, Vienna, Austria, 2004, paper GT2004-54310.

7 Rea, S., James, S., Goy, C., and Colechin, M. J. F. On-line

combustion monitoring on dry low NOxindustrial gas

turbines. Meas. Sci. Technol., 2003, 14, 1123–1130.

8 Roux,S.,Lartigue,G.,Poinsot,T.,Meier,U.,and Bérat,C.

Studies of mean and unsteady flow in a swirled combus-

tor using experiments, acoustic analysis, and large eddy

simulations. Combust. Flame, 2005, 141, 40–54.

9 Schuermans, B., Luebcke, H., Bajusz, D., and Flohr, P.

Thermoacoustic analysis of gas turbine combustion sys-

tems using unsteady CFD. In Proceedings of the ASME

Turbo Expo, Reno-Tahoe, USA, 2005, paper GT2005-

68393.

Modelling and analysis of complex thermoacoustic systems 667

10 Polifke, W. Combustion instabilities. In Advances in

aeroacoustics and applications, lecture series, VKI LS

2004-05, 2004 (Von Karman Institute, Brussels, Belgium).

11 Sattelmayer, T. and Polifke, W. Assessment of methods

for the computation of the linear stability of combustors.

Combust. Sci. Technol., 2003, 175, 453–476.

12 Sattelmayer, T. and Polifke, W. A novel method for

the computation of the linear stability of combustors.

Combust. Sci. Technol., 2003, 175, 477–497.

13 Evesque, S. and Polifke, W. Low-order modelling for

annular combustors: validation and inclusion of modal

coupling. In Proceedings of the ASME Turbo Expo, Ams-

terdam, The Netherlands, 2002, paper GT-2002-30064.

14 Stow, S. R. and Dowling, A. P. Thermoacoustic oscil-

lations in an annular combustor. In Proceedings of

the ASME Turbo Expo, New Orleans, USA, 2001, paper

2001-GT-0037.

15 Schuermans, B. B. H., Paschereit, C. O., Polifke,W., and

van der Linden, J. H. Prediction of acoustic pressure

spectra in combustion systems using swirl stabilized gas

turbine burners. In Proceedings of the IGTI Turbo Expo,

Munich, Germany, 2000, paper 2000-GT-0105.

16 Paschereit, C. O., Moeck, J. P., and Bothien, M. R. State-

space modelling of thermoacoustic systems for stability

analysis and time domain simulation. In Proceedings of

the 13th International Congress on Sound andVibration,

Vienna, Austria, 2006.

17 Dupère, I. D. J. and Dowling, A. P. The use of Helmholtz

resonators in a practical combustor. Trans. ASME, J. Eng.

Gas Turbines Power, 2005, 127, 268–275.

18 Paschereit, C. O. and Gutmark, E. The effectiveness of

passive combustion control methods. In Proceedings

of the ASME Turbo Expo, Vienna, Austria, 2004, paper

GT-2004-53587.

19 McManus, K. R., Poinsot, T., and Candel, S. M. A review

of active control of combustion instabilities. Prog. Energy

Combust. Sci., 1993, 19, 1–29.

20 Munjal,M. L. Acoustics of ducts and mufflers, 1976 (Wiley

& Sons, Inc., NewYork).

21 Peat, K. S. Evaluation of four-pole parameters for ducts

with flow by the finite element method. J. Sound Vib.,

1982, 84(3), 389–395.

22 Zhou, K. and Doyle, J. Essentials of robust control, 1997

(Prentice Hall, Upper Saddle River, NJ).

23 Schuermans, B. B. H., Polifke,W., and Paschereit, C. O.

Modeling transfer matrices of premixed flames and com-

parison with experimental results. In Proceedings of

the ASME Turbo Expo, Indianapolis, USA, 1999, paper

1999-GT-0132.

24 Paschereit, C. O., Schuermans, B., Polifke, W., and

Mattson, O. Measurement of transfer matrices and

source terms of premixed flames. Trans. ASME, J. Eng.

Gas Turbines Power, 2002, 124, 239–247.

25 Åbom, M. A note on the experimental determination of

acoustical two-port matrices. J. SoundVib., 1992, 155(1),

185–188.

26 Gustavsen, B. and Semlyen, A. Rational approximation

of frequency domain responses by vector fitting. IEEE

Trans. Power Deliv., 1999, 14(3), 1052–1061.

27 Kollár, I., Pintelon, R., and Schoukens, J. Frequency

domain system identification toolbox for matlab:

improvements and new possibilities. In Proceedings of

the XI. IFAC/IFORS International Symposium on Sys-

tem Identification and Parameter Estimation, Fukuoka,

Japan, 1997, pp. 991–994.

28 Antoulas, A. C. and Sorensen, D. C. Approximation of

large-scale dynamical systems: an overview. Int. J. Appl.

Math. Comput. Sci., 2001, 11(5), 1093–1121.

29 Schuermans, B., Bellucci, V., Guethe, F., Meili, F.,

Flohr, P., and Paschereit, C. O. A detailed analysis of

thermoacoustic interaction mechanisms in a turbulent

premixed flame. In Proceedings of the ASMETurbo Expo,

Vienna, Austria, 2004, paper GT 2004-53831.

30 Levine, H. and Schwinger, J. On the radiation of sound

from an unflanged circular pipe. Phys. Rev., 1948, 73(4),

383–406.

31 Polifke, W., Kopitz, J., and Serbanovic, A. Impact of

the fuel time lag distribution in elliptical premix noz-

zles on combustion stability. In Proceedings of the 7th

AIAA/CEAS Aeroacoustics Conference, Maastricht, The

Netherlands, 2001, paper 2001-2104.

32 Haber, L. C., Vandsburger, U., Saunders, W. R., and

Khanna, V. K. An examination of the relationship

between chemiluminescent light emission and heat

release rate under non-adiabatic conditions. In Proceed-

ings of the ASME Turbo Expo, Munich, Germany, 2000,

paper 2000-GT-0121.

33 Albrecht, P., Bauermeister, F., Bothien, M. R., Lacarelle,

A., Moeck, J. P., Paschereit, C. O., and Gutmark, E.

Characterization and control of lean blowout using

periodically generated flame balls. In Proceedings of

the ASME Turbo Expo, Barcelona, Spain, 2006, paper

GT2006-90340.

34 Moeck, J. P., Bothien, M. R., Guyot, D., and Paschereit,

C. O. Phase-shift control of combustion instability using

(combined) secondary fuel injection and acoustic forc-

ing. In Proceedings of the 1st Active Flow Control

Conference, Berlin, Germany, 2006.

35 Lacarelle, A., Moeck, J., Konle, H., Nayeri, C. N., and

Paschereit, C. O. Effect of the fuel/air mixing on NOx

emissions and stability in a gas premixed combustion

system. In Proceedings of the 45th AIAA Aerospace Sci-

ence Meeting, Reno-Tahoe, USA, 2007, paper 2007-1417.

36 Banaszuk,A.,Jacobson, C. A., Khibnik, A. I., and Mehta,

P. G . Linear and nonlinear analysis of controlled com-

bustion processes. Part I: linear analysis. In Proceed-

ings of the IEEE International Conference on Control

Applications, Hawaii, USA, 1999, pp. 199–205.

APPENDIX

Notation

A,B,C,Dstate-space matrices (−)

cspeed of sound (m/s)

fdownstream travelling wave (m/s)

gupstream travelling wave (m/s)

kwave number (1/m)

MMach number (−)

pacoustic pressure scaled with ρc(m/s)

rradius (m)

Rreflection coefficient (−)