Numerical study on intrinsic thermoacoustic instability of a laminar premixed flame

Combustion and Flame 162 (2015) 3370–3378

Contents lists available at ScienceDirect

Combustion and Flame j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m b u s t fl a m e

Numerical study on intrinsic thermoacoustic instability of a laminar premixed flame Camilo F. Silva ⇑, Thomas Emmert, Stefan Jaensch, Wolfgang Polifke Professur für Thermofluiddynamik, Technische Universität München, D-85747 Garching, Germany

a r t i c l e

i n f o

Article history: Received 24 November 2014 Received in revised form 25 May 2015 Accepted 1 June 2015 Available online 17 June 2015 Keywords: Intrinsic thermoacoustic feedback Finite Impulse Response (FIR) Flame Transfer Function (FTF) Combustion instability

a b s t r a c t A study on the velocity sensitivity and intrinsic thermoacoustic stability of a laminar, premixed, Bunsen-type flame is carried out. Direct numerical simulation (DNS) of the flame, placed in an acoustically anechoic environment and subjected to broad-band, low-amplitude acoustic forcing, generates time series of fluctuating heat release rate, velocities and pressure. The time series data is post-processed with system identification to estimate the impulse response and transfer function of the flame. The associated frequency response is validated against experiment with good accuracy. DNS results obtained with acoustic excitation from the inlet or outlet boundary, respectively, confirm that the flame responds predominantly to perturbations of velocity. The stability of eigenmodes related to intrinsic thermoacoustic feedback is investigated with a network model. Both stable and unstable intrinsic thermoacoustic modes are predicted, depending on details of the configuration. The predicted modes are directly observed in direct numerical simulations, with good agreement in frequencies and stability. Ó 2015 The Combustion Institute. Published by Elsevier Inc. All rights reserved.

1. Introduction Thermoacoustic combustion instabilities have been an important subject of combustion research for many decades, originating from the enormous problems encountered during the design and test of rocket engines [1,2]. The increased interest in recent years is due to the occurrence of thermoacoustic instabilities in gas turbine power-plants, aero engines and other low emission combustion systems, which often operate under lean premixed conditions [3,4]. The established understanding of self-excited thermoacoustic instabilities involves feedback between unsteady heat release by the flame and acoustic waves in fuel and air supply and the combustion chamber: The unsteady flame acts as a monopole source of sound and generates acoustic waves, which are reflected by the combustion system back towards the flame, where they perturb in turn the heat release rate. This feedback loop may go unstable provided phase lags are favorable [5] and the loss of perturbation energy not excessive. In the simplest case, premix flames are described as velocity sensitive: the global heat release rate Q_ 0 responds to perturbations u0 of velocity upstream of the flame. The link between these two ⇑ Corresponding author. E-mail addresses: [email protected] (C.F. Silva), [email protected] (T. Emmert), [email protected] (S. Jaensch), [email protected] (W. Polifke).

quantities is commonly given in terms of the flame frequency response or flame transfer function F ðxÞ, such that

^_ ^ ð xÞ u Qð xÞ ;  ¼ F ðxÞ u  Q_

ð1Þ

or the corresponding flame impulse response hðtÞ

Z Q_ 0 ðtÞ 1 T ¼ hðsÞ u0 ðt  sÞ ds:   0 u Q_

ð2Þ

In these equations, ð0 Þ denotes the deviation of a quantity from its mean value ðÞ and the ð^Þ stands for the Fourier transform. The symbols x and s represent angular frequency and time delay, respectively. T accounts here for the duration of the impulse response hðtÞ. The flame transfer function (FTF) as defined in Eq. (1) has been used extensively in the study of thermoacoustic instabilities of laminar and turbulent flames. Experimentally, the FTF is typically evaluated by imposing a sinusoidal excitation of the flow upstream of the flame and by measuring the corresponding response in heat ^ release Q_ ðxÞ [6–11]. Sinusoidal excitation also gives satisfactory results in numerical studies of flame dynamics [12,13], but only at considerable computational cost. The identification of the flame response through the corresponding impulse response hðsÞ – see Eq. (2) – offers an alternative, which is computationally more efficient. This method was introduced in the work of Polifke et al. [14],

http://dx.doi.org/10.1016/j.combustflame.2015.06.003 0010-2180/Ó 2015 The Combustion Institute. Published by Elsevier Inc. All rights reserved.

3371

C.F. Silva et al. / Combustion and Flame 162 (2015) 3370–3378

Zhu et al. [15] and Gentemann et al. [16], and subsequently used and refined in a number of studies [17–20]. Once the flame transfer function is known, it can be implemented in analytical dispersion relations [21–23], network models [21,24,25] or Helmholtz solvers [26–28] to study the thermoacoustic stability of combustion systems. The present paper originated from an attempt to estimate the dynamics of a laminar premix flame by direct numerical simulation (DNS) combined with system identification (see the Ref. [29] for a description of the approach and previous applications in aero- and thermoacoustics). A computational model (case B in Section 4.2) with acoustically non-reflecting boundary conditions was set up, because it is understood that low reflection coefficients foster accurate identification [30]. However, it was observed that this computational model responded to low amplitude broad-band excitation with very strong, resonant, low-frequency oscillations in flame shape and heat release rate, such that identification of the flame dynamics in the limit of linear behavior was not possible. This was not expected, since acoustically non-reflecting boundaries imply a significant loss of acoustic energy. Furthermore, the low frequencies of the exhibited thermoacoustic instability could not be explained. Due to the small size of the domain, any acoustic cavity resonances would occur at frequencies much higher than that of the combustion dynamics. Hoeijmakers et al. [31] observed in an independent study that low-order acoustic network models of a combustion system with zero acoustic reflection coefficients can indeed possess unstable eigenmodes. This might appear unphysical at first sight, because non-reflecting boundaries should break the feedback between the flame and the system acoustics introduced above. However, Hoeijmakers et al. [31] demonstrated analytically with a simple n-s model for the flame dynamics that under some conditions the flame scattering matrix may indeed be an ‘‘intrinsically unstable element’’. An interpretation of the physics of the instability was not developed, but instabilities observed in a combustion test rig with low acoustic reflection coefficients [32] lend some support to the argument that flame-acoustic coupling may be intrinsically unstable. A physics-based explanation of Intrinsic Thermo-Acoustic (ITA) modes in terms of flow – flame – acoustic interactions was developed by Bomberg et al. [33]. The interactions were analyzed in a representation that respects causality, i.e. a formulation that allows to distinguish input from output signals and thus cause from effect. In this framework a thermoacoustic feedback loop intrinsic to the thermo-acoustic coupling was identified, which may go unstable, or exhibit resonant amplification. As a result, the generation of acoustic energy by a perturbed flame may exhibit very strong peaks, which can be associated with poles of the flame scattering matrix ([34,31]). The ITA feedback loop was detailed by Bomberg et al. [33] and Emmert et al. [34] as follows: a flow perturbation u0 just upstream of the flame causes perturbations that are convected along the length of the flame and cause in turn a fluctuation in heat release rate Q_ 0 , see Eq. (1). The fluctuation Q_ 0 then generates acoustic waves that propagate away from the flame. The acoustic wave propagating in the upstream direction modulates the upstream velocity u0 , thus closing the feedback loop. A more detailed analysis of the ‘‘physics’’ of ITA feedback and the structure of the corresponding modes is given by Courtine et al. [35]. In the light of these findings, it was understood that the unexpected behavior of the laminar premix flame DNS described above should be interpreted as a result of resonance of the external excitation with the ITA feedback loop. Without persistent excitation, the DNS develops a strong instability, which can be identified as an unstable ITA mode (see Section 4.2). Further analysis of the case

