Three-dimensional nonlinear thermoacoustic instability analysis based

on Green’s function approach Weipeng Zhou 1 Research Institute of Aero-Engine, Beihang University 100191, Beijing, China Shuo Tian School of Energy and Power Engineering, Beihang University 100191, Beijing, China Guangyu Zhang Research Institute of Aero-Engine, Beihang University 100191, Beijing, China Xiaoyu Wang Research Institute of Aero-Engine, Beihang University 100191, Beijing, China Xiaofeng Sun School of Energy and Power Engineering, Beihang University 100191, Beijing, China

ABSTRACT In this work, nonlinear thermoacoustic instability analysis is made in a hard-wall box to investigate the limit cycle of thermoacoustic oscillation and the interplay of two modes. The flame is assumed to be compact, and described by a generic flame describing function, i.e. an amplitude-dependent flame transfer function. The acoustic field in the hard-wall box could be described by an integral equation using a Green’s function tailored to a three-dimensional (3- D) rectangular box with hard-wall boundary conditions. Two different methods are used to solve this equation. The first is an iteration method to obtain the time history of the acoustic velocity at the flame. The other method is done in the frequency domain to determine the thermoacoustic eigenfrequency and growth rate of thermoacoustic modes. We observed the phenomenon of modes “jump” after two different time windows of the time history operated with the Fourier transformation. This result reveals that when there are two modes in the system, the interference between modes will occur. In the early stage, the two modes work together. After reaching the limit cycle oscillation, the role of one mode is dominant and the other mode basically disappears. 1. INTRODUCTION ：

Thermoacoustic instability, driven by constructive coupling between unsteady heat release and acoustic waves at one or multiple frequencies [1], has long been an urgent problem in the field of combustion. It is widely observed in boilers, furnaces, ground gas turbines, aero-engines, rocket engines and other devices with flame burning continuously in the cavity [2, 3]. Thermoacoustic instabilities will trigger the combustion device even the whole system to vibrate strongly, resulting in larger noise and heat load, higher pollutant emission ratio, affecting the normal operation of the system and even causing the failure of the whole combustion system in serious cases. In order to understand the mechanism of thermoacoustic instabilities and develop suitable control method with affordable time and computation cost, analytical tools that are able to predict the thermoacoustic instabilities including the stability and limit cycle are needful. 1 weipengz@buaa.edu.cn

The present method is based on the tailored Green’s function (the acoustic field that is generated by an impulsive point source in a cavity). It’s firstly used by Heckl and Howe [4, 5] to model 1-D thermoacoustic systems. Wang and Heckl [6] then developed this method to 3-D thermoacoustic systems in order to study the 3-D effects, as is shown in figure 1. But the linear heat release model used in her model fails to capture the important thermoacoustic oscillation feature--limit cycle, which is important in studying thermoacoustic instabilities. Because we need only control the limit cycle amplitude in acceptable range rather than eliminate thermoacoustic instabilities completely when controlling thermoacoustic instabilities. Therefore, we use the nonlinear heat release model in order to capture the limit cycle in this paper as an extension of Wang and Heckl [6].

(a) «, VA t peat a — 4,—

Figure 1: Schematic of the rectangular combustion chamber; (a) 3-D geometry, (b) geometry for

case of no z-dependence [6].

2. THE INTEGRAL GOVERNING EQUATION AND SOLUTION

2.1 The integral governing equation

An integral governing equation giving the velocity at the flame in terms of the heat release is given

    −  −     = − −       

t

G r r t t G u t q t dt c r c r t  

2 0 2 2 0 0

1 ( , , ) ( ) ( )

(1)

q r r t r r r r

= =

q

  = =

q q

=

r r

q

with the two initial condition

r r t    =

 = = =  

(2)

0 , 0 at t q t

0 0

=

where c is speed of sound and  is specific heat ratio. 0  is the initial value of the velocity

potential at the flame. ( , , ) G r r t t   −

is the 3-D tailored Green function of a rectangular hard-wall cavity with side length x H , y H , z H . Its meaning is the response of an observer point ( , , ) r x y z =

in the cavity to an impulsive point source at ( , , ) r x y z    =

. t t  − is the time it takes the signal to travel from r 

to r

.

mnk mnk j t t j t t mnk m n k G r r t t g r r e e   

1 ( , , ) ( , ) [ ] 2

=    − = − 

  − − −

