Thermoacoustic Instability in a Rijke Tube: Impact of Linear Coupling and Stochastic Sources
A Project Report submitted by
KUSHAL SHARAD KEDIA in partial fulfilment of the requirements for the award of the degrees of
BACHELOR OF TECHNOLOGY and MASTER OF TECHNOLOGY (DUAL DEGREE)
DEPARTMENT OF AEROSPACE ENGINEERING INDIAN INSTITUTE OF TECHNOLOGY, MADRAS. MAY 2008
CERTIFICATE
This is to certify that the project report entitled Thermoacoustic Instability in a Rijke Tube: Impact of Linear Coupling and Stochastic Sources, submitted by Kushal Sharad Kedia, to the Indian Institute of Technology, Madras, for the award of the degrees of Bachelor of Technology and Master of Technology (Dual Degree), is a bonafide record of the research work done by him under our supervision. The contents of this thesis, in full or in parts, have not been submitted to any other Institute or University for the award of any degree or diploma.
————————Prof. R. I. Sujith Research Guide Professor Dept. of Aerospace Engineering IIT Madras Chennai 600 036 Place: Chennai, India Date: ___ May 2008
————————Prof. Job Kurian Head of the Department Professor Dept. of Aerospace Engineering IIT Madras Chennai 600 036
DEDICATION To my parents, who gave me freedom to venture out in the world, and my sister, Ashita, whose incessant support and encouragement redefined the upper bounds of my passion to pursue my interests.
i
ACKNOWLEDGEMENTS With a deep sense of gratitude, I wish to express my sincere thanks to my research guide, Prof. R. I. Sujith, for his immense help in planning and executing the works in time. The confidence and dynamism with which Prof. Sujith guided the work requires no elaboration. Profound knowledge and timely wit came as a boon under his guidance. My sincere thanks are due to Prof. J. Kurian, head of the department, for providing me the initial push and enthusiasm towards scientific research. I specially thank Dr. Sunetra Sarkar for the useful discussions on stochastic process. I am thankful to Prof. M. Ramakrishna for his highly motivating and enthusiastic way of teaching fluid mechanics and numerical methods. The cooperation I received from other faculty members of this department is also gratefully acknowledged. I wish I would never forget the company I had from my fellow students and research scholars of my department and some friends of my hostel. In particular, I am thankful to Sharath for the endlessly long brain storming sessions without which my research would have been incomplete. I would also like to mention Nishant, Priya, Rahul, Sathesh, Koushik, Bharath, Vikrant, Ramya and Vineet for the discussions that helped me strengthen my understanding in the challenging field of thermacoustics. I would like to specially thank Vainatha for her constant moral support and encouragement. I would also like to thank Jatin, Lal, Sharath, Nishant, Anshuman, Arnab, Pahwa, Somani, Abhinav, Manuj, Shyam, Shruthi and the endless list of my friends who have made my life at IIT worth cherishing forever. The acknowledgement cannot be complete without thanking my parents Shri. Sharad Kedia and Smt. Meena Kedia, who taught me the value of hard work by their own example and showered upon me their constant blessings. I would like to share this moment of happiness with my beloved sister Ashita, my friend Preeti and my cousins Pratik and Ankit. They rendered me enormous support during the whole tenure of my research.
ii
ABSTRACT KEYWORDS: Combustion, linear coupling, non-normality, Rijke Tube, thermoacoustic instabilities, Galerkin, pseudospectra, stochastic, singularvalue, eigenvalue.
There have been conflicting opinions in literature, about the consequences of ignoring the linear coupling terms in finite-dimensional thermoacoustic models. This thesis attempts to clarify this issue. The consequences of ignoring the linear coupling terms on the non-normal thermoacoustic interactions in a horizontal Rijke tube are studied. In the Galerkin formulation, pseudospectra and Henrici index of non-normality are used to study the impact of linear coupling between the Galerkin modes on the non-normality and the ensuing transient growth. This linear coupling is important for the energy transfer between the eigenmodes of the system (which are not the same as the natural acoustic modes of the duct). Although ignoring this linear coupling results in negligible shift in the eigenfrequencies, it may have serious consequences on the non-normality of the thermoacoustic system and its dynamics. Nonlinear phenomena such as triggering, boot-strapping and limit cycle oscillations may be completely missed out when this linear coupling is ignored. Ignoring the linear coupling between the Galerkin modes is justified only when the non-normality and the nonlinear dynamics of the system are not altered significantly. Stochastic sources in a thermoacoustic system may have serious consequences on its linear and nonlinear stability. The impact of such random sources on a linearized nonnormal thermoacoustic system is studied here. The stochastic sources are assumed to be multiplicative. Deterministic evolution equations of the ensemble mean state space vector and its covariance are obtained from the stochastic evolution equations. Optimal initial conditions for the maximum transient energy growth in the presence of uncertainties are obtained. The robustness of similar calculations made for the corresponding deterministic thermoacoustic system is verified. iii
TABLE OF CONTENTS DEDICATION
i
ACKNOWLEDGEMENTS
ii
ABSTRACT
iii
LIST OF FIGURES
vii
NOMENCLATURE
viii
1
2
3
INTRODUCTION
1
1.1
Thermoacoustic Instabilities . . . . . . . . . . . . . . . . . . . . .
1
1.2
Linear Coupling in Thermoacoustics . . . . . . . . . . . . . . . . .
2
1.3
Stochastic Processes in Thermoacoustics . . . . . . . . . . . . . . .
5
1.4
Outline of the Thesis . . . . . . . . . . . . . . . . . . . . . . . . .
7
NON-NORMALITY AND LINEAR COUPLING
8
2.1
Non-normality . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8
2.2
Apparent Linear Coupling and Physical Linear Coupling . . . . . .
9
IMPACT OF LINEAR COUPLING ON THERMOACOUSTIC INSTABILITIES IN A RIJKE TUBE
14
3.1
Governing Equations . . . . . . . . . . . . . . . . . . . . . . . . .
14
3.2
Linearized Rijke Tube Model . . . . . . . . . . . . . . . . . . . . .
17
3.2.1
Impact of Linear Coupling on Non-normality . . . . . . . .
20
3.2.2
Results and Discussions . . . . . . . . . . . . . . . . . . .
22
Nonlinear Rijke Tube Model . . . . . . . . . . . . . . . . . . . . .
26
3.3.1
Impact of Linear Coupling on Nonlinearity . . . . . . . . .
26
3.3.2
Results and Discussions . . . . . . . . . . . . . . . . . . .
28
Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
31
3.3
3.4
iv
4
5
IMPACT OF STOCHASTIC PROCESSES ON THERMOACOUSTIC INSTABILITIES IN A RIJKE TUBE
32
4.1
Elementary Definitions and Concepts . . . . . . . . . . . . . . . .
33
4.1.1
Probability Space . . . . . . . . . . . . . . . . . . . . . . .
33
4.1.2
Random Variable or Stochastic Variable . . . . . . . . . . .
33
4.1.3
Probability Distribution Function and Probability Density Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
34
4.1.4
Joint Distribution Function and Joint Density function . . .
34
4.1.5
Mean, Variance and Standard Deviation . . . . . . . . . . .
35
4.1.6
Auto-correlation and Auto-Covariance . . . . . . . . . . . .
35
4.1.7
Cumulants . . . . . . . . . . . . . . . . . . . . . . . . . .
36
4.1.8
Random Process or Stochastic Process . . . . . . . . . . . .
36
4.1.9
Gaussian Distribution . . . . . . . . . . . . . . . . . . . . .
37
4.2
Linearized Governing Equations . . . . . . . . . . . . . . . . . . .
38
4.3
Results and Discussions . . . . . . . . . . . . . . . . . . . . . . . .
42
4.4
Remarks and Further Scope . . . . . . . . . . . . . . . . . . . . . .
44
CONCLUSIONS
46
A PSEUDOSPECTRA
48
B EFFECT OF NON-IDEAL BOUNDARY CONDITIONS ON GALERKIN MODES 50 C MAXIMUM TRANSIENT ENERGY GROWTH AND OPTIMAL INITIAL CONDITIONS FOR DETERMINISTIC SYSTEMS
53
D DETERMINISTIC EVOLUTION EQUATION OF THE ENSEMBLE MEAN FROM STOCHASTIC EVOLUTION EQUATION
55
E MAXIMUM TRANSIENT ENERGY GROWTH AND OPTIMAL INITIAL CONDITIONS FOR STOCHASTIC SYSTEMS
59
REFERENCES
64
LIST OF PAPERS BASED ON THIS THESIS
65
LIST OF FIGURES 1.1
Schematic of Thermoacoustic Coupling . . . . . . . . . . . . . . .
2
1.2
The phase portrait of the limit cycle pressure oscillations slightly vary from cycle to cycle in a seemingly random fashion. From Lieuwen (2003) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6
The pressure oscillations and power spectrum of a combustor exhibiting combustion instability. From Burnley (1996) . . . . . . . . . . . .
7
(a) shows monotone decay of a normal system, (b) shows transient growth of a non-normal system. The initial state is φ(0) = d1 (0)e1 + d2 (0)e2 , and the final state φ(t) = d1 (t)e1 + d2 (t)e2 . The dashed lines denote the vectors at time t = 0 and the solid lines denote the vector at some time t. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9
1.3
2.1
2.2
(a) shows pseudospectra in the presence of linear coupling terms (b) shows pseudospectra in the absence of linear coupling terms (The grayscale bar indicates the values of log10 ε). . . . . . . . . . . . . . . . . . . 12
3.1
shows a schematic of the horizontal Rijke tube setup . . . . . . . .
14
3.2
shows the pseudospectra with the linear coupling of the Galerkin modes present (The grayscale bar indicates the values of log10 ε. . . . . . .
23
shows the pseudospectra with the linear coupling of the Galerkin modes absent (The grayscale bar indicates the values of log10 ε). . . . . . .
24
shows the linear plot of the norm of the matrix exponential of the stability operators multiplied by non dimensional time t, ||eLt || as a function of non dimensional time. . . . . . . . . . . . . . . . . . . . . . . .
25
shows the evolution of acoustic velocity at xf with non-dimensional time, for the case with the apparent linear coupling present. Initial conditions are ηi (0) = 0 for i = 2, 3, ..., N , η˙ i (0) = 0 for i = 1, 2, 3, ..., N and (a) η1 (0) = 0.1, (b) η1 (0) = 0.25 . . . . . . . . . . . . . . . .
28
shows the evolution of acoustic velocity at xf with non-dimensional time, for the case with the apparent linear coupling absent. Initial conditions are ηi (0) = 0 for i = 2, 3, ..., N , η˙ i (0) = 0 for i = 1, 2, 3, ..., N and (a) η1 (0) = 0.1, (b) η1 (0) = 0.25, (c) η1 (0) = 1.0 . . . . . . . .
29
3.3 3.4
3.5
3.6
vi
3.7
shows the evolution of acoustic velocity at xf with non-dimensional time. Initial conditions are ηi (0) = 0 for i = 2, 3, ..., N , η˙ i (0) = 0 for i = 1, 2, 3, ..., N and (a) for the case with apparent linear coupling present, η1 (0) = 0.02, (b) for the case with apparent linear coupling absent, η1 (0) = 0.02 . . . . . . . . . . . . . . . . . . . . . . . . .
30
4.1
Thermoacoustic System with input - output viewpoint . . . . . . . .
32
4.2
Gaussian Probability Density Function for zero mean and different variances . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
37
Gaussian Probability Distribution Function for zero mean and different variances . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
38
Optimal initial conditions , u0 (x, 0) and p0 (x, 0), for maximum transient energy growth: for randomness in the interaction index, n . . . . . .
42
Optimal initial conditions , u0 (x, 0) and p0 (x, 0), for maximum transient energy growth for high autocorrelation time: for randomness in the interaction index, n . . . . . . . . . . . . . . . . . . . . . . . . . . .
43
Optimal initial conditions , u0 (x, 0) and p0 (x, 0), for maximum transient energy growth: for randomness in the time lag , τ . . . . . . . . . .
44
4.3 4.4 4.5
4.6
vii
NOMENCLATURE Roman b
Defined in Equation 3.9
co
Mean speed of sound
Cv
Specific heat capacity per unit volume of air
dw
Diameter of the wire
F
Fluctuating component of the stability matrix
Gmax
Maximum growth factor in a given time interval
He
Henrici number
kj
Wave number of the jth Galerkin mode
L
Mean stability matrix
La
Duct Length
Lw
Equivalent length of the wire
M
Mach number
p
Pressure
Q˙
Heat release rate per unit volume
S
Cross-section area of the Rijke tube
To
Mean Temperature
Tw
Temperature of the wire
t
Time
tc
Noise autocorrelation time
u
Velocity
xf
Flame location
x¯
Dimensional mean variable x
x˜0
Dimensional perturbation variable x
x0
Non-dimensional perturbation variable x
A†
Adjoint or Hermitian transpose of matrix A
viii
Greek α
Amplitude of noise
γ
Specific heat ratio
ε
Measure of perturbation over non-normal matrix L. See Appendix A
ε0
Order of acoustic velocity perturbation
ηj
Velocity amplitude of jth Galerkin mode
λ
Heat conductivity of air
µ
Order of mean flow
ν
Defined in Equation 3.14
νc
1/tc
ρ
Density
κ
Kreiss Constant
τ
Non-dimensional time lag
ξj
Damping constant of the jth Galerkin mode
δ(.)
Dirac delta function
ζ(t)
Gaussian stochastic process
h.i
Notation for the ensemble mean quantity
ix
CHAPTER 1
INTRODUCTION
1.1
Thermoacoustic Instabilities
The coupling between acoustics and heat release dynamics leads to unsteady combustion in the combustors of rockets, aircraft engines and power generating gas turbines, making these thermoacoustic systems unstable. These instabilities are popularly known as combustion instabilities or thermoacoustic instabilities. These instabilities can give rise to high amplitude pressure oscillations which can lead to flame blowout, structural damage etc. (Zinn and Lieuwen, 2005). Thermoacoustics play a vital role in many technical applications (Matveev, 2003). Rocket motors (solid or liquid) and gas turbines can exhibit combustion instabilities. Oscillations inside the combustion chamber or a motor can lead to unacceptable levels of vibration and enhanced heat transfer that can degrade propulsive efficiency. This can damage or even destroy the system. Recent regulations for environmental protection require significant reduction in emission of pollutants. This requires leaner combustion in the chambers, which further increases the chances of thermoacoustic instabilities. Intense noise generation is another serious problem associated with unstable combustion chambers. Thermoacoustic instabilities is not always a negative phenomenon. Active research is going on in the direction of pulse combustors and thermoacoustic engines that can be used as prime movers, refrigerators, etc. An in-depth understanding of combustion instabilities is thus necessary for their accurate prediction and control. Every combustor possesses certain acoustic properties. If the heat released in the combustor depends on the pressure and velocity fluctuations, a feedback loop is established which might lead to instabilities in the system. This can be understood from the schematic shown in Figure 1.1. The phase between the pressure and heat release oscillations at the flame determines whether combustion instabilities can occur in the system
Figure 1.1: Schematic of Thermoacoustic Coupling
or not (Rayleigh, 1878). Acoustic pressure oscillations are amplified (until limited by some nonlinear characteristic) if the pressure oscillations and the heat release fluctuations are in phase; however they are attenuated when the pressure oscillations and the heat release fluctuations are out of phase. The closed-loop system, as shown in Figure 1.1, represents the interaction between the acoustics and the thermodynamics of a heat source. The dynamics of this mechanism is still not thoroughly understood by the thermoacoustic community. This study will further the understanding in the mechanism of combustion instability or thermoacoustic instabilities.
1.2
Linear Coupling in Thermoacoustics
Finite-dimensional or Galerkin models (Meirovitch, 1967) have been widely used to study thermoacoustic interactions (Culick, 1976; Dowling, 1995). Most of these models adopt the natural modes of the combustors as the basis expansion functions or the Galerkin modes. Extensive research on longitudinal thermoacoustic instabilities in solid rocket motors, with a distributed heat release, has been carried out assuming that these acoustic Galerkin modes are linearly uncoupled (Culick, 1976; Culick et al., 1995; Ananthkrishnan et al., 2005; Wicker et al., 1996). The primary purpose of this assumption is to obtain closed form solutions for frequency increments and damping factors consisting of sums of terms, each associated with an individual physical effect. This methodology allows new physical effects to be added to the analysis without the need to do the entire analysis over again.
2
Culick (1997) shows that ignoring the linear coupling between the Galerkin modes in a combustor with a localized heat release produces negligible change in the eigenfrequencies of the system. Hence, Culick (1997) claimed that the linear coupling terms can be ignored for analytical and numerical convenience. On the contrary, Annaswamy et al. (1997) show the importance of the linear coupling between the Galerkin modes in model-based active control design. They develop a thermoacoustic model for a bench top premixed combustor. They show that the response of the thermoacoustic dynamics to external actuation depend strongly on this linear coupling, directly and through the heat release. In this thesis, these conflicting views for a thermoacoustic system with a localized heat release source, in the context of a horizontal Rijke tube is clarified. Hwang and Ma (1993) and Park et al. (1994) studied the impact of ignoring the off-diagonal elements of a modal damping matrix in a non-classically damped linear system, referred by them as the decoupling approximation. They note that this decoupling approximation method is a common procedure adopted to substantially reduce the computational effort in large scale systems. However, this approximation might lead to significant errors which may make the model unsuitable for engineering applications. In thermoacoustics, neglecting the linear coupling between the Galerkin modes is equivalent to changing a non-classically damped system to a classically damped system. The non-normal nature of thermoacoustic interactions has received attention only recently. Balasubramanian and Sujith (2007a,b, 2008a,b) have highlighted the role of non-normality in combustion-acoustic interaction. Non-normality can lead to transient growth of a system even when the eigenvalues indicate linear stability. Transient growth can help the perturbations grow in amplitude to significant levels where the nonlinearities may get triggered. The ensuing nonlinear driving can make the system, which was thought to be classically linearly stable, saturate to a limit-cycle. Nicoud et al. (2007) has shown that the eigenvectors of a thermoacoustic system are non-orthogonal in the presence of heat release or in the presence of a general complex impedance boundary condition. Non-normality has been extensively studied in the context of atmospheric flows (Trefethen et al., 1993; Farrell and Ioannou, 1996b,a) and hydrodynamic stability (Schmid and Henningson, 2001; Reddy and Henningson, 1993; Baggett et al., 1995). Schmid 3
(2007) is an excellent review paper that can be referred to understand non-normality in the context of hydrodynamic stability. The impact of the linear coupling between the Galerkin modes on the non-normality of a thermoacoustic system, to the best of our knowledge, has never been taken into account in the literature. It is shown that the linear coupling between the Galerkin modes may have significant effects on the calculation of transient growth of acoustic oscillations. In Chapter 2, it is shown that non-normality arises due to the linear coupling between the eigenmodes and that the linear coupling between the Galerkin modes is only a projection of the coupling between the eigenmodes onto the space spanned by the basis functions. With the help of pseudospectra (see Appendix A), the change in non-normality of a dynamical system in the presence and absence of the linear coupling between Galerkin modes is shown . Owing to its simplicity, Rijke tube has been widely studied in literature to investigate the thermoacoustic instabilities (Yoon et al., 1998, 2001; Heckl, 1990; Matveev, 2003; Matveev and Culick, 2003; Kopitz and Prolifke, 2005). Heckl and Howe (2007) present a stability analysis of a Rijke tube using Green’s function as an alternative to the eigenvalue approach to the stability analysis. Balasubramanian and Sujith (2008b) show the consequences of non-normality in a horizontal Rijke tube. Non-normality of the linearized system plays a crucial role in the subcritical regime (regime of classical linear stability), where the system might show triggering for finite amplitude perturbations. In the supercritical regime (regime of linear instability), non-normality of the linearized system becomes important in the presence of control. Active control of combustion instabilities in the supercritical regime, including the effects of non-normality, has been addressed by Kulkarni et al. (2008) in the context of a horizontal Rijke tube. The scope of the present thesis is limited to the subcritical case alone. In this thesis, using a simple model for a Rijke tube, it is demonstrated that in the presence of the linear coupling between the Galerkin modes, there is a negligible shift in the eigenfrequencies, in agreement with the analysis of Culick (1997). However, because of the non-normal nature of the thermoacoustic interactions, acoustic oscillations may experience large transient growth, in the presence of this linear coupling. When these oscillations reach sufficiently high amplitude, nonlinear driving may occur lead4
ing to phenomena such as triggering and bootstrapping. Neglecting the linear coupling between the Galerkin modes may alter the system dynamics to such an extent that the fore-mentioned nonlinear phenomena might be missed.
1.3
Stochastic Processes in Thermoacoustics
Combustion takes place in a highly noisy, turbulent environment. There is often flow separation present in the combustion chamber. The resulting vorticity and entropy waves form sources of uncertainties in the thermoacoustic system (Burnley, 1996). Alternatively, they are the sources of stochastic or random processes. Combustion instability hence involves an interplay of a variety of complex physical phenomena, which may involve linear, non-normal, nonlinear and stochastic processes. Completely deterministic models do not capture several important characteristics of combustion instabilities. The practical combustor exhibits many stochastic features that cannot be characterized by deterministic models. Figure 1.2, taken from Lieuwen (2003), shows that the pressure oscillations in the limit cycle follow an orbit which differs from cycle to cycle in a seemingly random fashion. Lieuwen (2003) indicates that the variability shown in Fig 1.2 is uncorrelated from one cycle to the next. This suggests that there are high degree of freedom processes with short correlation times compared to the acoustic period. Figure 1.3, taken from Burnley (1996), shows the power spectrum of an experimental data of a combustor exhibiting instability. The power spectrum shows well defined peaks in addition to the background noise over the entire range of frequencies. Burnley (1996) and Clavin et al. (1994) suggest that the instability characteristics are random variables due to the presence of stochastic sources. Lieuwen (2003) carries out a statistical analysis of the features of combustion instability. There can be very serious practical consequences of stochastic processes. They can modify both the linear and nonlinear characteristics of the combustor. Stochastic processes can affect the linear stability boundaries of the combustor. This may even lead to instabilities or "noiseinduced transitions". Nonlinear characteristics can be altered by altering instability amplitudes. A motor which may be stable during its initial testing can suddenly become 5
Figure 1.2: The phase portrait of the limit cycle pressure oscillations slightly vary from cycle to cycle in a seemingly random fashion. From Lieuwen (2003)
unstable for no conspicuous reason. Stochastic processes can also limit the extent to which active control systems are capable to control combustion instabilities. A rigorous analysis of random processes and their effects on thermoacoustic instabilities is not yet carried out by engineers or scientists. This analysis is extremely complex and is very poorly understood by the combustion instability community. Lieuwen and Banaszuk (2005) studies the effect of noisy interaction index of heat release, frequency and time lag in a time delayed thermoacoustic system. They observe qualitative changes in the system dynamics and quantitative changes in the stability boundaries. They conclude that a detailed analysis of background turbulent fluctuations in the combustor is warranted. In the present thesis, only a linearized thermoacoustic system is focussed upon for the analysis of stochastic processes. Farrell and Ioannou (2002a,b) deal with characterizing uncertain perturbations in the context of atmospheric flows. Analytical dynamical system equations are obtained in the presence of stochastic sources and non-normal features are characterized. Maximum energy growth and the corresponding optimal initial conditions are obtained in the presence of uncertainties. Schmid (2007) presents similar analysis for a linearized Navier Stokes equation in the context of hydrodynamics stability. The present thesis aims to study the impact of uncertain perturbations in a linearized
6
Figure 1.3: The pressure oscillations and power spectrum of a combustor exhibiting combustion instability. From Burnley (1996)
thermoacoustic system, having a time delayed heat response. It lays the foundation for similar analysis applied to even more complicated thermoacoustic systems.
1.4
Outline of the Thesis
In this chapter, a general introduction and motivation of linear coupling and stochastic process in thermoacoustics is given. Relevant literature survey is also done. In chapter 2, a brief introduction to non-normality is given. The concept of linear coupling is clarified and the terms, apparent linear coupling and physical linear coupling, are introduced. Chapter 3 deals with one of the two central issues which this thesis focuses. It deals with the impact of linear coupling of the Galerkin modes on the non-normality and nonlinearity in a horizontal Rijke tube. Chapter 4 deals with a linearized Rijke tube model and studies the effects of uncertainties or stochastic processes in the thermoacoustic system. Analytical deterministic system equations are derived from stochastic differential equations. Finally, the work is summarized and conclusions, based on the discussions in the previous chapters, are drawn. Relevant derivations are added as appendices.
7
CHAPTER 2
NON-NORMALITY AND LINEAR COUPLING
2.1
Non-normality
A matrix operator A is said to be non-normal if it does not commute with its adjoint 1
; i.e. AA† 6= A† A where A† is the adjoint of the matrix (Schmid and Henningson,
2001). The eigenvectors of a non-normal operator are nonorthogonal. Balasubramanian and Sujith (2008b,a) have shown that a dynamical system such as the thermoacoustic system in a combustor is non-normal as it’s linearized evolution operator is non-normal.
Classical linear stability analysis becomes a poor predictor of the short term dynamics for non-normal systems because the nonorthogonality of the eigenvectors can lead to transient growth. This property is illustrated in Figures 2.1a and 2.1b, where e1 and e2 represent the direction of the eigenvectors and φ is the vector sum of the eigenvectors. Figure 2.1 shows that for a normal system, the eigenvectors are orthogonal. If the amplitude of the individual eigenvectors decays, then φ decreases monotonically. On the contrary, for the non-normal system shown in Figure 2.1a, the eigenvectors are nonorthogonal. The case shown is a special one wherein φ increases even when the amplitudes of individual eigenvectors decay. However, if in the process of transient growth the nonlinear effects continue to remain insignificant, then for a linearly stable system, φ decays after a sufficiently long time. On the contrary, during transient growth, if the amplitudes become sufficiently high to make the nonlinear processes significant, nonlinear driving may occur and the system which was classically linearly stable might exhibit instabilities. 1
The adjoint of a real matrix is the transpose of the matrix. For a differential operator T , the adjoint T is such that hu, T vi = hT † u, vi where h·, ·i is the inner product. u and v are two vectors of same length. †
Figure 2.1: (a) shows monotone decay of a normal system, (b) shows transient growth of a non-normal system. The initial state is φ(0) = d1 (0)e1 + d2 (0)e2 , and the final state φ(t) = d1 (t)e1 + d2 (t)e2 . The dashed lines denote the vectors at time t = 0 and the solid lines denote the vector at some time t.
2.2
Apparent Linear Coupling and Physical Linear Coupling
A linear dynamical system can be expressed as dχ = Aχ dt
(2.1)
where χ = [η1 , η2 ....ηn ]T is a vector space of n modes and A is the n × n linear matrix operator. When A has non-zero off-diagonal terms, it is preferred to refer to the coupling between the modes η1 , η2 ....ηn as apparent linear coupling. The space spanned by the basis χ can be transformed to the space spanned by the eigenvectors using the following transformation: E = V Tχ
(2.2)
where V is the eigenvector matrix of A and E is the space spanned by the eigenvectors of A. The non-normality of a system can be visualized using the angle between
9
the eigenvectors. The dot product of the eigenvectors, giving the cosine of the angles between them, can be represented as
e .e . . . e1 .en 1 1 .. . .. V † V = .. . . en .e1 · · · en .en
(2.3)
where V † is the adjoint of V and ej is the j th eigenvector. For a system with orthogonal eigenvectors, V † V will be a diagonal matrix (property of a normal matrix A with eigenvector matrix V ). In general V † V is a symmetric matrix, since the dot product is commutative. The off-diagonal terms in V † V represent the coupling between the nonorthogonal eigenvectors referred to as the physical linear coupling by the authors. The physical coupling provides a mechanism for energy transfer among the eigenvectors and results in the non-normality of the system. A non-normal system can experience transient growth as discussed by Balasubramanian and Sujith (2008b,a). However it must be emphasized that the presence of large apparent linear coupling terms does not guarantee higher transient growth. Transient growth depends on how non-normal the system is, and how fast the eigenvectors are decaying relative to each other. In classical linear stability, i.e., with all eigenvalues having negative real parts, if the real part of one eigenvalue is more negative as compared to the real part of another eigenvalue, then the corresponding eigenvectors decay at different rates, which may result in higher transient growth. The off-diagonal elements of V † V give us a measure of the non-normality of the system. A number of scalar measures have been used to quantify non-normality in literature. Trefethen (1999) provides an excellent introduction to interpret and compute some of these scalar measures of non-normality. The Henrici number, as a scalar measure of non-normality of a matrix (Henrici, 1969), is widely used in literature. It is defined as He =
° ° † °A A − AA† ° kA2 k
F
(2.4)
F
where A and A† are the matrix under study and its adjoint respectively. kAkF denotes
10
the Frobenius norm of A, which is defined as v u m n 2 uX X t kAkF = |aij |
(2.5)
i=1 j=1
where aij denotes an element of A. The Henrici number quantifies the extent to which the eigenvectors are coupled to each other. For a system with orthogonal eigenvectors, the Henrici number will be zero indicating that the eigenvectors are not coupled. Like all other scalar measures of non-normality, the Henrici number suffers from a basic limitation. Non-normality is far too complex to be completely described by a single scalar measure (Trefethen and Embree, 2005). In finite-dimensional Galerkin models, increasing the number of modes changes the dimensions of the stability operator or the matrix, thereby changing the Henrici number. It must be emphasized that, for the sake of comparison of the degree of non-normality of two different linear stability operators having the same dimension, the Henrici number may serve a purpose. However, it cannot be used alone to quantify and compare non-normality. To understand the role played by the apparent and physical linear coupling terms in determining the non-normality of a linear matrix operator, we consider two simple examples. For a 2D linear system expressed by 2.1 with χ = [η1 , η2 ]T , as a first example, we choose
A=
1
1000
1000
1
(2.6)
The matrix A is symmetric and hence is normal. For this matrix, V † V =
1 0
.V † V
0 1 is a diagonal matrix; i.e., the physical linear coupling is absent, verifying that the eigen-
vectors are orthogonal to each other. Hence matrix A is normal. However, there is strong apparent linear coupling between η1 and η2 as shown by the off-diagonal terms in the matrix A. This example demonstrates the fact that apparent linear coupling between the modes is not a measure of non-normality of the linear matrix operator.
11
As a second example, we choose A=
−0.1 −0.001 0.9
This matrix has V † V =
(2.7)
−0.1
1
0.9978
. It can be seen that the order of mag-
0.9978 1 nitude of one of the apparent linear coupling terms in A is much less compared to the order of magnitude of the diagonal terms. This indicates that the physical linear coupling and hence the non-normality is much higher than anticipated. If we ignore the
apparent linear coupling terms of A, the system will become normal (V † V will then be a diagonal matrix) and the physics of the problem would show drastic changes.
Figure 2.2: (a) shows pseudospectra in the presence of linear coupling terms (b) shows pseudospectra in the absence of linear coupling terms (The grayscale bar indicates the values of log10 ε).
ε-pseudospectrum is a useful tool in understanding and interpreting non-normality. The change in the non-normality of A =
−0.1 −0.001
in the presence and ab-
0.9 −0.1 sence of the coupling terms is clearly seen in the pseudospectra. The ε-pseudospectrum is obtained using Eigtool (Wright, 2002). In the presence of the apparent linear cou-
pling terms, the pseudospectra contours spill to the right of the imaginary axis as shown in Figure 2.2a. This indicates that the system may show transient growth (see Section 3.2.2 and Appendix A for detailed explanation). In the absence of the apparent linear 12
coupling terms, the system becomes normal. The pseudospectra contours are closed circles, and lie entirely to the left of the imaginary axis as shown in Figure 2.2b. In summary, if the apparent linear coupling terms, i.e., the off-diagonal terms of the matrix A are of the same order of magnitude, then the system is normal. In this case, the physical linear coupling terms, i.e., the off-diagonal terms of V † V are zero. The system will be non-normal even if the apparent linear coupling is ignored. On the contrary, if the apparent linear coupling terms are not of the same order, the physical coupling terms are non-zero. In this case ignoring the apparent linear coupling terms would wrongly render the system normal. Thus, the decision to ignore the apparent linear coupling terms should be made considering the presence or absence of the physical linear coupling terms. In Chapter 3, the impact of ignoring the apparent linear coupling terms in a thermoacoustic system in the context of a Rijke tube is discussed.
13
CHAPTER 3
IMPACT OF LINEAR COUPLING ON THERMOACOUSTIC INSTABILITIES IN A RIJKE TUBE
Balasubramanian and Sujith (2008b) discuss the consequences of non-normality and nonlinearity in the context of a horizontal Rijke tube. The current study adopts their horizontal Rijke tube model, a schematic of which is shown in Figure 3.1. A similar schematic setup was adopted by Matveev (2003) and Matveev and Culick (2003).
Figure 3.1: shows a schematic of the horizontal Rijke tube setup
3.1
Governing Equations
In this section, coupled ordinary differential equations that describe the thermoacoustic interactions in a horizontal Rijke tube are derived. A perfect, inviscid, and non-heatconducting gas is assumed. The one-dimensional, linearized, acoustic momentum and acoustic energy equation, in the presence of a heat source, for a negligible mean flow and negligible mean temperature gradient, can be respectively written as (Balasubramanian and Sujith, 2008a) ρ¯
∂ u˜0 ∂ p˜0 + =0 ∂ x˜ ∂ t˜
(3.1)
∂ p˜0 ∂ u˜0 ˜˙ 0 + γ p¯ = (γ − 1) Q ∂ x˜ ∂ t˜
(3.2)
and, c2o =
γ p¯ ρ¯
(3.3)
Equations 3.1 and 3.2 can be non-dimensionalized by using the following non-dimensional variables x=
x˜ c0 u˜0 p˜0 u¯ ˙ 0 La ; t = t˜ ; u0 = ; p0 = ; M = ; Q˙ 0 = ˜Q 3 La La u¯ p¯ c0 ρ¯c0
The mean flow is assumed to be steady. Assuming M << 1, the zero mean flow acoustic equations are used (Nicoud et al., 2007). The corresponding non-dimensional linear acoustic momentum and energy equations can hence be written as (Balasubramanian and Sujith, 2008b) γM
∂u0 ∂p0 + =0 ∂t ∂x
(3.4)
∂p0 ∂u0 + γM = γ(γ − 1)Q˙ 0 ∂t ∂x
(3.5)
The open-open boundary conditions for a duct are applied; i.e.,p0 = 0 at both x = 0 and x = 1. The effect of non-ideal boundary conditions on the Galerkin modes can be seen in Appendix B. The above set of partial differential equations can be reduced to a set of ordinary differential equations using the Galerkin technique (Meirovitch, 1967). The acoustic velocity and the acoustic pressure can be written in terms of the first N natural modes of the duct (Dowling, 1995; Dowling and Stow, 2003) as u0 (x, t) = ε0
N X
j=1
ηj (t) cos(kj x)
and p0 (x, t) = −ε0
N X γM
j=1
η˙ j (t) sin(kj x) kj
(3.6)
where kj = jπ is the wavenumber and ε0 is the order of the velocity perturbation over the mean flow 1 . This is used for book keeping, in order to be consistent with the orders of magnitudes. The basis functions are chosen in such a way that they satisfy the boundary conditions. It must be emphasized that the basis functions chosen are the 1 Conventionally, ε is used to denote the order of velocity perturbation (Culick, 1997). However in this thesis, ε is reserved for pseudospectra (see Appendix A), and hence to avoid ambiguity, ε0 is used as the order of velocity perturbation.
15
eigenmodes of the self-adjoint part of the linearized system, and not the eigenmodes of the linearized system; i.e., the basis functions are the natural modes of the duct in the absence of the heat release. A comprehensive study of linear and ad-hoc nonlinear Rijke tube velocity coupled oscillatory heat release models is done by Yoon et al. (1998, 2001). To qualitatively study the effect of linear coupling, a modified form of King’s law is used, proposed by Heckl (1990). ˜˙ 0 = µ Q
(
2Lw (Tw − T0 ) √ S 3
r
dw πλCv ρ¯u ¯ 2
"s¯ ) r # ¯ ¯1 ¯ 1 ¯ + u0 (t − τ )¯ − δ(˜ x−x ˜f ) ¯3 ¯ 3
(3.7)
The interaction between the combustion process and the mean flow is of the order of the mean flow µ (Culick, 1997). Thus, the order of the oscillatory heat release is µ times over the order of the perturbation ε0 . The acoustic energy, given by Eq. 3.5, can thus be written as ( "s¯ ) ¯ r # ¯1 ¯ ∂p ∂u 1 ¯ + u0 (t − τ )¯ − + γM =µ b δ (x − xf ) ¯3 ¯ ∂t ∂x 3 0
0
(3.8)
where the constant b is defined as 2La Lw (Tw − T0 ) √ b = γ (γ − 1) ρ¯c30 S 3
r πλCv ρ¯u¯
dw 2
(3.9)
Substituting the Galerkin expansions given by Eq. 3.6 in Eqs. 3.3 and 3.4 and projecting them along the basis functions 2 , we get the following set of evolution equations: 2
d ηj µ + kj2 ηj = − 0 2 dt ε
(
"s¯ ) ¯ r # ¯ ¯ 2b ¯ 1 + u0 (t − τ )¯ − 1 sin(kj xf ) kj f ¯3 ¯ γM 3
(3.10)
for j = 1, 2, ..., N where j indicates the j th Galerkin mode. In the presence of damping, the set of the equations given by Eq. 3.10 can be modified as (Matveev, 2003) µ d2 ηj + 2ξj kj η˙ j + kj2 ηj = − 0 2 dt ε
(
2bkj γM
"s¯ ) ¯ r # ¯ ¯1 1 ¯ + u0 (t − τ )¯ − sin(kj xf ) (3.11) f ¯ ¯3 3
2 The projection of Ra function f along a basis function ψ is given by their inner product < f |ψn > which is defined as f (x)ψn (x)dx. domain
16
where the damping constant is given as s # " kj 1 k1 ξj = c1 + c2 2π k1 kj
(3.12)
The higher frequency waves are damped at a higher rate when compared to the lower frequency waves. Acoustic boundary layer losses, sound radiation losses at the ends and the convection of sound by the mean flow are a few damping mechanisms. It must be noted that Eq. 3.11 is derived by assuming that the nonlinearities present are due to the heat transfer rate alone, and the nonlinear gas dynamics terms are ignored.
3.2
Linearized Rijke Tube Model
For small values of u0 , Eq. 3.11 can be linearized as ª d2 η˙ j µ© 2 0 + 2ξ k η ˙ + k η = − νk u (t − τ ) sin(k x ) j j j j j j f j f dt2 ε0 where
(3.13)
√ ν=
3b γM
(3.14)
Here τ is the non-dimensional time lag in the heat release rate fluctuation due to the velocity perturbations. In this subsection, we choose τ = 0 for mathematical simplicity. Following Culick (1997), we assume a two mode expansion for analytical convenience. We can write Eq. 3.13 as µ d2 ηj + 2ξj kj η˙ j + kj2 ηj = − 0 ε0 2 dt ε
( νkj sin(kj xf )
2 X
) ηm cos(km xf )
(3.15)
m=1
where j = 1, 2. Ignoring the damping for the time being, the evolution of the individual Galerkin modes can be expressed as d2 η1 + k12 η1 = −µ{νk1 sin(k1 xf ) cos(k1 xf )η1 + νk1 sin(k1 πxf ) cos(k2 xf )η2 } (3.16) dt2
17
d2 η2 + k22 η2 = −µ{νk2 sin(k2 xf ) cos(k1 xf )η1 + νk2 sin(k2 πxf ) cos(k2 xf )η2 } (3.17) dt2 Equations 3.16 and 3.17 can be written as
d2 dt2
+
(k12
+ µA11 ) d2 dt2
µA21
η1 µA12 =0 + (k22 + µA22 ) η2
A11 = νk1 cos(k1 xf ) sin(k1 xf ), A12 = νk1 cos(k2 xf ) sin(k1 xf )
(3.18)
(3.19)
A21 = νk2 cos(k1 xf ) sin(k2 xf ), A22 = νk2 cos(k2 xf ) sin(k2 xf ) Equation 3.18 is deliberately written into a form expressed in Culick (1997). Following the analysis in Culick (1997), the eigenfrequencies and the eigenmodes are found out by assuming that the amplitudes have the same time dependence, η1 = ηˆ1 eλt and η2 = ηˆ2 eλt
(3.20)
The above step attempts to give the eigenvalue formulation of the system under consideration. Here λ and ηˆ represent the eigenvalue and the eigenvectors respectively. The equations for ηˆ1 and ηˆ2 are then, ηˆ1 λ + + µA11 ) µA12 =0 µA21 λ2 + (k22 + µA22 ) ηˆ2
2
(k12
(3.21)
The roots of the characteristic equation formed give the values of λ2 λ21 = −k12 (1 + µ2A1 − µ2 2B1 )
(3.22)
λ22 = −k22 (1 + µ2A2 − µ2 2B2 )
(3.23)
A11 A22 A12 A21 A12 A21 , A2 = 2 , B1 = 2 2 and B2 = 2 2 2 2 2k1 2k2 2k1 (k2 − k1 ) 2k2 (k2 − k12 )
(3.24)
where the constants are A1 =
18
To second order in µ, Eqs. 3.22 and 3.23 give the perturbed eigenfrequencies as ·
µ
kˆ1 = k1 1 + µA1 − µ
2
A2 B1 + 1 2
¶¸
· µ ¶¸ A22 2 ˆ k2 = k2 1 + µA2 − µ B2 + 2
(3.25)
(3.26)
The steps from Eq. 3.18 to Eq. 3.26 are originally from Culick (1997). It is clear from Eqs. 3.25 and 3.26 that the shift in the eigenfrequencies is of the order µ2 . Even if the damping is not ignored, it can still be analytically worked out that the perturbed eigenfrequencies will be of the order of µ2 . The terms providing the apparent linear coupling indeed contribute very little to the shift in the eigenfrequencies; however, ignoring the apparent linear coupling terms may have a serious impact on the non-normality of the system as seen in chapter 2. As shown by Balasubramanian and Sujith (2008a,b), if a system is non-normal, effects such as transient growth of a classically linearly stable system might trigger nonlinearities leading to thermoacoustic instabilities. Equation 3.15 can be written as a set of first order ordinary differential equations (ODE) as dηj = η˙ j dt
(3.27)
2 X dη˙ j 2 + 2ξj kj η˙ j + kj ηj = −µνkj sin(kj xf ) ηm cos(km xf ) dt m=1
(3.28)
Equations 3.27 and 3.28 can be written in a general form as dχ = Lχ dt
(3.29)
where χ = [η1
η˙ 2 T η˙ 1 η2 ] k1 k2
19
(3.30)
0
k1
−k1 − µν cos(k1 xf ) sin(k1 xf ) −2ξ1 k1 L= 0 0 −µν cos(k1 xf ) sin(k2 xf ) 0
0
0
−µν cos(k2 xf ) sin(k1 xf )
0
0
k2
−k2 − µν cos(k2 xf ) sin(k2 xf ) −2ξ2 k2 (3.31)
or
0 k1 0 0 −k1 − µ Ak11 −2ξ1 k1 −µ Ak12 0 1 1 L= 0 0 0 k2 −µ Ak21 0 −k2 − µ Ak22 −2ξ2 k2 2 2
(3.32)
It can be clearly seen that L does not commute with its adjoint and hence the system given by Eq. 3.29 is non-normal. For a thermoacoustic system, the apparent linear coupling terms are the linear coupling terms between the Galerkin modes (which in this case are the natural duct acoustic modes). It is emphasized again that the apparent linear coupling does not always imply the physical linear coupling between the eigenmodes of the system. The nonnormality in the system may not necessarily be due to the apparent linear coupling but it is necessarily due to the physical linear coupling. The following sections will show the importance of retaining the apparent linear coupling terms for the calculation of non-normal effects such as transient growth, even when the eigenfrequencies do not shift significantly as seen by Eqs. 3.25 and 3.26.
3.2.1
Impact of Linear Coupling on Non-normality
Non-normality in a system may cause transient growth of oscillations which may result in triggering of the nonlinearities. Balasubramanian and Sujith (2008a,b) discuss a method to analyze transient growth of oscillations in a non-normal thermoacoustic system in the context of a Rijke tube and diffusion flames. They show that the non-normal system given by Eq. 3.29 has the following solution: χ(t) = exp(Lt)χ(0)
20
(3.33)
where L is the linear stability operator. The transient growth is quantified by the maximum growth factor, which, over a time interval [0, t], is defined as (Schmid and Henningson, 2001) follows: (See Appendix C for derivation) Ã
Gmax
kχ(t)k2 = max max 2 t χ(0) kχ(0)k
!
¢ ¡ = max kexp(Lt)k2 t
(3.34)
The growth rate at any instant of time is maximized over all the initial conditions. This is achieved by the singular value decomposition. Maximizing it in the time interval [0, t] gives Gmax , which is a measure of the maximum amplification of the energy, for a certain initial condition. As discussed by Balasubramanian and Sujith (2008a,b) and Schmid and Henningson (2001), Gmax = 1 indicates a system not exhibiting any transient growth in the given time interval. A system exhibiting transient growth has Gmax > 1. Extending Eqs. 3.29-3.31 for N Galerkin modes, the linear stability operator L can be written as
0
−k1 − ν C˜1 S˜1 0 −ν C˜1 S˜2 L= ......... ......... 0 ˜ −ν C1 S˜N
k1
0
0
..... .....
0
0
−2ξ1 k1
−ν C˜2 S˜1
0
..... .....
−ν C˜N S˜1
0
0
0
k2
..... .....
0
0
0
−k2 − ν C˜2 S˜2
−2ξ2 k2
..... .....
−ν C˜N S˜2
0
.........
.........
.........
..... .....
....
.........
.........
.........
.........
..... .....
....
.........
0
.........
.........
.....
0
kN
0
−ν C˜2 S˜N
0
.....
..... .....
−kN − ν C˜N S˜N
−2ξN kN (3.35)
where C˜j = cos(kj xf ) and S˜j = sin(kj xf ). Note that L is the linear stability operator with the apparent linear coupling terms present. We get the linear stability operator L0 for the case where the apparent linear
21
coupling is absent as
0
−k1 − ν C˜1 S˜1 0 0 L0 = ......... ......... 0 0
k1
0
0
..... .....
0
0
−2ξ1 k1
0
0
..... .....
0
0
0
0
k2
..... .....
0
0
0
−k2 − ν C˜2 S˜2
−2ξ2 k2
..... .....
0
0
.........
.........
.........
.....
.....
....
.........
.........
.........
.........
.....
.....
....
.........
0
.........
.........
..... .....
0
kN
0
0
0
..... .....
−kN − ν C˜N S˜N
−2ξN kN (3.36)
where C˜j = cos(kj xf ) and S˜j = sin(kj xf ). Note: µ was used in Eqs. 3.25 and 3.26 to infer that the shift in the eigenfrequencies is negligible (Culick, 1997). µ is used only for the book keeping of the order of magnitudes. It does not have any specific value and hence it is dropped while writing the equation for the stability operators L and L0 .
3.2.2
Results and Discussions
The following example of a linearized Rijke tube model demonstrates the effect of the apparent linear coupling terms on the non-normality of the system and the maximum transient growth. The following numerical values are used for computing the stability operators L and L0 from Eqs. 3.35 and 3.36. The heating element is located at the quarter location of the duct or xf = 0.25 from the inlet. The number of Galerkin modes chosen is N = 10 for all the results presented in this section, so that the change in the solution with an increase in the number of modes is less than 5%. The other parameters chosen are Lw = 2.5m, Tw = 900K, u¯ = 0.05m/s, c1 = 0.135, c2 = 0.015, S = 1.0 × 10−3 m2 , γ = 1.4, c0 = 399.6m/s, ρ¯ = 1.025kg/m3 , Cv = 719J/kgK, λ = 0.0328W/mK, dw = 0.0005m and La = 1.0m. Figure 3.2 shows the ε-pseudospectrum of the stability operator L, in which the apparent linear coupling is present. The range of is chosen such that ε << kLk (Note kLk ≈ 850 for the current case). The non-circular nature of the ε-pseudospectrum contours confirm that the system is non-normal. All the eigenvalues of the system lie in the left half plane (or the real part of all the eigenvalues is negative), indicating a 22
classically linearly stable system. On a contour corresponding to an ε, the real part of the point on the contour which is farthest to the imaginary axis, in the right half plane, is the pseudospectral abscissa (see Appendix A for exact mathematical definition). The ratio of the pseudospectral abscissa and its corresponding ε forms a lower bound for the growth factor (see Appendix A). If this ratio is greater than 1, then there is transient growth in the system for certain initial conditions. In Figure 3.2, the pseudospectral abscissa is computed to be 2.08 for ε = 100.2 . The ratio 2.08/100.2 = 1.31, being the lower bound for the growth factor, indicates that there is a transient growth.
Figure 3.2: shows the pseudospectra with the linear coupling of the Galerkin modes present (The grayscale bar indicates the values of log10 ε.
Similar features are also observed in Figure 3.3 which is the ε-pseudospectrum of the stability operator L0 , in which the apparent linear coupling is absent. It can be seen that the ε-pseudospectrum contours in Figure 3.3 tend to protrude in the right half plane to a lower extent as compared to Figure 3.2 for the same values of parameter ε. This clearly indicates that the system is less non-normal in the case without the apparent linear coupling as compared to the case with the apparent linear coupling. In Figure 3.3, the pseudospectral abscissa is computed to be 1.19 for ε = 100.2 . This is much less than the corresponding pseudospectral radius for the same ε, computed in Figure 3.2, further strengthening our argument of significant change in non-normality. The lower 23
Figure 3.3: shows the pseudospectra with the linear coupling of the Galerkin modes absent (The grayscale bar indicates the values of log10 ε).
bound for the growth factor here is 1.19/100.2 = 0.75. It must be noted that although the ratio is less than 1, transient growth is possible in this system, as this ratio is only a lower bound of the growth factor. The spilling of pseudospectra in the right half plane is a necessary but not sufficient condition for transient growth to occur. The ratio of the pseudospectral abscissa to ε, when maximized over all ε, is called the Kreiss constant, κ (see Appendix A for an exact mathematical definition). Kreiss constant provides a tight lower bound for the energy growth factor (Trefethen and Embree, 2005). The non-normality index, given by the Henrici number (Eq. 2.4), also shows a small difference for the two cases. However, the maximum transient energy growth differs considerably, as it will be evident from the following paragraph. The Henrici number for L (Figure 3.2) is 0.119 and the Henrici number for L0 (Figure 3.3) is 0.116. The eigenvalues in these two cases do not change significantly (approximately 2 % change). This can be seen by observing the change in the position of the corresponding eigenvalues by comparing Figure 3.3 and Figure 3.3. Thus, the impact of the apparent linear coupling on the shift in the eigenfrequencies is indeed negligible, as noted by Culick (1997), reproduced in section 3.2, Eqs. 3.25 and 3.26. 24
Figure 3.4: shows the linear plot of the norm of the matrix exponential of the stability operators multiplied by non dimensional time t, ||eLt || as a function of non dimensional time. ° ° ° 0° Figure 3.4 shows the plot of the norm of the matrix exponentials °eLt ° and °eL t ° versus t. The evolution of the norm of the matrix exponential, in the presence of the apparent linear coupling (L), is much higher as compared with the case where the apparent linear coupling is absent (L0 ). The maximum growth factor in the energy of the system can be calculated if the evolution of the norm of the matrix exponential is known, Eq. 3.34. For the case with the apparent linear coupling present (L), the value of Gmax is 4.66 and the non-dimensional time taken to achieve this maximum growth is 0.36. For the case with the apparent linear coupling absent (L0 ), the value of Gmax is 1.41 and the non-dimensional time taken to achieve this maximum growth is 0.28. Thus the maximum energy amplification and the time taken to achieve this maximum amplification are seen to decrease in this example when the apparent linear coupling is ignored. The Galerkin modes are basis functions used to describe the original state space of the system dynamics, and have limited physical significance. The linear coupling between the acoustic Galerkin modes, mentioned in the literature, is thus apparent and not entirely physical. The ε-pseudospectrum, index of non-normality through the Henrici ° ° number, plot of the norm °eLt ° and the maximum amplification of the energy Gmax
25
clearly show the change in non-normality of the system by the exclusion of the apparent linear coupling terms. As shown in Chapter 2, the non-normality in the system is because of the physical linear coupling of the eigenvectors. So, if the apparent linear coupling is ignored, the physical linear coupling between the eigenvectors may change, as seen in this case. This may have serious impact on the nonlinear effects such as triggering, bootstrapping, etc. (as will be seen in the following section), since they depend on the energy transfer between the eigenmodes of the system. The Gmax for the system in the presence of the apparent linear coupling terms is shown to be higher in the above example than the Gmax for the system with no apparent linear coupling. This might lead to a situation where the nonlinearities get triggered for a system with the apparent linear coupling terms present; however, the system with the apparent linear coupling neglected might not show triggering for similar initial conditions.
3.3
Nonlinear Rijke Tube Model
3.3.1
Impact of Linear Coupling on Nonlinearity
As noted in Section 3.1, Eq. 3.11 is linear in gas dynamics and nonlinear in heat release rate response. This second order ODE can be rewritten as two first order ODEs as follows 3 : dηj = η˙ j dt
2bkj dη˙ j = −2ξj kj η˙ j − kj2 ηj − dt γM
"s¯ ¯ r # ¯ ¯1 ¯ + u0 (t − τ )¯ − 1 sin(kj xf ) f ¯ ¯3 3
(3.37)
(3.38)
As discussed for Eqs. 3.35 and 3.36, ε0 and µ are used in the governing equations only to indicate the order of magnitudes. They are dropped while writing the evolution Eqs. 3.37 and 3.38. 3
26
The apparent linear coupling is due to the square root term;
q¯ ¯ ¯ 1 + u0 (t − τ )¯. For f 3
small values of u0f (t − τ ), the square root can be linearized as, s¯ ¯ ¯ ¯1 ¯ + u0 (t − τ )¯ = √1 (1 + 3 u0f (t − τ ) + higher order terms) f ¯ ¯3 2 3
(3.39)
It can be clearly seen from Eqs. 3.6 and 3.39 that the apparent linear coupling is provided by the 32 u0f (t − τ ) term. Hence, after the apparent linear coupling is removed, the nonlinear left hand side term in Eq. 3.39 can be expressed as "s¯ ¯# ¯1 ¯ 0 ¯ + u (t − τ )¯ f ¯3 ¯
no linear coupling
s¯ √ ¯ √ ¯1 ¯ 3 3 0 0 0 = ¯¯ + uf (t − τ )¯¯ − u (t − τ ) + u (t − τ ) (3.40) 3 2 f 2 jf
where u0jf (t−τ ) is the velocity contribution of the j th mode to the total acoustic velocity, or u0jf (t − τ ) = ηj (t − τ ) cos(kj xf )
(3.41)
Substituting Eqs. 3.40 and 3.41 in Eq. 3.11, we get the following oscillator equations for the case with the apparent linear coupling is absent: dηj = η˙ j dt
(3.42)
dη˙ j = −2ξj kj η˙ j − kj2 ηj dt "s¯ r # √ ¯ √ ¯1 ¯ 2bkj 3 3 1 0 0 ¯ + u0 (t − τ )¯ − u (t − τ ) + u (t − τ ) − sin(kj xf ) − f jf f ¯3 ¯ γM 2 2 3 (3.43) The nonlinear equations 3.37), 3.38 and 3.42, 3.43 are integrated numerically using a standard four step Runge-Kutta Method. The solution along with the Galerkin expansion given by Eq. 3.6 gives the evolution of the non-dimensional acoustic velocity. The evolution is obtained by exciting only the first Galerkin mode by initializing η1 at time t = 0, denoted by η1 (0). In all the cases, ηi=2...10 = 0 and η˙ i=1...10 = 0 .
27
3.3.2
Results and Discussions
The parameters that are kept fixed while presenting the results are Lw = 2.2m, Tw = 1000K, u¯ = 0.05m/s, S = 1.56 × 10−3 m2 , γ = 1.4, c0 = 399.6m/s, ρ¯ = 1.025kg/m3 , Cv = 719J/kgK, λ = 0.0328W/mK, dw = 0.0005m and La = 1.0m. The number of Galerkin The number of Galerkin modes used is N = 10 for all the results presented in this section, so that the change in solution with increase in number of modes is less than 5%.
Figure 3.5: shows the evolution of acoustic velocity at xf with non-dimensional time, for the case with the apparent linear coupling present. Initial conditions are ηi (0) = 0 for i = 2, 3, ..., N , η˙ i (0) = 0 for i = 1, 2, 3, ..., N and (a) η1 (0) = 0.1, (b) η1 (0) = 0.25
Similar to the results obtained by Balasubramanian and Sujith (2008b), Figure 3.5 shows the presence of triggering for the case with apparent linear coupling present. The damping coefficients used are c1 = 0.135 and c2 = 0.015. The localized heat release location is chosen as xf = 0.3. For η1 (0) = 0.1 (Figure 3.5a), the velocity shows a transient growth but eventually goes to zero asymptotically. For a duct with heat source, the natural modes are not the same as the eigenmodes of the system. Hence this initial excitation is not in the direction of a specific eigenmode and will have components along all the eigenmodes. Due to this non-normal phenomenon, a transient growth is seen in the evolution, but the amplitude of the initial excitation is not large enough to excite the nonlinearities and cause triggering. However, for η1 (0) = 0.25 (Figure 3.5b), the transient growth triggers the nonlinearity and eventually a limit cycle is attained. The 28
Figure 3.6: shows the evolution of acoustic velocity at xf with non-dimensional time, for the case with the apparent linear coupling absent. Initial conditions are ηi (0) = 0 for i = 2, 3, ..., N , η˙ i (0) = 0 for i = 1, 2, 3, ..., N and (a) η1 (0) = 0.1, (b) η1 (0) = 0.25, (c) η1 (0) = 1.0
phase difference between the heat release and the acoustic pressure evolves to a value such that the driving exactly balances out the damping in the system and hence a saturation of the amplitude is seen. Figure 8 shows that this triggering phenomenon is not observed for the case where the apparent linear coupling is ignored. Triggering is not observed for either of the initial conditions η1 (0) = 0.1 (Figure 3.6a) and η1 (0) = 0.25 (Figure 3.6b). Even if the initial amplitude is increased to as high as η1 (0) = 0.1 (Figure 3.6c), triggering is not observed. The velocity perturbations die down asymptotically. It can be inferred that, for this example, inclusion of the apparent linear coupling is critically important to accurately predict the triggering phenomena. As seen in Chapter 3.2, ignoring the apparent linear coupling may result in altering the physical linear coupling among the eigenmodes. This physical linear coupling is crucial for the energy exchange among the modes, which may play an important role in triggering. The acoustic velocity projected on some modes may be decaying whereas the projection on some other modes may be growing. When sufficient energy is projected on the growing modes, they may start transferring the energy back to the decaying modes. This phenomenon was referred to as bootstrapping by Yoon et al. (2001). Bootstrap29
Figure 3.7: shows the evolution of acoustic velocity at xf with non-dimensional time. Initial conditions are ηi (0) = 0 for i = 2, 3, ..., N , η˙ i (0) = 0 for i = 1, 2, 3, ..., N and (a) for the case with apparent linear coupling present, η1 (0) = 0.02, (b) for the case with apparent linear coupling absent, η1 (0) = 0.02
ping results in a shift in the dominant frequency during the evolution. Figure 3.7 shows the nonlinear phenomenon of bootstrapping in a Rijke tube. The damping coefficients used are c1 = 0.05 and c2 = 0.005. The localized heat release location is chosen as xf = 0.25. The initial condition chosen for both Figure 3.7a and Figure 3.7b is η1 (0) = 0.02. Figure 3.7a shows the evolution of velocity in the presence of the apparent linear coupling and Figure 3.7b shows the evolution of velocity in the absence of the apparent linear coupling. It is clearly seen that the bootstrapping phenomenon is not predicted if the apparent linear coupling between the Galerkin modes is ignored. The energy from one mode is transferred to another mode due to the existence of the physical coupling among the eigenmodes. The velocity projected on the first Galerkin mode is decaying in the beginning of the evolution. However, once the energy has been transferred from one mode to another, there is again a growth of oscillations. This triggers the nonlinearities and eventually saturation is attained as the phase difference between the heat release and the acoustic pressure evolves to a point where the driving and the damping cancel exactly (Figure 3.7a). However, as argued in the case of triggering, neglecting the apparent linear coupling implies that the pathways of energy exchange among the eigenmodes are altered and phenomenon such as bootstrapping is not observed.
30
The nonlinearity in King’s law is square root in nature. The perturbations are very small and hence the physical linear coupling along with the nonlinearities together may provide a pathway for energy exchange among the eigenmodes of the system, which may result in phenomena such as triggering and bootstrapping. If the apparent linear coupling is ignored, pathway of the energy exchange due to the physical linear coupling may get altered and the energy exchange brought about by nonlinearities due to higher order coupling alone might not be sufficient to exhibit triggering and bootstrapping effects. Thus, the inclusion of the apparent linear coupling may be critical for accurately capturing the nonlinear effects of the system as shown in the examples. Ignoring it may result in a failure to capture the essential dynamics of the system.
3.4
Remarks
The approximation of ignoring the apparent linear coupling provides analytical simplicity, especially in the case of solid rocket motors where it is tedious to get closed form solutions with the coupling terms. Ad-hoc combustion nonlinearities are then added to the thermoacoustic system, without the apparent linear coupling terms, to obtain bifurcation plots and to study the nonlinear system dynamics (Ananthkrishnan et al., 2005; Wicker et al., 1996). The effects of this assumption, in the case of a distributed heat release (solid rocket motor), is not investigated in this thesis. A natural extension of the present work is to investigate the effect of this approximation in the presence of combustion nonlinearities in a solid rocket motor. In the present thesis, it is shown that the approximation may not be valid in the presence of combustion nonlinearities with a localized heat release. In the heat release rate model used in this thesis, the linear coupling and the combustion nonlinearities are both present in King’s law and hence it may seem tedious and impertinent to ignore the apparent linear coupling terms. However, this may not be the case in other thermoacoustic systems.
31
CHAPTER 4
IMPACT OF STOCHASTIC PROCESSES ON THERMOACOUSTIC INSTABILITIES IN A RIJKE TUBE
Figure 4.1: Thermoacoustic System with input - output viewpoint
A non-normal thermoacoustic system can be visualized as a system with some transfer function that transforms the input perturbation into output perturbations. There may be external forcing present in the system too. Figure 4.1 shows a schematic of such a system. The input, external forcing and the system itself can be deterministic or random (stochastic). The thermoacoustic system, considered in Chapter 3, is deterministic with no external forcing and a deterministic input or initial condition. In this chapter, the thermoacoustic system considered is linear and stochastic, with no external forcing and with a deterministic input or initial condition.
4.1
Elementary Definitions and Concepts
In this section, relevant basic definitions in Statistics and Random theory are reviewed. Any introductory book on this subject can be read for detailed explanation and examples of these concepts (Papoulis and Pillai, 2002; Chung, 1975).
4.1.1
Probability Space
A triple (S, E, P ) is called a probability space such that S is a non-empty set whose elements are known as outcomes. A set of all possible outcomes is called as sample space. E is a collection of subsets of S or the σ-algebra of the subsets of S and P is the probability measure on E. The elements of E are known as events. For example, consider a probability space defined by a flip of a fair coin. The possible outcomes are S = {H, T } , where H denotes heads and T denotes tails. E = {{H}, {T }, {}, {H, T }} and the corresponding probability for the elements of E is, P ({H}) = P ({T }) = 1/2, P ({}) = 0, i.e., the probability of getting null outcome after tossing the coin is zero and P ({H, T }) = 1, i.e.,the probability of getting either a heads or a tails after tossing the coin is 1. Here, the set of possible events or the sample space is S = {H, T }.
4.1.2
Random Variable or Stochastic Variable
It is a mapping from a sample space (called probability space) into a real line (random variable or stochastic variable). Complex valued random variable is also defined similarly. Let X : S → R be the function mapping S into the real line such that for each s ∈ S, there exists a unique X(s) ∈ R. Then, X is called a random variable or a stochastic variable. Consider a sample space defined by outcomes of tossing 2 coins, i.e., S = {HH, HT, T H, T T }. Let the mapping be such that the corresponding values of the random variables are say X = {0.1, 0.2, 0.3, 0.4}. The corresponding probabilities will then be P = {1/4, 1/4, 1/4, 1/4} since all the four possible outcomes in set S are 33
equally probable. Alternatively, X(HH) = 0.1 and P (X = 0.1) = 1/4 and so on. This is an example of a discrete random variable. Continuous random variable can be understood analogously.
4.1.3
Probability Distribution Function and Probability Density Function
The cumulative probability is called the probability distribution function, FX (x) = P (X ≤ x) For the above example, FX (0.2) = P (X ≤ 0.2) = 1/4 + 1/4 = 1/2. Similarly, FX (0.1) = 1/4, FX (0.3) = 3/4 and FX (0.4) = 1 The distribution of the probability associated with a random variable P (X = x) versus the random variable X = x is called the probability density function or the PDF, fX (x). Alternatively, the derivative of the probability distribution function is the PDF, provided it exists. Mathematically, Z
x
FX (x) =
fX (x)dx −∞
or fX (x) =
d FX (x) dx
The PDF satisfies a basic property of probability that
4.1.4
R∞ −∞
fX (x)dx = 1.
Joint Distribution Function and Joint Density function
If X and Y are two random variables defined on the same sample space S, then the joint distribution function is defined as, FX,Y (x, y) = P (X ≤ x, Y ≤ y)
34
The joint density function is defined as ∂2 FX,Y (x, y) ∂x∂y
fX,Y (x, y) = provided it exists. Alternatively, Z
Z
x
y
FX,Y (x, y) =
fX,Y (x, y)dxdy −∞
4.1.5
−∞
Mean, Variance and Standard Deviation
The mean or the expected value of a distribution of random variables, E(X), is the average value of the random variables. Mathematically, Z
∞
E(X) =
xfX (x)dx −∞
The variance, V ar(X), is one measure of statistical dispersion, averaging the squared distance of the values of a random variable from the expected value (mean). Mathematically, V ar(X) = E((x − E(X))2 ) or
Z
∞
V ar(X) =
(x − E(X))2 fX (x)dx
−∞
The variance is a way to capture the scale or degree of being spread out whereas the mean is a way to describe the location of a distribution. The variance is often denoted as σ(X)2 also. The Standard Deviation, σ(X), is the positive square root of the variance.
4.1.6
Auto-correlation and Auto-Covariance
The auto-correlation between two random variables defined on the same sample space, S, is measured by Cor(X, Y ) = E(xy) 35
Alternatively,
Z
∞
Z
∞
Cor(X, Y ) =
xyfX,Y (x, y)dxdy −∞
∞
The auto-covariance between two random variables defined on the same sample space, S, is measured by Cov(X, Y ) = E(x − E(X))(y − E(Y )) Alternatively, Z
∞
Z
∞
Cov(X, Y ) =
(x − E(X))(y − E(Y ))fX,Y (x, y)dxdy −∞
∞
The auto-covariance of uncorrelated random variables X and Y is zero.
4.1.7
Cumulants
The following cumulant generating function g(t) is used for defining cumulants. tX
g(t) = log(E(e )) =
∞ X n=1
cn
tn n!
The coefficient cn is the cumulant of nth order. For a distribution, the mean and the variance are the first and the second cumulants respectively.
4.1.8
Random Process or Stochastic Process
A real random variable maps each point in the sample space to a point on the real line. Similarly, a Random Process or Stochastic Process maps each point in the sample space to a waveform. A random process is an indexed family of random variables denoted by X(s, t) where s is a sample point and t is an index set usually denoting time. The sample space is constant and does not vary with time. For a fixed t, X(s) is a random variable and for a fixed point in sample space s, X(t) is a single realization of the random process and is a deterministic function. Random process can be defined in a variety of ways, one of them being the specification of the autocorrelation function. 36
4.1.9
Gaussian Distribution
The Gaussian distribution or the normal distribution is one of the most commonly used distributions by engineers. The PDF is given by ¶ µ 1 (x − E(X))2 √ exp − fX (x) = 2σ(X)2 σ(X) 2π The Probability distribution function is given by · µ ¶¸ 1 x − E(X) √ FX (x) = 1 + erf 2 σ(X) 2 where 2 erf (y) = √ π
Zy 2
e−t dt 0
Figure 4.2: Gaussian Probability Density Function for zero mean and different variances The Gaussian PDFs are shown in Figure 4.2 for zero mean and different variances. The Gaussian PDF is also known as a bell curve due to its typical shape. The probability distribution functions corresponding to Figure 4.2 are shown in Figure 4.3. A Gaussian distribution with zero mean and unit variance is also known as the "standard normal distribution". The cumulants of order greater than 2, of a gaussian distribution, are 0.
37
Figure 4.3: Gaussian Probability Distribution Function for zero mean and different variances
4.2
Linearized Governing Equations
The linearized equations developed for a horizontal Rijke tube in chapter 3 are used in these chapter. The fluctuating heat release given in Eq. 3.7 can be linearized and written in the time lag , n − τ , form as (γ − 1)Q˙ 0 = nu0 (x, t − τ )δ(˜ x − x˜f )
(4.1)
where the interaction index n is given by √ b 3 n= 2
(4.2)
where b is defined in Eq. 3.9. The popular time lag n − τ model, proposed by Crocco (1952), has been useful in quantitative prediction of combustion instabilities in a variety of cases. This linearized heat release model is used in this chapter to study the impact of stochastic sources in thermoacoustic system. Similar analysis can be extended to other linear thermoacoustic systems as well. The values of n chosen here are obtained using representative values used in chapter 3. For the open - open boundary conditions, using the Galerkin modes (Eq. 3.6) and
38
the procedure shown in chapter 3, the oscillator equations can be easily obtained as dηj d2 ηj 2nkj 0 2 + ξ + k η = − u (xf , t − τ ) sin(jπxf ) j j j dt2 dt γM
(4.3)
where ξj is defined in Eq. 3.12. Splitting each oscillator equation into two first order ordinary differential equation and using Taylor series expansion, Eq. 4.3 can be written in the dynamical system equation of the form of Eq. 3.29 with χ defined in Eq. 3.30 and the linear stability operator L as
0
−k1
k1 + nβ1 C˜1 ξ1 − τ k1 nβ1 C˜1 0 0 nβ2 C˜1 −τ k2 nβ2 C˜1 L= ......... ......... ......... ......... 0 0 nβN C˜1 −τ kN nβN C˜1
..... .....
0
0
..... .....
nβ1 C˜N
−τ k1 nβ1 C˜N
..... .....
0
0
..... .....
nβ2 C˜N
−τ k2 nβ2 C˜N
..... .....
....
.........
..... .....
....
.........
..... .....
0
−kN
..... ..... kN + nβN C˜N
ξN − τ kN nβN C˜N (4.4)
where C˜j = cos(kj xf ) and βj = (2/γM ) sin(jπxf ). L is a deterministic or a certain operator. However, in combustion instabilities, complete knowledge of the deterministic system is often lacking. Also, from experimental observations, it has been shown that the operator governing the thermoacoustic system dynamics are often stochastic or uncertain. The uncertainties can be additive or multiplicative in nature. dχ = M 0 χ + A0 dt
(4.5)
where M 0 denotes the multiplicative stochastic operator and A0 denotes the additive stochastic operator. For simplicity, we assume that our system is multiplicatively uncertain but additively certain, i.e. , A0 = 0 here. We also assume that M 0 can be split into a mean certain operator and a fluctuating operator which is noisy or M 0 = Lχ+αζ(t)F . Thus consider the uncertain linear system for the state space vector χ dχ = Lχ + αζ(t)F χ dt 39
(4.6)
where L is the mean stability operator, α is the amplitude of noise level, F is the matrix of fluctuating structure. ζ(t) is a Gaussian stochastic process with zero mean and unit variance and autocorrelation time tc = 1/νc such that hζ(t1 )ζ(t2 )i = exp(−νc |t1 − t1 |)
(4.7)
F can be obtained by assuming noise in parameters like the interaction index n, time lag τ or the damping factor ξ. For mathematical convenience it is assumed that only one parameter is noisy or stochastic so as to avoid dealing with the interaction among multiple noises in the system. The fluctuating matrix F for noise in n and τ can be represented as follows:
Fn =
0
0
..... .....
0
nβ1 C˜1
ξ1 − τ k1 nβ1 C˜1
0
0
nβ2 C˜1
−τ k2 nβ2 C˜1
.........
.........
..... .....
....
.........
.........
.........
..... .....
....
.........
0
0
..... .....
0
0
nβN C˜1
−τ kN nβN C˜1
..... ..... nβ1 C˜N ..... .....
0 −τ k1 nβ1 C˜N
0
..... ..... nβ2 C˜N
0 −τ k2 nβ2 C˜N
..... ..... nβN C˜N ξN − τ kN nβN C˜N
(4.8)
0
0
0 −τ k1 nβ1 C˜1 0 0 0 −τ k2 nβ2 C˜1 Fτ = ......... ......... ......... ......... 0 0 0 −τ kN nβN C˜1
..... .....
0
..... .....
0
0
−τ k1 nβ1 C˜N ..... ..... 0 0 ..... ..... 0 −τ k2 nβ2 C˜N ..... ..... .... ......... ..... ..... .... ......... ..... ..... 0 0 ˜ ..... ..... 0 −τ kN nβN CN
(4.9)
A system of deterministic ensemble mean evolution equations can be obtained from
40
the system of stochastic evolution equation given in Eq. 4.6. See Appendix D for complete derivation.
¢ dhχi ¡ = L + α2 F Θ(t) hχi dt
where
(4.10)
Zt eLs F e−Ls e−νc s ds
Θ(t) =
(4.11)
0
For short autocorrelation times tc << 1, Eq. 4.10 can be approximated as (see Appendix D)
· ¸ dhχi α2 2 = L + F hχi dt νc
(4.12)
From Eq. 4.12, it can be seen that asymptotically, the mean operator L has been modified by
α2 2 F . νc
For imaginary eigenvalues of F , this term contributes to additional
damping asymptotically and for real eigenvalues of F it contributes to additional driving or growth asymptotically. Appendix C deals with the procedure to determine the maximum transient energy growth and the optimal initial conditions for maximum energy growth for deterministic, autonomous, linear dynamical systems governed by non-normal operators. The procedure to obtain similar quantities for uncertain systems is derived in Appendix E. The maximum of the ensemble evolution of the perturbation square amplitude hχ† (t)χ(t)i is corresponds to the maximum transient energy growth. The method described in Appendix C cannot be used on the ensemble mean equation obtained in Appendix D (Farrell and Ioannou, 2002a,b). The covariance dynamics has to be used to obtain the required optimals. The optimal initial conditions of the stochastic dynamical system given by Eq. 4.6 at time t is the eigenvector associated with the largest value of the Hermitian ensemble matrix hS(t)i which is obtained at time t by integrating the integral-differential equation given by (see Appendix E for complete derivation) dhSi = (L+α2 Ξ(t)F )† hSi+hSi(L+α2 Ξ(t)F )+α2 (Ξ(t)† hSiF +F † hSiΞ(t)) (4.13) dt
41
Figure 4.4: Optimal initial conditions , u0 (x, 0) and p0 (x, 0), for maximum transient energy growth: for randomness in the interaction index, n
where
Zt e−Lx F eLx e−νc x dx
Ξ(t) =
(4.14)
0
with initial condition as hS(0)i = I. For short autocorrelation times tc << 1, Eq. 4.13 can be approximated as (see Appendix E) dhSi = dt
µ
α2 L+ F νc
¶†
µ
α2 hSi + hSi L + F νc
¶
2α2 † + F hSiF νc
(4.15)
The short autocorrelation approximation is known as equivalent white noise approximation. The equations with no approximation are known as color noise equations.
4.3
Results and Discussions
This section discusses few results obtained using the equations derived above. The results of equivalent white noise approximation, color noise and a completely deterministic system will be compared. Equations 4.10 , 4.12, 4.13 and 4.15 are numerically solved using simple forward difference schemes and trapezoidal integration technique. Figures 4.3, 4.3 and 4.3 show the optimal initial conditions projected on the spatial velocity and pressure distribution. These plots compare the optimal initial conditions obtained for a deterministic (Appendix C) and stochastic system (Appendix E). The optimals obtained by using the equivalent white noise approximation is also compared to cross check the validity of the approximation. The parameters that are fixed for all 42
Figure 4.5: Optimal initial conditions , u0 (x, 0) and p0 (x, 0), for maximum transient energy growth for high autocorrelation time: for randomness in the interaction index, n
the plots in this section are N = 8; xf = 0.3; c2 = 0; n = 0.0003; c1 = 0.05 and α = 0.3 Figure 4.3 is obtained for a case with randomness in the interaction index n. τ = 0.002; and tc = 0.1. The optimal initial conditions are almost identical for all the 3 cases: deterministic, equivalent white noise stochastic and stochastic. The maximum growth factors are 1.81 for deterministic system, 1.86 for equivalent white noise case and 1.82 for the stochastic system. The equivalent white noise approximation clearly holds good. However, it can be seen in Figure 4.3 that this approximation fails for same conditions as Figure 4.3, but with high autocorrelation time tc = 2. The maximum growth factors corresponding to Figure 4.3 are 1.78 for deterministic system, 2.78 for equivalent white noise case and 1.78 for the stochastic system. Figure 4.3is obtained for the case of randomness in the time lag τ . Here, τ = 0.01 and the remaining parameters are same as Figure 4.3. Even here, the optimal initial conditions are almost identical for the 3 cases as seen in Figure 4.3. The maximum growth factors corresponding to Figure 4.3 are approximately 1.63 for all the 3 cases. Similar plots were obtained for a large combination of parameters. In all the computations, the optimal initial conditions were almost identical like in the other figures shown in this section. The thermoacoustic system considered here is a very simplified model. Also, the nature of the noise is highly simplified, with only multiplicative noise present and only one of the parameters is assumed to be noisy. The identical optimal 43
Figure 4.6: Optimal initial conditions , u0 (x, 0) and p0 (x, 0), for maximum transient energy growth: for randomness in the time lag , τ
initial conditions for the deterministic and stochastic cases may be attributed to these reasons. However, this may not be the case for other thermoacoustic systems. These optimals vary quite significantly for noisy atmospheric flows (Farrell and Ioannou, 2002b). The similarity in the maximum growth factors for the deterministic and stochastic case ensures the robustness of the computations in obtaining maximum growth factors for the deterministic system. These results should not be too sensitive to noise in the system, which is evident from the numbers shown above.
4.4
Remarks and Further Scope
Ergodicity Each realization of the stochastic differential equation is the evolution of the system with a random process corresponding to one sample point in the sample space. The ensemble mean is the average of a large number of such realizations. The ensemble mean is not defined for a single realization. But it is of theoretical and practical interest to have an understanding of the ensemble mean with reference to a single realization. If the system reaches a stationary state, the time mean is same as the ensemble mean or ¯ = hχ(t)i. This is known as the ergodic assumption and it has been verified using χ(t) Monte Carlo simulations for uncertain systems excited by both additive and stochastic forcing (Farrell and Ioannou, 2002a).
44
Further Scope The interaction between stochastic, linear and nonlinear processes is very complex and is not well understood by the combustion instability community at present. The optimal initial conditions for maximum energy growth were not too different from the corresponding optimals of the deterministic thermoacoustic system for the case studied here. However, it may not be the case with other linearized thermoacoustic system having the unsteady combustion due to premixed flames and diffusion flames. The impact of stochastic sources on such thermoacoustic system needs to be addressed. The effect of additive noise may also be crucial. In this thesis, the assumption of only a single stochastic process was made. However, in practical systems, there may be many stochastic processes simultaneously present. This may significantly affect the stability of the system. Extensive numerical Monte Carlo simulations are however required to carry out such computations. Realistic values of the order of autocorrelation times and the amplitude of the noise needs to be experimentally obtained. Effects of stochastic processes in the presence of nonlinearities needs to be investigated. Getting analytical deterministic expressions may not be possible in the presence of nonlinearities. Monte Carlo simulations using Stochastic Runge-Kutta methods or Polynomial Chaos techniques may be used to solve the stochastic differential oscillator equations. Noise induced triggering needs to be theoretically and numerically studied. The impact of stochastic sources on Bifurcation plots may also be interesting.
45
CHAPTER 5
CONCLUSIONS
Non-normality of a thermoacoustic system is due to the presence of the coupling between its eigenmodes. The term "linear coupling" used in literature refers to the linear coupling between the acoustic Galerkin modes (referred to as apparent linear coupling in this thesis). It may not represent the linear coupling between the eigenmodes (referred to as physical linear coupling in this thesis). Ignoring the apparent linear coupling does not result in a significant shift in the eigenfrequencies of the system. However, this approximation can alter the physical linear coupling, and hence alter the non-normality of the system. The physical linear coupling provides a pathway for energy exchange between the eigenmodes. This may result in significant change in the linear and nonlinear system dynamics.
The Henrici number measures the departure from normality of a dynamical system. However, it alone is insufficient to quantify and compare nonnormal systems. In the numerical example of a linearized Rijke tube, Henrici number decreases, by a very small amount, on ignoring the apparent linear coupling. The approximate system with no apparent linear coupling shows considerably lesser maximum energy amplification than the actual system with the apparent linear coupling present. As shown by a numerical example, this can lead to different consequences in the presence of nonlinearities. The nonlinearity in the Rijke tube is modeled here by King’s law, which is square root in nature. Under the assumption of small velocity perturbations, the Taylor series expansion provides the apparent linear coupling terms. Strongly nonlinear phenomena such as bootstrapping and triggering may be completely missed out if these terms are ignored. Ignoring the apparent linear coupling is justified only when the change in the system dynamics is insignificant.
With appropriate approximations for the type of noise in the system, it is possible to derive analytical expressions for the deterministic ensemble mean equations from the stochastic evolution equations. The maximum transient energy growth for deterministic and stochastic systems are close to each other thus proving the robustness of the calculations of the deterministic system. The equivalent white noise approximation can be used for short autocorrelation times. The optimal initial conditions for maximum transient energy growth, for the linearized n − τ system considered here, do not significantly differ from its corresponding optimal deterministic initial conditions. However, this result cannot be generalized to all thermoacoustic systems since these results may be system specific. Hence extension of this analysis to other thermoacoustic systems is required. The study of the interaction between the stochastic sources and system nonlinearities is also warranted.
47
APPENDIX A PSEUDOSPECTRA For a linear operator L which is normal or close to being normal, k(ωI − L)−1 k (the resolvent of L) is large only when ω is close to an eigenvalue of the system. In such cases, the resonant frequencies can be found from the spectrum. However for a non-normal system, k(ωI − L)−1 k can be large even when ω is far from any eigenvalue. This is known as pseudoresonance (Trefethen et al., 1993).Pseudospectra are plot of contours with equal resonant magnitudes. Pseudospectra can be used to study the evolution of dynamical systems governed by non-normal operators (Trefethen and Embree, 2005). The ε-pseudospectrum can be defined as ° ° z is an ε - pseudoeigenvalue of L if it satisfies °(zI − L)−1 ° ≥ ε−1
(A.1)
Alternatively, z is an eigenvalue of the matrix L + E, where E is a matrix of small perturbations such that kEk ≤ ε (see Trefethen and Embree (2005) for other equivalent definitions of pseudoeigenvalues). ε is a measure of the perturbation over the nonnormal operator L , or, ε ≤ kLk. The ε-pseudospectrum consists of disks of radius ε centered at each eigenvalue. It can be used to study the transient growth (Trefethen et al., 1993). The geometry of the pseudospectra reveals the non-normality and the transient growth factor. The ε-pseudospectrum of normal operators are closed circles. A dynamical system may exhibit transient growth when a contour corresponding to some ε value does not lie entirely in the left-half plane. The spilling of the pseudospectra in the right half plane is a necessary but not a sufficient condition for transient growth. The transient growth is a function of time and is different for different initial conditions. The maximum energy growth rate, maximized over all initial conditions and all times is given by
¡ ¢ Gmax = max kexp(Lt)k2 t
(A.2)
A system is stable if the real parts of all the eigenvalues are negative and if the maximum growth factor Gmax is 1. Such a system will not grow for any initial condition. Gmax > 1
is a necessary and a sufficient condition for a system to exhibit transient growth for certain initial conditions. Lower bounds for the maximum energy growth can be obtained using pseudospectra geometry (Trefethen and Embree, 2005). The maximum of the real parts of the ε-pseudoeigenvalues of L is known as the pseudospectral abscissa, denoted mathematically by γε = max (Re z). Here, Λε (L) denotes the ε-pseudospectrum of L. The z∈Λε (L)
ratio of the pseudospectral abscissa to the corresponding ε forms a lower bound for the maximum energy growth (Trefethen and Embree, 2005). p
Gmax = max (kexp(Lt)k) ≥ γε /ε
∀ε > 0
t
(A.3)
Thus if an ε- pseudospectrum of L has the ratio γε /ε > 1, there will be a transient growth in the system for certain initial conditions. It must be noted that the ratio γε /ε is a lower bound and not the upper bound to the maximum energy growth. Hence γε /ε > 1 is a sufficient, but not a necessary condition for the transient growth to occur. Maximizing this ratio over all the ε provides the Kreiss constant, κ. κ(L) = max(γε /ε)
(A.4)
ε
Kreiss constant thus forms a stronger lower bound to the maximum energy growth compared to the other ratios since it is maximized over all the ε values. Alternatively, mathematically the lower bound inequality can be expressed as p Gmax ≥ κ(L) ≥ γε /ε
∀ε > 0
(A.5)
Note that upper bounds to the maximum energy growth can also be obtained. For these definitions, see Trefethen and Embree (2005). In summary, the Kreiss constant being greater than one is a sufficient condition, spilling of the pseudospectra contours to the right half plane is a necessary condition and Gmax > 1 is both a necessary and a sufficient condition for the transient energy growth to occur in the system, for certain initial conditions.
49
APPENDIX B
EFFECT OF NON-IDEAL BOUNDARY CONDITIONS ON GALERKIN MODES
The mode shapes are the eigenfunctions of the homogenous wave equation. For nonideal boundary conditions, the mode shapes can be derived. The following is the derivation of the mode shapes for purely complex impedances at the boundary. The mode shape obtained here appropriately changes the L matrix of the non-normal system. Nonnormal analysis can then be done on the new L matrix. The unperturbed wave equation 2 0 ∂ 2 p0 2∂ p − c =0 0 ∂t2 ∂x2
(B.1)
The solution p0 (x, t) to Eq. B.1 is of the form p0 (x, t) = Re(ˆ p(x)eiωt )
(B.2)
¯ p0 ¯¯ = iZj (ω) i = 1, 2 u0 ¯xj
(B.3)
The non-ideal BC is given by,
where 1 denotes the inlet and 2 denote the outlet. Z1 and Z2 are assumed to be real functions of ω. The impedance is non-dimensionalized by ρ¯c0 . Hence the boundary conditions are purely complex impedance conditions. The wave equation, Eq. B.1, is simplified to d2 pˆ + k 2 pˆ = 0 2 dx where k =
ω . c0
(B.4)
The solution of Eq. B.4 is of the form pˆ = aeikx + be−ikx
(B.5)
The unperturbed momentum equation is ρ¯
∂u0 ∂p0 + =0 ∂t ∂x
(B.6)
Equation B.6 can be integrated to obtain u0 =
i dˆ p iωt e ρ¯ω dx
(B.7)
Combining Equations B.7 and B.3 ρ¯ω dˆ p (xj ) = − pˆ(xj ) j = 1, 2 dx Zj (ω)
(B.8)
Equation B.8 is substituted in Eq. B.5 and a system of linear equations is obtained as ikeikx1 a − ike−ikx1 b = −
¤ ρ¯ω £ ikx1 e a + e−ikx1 b Z1 (ω)
(B.9)
ikeikx2 a − ike−ikx2 b = −
¤ ρ¯ω £ ikx2 e a + e−ikx2 b Z2 (ω)
(B.10)
From Eqs. B.9 and B.10, nonzero values of a and b can be found if, µ
ρ¯2 c20 1+ Z1 (ω)Z2 (ω)
¶
ωL sin( )+ c0
µ
ρ¯c0 ρ¯c0 − Z1 (ω) Z2 (ω)
¶ cos(
ωL )=0 c0
(B.11)
Equation B.11 can be solved to obtain admissible values of ω. Equation B.11 will give only real values of ω. For a purely complex impedance, the eigenfrequencies of the homogenous equation (ω) are real. (Similar conclusion was arrived at in (Nicoud et al., 2007)). The frequency being real, there is no dissipation for such boundary conditions. Energy absorption or dissipation is a function of the real part of the impedance alone (Dowling and Ffowcs-Williams, 1983). Equation B.10 is used to obtain the relation between a and b. We obtain the following equation, 1 ³
pˆ(x) = r k2 +
· ρ¯ω Z1 (ω)
´2
¸ ρ¯ω sin(kx) k cos(kx) + Z1 (ω)
51
(B.12)
It can be seen that pˆ(x) is real. pˆ(x) can also be written as pˆ(x) = sin(kx + φ) where k =
ω c0
and tan φ =
(B.13)
k Z (ω). ρ¯ω 1
Thus, for purely imaginary impedance boundary conditions, the mode shape is of the form of Eq. B.13. For an open-open boundary condition, Z1 = Z2 = 0, giving the standard pressure mode shape sin(kx). For a closed-closed boundary condition, Z1 = Z2 → ∞, giving the standard pressure mode shape cos(kx). It must be noted that Z2 affects the solution of Eq. B.11; hence the wave number k of the mode shape. As shown by Nicoud et al. (2007), only for purely complex impedance boundary conditions, the Galerkin modes (eigenvectors of the homogenous wave equation) are orthogonal. Hence Galerkin projection would form a complete basis and convenient to work with. There is no dissipation by the purely complex impedance boundary conditions since the eigenfrequencies are real. For general complex impedance boundary conditions, dissipation will occur (Dowling and Ffowcs-Williams, 1983; Nicoud et al., 2007). For such boundary conditions, the eigenvectors of the homogenous wave equations will not be orthogonal, (Nicoud et al., 2007), and hence the Galerkin method can hardly be generalized (Nicoud et al., 2007). This is one of the disadvantages of the Galerkin technique. The internal dissipation due to viscosity (acoustic boundary layer) is neglected. Thus, the assumption of neglecting the dissipation due to complex impedance boundary condition for near open-open boundary conditions and for near purely imaginary complex impedance conditions is reasonable. To take into account the effect of damping due to non-ideal boundary conditions and viscosity in the Galerkin modal analysis, damping models may be introduced. In this way, the simplicity of the Galerkin technique can be kept intact and the dissipation effects can be accounted for. The constants for such damping models may be found experimentally. One such damping model is proposed by Matveev and Culick (2003). They ignore the non-orthogonality of the eigenmodes of the homogenous wave equation and propose a linear mode damping model.
52
APPENDIX C
MAXIMUM TRANSIENT ENERGY GROWTH AND OPTIMAL INITIAL CONDITIONS FOR DETERMINISTIC SYSTEMS
Consider the non-normal dynamical system given by Eq. 3.29, where the state space vector, χ(t) is given by Eq. 3.30. The 2-norm of the state space vector, kχ(t)k, can be expressed as kχ(t)k2 =
N X
Ã
i=1
η˙ i (t)2 ηi (t)2 + ki2
! (C.1)
The non-dimensional acoustic energy at any time t can be written as Z1 E(t) =
1 0 1 p (x, t)2 + γM u0 (x, t)2 dx 2 2
(C.2)
0
On substituting the Galerkin expansions from Eq. 3.6 and the norm from Eq. C.1 in Eq. C.2, we get
Z1 µ E(t) =
γM 2
¶2 kχ(t)k2
(C.3)
0
The solution χ(t) is given by Eq. 3.33. Hence the greatest possible energy growth at time t, maximized over all possible initial perturbations is measured by the growth function defined as follows (Reddy and Henningson, 1993; Schmid and Henningson, 2001)
kχ(t)k2 2 G(t) = max 2 = kexp(Lt)k χ(0) kχ(0)k
(C.4)
For any general matrix A, the 2-norm can be geometrically interpreted as the radius of the smallest circle that contains the image of a unit disk under the map A. Mathematically, the 2-norm of any matrix is its principal singularvalue. It can be computed by
singularvalue decomposition (SVD), using standard subroutines available in most software libraries. The corresponding right singularvector gives the most sensitive or the optimal initial condition for maximum energy growth. It must be emphasized that the maximum amplification obtained by the above procedure is for a particular time instant t maximized over all initial conditions (Eq. C.4). Thus, in a given time interval [0, t], the global maximum amplification factor is given as follows: ¡ ¢ Gmax = max (G(t)) = max kexp(Lt)k2 t
t
(C.5)
An excellent geometric interpretation of the optimal initial condition required for maximum growth is provided by Farrell and Ioannou (1996a) in the context of atmospheric flows. Note that the above definitions are true for autonomous stability operators only. For non-autonomous operators, the expressions to calculate the maximum transient energy growth and the corresponding optimal initial conditions are different and can be referred in Schmid and Henningson (2001). Appendix E lays out a procedure to obtain these growth factors and optimal initial conditions in the presence of multiplicative stochastic sources in the system.
54
APPENDIX D
DETERMINISTIC EVOLUTION EQUATION OF THE ENSEMBLE MEAN FROM STOCHASTIC EVOLUTION EQUATION
Approximate deterministic dynamical system equations for the ensemble mean evolution of an uncertain system is obtained by van Kampen (1992), Farrell and Ioannou (2002a) and Farrell and Ioannou (2002b). This is not restricted to white noise stochastic processes and is valid in general for colored noise. In this section, the evolution equation of the ensemble mean field are derived. Consider the uncertain linear system for the state space vector χ dχ = Lχ + αζ(t)F χ dt
(D.1)
where L is the mean stability operator, α is the amplitude of noise level, F is the matrix of fluctuating structure. ζ(t) is a Gaussian stochastic process with zero mean and unit variance and autocorrelation time tc = 1/νc such that hζ(t1 )ζ(t2 )i = exp(−νc |t1 − t1 |)
(D.2)
Consider the interaction perturbation φ(t) defined by χ(t) = eLt φ(t)
(D.3)
dφ = αζ(t)e−Lt F eLt φ = αζ(t)H(t) dt
(D.4)
H(t) = e−Lt F eLt
(D.5)
Equation D.1 is then given by
where
For a sure (certain) initial condition, φ(0) = χ(0) and φ(t) is then given by φ(t) = G(t)χ(0)
(D.6)
with the fundamental matrix G(t) given by the time-ordered exponential
Zt
G(t) = expo α
ζ(s)H(s)ds
(D.7)
0
or in expansion series as Zt
Zt ζ(t1 )H(t1 )dt1 + α2
G(t) = I+α 0
Zt1 dt1
0
ζ(t1 )ζ(t2 )H(t1 )H(t2 )dt2 + ...... (D.8) 0
The time ordered exponential is the exponential with the convention that the terms in the series expansion are grouped in ascending time order. Kubo (1962) has shown that the ensemble average of G(t) can be expanded in the cumulants of ζ(t) in the same way that the exponents of ordinary stochastic processes are expanded in cumulants of their arguments
α2 hG(t)i = expo ( 2!
Z t Zt1 hhζ(t1 )ζ(t2 )iiH(t1 )H(t2 )dt1 dt2 + .... 0
2n
+
α n!
Zt
0
... 0
(D.9)
Ztn hhζ(t1 )...ζ(tn )iiH(t1 )...H(tn )dt1 ...dtn + ......) 0
where hh . ii denotes the cumulants of ζ. For a Gaussian stochastic process, cumulants of order higher than 2 vanish and hhζ(t1 )ζ(t2 )ii = hζ(t1 )ζ(t2 )i
(D.10)
From eqs. D.2, D.9 and D.10, for a Gaussian stochastic random process, we obtain hG(t)i = expo
2
α 2!
Z t Zt1
H(t1 )H(t2 )e−νc (t1 −t2 ) dt1 dt2 0
0
56
(D.11)
From Eq. D.6, we can easily get dhφ(t)i dhG(t)i = χ(0) dt dt
(D.12)
From Eq. D.11, dhG(t)i 2 = α H(t) dt
Zt
H(s)e−νc (s−t) ds hG(t)i
(D.13)
0
Thus, from Eqs. D.6, D.12 and D.13, the ensemble average of the interaction perturbation obeys
dhφi 2 = α H(t) dt
Zt
H(s)e−νc (s−t) ds hφi
(D.14)
0
Differentiating Eq. D.3, dhχi dhφi = Lhχi + eLt dt dt
(D.15)
From Eq. D.15, the exact deterministic evolution equation for the ensemble average perturbation χ(t) can be written as ¢ dhχi ¡ = L + α2 F Θ(t) hχi dt where
(D.16)
Zt eLs F e−Ls e−νc s ds
Θ(t) =
(D.17)
0
Equivalent White Noise Approximation Lieuwen (2003) mentions that in practical thermoacoustic systems, there are high degree of freedom processes with short correlation times compared to the acoustic period. For short autocorrelation times tc << 1, Eq. D.16 can be approximated. The integrand in Eq. D.17 can be expanded as
Ls
e Fe
−Ls −νc s
e
=e
−νc s
¶ µ s2 L + s[L, F ] + [L, [L, F ]] + ... 2!
57
(D.18)
where [L, F ] = LF − F L is the commutator. Equation D.16 can then be written in the exact form as · µ ¶¸ dhχi I2 2 = L + α F I0 F + I1 [L, F ] + [L, [L, F ] + ...] hχi dt 2! where In =
Rt 0
(D.19)
sn e−νc s ds. Since In+1 = O(1/νcn+1 ), for tc << 1, the ensemble mean
accurately evolves retaining only the first power of 1/νc in Eq. D.19. ¸ · dhχi α2 νc t 2 = L + (1 − e )F hχi dt νc
(D.20)
In the limit, Eq. D.20 becomes · ¸ α2 2 dhχi = L + F hχi dt νc
(D.21)
Equation D.21 is identical to the ensemble mean equations obtained for a white noise process with a variance of 2α2 /ν in the physically relevant Stratonovich limit (Arnold, 1992).
58
APPENDIX E
MAXIMUM TRANSIENT ENERGY GROWTH AND OPTIMAL INITIAL CONDITIONS FOR STOCHASTIC SYSTEMS
Farrell and Ioannou (2002b) derives the equations to find the maximum energy growth and optimal initial conditions for stochastic systems. Let ϕ(t, 0) be the fundamental matrix associated with each realization of the operator L + αζ(t)F . At time t, the perturbation square amplitude for each realization of the fluctuation is χ† (t)χ(t) = χ† (0)ϕ† (t, 0)ϕ(t, 0)χ(t)
(E.1)
As seen in Appendix C, the maximum of the largest eigenvalues of the ensemble mean of the matrix S(t) = ϕ† (t, 0)ϕ(t, 0)
(E.2)
gives the maximum transient energy growth, the corresponding eigenvector gives the initial condition for maximum energy growth, and the time at which this maximum of the largest eigenvalues occurs gives the time at which the energy growth maximizes. The other eigenvectors complete the set of mutually orthogonal initial conditions ordered according to their growth at time t. The matrix S(t) can be determined by integrating dS = (L + αζ(t)F )† S + S(L + αζ(t)F ) dt
(E.3)
The ensemble mean evolution equation corresponding to Eq. E.3 needs to be obtained. First, Eq. E.3 is expressed as a vector equation using tensor products. This is achieved by consecutively stacking the columns of the n × n matrix S(t) to for a n2 column
vector s(t). In tensor form, Eq. E.3 can be written as ds = (l + αζ(t)f )s dt
(E.4)
where l = I ⊗ L† + LT ⊗ I f = I ⊗ F† + FT ⊗ I where, T denotes transpose, † denotes Hermitian conjugate and I is the Identity matrix. The ensemble evolution equation for s can be obtained as (see Appendix D) ¢ dhsi ¡ = l + α2 f θ(t) hsi dt
(E.5)
Zt elx f e−lx e−νc x dx
θ(t) =
(E.6)
0 Tt
Since elt = eA
†
⊗ eA t , repeated application of tensor product properties gives θ(t) = I ⊗ Ξ(t)† + Ξ(t)T ⊗ I
where
(E.7)
Zt e−Lx F eLx e−νc x dx
Ξ(t) =
(E.8)
0
The ensemble average equation for the vector s can be written as dhsi = (I ⊗ L† + LT ⊗ I + α2 (I ⊗ F † + F T ⊗ I)(I ⊗ Ξ(t)† + Ξ(t)T ⊗ I))hsi dt (E.9) In the matrix notation, the deterministic ensemble mean evolution equation can be written as dhSi = (L+α2 Ξ(t)F )† hSi+hSi(L+α2 Ξ(t)F )+α2 (Ξ(t)† hSiF +F † hSiΞ(t)) (E.10) dt
60
Equivalent White Noise Approximation For short autocorrelation times tc << 1, Eq. E.10 can be approximated, to obtain the equivalent white noise approximation, as dhSi = dt
µ
α2 L+ F νc
¶†
µ
α2 hSi + hSi L + F νc
¶ +
2α2 † F hSiF νc
(E.11)
The ensemble mean evolution thus obtained can be solved with initial condition hS(0)i = I to obtain hS(t)i. The eigen-analysis of hS(t)i will give the maximum energy growth factor, corresponding optimal initial condition and the time to attain this maximum growth.
61
REFERENCES 1. Ananthkrishnan, N., S. Deo, and F. E. C. Culick (2005). Reduced-order modeling and dynamics of nonlinear acoustic waves in a combustion chamber. Combustion Science and Technology, 177, 221–247. 2. Annaswamy, A. M., M. Fleifel, J. P. Hathout, and A. F. Ghoniem (1997). Impact of linear coupling on design of active controllers for the thermoacoustic instability. Combustion Science and Technology, 128, 131–180. 3. Arnold, L., Stochastic Differential Equations: Theory and Applications. Kreiger, 1992. 4. Baggett, J., T. Driscoll, and L. Trefethen (1995). Transition in shear flows:non-linear normality versus non-normal linearity. Physics of Fluids, 7, 833–838. 5. Balasubramanian, K. and R. I. Sujith (2007a). The role of nonnormality in combustion interactions in diffusion flames. AIAA 2007-0567, Reno, Nevada, USA. 6. Balasubramanian, K. and R. I. Sujith (2007b). Thermoacoustic instability in a Rijke tube: nonnormality and nonlinearity. AIAA 2007-3428, Rome, Italy. 7. Balasubramanian, K. and R. I. Sujith (2008a). Non-normality and nonlinearity in combustion-acoustic interaction in diffusion flames. Journal of Fluid Mechanics, 594, 29–57. 8. Balasubramanian, K. and R. I. Sujith (2008b). Thermoacoustic instability in a rijke tube: nonnormality and nonlinearity. Physics of Fluids, 20, 044103. 9. Burnley, V. S. (1996). Nonlinear Combustion Instabilities and Stochastic Sources. Ph.D. thesis, California Institute of Technology. 10. Chung, K., Elementary Probability Theory with Stochastic Processes. Springer, 1975. 11. Clavin, P., J. S. Kim, and F. A. Williams (1994). Turbulence-induced noise effects on high-frequency combustion instabilities. Combustion Science and Technology, 96, 61–84. 12. Crocco, L. (1952). Aspects of combustion stability in liquid propellant rocket motors. Journal of American Rocket Society, 21, 163–178. 13. Culick, F. E. C. (1976). Nonlinear behavior of acoustic waves in combustion chambers, parts i and ii. Acta Astronautica, 3, 715–756. 14. Culick, F. E. C. (1997). A note on ordering perturbations and insignificance of linear coupling in combustion instabilities. Combustion Science and Technology, 126, 359– 379. 15. Culick, F. E. C., V. Burnley, and G. Swenson (1995). Pulsed instabilities in solid propellant rockets. Journal of Propulsion and Power, 11, 657–665. 62
16. Dowling, A. P. (1995). The calculation of thermoacoustic oscillations. Journal of Sound and Vibration, 180, 557–581. 17. Dowling, A. P. and S. R. Stow (2003). Acoustic analysis of gas turbine combustors. Journal of Propulsion and Power, 19, 751–764. 18. Dowling, P. and J. E. Ffowcs-Williams, Sound and Sources of Sound. Ellis Horwood, West Sussex, U.K., 1983. 19. Farrell, B. and P. Ioannou (1996a). Generalized stability theory. part 1: Autonomous operators. Journal of Atmospheric Sciences, 53, 2025–2040. 20. Farrell, B. and P. Ioannou (1996b). Generalized stability theory. part 2: Nonautonomous operators. Journal of Atmospheric Sciences, 53, 2041–2053. 21. Farrell, B. F. and P. J. Ioannou (2002a). Perturbation growth and structure in uncertain flows. part 1. Journal of Atmospheric Sciences, 59, 2629–2646. 22. Farrell, B. F. and P. J. Ioannou (2002b). Perturbation growth and structure in uncertain flows. part 2. Journal of Atmospheric Sciences, 59, 2647–2664. 23. Heckl, M. A. (1990). Nonlinear acoustic effects in the rijke tube. Acoustica, 72, 63–71. 24. Heckl, M. A. and M. S. Howe (2007). Stability analysis of a rijke tube with a green’s function approach. Journal of Sound and Vibration, 305, 672–688. 25. Henrici, P. (1969). Bounds for iterates, inverses, spectral variation and field of values of nonnormal matricies. Numerical Mathematics, 4, 24–40. 26. Hwang, J. H. and F. Ma (1993). On the approximate solution of nonclassically damped linear systems. ASME Journal of Applied Physics, 60, 695–701. 27. Kopitz, J. and W. Prolifke (2005). CFD based analysis of thermoacoustic instabilities by determination of open-loop-gain. The International Institute of Acoustics and Vibration (IIAV), Paper no. 389, Lisbon. 28. Kubo, R. (1962). Generalized cumulant expansion method. Journal of the Physical Society of Japan, 17, 1100–1120. 29. Kulkarni, R. R., K. Balasubramanian, and R. I. Sujith (2008). The role of nonnormality in active control of combustion instability. 30. Lieuwen, T. C. (2003). Statistical characteristics of pressure oscillations in a premixed combustor. Journal of Sound and Vibration, 260, 3–17. 31. Lieuwen, T. C. and A. Banaszuk (2005). Background noise effects on combustor stability. Journal of Propulsion and Power, 21, 25–31. 32. Matveev, K. I. (2003). Thermoacoustic instabilities in the Rijke tube: experiments and modeling. Ph.D. thesis, California Institute of Technology. 33. Matveev, K. I. and F. E. C. Culick (2003). A model for combustion instability involving vortex shedding. Combustion Science and Technology, 175, 1059–1083. 63
34. Meirovitch, L., Analytical Methods in Vibrations. Macmillan Publishing Co., Inc., New York, 1967. 35. Nicoud, F., L. Benoit, C. Sensiau, and T. Poinsot (2007). Acoustic modes in combustors with complex impedances and multidimensional active flames. AIAA Journal, 45, 426–441. 36. Papoulis, A. and S. U. Pillai, Probability, Random Variables and Stochastic Processes. McGraw-Hill, 2002. 37. Park, I. W., J. S. Kim, and F. Ma (1994). Characteristics of modal coupling in nonclassically damped systems under harmonic excitation. ASME Journal of Applied Mechanics, 61, 77–83. 38. Rayleigh, J. W. S. (1878). The explanation of certain acoustic phenomena. Nature, 18, 319–321. 39. Reddy, S. and D. Henningson (1993). Energy growth in viscous channel flows. Journal of Fluid Mechanics, 252, 209–238. 40. Schmid, P. J. (2007). Nonmodal stability theory. Annual Review of Fluid Mechanics, 39, 129–162. 41. Schmid, P. J. and D. S. Henningson, Stability and Transition in Shear Flows. SpringerVerlag, New York, Berlin, Heidelberg, 2001. 42. Trefethen, L. N. (1999). Computation of pseudospectra. Acta Numerica, 8, 247–295. 43. Trefethen, L. N. and M. Embree, Spectra and Pseudospectra. Princeton University Press, Princeton and Oxford, 2005. 44. Trefethen, L. N., A. E. Trefethen, S. C. Reddy, and T. A. Driscoll (1993). Hydrodynamic stability without eigenvalues. Science, 261, 578–584. 45. van Kampen, N. G., Stochastic Processes in Physics and Chemistry. Elsevier, 1992. 46. Wicker, J. M., W. D. Greene, S. I. Kim, and V. Yang (1996). Triggering of longitudinal combustion instabilities in rocket motors: nonlinear combustion response. Journal of Propulsion and Power, 12, 1148–1158. 47. Wright, T. G. (2002). Eigtool. Oxford University. URL http://www.comlab. ox.ac.uk/pseudospectra/eigtool/. 48. Yoon, H. G., J. Peddieson, and K. R. Purdy (1998). Mathematical modeling of a generalized rijke tube. International Journal of Engineering Science, 36, 1235–1264. 49. Yoon, H. G., J. Peddieson, and K. R. Purdy (2001). Nonlinear response of a generalized rijke tube. International Journal of Engineering Science, 39, 1707–1723. 50. Zinn, B. T. and T. C. Lieuwen, Combustion instabilities: Basic Concepts, Chapter 1, Combustion instabilities in gas turbine engines: operational experience, fundamental mechanisms, and modeling. 210 Progress in Aeronautics and Astronautics, AIAA, Inc., edited by T. C. Lieuwen and V. Yang, 2005. 64
LIST OF PAPERS BASED ON THESIS 1. Kedia, K. S., Nagaraja, S. B. and Sujith, R. I., Impact of linear coupling on thermoacoustic instabilities in a Rijke tube. Combustion Science and Technology. (in press). 2. Nagaraja, S. B.,Kedia, K. S. and Sujith, R. I., Characterizing energy growth in thermoacoustic instabilities: singularvalues or eigenvalues?. International Combustion Symposium, Montreal, Canada, August 2008.
65