with low-order models helped to identify a configuration where the intrinsic mode is stable, such that identification of flame dynamics is possible. Before completion of the present study, several questions, which indeed are strongly related to each other, were pending. It was questioned whether the thermoacoustic intrinsic instability is a real physical phenomenon or merely a spurious result of simplistic network models. Furthermore, it was debated whether the response of a compact premixed flame is actually velocity sensitive, or whether flame dynamics should be described as a response also to pressure. In addition, more direct evidence for the validity of the ITA feedback as described by Bomberg et al. [33] was sought. Accordingly, the present study has the following objectives: First, to demonstrate, by Direct Numerical Simulation, that the ITA feedback is an authentic physical phenomenon that is present in premixed combustion systems. Second, to give evidence, by an adequate identification of the FTF, that compact laminar premixed flames indeed respond predominantly to perturbations of upstream velocity and thus may be treated as velocity sensitive elements. Furthermore, we confirm that thermoacoustic network models in conjunction with the identified FTF are indeed capable of analyzing the stability of ITA modes. This article is organized as follows. In the next section, a brief overview of system identification is given. In the third section the flame transfer function of the laminar flame under study is estimated with DNS/SI, by imposing upstream as well as downstream broadband acoustic forcing, and validated against experimental and computational results published previously [36,12]. Results are analyzed under the description of a laminar premix flame considered as a velocity sensitive element. In the fourth section, low-order models of a burner-stabilized laminar flame are developed to compute the frequencies and growth rates of acoustic modes for various acoustic reflection coefficients of the inlet–outlet boundaries, including the limiting cases of zero reflection coefficients. Stable as well as unstable intrinsic modes are found, depending on the details of the configuration. Corresponding results from DNS of laminar premix flame corroborate the results of the low order models, and some properties of the ITA mode are studied. In concluding this study, the reader is cautioned that intrinsic thermoacoustic instabilities are not to be mistaken for intrinsic thermo-diffusive or hydrodynamic instabilities of premix flames, as reviewed e.g., by Clavin [37]. 2. Overview of system identification In order to infer a model of the relation between inputs and outputs of a system, system identification uses the information contained in corresponding time series data [38]. If a Single Input Single Output (SISO) system, considered as a black box, is assumed to be Linear Time Invariant (LTI), it is then completely characterized by its Impulse Response (IR). Be aware that the IR is a property of a such a system, not a model. When the system does not exhibit internal feedback, i.e. when the outputs of the system can be assumed to be related to the input by merely time delays, a model for the system can be written as

rt ¼ b0 st þ b1 st1 þ    þ bnb stnb ;

ð3Þ

where r t and st denote the output and the input of the system at discrete time t, respectively. In this case, the IR of the system can be modeled as a polynomial

GðqÞ ¼ b0 þ b1 q1 þ    þ bnb qnb ;

ð4Þ k

where the shift operator is defined by the property q st ¼ stk . The coefficients bk describe the response of the system to a unit impulse

3372

C.F. Silva et al. / Combustion and Flame 162 (2015) 3370–3378

st ¼ dt0 . In this situation we consider GðqÞ as a Finite Impulse Response (FIR), because after a certain (discrete) time the contribution of the coefficients bk , for k > nb , is small and therefore negligible. The FIR model has been used in several studies that combine numerical simulations and system identification (see review of Ref. [29]) when studying the dynamic response of flames. In addition to the discrete time representation of LTI systems, the corresponding response in the frequency domain is very useful, since it characterizes the associated filter behavior of the system. The flame transfer function (FTF) is derived as the z-transform of the impulse response GðqÞ. It reads:

F ð xÞ ¼

ng X

bk eixkDt :

Table 1 Operating and thermal conditions. Equivalence ratio

Inlet velocity uin

Premixture temperature Ta

Duct temperature Td

Combustor wall temperature T w

0.8

0.4 m s1

293 K

373 K

373 K

ð5Þ

k¼0