( ) ( )

(3)

, , 0

    =   

2 ( , ) ( ) ( ) m n k mnk mnk mnk mnk mnk g r r c r r j

(4)

w here

  

m n k r k x k y k z k k k H H H

 = = = =

( ) cos( )cos( )cos( ); , , ;

mnk x y z x y z x y z

(5)

=  = + + = =  =  

 

1 if 0 ; ; ; . 2 if m >0

m k k k k k H H H c

2 2 2

0

mnk x y z m mnk x y z

, , x y z k k k are the allowed wave numbers of the cavity in the , , x y z direction respectively, and

, , m n k are the corresponding mode numbers. The detailed deviation about the integral equation and the Green function of a rectangular hard-wall cavity can be found in Wang and Heckl [6].

The flame is assumed to be compact and locates at q r r =

and the flame is flat and its surface

normal is described by the angle  and  , as is shown in figure 1. We adopt the heat release model fitted by Heckl [7] according to experimental data.

1 0 ( ) [ ( ) ( )] ( ) q r r q t K n u t n u t r r      = − − −

(6)

where Q K US  = is the heater power per mass flow, Q is the mean value of global heat release

rate, U is the mean flow velocity,  is the mean density. The coupling coefficients and delay time is described as follows.

2

0 1 0 1 0 1 1 1 ; ; ; 2 2 A A g g g g g n n U U    − +   = + = − = =    

(7)

where 3 3 0 1 0 1 5 10 , 4.4 10 , 1.4, 0.3 g g   − − =  =  = = that are measured experimentally. A is

acoustic velocity amplitude. U is the background mean velocity.

2.2 Solution of the integral governing equation

2.2.1 Numerical iteration method

Equation (1) can be change to

      − − − 

t j t t j t t

  −     = − − −       

( ) ( ) 2 0 2 2 0 , , 0 0

1 ( ) [ ] ( ) 2 2

e e G u t G q t dt c c r t

mnk mnk

(8)

q mnk t m n k r r r r

= = =

q

=

q

where

  = 

( , )

g r r G r = =

mnk mnk r r r r

(9)

q

q

We adopt the numerical iteration method described by Heckl [8].

   − − −    = =  =

       = + 

  

j t j t j t t t t t t t t

e e e I t q t dt q t dt q t dt

mnk mnk mnk

-

  

( ) ( ) = ( ) ( ) 2 2 2



1 0 0 -

n t t t t t

   = = = 

(10)

      = =  =

  

j t j t j t t t t t t t t

e e e I t q t dt q t dt q t dt

 

mnk mnk mnk

-

  

      = +

( ) = ( ) ( ) ( ) 2 2 2

2 0 0 -

n t t t t t

   = = = 

( ) q t  is assumed to be a constant in the time interval t  , so

( ) ( )

 −   =  +  −  −     + 

 

− 

j t j t

1 - ( ) ( - ) 2

e e q t t I t I t t j

mnk mnk

1 1

n n

  

mnk

(11)

( ) ( )

− 

j t j t

1 - ( ) = ( - ) 2

e e q t t I t I t t j

mnk mnk

2 2

n n



mnk

and

   = − − = −       

  −

2

1 ( ) cos( ) 2

G G j e e G j t r t

    

j t j t mnk mnk mnk mnk mnk t m n k m n k r r r r

(12)

mnk mnk

= = = =

0 , , 0 , , 0

q

=

q

substitute the above two equations into equation (8)

      

−  = − − − − −   (13)

  − −

1 1 ( ) [ ( ) ( )] ( ) 2

mnk mnk mnk mnk j t j t j t j t q mnk n n mnk mnk m n k m n k u t G e I t e I t G j e e c c

0 1 2 2 2 , , 0 , , 0

= =

Later, we can obtain the time history of the acoustic velocity by the iteration method.

2.2.2 frequency domain method

Laplace transform is performed on equation (1). The result is



1 1 1 ( ) ( )( ) ( ) 2



  = − − − − +

  

−

s q q mnk q m n k mnk mnk

u s c G j n e n u s j s j s j

1 0 , , 0

=

(14)

   



1 1 1 ( ) 2



+ + − +

G j c s j s j

0 2

mnk mnk m n k mnk mnk

=

, , 0 0

q c K c − = .

1

where 2

The equation of ( ) q u s  can be solved from it.

   



1 1 1 ( ) 2 ( ) 1 1 1 [1 ( )( )] 2



+ − +  =

G j c s j s j u s

0 2

mnk mnk m n k mnk mnk q

= 

, , 0 0

(15)





−

+ − − − +

s q mnk m n k mnk mnk

c G n e n s j s j

1 0 , , 0

 

=

The equation of ( ) q u t  will be obtained after applying inversion of Laplace transform to ( ) q u s  . It

can be solved with residue theorem. In the process of solving, we can obtain an equation about complex eigenfrequency.



1 1 1 1 ( )( ) 2

= = − − − − +  (16)

  

−

s q mnk m n k mnk mnk c G n e n s j s j

1 0 , , 0

where s i =  . It can be solved by Newton iterative method. Ω is the complex eigenfrequency of the system. The real part of Ω denotes thermoacoustic oscillation frequency; the imaginary part of Ω denotes growth rate. When Im(Ω) is positive, the system is stable; when Im(Ω) is negative, the system is unstable with the time dependence

i t e  . 3. MODEL VALIDATION

3.1 Comparison with a 1-D network model

We reduce our model to 1-D situation with 1 x H m = , which can be analyzed by a simple 1-D network model. The 1-D network model is described in Wang and Heckl [6]. The difference is that the n  - law is replaced by the nonlinear heat release model described above. An equation about complex eigenfrequency obtained from the network model is

    − − − − − −     − −+ − −            

( ) 0 0 0 0 2 2 2 2

( )

q q i x i L i L x i L c c c c i q c n e n e e e e 

1 0 1 +2 1 =0

(17)

It can be solved by Newton iterative method to obtain complex eigenfrequency which is compared to that obtained by Green function method.

Figure 2 shows Ω as a function of A/U (ratio of the acoustic velocity amplitude to background velocity): Re(Ω) is shown in sub-figure(a) and Im(Ω) is shown in sub-figure(b). The black curves denote the results of Green’s function method; the red curves denote the results of network model. We could conclude from figure 2 that the results of two models are basically consistent, so the present model is effective to predict thermoacoustic instability. From sub-figure(b) of figure 2, we can conclude that system is unstable for 0 / 0.46 A U   , stable for 0.46 / 0.92 A U   , unstable for 0.92 / 1.25 A U   , stable for 1.25 / 1.48 A U   . The magnitude of velocity disturbance at the flame is called the limit cycle when the system transitions from unstable state to stable state. As the name implies, it is the state oscillating with a fixed amplitude. When the system transitions from stable state to unstable state, the magnitude of velocity disturbance at this time isn’t the limit cycle but the cut- off point of the next limit cycle. So, the red numbers 0.46, 1.25, 1.70 in figure 2(b) denote the limit

cycle. Next, we will give the time history of velocity fluctuations using the numerical iteration method to verify our prediction.

Green function method Network model

1080

Green function method Network model

10

1075

5

Re  [s -1 ]

Im  [s -1 ]

1070

1.48

0.92

0

1.70

1.25

0.46

1065

-5

1060

0.0 0.5 1.0 1.5 2.0 -10

0.0 0.5 1.0 1.5 2.0

A/U

A/U

(a) Oscillation frequency (b) Growth rate

Figure 2: Comparison of the complex eigenfrequency at different velocity disturbance calculated from the Green’s function method (black curve) and from network model (red curve). The flame is located at 0.25m, 0.05 q c = ; (a) real part (oscillation frequency), (b) imaginary part (growth rate).

3.2 Comparison of the two solving methods

In linear situation, we can compare the time history obtained by the two methods. But in nonlinear situation, the method done in frequency domain can’t give the time history of velocity disturbance because the thermoacoustic eigenfrequency is changed simultaneously with A/U. So, we compare the limit cycle obtained by the two methods.

Figure 3 ： the time history of acoustic velocity at the flame. The initial perturbation amplitude is (a) 0.01, (b) 0.62, (c) 0.94 and (d) 1.35.

(a) ous 0) 03 os 0.0. S00. 03 04 os! . 08 0 os To 00 os Yo 13) 108) 810 Ls 10 00 os ro 16) 1s)

Figure 3 shows the time history of velocity disturbance with different initial velocity disturbance

0 / q t u U =  . As can be seen from the figure 3, when the initial velocity disturbance is 0.01 and 0.94,