where i denotes the imaginary unit and Dt is the sample interval in the discrete time domain. Note that the angular frequency x is considered here as a complex-valued variable. Identifying or ‘‘estimating’’ a FIR model means finding the optimal number and values of the corresponding coefficients. The coefficients nb and bk in a FIR model are commonly obtained by correlation analysis in terms of the Wiener–Hopf equation [38]. In the present study, the system identification tool box of Matlab 2013b was used. 3. Numerical assessment of the flame impulse response and flame transfer function This section describes the methodology to obtain by direct numerical simulation (DNS) the flame impulse response and the corresponding flame transfer function of the laminar flame studied experimentally by [36]. Two different cases are evaluated with acoustic excitation from the inlet or outlet boundary, respectively. This procedure is based on System Identification (SI) and is similar to the methods used by Kaess et al. [17], Föller and Polifke [39] and Tay-Wo-Chong et al. [19]. Discussion on the velocity sensitivity of the flame under study is presented. 3.1. Experimental set-up The test rig used by Kornilov et al. [36] consists of a vessel with 8 laminar Bunsen-type flames stabilized at a flat perforated plate. This plate is of 1 mm thickness and perforated by 8 rectangular slits of dimensions 12  2 mm, separated from each other by 3 mm. The premixture is composed of air and methane at an equivalence ratio equal to 0.8. The operating and thermal conditions are found in Table 1 and illustrated in Fig. 1. Following the study of Duchaine et al. [12], the temperatures T d and T w are assumed to be around 373 K, even though it is known by experiments [36] that this temperature lays between 373 K ± 50 K during steady combustion.

Fig. 1. Test rig configuration.

needs as an input the plane acoustic waves that are leaving the computational domain, a Characteristic Based Filter [42] is also considered. A two-step chemistry (2S-CM2) [43] is used since full chemistry schemes, for the flame exposed in the present work, do not present considerable differences in results. Detailed explanations on the subject can be found in the work of Duchaine et al. [12]. The flame front is resolved by 8 grid points on a 3D mesh, which leads to a tetrahedral computational grid containing around 100,000 nodes. 3.3. Numerical identification of flame dynamics In the time domain, Eq. (2) relates the response Q_ 0 ðtÞ to incomu ing velocity perturbations u0 ðtÞ. Let us define the input as s ¼ u0u =u  _  u and Q_ 0 ¼ Q_  Q and the output as r ¼ Q_ 0 =Q_ , where u0 ¼ uu  u . u

The index ‘u’ denotes the flow at a location immediately upstream of the flame (see Fig. 2). Eq. (2) is written in discrete time as an FIR model (see Eq. (3)) and the corresponding FTF is computed as the z-transform of the impulse response based on Eq. (5). When

3.2. Numerical solver The DNS/LES (Large Eddy Simulation) numerical tool AVBP1 developed by CERFACS and IFP is used in this work to solve the compressible set of Navier–Stokes equations. In the present study, AVBP discretizes the spatial terms of the equations considering a cell-vertex finite volume formulation based on the second-order accurate Lax-Wendroff scheme. Temporal evolution is computed through a second order accurate Runge–Kutta approach. The Navier Stokes Characteristic Boundary Conditions (NSCBC) [40] are implemented in AVBP to treat waves crossing the boundaries. Totally non reflecting acoustic boundary conditions are imposed based on a technique known as plane wave masking [41]. Since this approach 1

http://www.cerfacs.fr/4-26334-The-AVBP-code.php.

Fig. 2. DNS set-ups for identification of flame dynamics with excitation signal f in at the upstream boundary (left), and signal g out at the downstream boundary (right). Plate thickness l1 ¼ 1 mm, length l2 ¼ 8 mm.

C.F. Silva et al. / Combustion and Flame 162 (2015) 3370–3378

determined in this way, the FTF is well-defined not only for purely real frequencies, but throughout the complex frequency plane. This is important for the accurate determination of growth rates of thermoacoustic instabilities [44]. The velocity perturbation u0u just upstream of the flame – see location (u) in Fig. 2 (right) – is regarded as the input signal of the FIR model. This signal cannot be imposed directly in the numerical simulations, but results from acoustic perturbations imposed at the boundaries of the computational domain. In this context, recall the definition of characteristic wave amplitudes f ; g

f 

  1 p0 þ u0 ;  c 2 q

g

  1 p0  u0 ;  c 2 q

ð6Þ

where q denotes density and c the speed of sound. The characteristic wave f is propagating in the downstream direction, while g is traveling upstream. In this study, two different set-ups are used to identify the impulse response and subsequently the corresponding FTF. In the first case, called DNSu and illustrated in Fig. 2 (left), the acoustic excitation is imposed at the inlet boundary of the computational domain as a wave f in traveling in the downstream direction. Due to the non-reflecting boundary conditions, any upstream traveling signal g in leaves the domain without reflection at the inlet, while g out is equal to zero by design. In the second case, called DNSd and depicted in Fig. 2 (right), the excitation signal is imposed as g out at the outlet boundary. The anechoic conditions at the boundaries assure that f in is zero and f out leaves the domain without being reflected at the outlet. In order to obtain information suitable for identification of flame dynamics from the DNS, the signals f in and g out are designed to exhibit the following properties: First, the spectral energy should be uniformly distributed over the frequency range of interest. This assures that all the frequencies of interest are excited. Second, it should demonstrate fast decorrelation with itself, so that the length of the time series, a key parameter for a reliable estimation of the coefficients bk under correlation analysis, remains reasonably short. Third, the signal should exhibit a crest factor close to one. This is done with the purpose of imposing high energy signals with the lowest amplitude possible so that non-linearities are avoided [38]. In contrast to the works of Kaess et al. [17], Föller and Polifke [39] and Tay-Wo-Chong et al. [19], the signal preferred in the present study is not a Discrete Random Binary Signal (DRBS), but instead is constructed based on Daubechies Wavelets [45]. Although this type of signal has a slightly worse (larger) crest factor, the spectral energy distribution is closer to the desired one. The velocity perturbation signal u0u is evaluated in the post-processing part of the simulations from the time series of

3373

volume flow rate at location (u). The relation of the perturbation signal to the acoustic waves at this location is simply

u0u ¼ f u  g u ;

ð7Þ