the system is unstable, and the velocity disturbance increases, and finally reaches the limit cycle 0.46 and 1.25 respectively. When the initial velocity disturbance is 0.62 and 1.35, the system is stable, the velocity disturbance decreases and reaches the limit cycle 0.46 and 1.25 respectively. The amplitudes of the limit cycle obtained from the numerical iteration method agree with that predicted by growth rate. 4. RESULTS AND DISCUSSION FOR A TWO-FREQUENCY RESONATOR

We now consider the case where two frequencies are present in the resonator so that interference effects can be expected to occur. Because the nature of the interference depends on whether or not the frequencies are close to each other [6], we therefore treat the two cases in separate sections.

4.1 Frequencies not close to each other

The box is assumed to be rectangular with different side lengths: 1 x H m = , 1.5 y H m = . The mode

(1,0,0) and (0,1,0) then have the frequencies 1 100 1068 s  − = and 1 010 712 s  − = .

4.1.1 Complex eigenfrequency calculations

(a) Oscillation frequency (b) Growth rate Figure 4: complex eigenfrequency as a function of velocity disturbance. 0.1 q q x y m = = , 45 =

,

0.05 q c = .

In this situation, Re(Ω 100 ) fluctuates within the range [1065s -1 ,1071s -1 ]; Re(Ω 010 ) fluctuates within the range [710s -1 ,715s -1 ]. It can be observed from figure 4(b) that two growth rate curves are interlaced and the amplitude of the limit cycle can no longer be obtained by this method when there are two modes simultaneously. Therefore, the amplitude of the limit cycle can only be obtained by numerical integration method.

1072. 1070. Foes. g = 1066 =a] oo 00S) ns: 74 B73 $n2. mM ‘no. ‘00 os 10 S20 AU

4.1.2 Time history calculation

We could conclude from figure(5) that the system is unstable and velocity disturbance increases to the limit cycle of 0.93 when the initial velocity disturbance is 0.85. We take the fast Fourier transform (FFT) to two different time windows of time history, [0s-1s] and [9s-10s], corresponding to the initial oscillation and the limit cycle oscillation respectively, and the result is shown in figure 5(b). It is obvious that the FFT of the time window [0s-1s] has two different frequencies, 710s -1 and 1068.2s -1 , corresponding to mode (0,1,0) and mode (1,0,0) respectively. While the FFT of the time window [9s- 10s] has only one frequency 710s -1 , corresponding to mode (0,1,0). Therefore, there will be interference between two modes when there are two modes in the box at the same time. In the early stage, the two modes work together. After reaching the limit cycle oscillation, the role of one mode is dominant and the role of the other mode basically disappears. This phenomenon is called mode “jump”.

1.0

0.93

0.5

u  q / U

0.0

-0.5

-1.0

0 2 4 6 8 10

2100, 0 0S OS AU 2.0

t (s)

(a) Time history (b) Frequency spectrum Figure 5 ： the time history of acoustic velocity and initial condition 0.85. 0.1 q q x y m = = , 45 =

,

0.05 q c = .

4.2 Frequencies close to each other

The box is assumed to have very similar side length: 1 x H m = , 1.01 y H m = . The mode (1,0,0)

and (0,1,0) therefore have similar frequencies: 1 100 1068 s  − = and 1 010 1057 s  − = .

4.2.1 Complex eigenfrequency calculations

(a) Oscillation frequency (b) Growth rate Figure 6 ： thermoacoustic eigenfrequency as a function of velocity disturbance. 0.1 q q x y m = = ,

45 =

, 0.05 q c = .

We could find from figure 6(a) that Re (Ω 100 ) fluctuates within the range [1065s -1 ,1075s -1 ], while Re (Ω 010 ) fluctuates within the range [1052s -1 ,1060s -1 ]. About figure 6(b), we could draw a same conclusion as the rectangular box with definitely different side lengths.

0.4 1068.26" 03 B02 ‘Tos! 201 oo! 500 10001500 2000 10 70s" =o 3 2 Bos 2 00 0 500 1000 1300 2000 als")

4.2.2 Time history calculations

We can see from figure 7(a) that the system is stable and velocity disturbance decreases to the limit cycle of 0.47 when the initial velocity disturbance is 0.85. In a nearly square box, there will be a phenomenon of “beating” in linear situation when frequencies of two modes is close to each other [6]. But as can be seen from figure 7(a), the phenomenon of beating occurs at the beginning, but disappears at the limit cycle oscillations. This is because when there are two modes in the box in the nonlinear situation interference between two modes will occur. In the early stage, the two modes work

together, so the phenomenon of beating occurs. In the stage of limit cycle oscillations, the role of one mode is dominant and the role of the other mode basically disappears, so the phenomenon of “beating” disappears.

1.0

0.47

0.5

u  q /U

0.0

-0.5

0 2 4 6 8 10 -1.0

t (s)

(a) Time history (b) Frequency spectrum Figure 7 ： the time history of acoustic velocity and initial condition 0.85. 0.1 q q x y m = = , 45 =

,

0.05 q c = .

5. CONCLUSION

In this paper, the three-dimensional theoretical model based on the Green's function is developed to consider nonlinear effects by using the nonlinear heat release model. The differences from the linear heat release model are that it contains the instantaneous acoustic velocity except the acoustic velocity at an earlier time  , and the delay time  and coupling coefficients are dependent on the acoustic velocity rather than constant. So, it can capture the characteristic of limit cycle. The cavity is described by its tailored Green’s function. By combining it with the acoustic analogy equation of velocity potential and through a series of algebraic transformations, the expression of the acoustic velocity at the flame can be obtained. We solve the equation in two ways:

1085684 Sone gos “500-750-1009 12501500, ofs'] 04 1088.66" = o.i05 FFT] $03 202 Zo 6. 1300 zg ss 23

(1) Numerical iteration method. The integration equation was solved directly with a numerical time stepping approach to obtain the time history of the acoustic velocity and from it we could obtain the stability and limit cycle of the system.

(2) Application of Laplace transform to obtain a function about complex eigenfrequency. And it can be solved using Newton iteration method. The real part of complex eigenfrequency denotes oscillation frequency and the imaginary part of complex eigenfrequency denotes the growth rate (the time dependence assumption

 is adopted, so the system is stable when imaginary part is positive and the system is unstable when imaginary part is negative). By performing case studies, we could draw the following conclusions.

j t e

(1) When there is only one mode in the system, both methods can tell the stability of the system and obtain the limit cycle of the system, and the limit cycle obtained by the two methods agrees well.

(2) Interference between modes occurs when there are two modes in the system. The method by solving complex eigenfrequency can no longer obtain the limit cycle of the system. It can only be obtained by numerical integration method. Fast Fourier transform is performed for two different time windows of the time history obtained by numerical integration method. It can be found that, due to the interference between modes, the two modes work together in the early stage, and gradually evolve into the dominant role of one mode in the later stage, while the role of the other mode basically disappears.

(3) When the frequencies of the two modes in the cavity are very close, the phenomenon of "beating" will appear in the linear case. But in the nonlinear case, due to the interference between the modes, the phenomenon of "beating" will appear in the early stage, and it will disappear when the

limit cycle oscillation is reached. 6. ACKNOWLEDGEMENTS

The study was supported by the National Natural Science Foundation of China (Grant Nos.51790514 and 52106038). 7. REFERENCE

1. Rayleigh J S W. The Theory of Sound[M]. NewYork: Dover Publications, 1945.

2. Crocco, L., and Cheng, S., Theory of Combustion Instability in Liquid Propellant Rocket Motors,

Butterworths Scientic, London, 1956.

3. Keller, J. J., “Thermoacoustic Oscillations in Combustion Chambers of Gas Turbines,” AIAA

Journal, Vol. 33, No. 12, 1995, pp. 2280–2287.

4. Heckl, M.A. & Howe, M.S. (2007b) Stability analysis of the Rijke tube with a Green's function

approach. J. Sound Vib. 305, 672-688.

5. Heckl, M.A. & Howe, M.S. (2007a) The Rijke tube: Green's function approach in the time and

frequency domain. Proceedings of the 14th International Congress on Sound and Vibration, Cairns, Australia, 9-12 July 2007.

6. Wang X, Heckl M. 3-D thermoacoustic instability analysis based on Green's function approach[J].

Journal of Sound and Vibration, 2022: 116816.

7. Heckl, Maria A. Analytical model of nonlinear thermo-acoustic effects in a matrix burner[J].

Journal of Sound & Vibration, 2013, 332(17):4021-4036.

8. Heckl, Maria, A, et al. A Green's function approach to the rapid prediction of thermoacoustic

instabilities in combustors[J]. Journal of Fluid Mechanics, 2016.