which follows immediately from Eq. (6). In order to evaluate a satisfactory impulse response from data using correlation analysis, as in the present study, experience suggests that the length of data series should extend at least over ten times the length of the impulse response [38,46]. In our case, data captured over 100 ms is considered enough to evaluate the impulse response and subsequently the FTF of the flame, since the length of the IR is around 10 ms (see Fig. 3(a)). Note that only one simulation is needed to retrieve the frequency response for the frequency band of interest. If, in contrast, harmonic excitation is performed, the frequency response for only one frequency can be computed in one simulation. Moreover, if data during 100 ms is captured with such an input signal, good results in the frequency response would be expected only for frequencies equal or higher than 10 Hz (in the most favorable situation when data over one cycle is considered sufficient). Contrary to broadband excitation, information at lower frequencies would not be well captured. The impulse response of the flame, modeled by a FIR structure and obtained after solving Eq. (3) (with na ¼ 0), is shown in Fig. 3(a). It exhibits a maximum at s1 ¼ 3:9 ms and a minimum s2 ¼ 6:4 ms, which is indeed a dynamics that commonly characterizes laminar flames. The maximum (minimum) impulse response corresponds to the moment of maximum (minimum) flame elongation and surface area. Figure 3(b) shows the flame frequency response computed from Eq. (5) and compares it with the numerical results of Duchaine et al. [12] and the experimental measurements of Kornilov et al. [36], respectively. Good agreement is observed for upstream as well as downstream forcing. The overshoot of the FTF at 107 Hz results from constructive superposition of contributions of the positive and negative coefficients bk of the impulse response. If all coefficients bk were of the same sign, an overshoot could not exist [47]. It should be pointed out that in the study of Duchaine et al. [12], which was based on harmonic excitation, seven DNS runs were required to obtain the values of the flame frequency response at seven discrete points in the frequency spectrum. In the present study, on the other hand, a single DNS run was sufficient to obtain the FTF over the entire frequency band of interest. Now we turn to the question whether the flame is indeed velocity sensitive, which has been questioned in the context of ITA instabilities. Let us recall that a flame is considered velocity sensitive when the corresponding fluctuations of heat release rate are related only to upstream velocity perturbations. Note, however, that a flame in general can also respond vigorously to pressure

Fig. 3. Characterization of laminar premix flame dynamics. (a) Impulse response estimated from DNSu (dark circle) and DNSd (gray circle). (b) Flame frequency response from DNSu (dark line) and DNSd (gray line), Experiments of Kornilov et al. [36] (dashed line) and DNS of Duchaine et al. [12] (triangles).

3374

C.F. Silva et al. / Combustion and Flame 162 (2015) 3370–3378

fluctuations. Examples are partially premixed flames or non-compact flames, i.e. flames whose length is of the order of the acoustic wavelength. In what follows we will show, based on our results, that the flame under investigation, which is compact and perfectly premixed, is predominantly velocity sensitive. For the sake of our argument, we will perform a proof by contradiction. Let us assume that the flame response to upstream velocity perturbations is different depending whether the flame is acoustically excited from the upstream or downstream side. Subsequently ^ ^ u and p ^u as we assume that Q_ is linked to both u

^ ^u ^u Q_ u p þ F pres ðxÞ _ ¼ F vel ðxÞ u   u q cu u u u Q

ð8Þ

and suppose that pressure fluctuations are related to velocity fluctuations through an acoustic impedance Z

Z¼

^ p f þg ¼ :  qcu^ f  g

ð9Þ

It follows

^  u ^u Q_ F vel ðxÞ þ Z u F pres ðxÞ : _ ¼ |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} u u Q

ð10Þ

F ðx Þ

Recalling that both upstream and downstream boundary conditions are acoustically non reflective, we consider now the same two cases previously discussed. In the first case, when the flame is acoustically excited from downstream, the impedance Z u is equal to 1 for all frequencies since f u is zero for all frequencies. In the second case, when the flame is acoustically excited from upstream, the impedance Z u is not anymore 1, as illustrated in Fig. 4. Fol^ ^ u and p ^u lowing that idea, it is clear that if Q_ is related to both u as suggested in Eqs. (8) and (10), then F ðxÞ must be different for these two cases if F pres ðxÞ – 0. The fact that F ðxÞ is independent of the excitation position (upstream or downstream as shown in Fig. 3(b)) indicates that F pres ðxÞ ¼ 0. We conclude therefore that the premix flame investigated in this work is velocity sensitive, while it does not (or only very weakly) respond to pressure perturbations. In other words, the SISO model structure implied by Eq. (2) with u0u as input signal indeed respects the ‘‘physics’’ of the flame dynamics and in particular the causality of acoustic – flow – flame interactions: For the case DNSu with excitation from the inlet boundary, the resulting acoustic signal f u at location (u) generates a velocity perturbation according to Eq. (7), which then causes a modulation of heat release rate Q_ 0 . On the other hand, if the acoustic excitation is imposed at the outlet boundary for the case DNSd, the resulting characteristic wave g travels towards and through the flame, such that it eventually causes an acoustic perturbation g u in the slit. In the absence of a downstream traveling wave f u (see Fig. 2, right), the velocity signal is reduced to u0u ¼ g u according to Eq. (7). For this case, which is slightly more involved than DNSu, the causality of acoustic – flow – flame interactions can be represented as g ! g ! u0 ! Q_ 0 . Note that in the two cases out

u

u

discussed the fluctuations of heat release Q_ 0 will generate acoustic

perturbations that propagate away from the flame in both directions. As explained in the introduction, the upstream traveling perturbations will contribute to g u at the upstream location, resulting in ITA feedback. More detailed explanations of the physics involved can be found in the works of Courtine et al. [35], Bomberg et al. [33], Emmert et al. [34] and Silva et al. [48]. 4. DNS of intrinsic thermoacoustic instability This section presents results on ITA instability observed directly in DNS of a laminar premix flame. Two configurations are considered, both with a flame and flame holder that corresponds to the experimental set-up of Kornilov et al. [36]. The first case, here denoted ‘‘A’’, is the one already studied in Section 3. It consists of two cavities of different cross sections as illustrated in Fig. 5(a). One cavity represents the slit of the flame-holding plate, while the adjacent cavity represents the combustion chamber. The second case ‘‘B’’ considered includes an additional plenum cavity upstream of the plate, see Fig. 5(b). Before presenting DNS results for these two configurations, results of linear stability analysis obtained with a network model will be introduced. This network model uses as input the Flame Transfer Function F ðxÞ obtained in Section 3. Note that F ðxÞ is defined throughout the complex frequency plane (not only for purely real frequencies), since it was computed by the z-transform of the impulse response (Eq. (5)). The eigenmodes identified with the network model are useful for the interpretation of the simulation results. DNS results are presented subsequently, with acoustically non-reflecting boundary conditions imposed at inlet and outlet [41]. Both stable as well as unstable modes are observed, which cannot be related to the acoustic modes of the configuration. These acoustic modes, called AC modes in this work, are considered as directly related to the geometry, boundary conditions and mean flow of the system under study, but being independent of the flame dynamics. 4.1. Stability analysis with low order network model Network models are useful tools to determine and interpret eigenmodes of (thermo-)acoustic systems. The models can be

Fig. 5. Computational domains of DNS of laminar premix flame stabilized on a flat plat with and without upstream plenum cavity.

Fig. 4. Impedance Z u for the upstream excitation case (DNSu).

C.F. Silva et al. / Combustion and Flame 162 (2015) 3370–3378

conveniently formulated in terms of characteristic wave amplitudes f and g at the nodes of network elements that comprise the acoustic system. Each element is described by its transfer matrix. In this section, conservation principles are combined with flame transfer function obtained numerically to produce a semi-analytical formulation. Following Polifke et al. [25] or Paschereit et al. [49], the acoustic transfer matrix T for the characteristic waves f ; g between an upstream plane (in) and a downstream plane (out) is defined as



f out



 ¼

g out

T 11

T 12

T 21

T 22



f in

 ð11Þ

g in

Let us consider now the locations (u) and (d) immediately upand downstream of an infinitesimally thin, planar zone of heat release, as illustrated in Fig. 5. Linearization of the quasi-1D Euler equations in the low Mach number regime yields [22]:

^d ¼ p ^u p

ðmomentum equationÞ;

ðc  1Þ ^_ ^ d ¼ A1 u ^u þ A2 u Q cp

ð12Þ

ðmass and entropy equationsÞ

ð13Þ

where A denotes the cross-sectional area. The perturbations of the ^ global heat release rate Q_ is now described by the Flame Transfer

F ð xÞ 

^  Q_ ðxÞ=Q_ ^ u ðxÞ=u u u

ð14Þ

Defining now the area change a ¼ A1 =A2 and relative temperature increment h ¼ ðT d =T u  1Þ, where T denotes here temperature, we obtain a simple expression of Eq. (13):

^d ¼ u ^ u að1 þ hF ðxÞÞ; u

ð15Þ

Recalling the relation (6) between Riemann invariants and primitive acoustic variables

^ ¼ ðf þ gÞq  c; and p

^ ¼f g u

ð16Þ

and by considering Eqs. (12) and (15), we obtain the coefficients of the flame transfer matrix T ðFÞ between planes (u) and (d)

1 ðn þ a þ ahF ðxÞÞ; 2 1 ðFÞ ðFÞ T 12 ðxÞ ¼ T 21 ðxÞ ¼ ðn  a  ahF ðxÞÞ; 2 ðFÞ

ðFÞ

T 11 ðxÞ ¼ T 22 ðxÞ ¼

ð17Þ

where n ¼ qu cu =qd ud denotes the specific acoustic impedance. In ðAÞ

(see order to derive an expression of the transfer matrix T Fig. 5(a)) with respect to the inlet and outlet planes (in) and (out), the phase change due to wave propagation in the adjacent ducts must be included:

" T ðAÞ ¼

#"

e|xl2 =cd

0

0

e|xl2 =cd

T 11

ðFÞ

T 12

ðFÞ

T 22

T 21

ðFÞ

#"

e|xl1 =cu

0

0

e|xl1 =cu

ðFÞ

# :

ð18Þ

The transfer matrix of configuration B (see Fig. 5(b)) is computed as:

T ðBÞ ¼ T ðAÞ

" 1 ð1 þ 1=aÞ e|xl0 =cu 2 ð1  1=aÞ e|xl0 =cu

ð1  1=aÞ e|xl0 =cu ð1 þ 1=aÞ e|xl0 =cu

# ð19Þ

An acoustic network model representation of both configurations A, B result in a 4  4 linear system of equations

2

1

6 6 0 6 6 4 T 11

where T may be either T ðAÞ or T ðBÞ , and Rin and Rout stand for the acoustic reflection coefficients of inlet and outlet, respectively. In order to find a non-trivial solution of the linear system, it is necessary to assure that the solution satisfies detðMÞ ¼ 0. For the systems described in Eq. (20), the determinant of M is easily found, and consequently the characteristic equation to be solved reads

detðMÞ ¼ 0

Rin

0

0

Rout

T 12

1

2 3 0 76 7 6 7 1 7 6 g in 7 6 0 7 76 7¼6 7 76 7 6 7 0 5 4 f out 5 4 0 5

0

32

f in

0 1 T 21 T 22 g out |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} M

3

0

ð20Þ

)

T 22  Rout T 12 þ Rin T 21  Rin Rout T 11 ¼ 0:

ð21Þ

A case similar to configuration A was studied by Emmert et al. [34]. That study describes the acoustics of the system by an open loop scattering matrix SðbcÞ SðFÞ , where SðbcÞ ¼ ½Rin ; 0; 0; Rout  is a 2  2 matrix containing the information on the boundary conditions. The corresponding eigenfrequencies are found by solving the open loop characteristic equation detðSðbcÞ SðFÞ  IÞ ¼ 0. Identical results are obtained for cases Rin ; Rout – 0. However, the formulation used in the present study based on the transfer matrix – which is also used by Hoeijmakers et al. [31] – is more convenient in the limit of vanishing reflection coefficients Rin ; Rout ¼ 0. In that case the characteristic equation resulting from defined by Eq. (21) reduces for case A simply to ðAÞ

T 22 ¼ 0

Function:

3375

()

e|xl1 =cu e|xl2 =cd ðn þ a þ ahF ðxÞÞ ¼ 0:

ð22Þ

The first two terms cannot satisfy e|xl1 =cu ¼ 0 or e|xl2 =cd ¼ 0 and therefore they do not need to be considered in the analysis. The solution of n þ a þ ahF ðxÞ ¼ 0 can be realized and is related to the ITA eigenmodes in complete agreement with the results of Hoeijmakers et al. [31], Emmert et al. [34] and Silva et al. [48]. The characteristic equation for case B with non-reflecting boundary conditions reads ðBÞ

|xl =cu e 1 ð1  1=aÞðn  a  ahF Þ

þ e|xl1 =cu ð1 þ 1=aÞðn þ a þ ahF Þ ¼ 0:

T 22 ¼ 0

()

ð23Þ

In the general case, solutions of Eq. (23) are not easily separable into AC and ITA modes. The AC modes (see definition of AC at the beginning of this section) of the system can be found by setting the flame insensitive to acoustic perturbations (F ¼ 0). Eq. (23) then reduces to

e|xl1 =cu ð1  1=aÞðn  aÞ þ e|xl1 =cu ð1 þ 1=aÞðn þ aÞ ¼ 0:

ð24Þ

For cu ¼ 324 l1 ¼ 1 mm, a ¼ 0:4 and n ¼ 2:81, the first AC eigenfrequency is around 162,000 Hz, which is close to the frequency corresponding to the half-wavelength of the slit. This frequency is very high if compared with the solutions once the flame is present (F – 0), whereas the corresponding growth rate is highly negative. It should be emphasized that if reflecting acoustic boundary conditions are considered, solutions of Eq. (21) for case B, while considering F ¼ 0, still lie above 8000 Hz. Let us consider now Eq. (21) for both A and B systems with an active flame, where F ðxÞ is the one identified by DNSu (see Section 3). An instability map based on different values of Rin and Rout is constructed and shown in Fig. 6. The eigenmodes shown are linked to the ITA feedback, since only frequencies around 150 Hz are present, which cannot be related to any AC mode. For both cases A and B, it is observed that a positive growth rate is predicted for several combinations of Rin and Rout . Indeed, the growth rate of the ITA mode is enhanced when an open-end (Rin ¼ 1) is imposed at the inlet and a closed-end (Rout ¼ 1) is imposed at the outlet. On the contrary, the growth rate is damped when a closed-end and an open-end are considered at the inlet and outlet, respectively. These results exhibit a different trend than the stability analysis of Hoeijmakers et al. [31] for a simple duct-flame-duct configuration. In that study, the stability map shows maximum instability when a closed inlet is combined

3376

C.F. Silva et al. / Combustion and Flame 162 (2015) 3370–3378

Fig. 6. Eigenfrequencies derived from Eq. (21). Rout varies from 1 (black) to 1 (light gray). Symbols in red specify results for Rin , Rout ¼ 0. The dashed line indicates the frequency (147 Hz) which is the solution of n þ a þ ahF ðxÞ ¼ 0. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

with an open outlet. The description of the heat source by a simple n–s model may account for the differences in stability prediction. One particular case is of great interest. It is observed that when (Rin ; Rout ¼ 0) a positive growth rate is predicted for Case B. Although these results are expected in the light of the results of Hoeijmakers et al. [31] and Emmert et al. [34], it is still counter intuitive that instability occurs in this particular situation when a maximum of acoustic energy is evacuated through the boundaries and, accordingly, acoustic losses are at the highest level. 4.2. ITA modes in direct numerical simulation In order to validate the low order stability predictions for cases A and B, DNS are performed. Configuration A is considered first. Since the ITA mode of this case is predicted to be stable, one way to observe this mode in DNS is by exciting the system, which is seen as a stable acoustic oscillator, with a broadband signal of uniform spectral energy. If an acoustic mode of the system lies between the frequency range of forcing, it will be excited and a

resonance peak will be identified. The resulting resonance peak indicates the eigenfrequency of the stable mode. Figure 7(a) illustrates both the input signal f in and output signal f out during a period of 60 ms. It is observed that the output signal is oscillating with a period of roughly 6 ms. The power spectral density (PSD) of both signals is shown in Fig. 8(a): whereas the energy of the signal is uniformly distributed for the input, the spectrum of the output signal shows a resonance peak at 151 Hz, in good agreement with the mode predicted at 147 Hz by the network model analysis. Figure 7(b) illustrates both input and output signals for case B. Here the input signal f in contains very low energy and is only applied during 0:005 s at the beginning of the simulation in order to trigger the instability. The output signal f out exhibits clearly a high amplitude oscillation at a frequency of 164 Hz, which is in fair agreement with the ITA flame eigenfrequency of 152 Hz predicted by the network model analysis. Figure 8(b) shows the PSD for input and output signals of case B. It is interesting to remark that network model predictions are more accurate for Case A. One possible

Fig. 7. Time series of input signal f in (black) and output signal f out (gray). (a) Case A. (b) Case B.

C.F. Silva et al. / Combustion and Flame 162 (2015) 3370–3378

3377

Fig. 8. Power Spectral Density of Input signal f in (black) and output signal f out (gray). (a) Case A. (b) Case B.

 Fig. 9. Signals of interest. (black) u0in , (gray) Q_ 0 =Q_ , (dashed) p0out =ðqout cout Þ. Note that since f in ¼ 0 and g out ¼ 0 (totally non-reflecting acoustic boundary conditions), then uin ¼ g in and p0out =ðqout cout Þ ¼ f out .

explanation relies on the linearity assumption in which the network model is based. Whereas acoustics in Case A still can be described by a linear model, an accurate prediction of the limit cycle frequency of Case B might need the inclusion of non-linearities in the description of the flame dynamics. In order to verify that the thermal source of acoustic energy production is positive, according to the Rayleigh criterion [5], we superpose now the normalized signals of Q_ 0 and p0 in Fig. 9. Note out

p0out

is indeed very close to the fluctuating pressure throughout that the domain due to the large value of the wavelength with respect to the length of the system. In addition, the signal u0in , which exhibits the same phase of u0u , is also present in order to visualize the phase shift between u0 and Q_ 0 . It is observed that Q_ 0 and p0 are in

simplistic models, but a type of combustion instability that until recently has escaped the attention of the combustion dynamics research community and warrants further study. Acknowledgments The authors gratefully acknowledge CERFACS for making available the CFD solver AVBP, and the Research Association for Combustion Engines (Forschungsvereinigung Verbrennungskraftmaschinen e.V. FVV) for the Financial support to Stefan Jaensch. We are also indebted to Leibniz Rechenzentrum LRZ for providing access to HPC resources (SuperMUC).

out

perfectly in phase and therefore the Rayleigh criterion for thermoacoustic energy production is verified. Furthermore, it is seen that u0 and Q_ 0 are completely out of phase (/ ¼ p) as expected from in

the analysis in the previous section. These results agree with the observations of Courtine et al. [35], Emmert et al. [34], Bomberg et al. [33] and Hoeijmakers et al. [31] and confirm that the instability is thermo-acoustic in nature. 4.3. Conclusions We have demonstrated, by Direct Numerical Simulation, that intrinsic thermoacoustic feedback is an authentic physical phenomenon that is present in premixed combustion systems. The corresponding ITA modes are in the present case independent of acoustic (AC) modes. Both stable and unstable modes have been observed. We have given evidence, by an adequate identification of the FTF, that compact laminar premixed flames are velocity sensitive. This corroborates the structure of the flow-flame-acoustic interaction with ITA feedback suggested by Bomberg et al. [33]. Finally, we have shown that thermoacoustic network models in conjunction with an FTF identified from DNS are capable of predicting stability/instability of this kind of acoustic mode. This notion might seem obvious, but has been questioned in the context of intrinsic thermoacoustic instabilities. In summary, the present study gives strong evidence that intrinsic thermoacoustic instabilities are not an artifact of

References [1] H.S. Tsien, The transfer functions of rocket nozzles, J. Am. Rocket Soc. 22 (3) (1952) 139–143. [2] L. Crocco, S.I. Cheng, Theory of Combustion Instability in Liquid Propellant Rocket Motors, Cambridge Univ Press, 1956. [3] S. Candel, Combustion dynamics and control: progress and challenges, Proc. Combust. Inst. 29 (2002) 1–28. [4] T. Lieuwen, V. Yang (Eds.), Combustion Instabilities in Gas Turbine Engines: Operational Experience, Fundamental Mechanisms and Modeling, Prog. in Astronautics and Aeronautics AIAA 210. [5] L. Rayleigh, The explanation of certain acoustic phenomena, Nature 18 (1878) 319–321. July. [6] H. Büchner, C. Hirsch, W. Leuckel, Experimental investigations on the dynamics of pulsated premixed axial jet flames, Combust. Sci. Technol. 94 (1–6) (1993) 219–228. [7] S. Ducruix, D. Durox, S. Candel, Theoretical and experimental determinations of the transfer function of a laminar premixed flame, Proc. Combust. Inst. 28 (2000) 765–773. [8] K. Kunze, C. Hirsch, T. Sattelmayer, Transfer function measurements on a swirl stabilised premix burner in an annular combustion chamber, in: Int’l Gas Turbine and Aeroengine Congress & Exposition, No. ASME GT-2004-54106, Vienna, Austria, 2004. [9] B. Schuermans, F. Bellucci, V. Guethe, F. Meili, P. Flohr, C.O. Paschereit, A detailed analysis of thermoacoustic interaction mechanisms in a turbulent premixed flame, in: Int’l Gas Turbine and Aeroengine Congress & Exposition, No. ASME GT2004-53831, Atlanta, GA, U.S.A., 2004. [10] D. Durox, T. Schuller, N. Noiray, S. Candel, Experimental analysis of nonlinear flame transfer functions for different flame geometries, Proc. Combust. Inst. 32 (2009) 1391–1398. [11] K. Kim, J. Lee, B. Quay, D. Santavicca, Spatially distributed flame transfer functions for predicting combustion dynamics in lean premixed gas turbine combustors, Combust. Flame 157 (2010) 1718–1730. [12] F. Duchaine, F. Boudy, D. Durox, T. Poinsot, Sensitivity analysis of transfer functions of laminar flames, Combust. Flame 158 (12) (2011) 2384–2394.

3378

C.F. Silva et al. / Combustion and Flame 162 (2015) 3370–3378

[13] S. Hermeth, G. Staffelbach, L.Y.M. Gicquel, V. Anisimov, C. Cirigliano, T. Poinsot, Bistable swirled flames and influence on flame transfer functions, Combust. Flame 161 (1) (2014) 184–196. [14] W. Polifke, A. Poncet, C.O. Paschereit, K. Döbbeling, Reconstruction of acoustic transfer matrices by instationary computational fluid dynamics, J. Sound Vib. 245 (3) (2001) 483–510. [15] M. Zhu, A.P. Dowling, K.N.C. Bray, Transfer function calculations for aeroengine combustion oscillations, J. Eng. Gas Turb. Power 127 (2005) 18–26. [16] A.M.G. Gentemann, C. Hirsch, K. Kunze, F. Kiesewetter, T. Sattelmayer, W. Polifke, Validation of flame transfer function reconstruction for perfectly premixed swirl flames, in: Int’l Gas Turbine and Aeroengine Congress & Exposition, No. ASME GT-2004-53776, Vienna, Austria, 2004. [17] R. Kaess, W. Polifke, T. Poinsot, N. Noiray, D. Durox, T. Schuller, S. Candel, Cfdbased mapping of the thermo-acoustic stability of a laminar premix burner, in: Proc. of the Summer Program, Center for Turbulence Research, NASA AMES, Stanford University, USA, 2008, pp. 289–302. [18] A. Cuquel, D. Durox, T. Schuller, Experimental determination of flame transfer function using random velocity perturbations, Proceedings of ASME Turbo Expo 2011, vol. GT2011-54624, ASME, 2011, pp. 793–802. [19] L. Tay-Wo-Chong, S. Bomberg, A. Ulhaq, T. Komarek, W. Polifke, Comparative validation study on identification of premixed flame transfer function, J. Eng. Gas Turb. Power 134 (2) (2012) 021502–1–8. [20] R.S. Blumenthal, P. Subramanian, R. Sujith, W. Polifke, Novel perspectives on the dynamics of premixed flames, Combust. Flame 160 (7) (2013) 1215–1224. [21] A.P. Dowling, The calculation of thermoacoustic oscillations, J. Sound Vib. 180 (4) (1995) 557–581. [22] T. Poinsot, D. Veynante, Theoretical and numerical combustion, R.T. Edwards, 2005. [23] N. Noiray, D. Durox, T. Schuller, S. Candel, A unified framework for nonlinear combustion instability analysis based on the flame describing function, J. Fluid Mech. 615 (2008) 139–167. [24] J.J. Keller, Thermoacoustic oscillations in combustion chambers of gas turbines, AIAA J. 33 (12) (1995) 2280–2287. [25] W. Polifke, C.O. Paschereit, K. Döbbeling, Constructive and destructive interference of acoustic and entropy waves in a premixed combustor with a choked exit, Int. J. Acoust. Vib. 6 (3) (2001) 135–146. [26] C. Pankiewitz, T. Sattelmayer, Time domain simulation of combustion instabilities in annular combustors, Trans. ASME, J. Eng. Gas Turb. Power (125) (2003) 677–685. [27] F. Nicoud, L. Benoit, C. Sensiau, Acoustic modes in combustors with complex impedances and multidimensional active flames, AIAA J. 45 (2007) 426–441. [28] C.F. Silva, F. Nicoud, T. Schuller, D. Durox, S. Candel, Combining a Helmholtz solver with the flame describing function to assess combustion instability in a premixed swirled combustor, Combust. Flame 160 (9) (2013) 1743–1754. [29] W. Polifke, Black-box system identification for reduced order model construction, Ann. Nucl. Energy 67 (2014) 109–128. [30] S.W. Yuen, A.M.G. Gentemann, W. Polifke, Investigation of the influence of boundary conditions on system identifiability using real time system modeling, in: 11th Int. Congress on Sound and Vibration (ICSV11), IIAV, Saint-Petersburg, Russia, 2004, pp. 3501–3508. [31] M. Hoeijmakers, V. Kornilov, D.P.D.G.I. Lopez Arteaga a, H. Nijmeijer, Intrinsic instability of flame-acoustic coupling, Combust. Flame 161 (11) (2014) 2860– 2867.

[32] M. Hoeijmakers, I. Lopez Arteaga, V. Kornilov, H. Nijmeijer, P. de Goey, Experimental investigation of intrinsic flame stability, in: Proceedings of the European Combustion Meeting, Lund, Sweden, 2013. [33] S. Bomberg, T. Emmert, W. Polifke, Thermal versus acoustic response of velocity sensitive premixed flames, Proc. Combust. Inst. 35 (3) (2015) 3185– 3192. [34] T. Emmert, S. Bomberg, W. Polifke, Intrinsic thermoacoustic instability of premixed flames, Combust. Flame 162 (1) (2015) 75–85. [35] E. Courtine, L. Selle, F. Nicoud, W. Polifke, C.F. Silva, M. Bauerheim, T. Poinsot, Causality and intrinsic thermoacoustic instability modes, in: Proc. of the Summer Program, Center for Turbulence Research, Stanford University/NASA Ames, USA, 2014. [36] V. Kornilov, R. Rook, J. ten Thije Boonkkamp, L. de Goey, Experimental and numerical investigation of the acoustic response of multi-slit Bunsen burners, Combust. Flame 156 (10) (2009) 1957–1970. [37] P. Clavin, Dynamics of combustion fronts in premixed gases: from flames to detonations, Proc. Combust. Inst. 28 (1) (2000) 569–585. [38] L. Ljung, System Identification – Theory for the User, second ed., Prentice Hall, 1999. [39] S. Föller, W. Polifke, Identification of aero-acoustic scattering matrices from large eddy simulation: application to a sudden area expansion of a duct, J. Sound Vib. 331 (13) (2012) 3096–3113. [40] T. Poinsot, S. Lele, Boundary conditions for direct simulation of compressible viscous flows, J. Comput. Phys. 101 (1) (1992) 104–129. [41] W. Polifke, C. Wall, P. Moin, Partially reflecting and non-reflecting boundary conditions for simulation of compressible viscous flow, J. Comput. Phys. 213 (1) (2006) 437–449. [42] J. Kopitz, E. Bröcker, W. Polifke, Characteristic-based filter to identification of acoustic waves in numerical simulation of turbulent compressible flow, in: 12th Int. Congress on Sound and Vibration (ICSV12), No. 389, Lisbon, Portugal, 2005. [43] L. Selle, G. Lartigue, T. Poinsot, R. Koch, K.U. Schildmacher, W. Krebs, B. Prade, P. Kaufmann, D. Veynante, Compressible large-eddy simulation of turbulent combustion in complex geometry on unstructured meshes, Combust. Flame 137 (4) (2004) 489–505. [44] M. Schmid, R. Blumenthal, M. Schulze, W. Polifke, T. Sattelmayer, Quantitative stability analysis using real frequency response data, J. Eng. Gas Turb. Power 135 (12) (2013) 121601. [45] S. Föller, W. Polifke, Advances in identification techniques for aero-acoustic scattering coefficients from large eddy simulation, in: International Congress on Sound and Vibration, ICSV18, 2011. [46] A.K. Tangirala, Principles of System Identification: Theory and Practice, Taylor & Francis, 2014. [47] W. Polifke, System identification for aero- and thermo-acoustic applications, in: V.L. 2011-01 (Ed.), Advances in Aero-acoustics and Thermo-acoustics, Von Karman Institute, Brussels, BE, 2010. [48] C.F. Silva, S. Jaensch, T. Emmert, W. Polifke, On the autoregressive behavior of the intrinsic thermoacoustic feedback loop observed in premixed flames, in: 22st International Congress of Sound and Vibration (ICSV22), Florence, Italy, 2015. [49] C.O. Paschereit, B. Schuermans, W. Polifke, O. Mattson, Measurement of transfer matrices and source terms of premixed flames, J. Eng. Gas Turb. Power 124 (2002) 239–247.