The calculation of thermoacoustic oscillations

Journal of Sound and Vibration (1995) 180(4), 557–581

THE CALCULATION OF THERMOACOUSTIC OSCILLATIONS A. P. D Department of Engineering, University of Cambridge, Trumpington Street, Cambridge CB2 1PZ, England (Received 1 September 1992, and in final form 23 November 1993) Thermoacoustic oscillations occur in a wide variety of practical applications in which heat is supplied to an acoustic resonator. A simple geometry is investigated systematically to determine the importance of various flow parameters on the frequency of the oscillations. Detailed consideration of elementary examples shows that the form of the coupling between the heat input and the unsteady flow has a crucial effect on the frequency of oscillation. The same elementary examples are used to compare how well (if at all) different calculation methods in the literature account for this influence. A mean flow and a distributed region of heat input significantly complicate the analysis of thermoacoustic oscillations and are often neglected. Model problems are used to illustrate that mean flow effects can become significant even at modest inlet Mach numbers, and to indicate circumstances under which a distributed heat input may be treated as concentrated.

1. INTRODUCTION

There is a possibility of thermoacoustic oscillations whenever combustion or heat exchange take place within an acoustic resonator. They occur because unsteady heating generates sound waves producing pressure and velocity fluctuations. However, within a resonator, these in turn perturb the rate of heat input. Instability is possible if the phase relationship is suitable because, while the acoustic waves perturb the heat input, the unsteady heat input generates yet more sound. Rayleigh [1] gave a clear physical description of this phenomenon. Acoustic waves gain energy when the unsteady rate of heat input is in phase with the pressure perturbations. If this energy gain exceeds that lost on reflection from the boundaries of the resonator, linear acoustic waves grow in strength and the system is unstable. In practice, the amplitudes of the pressure waves are limited by non-linear effects. Nevertheless, they can become so intense that structural damage is done. Many practical devices are susceptible to destructive thermoacoustic oscillations. They occur in rockets, ramjets, gas boilers and aeroengines. Emphasis has been placed on investigating the frequency of the oscillations and the mean heat input level for the onset of the instability. Two excellent reviews of the literature have been given by Candel and Poinsot [2] and Culick [3]. One particular thermoacoustic oscillation, involving the propagation of low frequency longitudinal pressure waves, occurs in the afterburners of jet aeroengines. It is commonly called ‘‘reheat buzz’’. An experimental and theoretical investigation of this problem has involved a series of tests on a laboratory-scale rig, designed to model the essential features of an afterburner [4]. The rig geometry involves a flame burning in a duct in the wake of a bluff body. There is a mean flow and the heat release is both distributed and unsteady. It was found to be necessary to include a proper description of all these effects in the 557 0022–460X/95/090557 + 25 $08.00/0

7 1995 Academic Press Limited

558

. . 

development of a theory to explain the experimental results [5]. This is at variance with views and alternative approaches reported in the literature, where one or more of these effects are neglected. The aim of this paper is to investigate systematically, for a very simple geometry, the importance of various flow effects on the frequency of the thermoacoustic oscillations. This is followed by a discussion on how well these effects are described by different calculation methods. Section 2 is concerned with the effect of the form of the coupling between the heat input and the unsteady flow on the frequency of the combustion oscillation. In work on the ‘‘buzz’’ oscillation [4, 5], this coupling was found to have a crucial effect on the frequency. However, the complicated details of the experimental geometry and flow conditions prevented any general conclusions from being drawn from this observation. A very simply geometry is considered in section 2 to avoid this difficulty. One-dimensional disturbances are investigated in a duct with no mean flow and one closed and one open end. Unsteady heat is added at a single plane, across which the mean temperature changes from T1 to T2 . Two particular forms for the relationship between the rate of heat input and the flow are investigated in detail. Case I is that of no unsteady rate of heat input. This represents a limiting case when the reaction time of the heating process is much longer than the period of oscillation. Then, even though the fluid might enter the heating zone unsteadily, the rate of heat input would be unable to adjust sufficiently rapidly and would remain constant. In Case II the instantaneous rate of heat input is directly proportional to the instantaneous mass flow rate into the heating zone: i.e., the heat input per unit mass is constant. Such a heat input would be appropriate as a limiting case when the reaction time of the heating process is so much faster than the period of oscillation that the rate of heat input responds quasi-steadily to changes of flow rate into the heating zone. Both forms of heat input are, of course, idealizations. However, it is worth noting that the heat released by a non-premixed flame could approximate to Case I, when fuel is forced into the combustion region at a constant rate and the combustion is directly proportional to the fuel injection. Observations [5] of the rate of combustion for a premixed flame burning in a duct are of the form of Case II at low frequencies, but with an additional time delay. The frequencies of the thermoacoustic oscillations are calculated in section 2.1 for these two forms of the unsteady heat release rate. The relationship between the unsteady heat input and the flow is found to crucially affect the frequency of oscillation. Indeed, by a mean temperature ratio of six, there is nearly a 60% difference between the lowest frequency predicted for Case II (when the instantaneous rate of heat input is proportional to the instantaneous mass flow rate) and for Case I (when there is no unsteady rate of heat input). Some calculation methods commonly used in the literature take no account of the form of the coupling between the unsteady heat input and the flow. In sections 2.2 and 2.3 respectively, it is shown that the calculation of the duct acoustic frequency and the Green function technique recommended by Hedge et al. [6] are appropriate only when the rate of heat input is not altered by the unsteady flow. The linearized Galerkin method [3] acknowledges that the form of the coupling between the unsteady heat input and the flow affects the frequency, but treats the shift in frequency, from that for no unsteady heat release rate, as linear. This method is tested in section 2.4, by applying it to Case II, with the instantaneous rate of heat input proportional to the instantaneous mass flow rate. The predicted frequency and the true frequencies agree exactly for T2 = T1 , but diverge significantly for larger temperature ratios. By a temperature ratio of six, the Galerkin method predicts a frequency which is less that half the exact value.

 

559

The presence of a mean flow significantly complicates the analysis of a thermoacoustic oscillation. An entropy wave, or convected hot spot, becomes coupled to the acoustic waves and a simultaneous solution of equations of conservation of mass, momentum and energy is required. Many authors choose to neglect the mean flow rather than deal with this complication, their justification being that the inlet Mach number is low. The errors introduced by such an approximation are investigated in section 3. It is striking how significant Mach number effects can be. For example, one might be tempted to think that an inlet Mach number of 0·15 is sufficiently small for the mean flow to be neglected. But at this Mach number, the frequency of the thermoacoustic oscillation, for Case II, T01 /T02 = 6, is reduced to half its no flow value. Most practical ways of adding heat also exert a drag force on the flow. In work on the ‘‘buzz’’ oscillation, the drag force exerted by the bluff-body flame holder was found not to have a significant effect [7]. A systematic investigation into the effects of drag is given in section 4. The drag force exerted by a grid or flame holder with a blockage ratio of 25% or less is found to have a negligible effect on the frequency of thermoacoustic oscillations  1 E 0·15. A blockage ratio of 50% for the range of inlet Mach numbers considered, 0 E M alters the frequency for inlet Mach numbers greater than 0·10. When the heat release is distributed over an axial extent d rather than concentrated at a single plane, the strengths of the acoustic and entropy waves vary continuously with position. It is then easier to determine the unsteady flow directly by integrating the equations of motion. That is done in section 5 to investigate the effects of an extended region of heat input. It is sometimes said that distributed heat input can be considered as concentrated provided that its axial extent is very much smaller than the acoustic wavelength: i.e., provided that vd/c¯ W 1, where v is the frequency and c¯ is the speed of sound. However, this is an oversimplification when there is a mean flow. Then there is also an entropy wave, with the much shorter length scale u¯ /v, where u¯ is the mean flow speed. The entropy wave is affected by even a modest spatial distribution of the heat input, and the frequency of oscillation is altered significantly when this entropy wave is coupled to the acoustic waves. For example, results in section 5 show that for an inlet Mach number of 0·10 and a stagnation temperature ratio of nine, distributing the heat input over a length as short as 5% of the duct can lead to a 25% change in the lowest frequency of oscillation. To guarantee, in general, that heat input over a region of length d can be treated as concentrated requires not only vd/c¯ W 1, but also vd/u¯ W 1. The second inequality is highly restrictive and circumstances under which it can be relaxed are discussed in section 5.

2. THE EFFECT OF UNSTEADY HEAT INPUT

2.1.      In a region with heat input, the density r varies through changes in both the pressure p and specific entropy s. The chain rule of differentiation shows that

b

Dr 1 Dp 1r Ds = + . Dt c 2 Dt 1s p Dt

(2.1)

When viscous and heat conduction effects are neglected, rT Ds/Dt = q(x, t), where q is the heat input/unit volume and T is the absolute temperature. Moreover for a perfect gas, 1r/1s= p = −r/cp = −rT(g − 1)/c 2, where c is the speed of sound, cp is the specific heat at

. . 

560

constant pressure and g is the ratio of specific heats. Substitution into equation (2.1) leads to

0

1

Dr 1 Dp = − (g − 1)q . Dt c 2 Dt

(2.2)

Equation (2.2) may be applied to a combusting gas, provided that the reactants and products behave as perfect gases, and there is no molecular weight change during the chemical reaction [8]. For linear perturbations in a region where there is no mean heat input or mean velocity, equation (2.2) can be linearized:

0

1

Dr 1 1p' = − (g − 1)q' . Dt c¯ 2 1t

(2.3)

The overbar denotes a mean value and the prime a perturbation. A combination of equation (2.3) with the linearized equations of mass and momentum conservation, Dr/Dt + r¯ 9·u' = 0 and 1u'/1t + 9p'/r = 0, leads to

0 1

1 (g − 1) 1q' 1 1 2p' . − r¯ 9· 9p' = r¯ 1t c¯ 2 c¯ 2 1t 2

(2.4)

This inhomogeneous wave equation describes the pressure perturbation generated by the unsteady heat input q'(x, t). The solution of this wave equation can be readily calculated for one-dimensional disturbances in the geometry illustrated in Figure 1. A duct of length l has one closed and one open end. The heat input is harmonic of frequency v, concentrated at the plane x = b, q'(x, t) = Q'(t)d(x − b), where Q'(t) = Q eivt and d denotes the Dirac d-function. In x Q b, the mean temperature is T1 , the density r¯ 1 and the speed of sound c¯1 . The solution of equation (2.4), that satisfies the rigid end boundary condition at x = 0, is p'(x, t) = A eivt(e−ivx/c¯1 + eivx/c¯1 ), (2.5) where the complex constant A has yet to be determined. Similarly, in x q b, the mean fluid properties have values T2 , r¯ 2 and c¯2 , and the solution of the wave equation can be written in the form p'(x, t) = B eivt(e−iv(x − l)/c¯2 − eiv(x − l)/c¯2 ), (2.6) after application of an open end boundary condition, p'(l, t) = 0. Integration of equation (2.4), with q'(x, t) = Q'(t)d(x − b), across the region x = b shows that +

b [p']xx = = b− = 0

and

$ % 1 1p' r¯ 1x

x = b+

x = b−

(g − 1) dQ' , r¯ 1 c¯12 dt

=−

since for a perfect gas, r¯ 1 c¯12 = gp¯ = r¯ 2 c¯22. Equation (2.7b) is equivalent to (g − 1) b+ [u']xx = Q'(t), = b− = r¯ 1 c¯12

Figure 1. One-dimensional disturbances in a duct.

(2.7a, b)

(2.8)

 

561

a relationship between the volumetric expansion and the rate of heat input that has been frequently used in the literature. When Q'(t) is specified, the constants A and B may be calculated by substitution for p'(x, t) from equations (2.5) and (2.6) into the jump conditions (2.7). The pressure perturbation within the duct is then determined. However, thermoacoustic instabilities occur when the heat input and the flow perturbations are coupled. As well as the unsteady heating generating sound as described by equation (2.4), the flow perturbations in turn affect the rate of heat input. The details of this feedback depend on the precise form of the heating but, in all cases, lead to a relationship between Q and the complex constants A and B. This relationship, together with the jump conditions (2.7), constitute three homogeneous equations for the three unknowns A, B and Q , which can be written symbolically in the form A X B = 0, (2.9) Q

23

where X is a 3 × 3 matrix. Such an equation will have a non-zero solution only if the determinant of X vanishes, and this condition determines the frequency of the thermoacoustic oscillation. In this paper, two particular forms for the unsteady rate of heat input will be considered as illustrative examples. The first case is that of no unsteady rate of heat input: Case I,

Q'(t) = 0,

no unsteady rate of heat input.

(2.10)

Such a heat input would typically occur when the reaction time of the heating process is much longer than the period of oscillation. Then, even though the fluid might enter the heating zone unsteadily, the rate of heat input would be unable to adjust sufficiently rapidly and would remain constant. In the second case, the instantaneous rate of heat input is taken to be directly proportional to the instantaneous mass flow rate into the heating zone, the heat input per unit mass being that required to raise the mean temperature from T1 to T2 : i.e., cp (T2 − T1 ). Since the heat input per unit mass is constant, this case is referred to as that of no unsteady heat input per unit mass. The rate of heat input is then proportional to the mass flow rate into the heating zone: Case II,

Q'(t) = cp (T2 − T1 )r¯ 1 u'1 , no unsteady heat input per unit mass,

(2.11)

where u'1 (t) denotes u'(b−, t). Such a heat input would be appropriate when the reaction time of the heating process is so much faster than the period of oscillation that the rate of heat input responds quasi-steadily to changes of flow rate into the heating zone. Of course, the two forms for the rate of heat input in (2.10) and (2.11) are idealizations. However, it is worth noting that the heat released by a non-premixed flame could approximate to Case I (expression (2.10)), if the fuel was forced into the combustion region at a constant rate and the rate of combustion was directly proportional to the rate of fuel injection. From a series of experiments on a confined premixed flame, Bloxsidge et al. [5] postulated a rate of heat input which, for harmonic disturbances, can be expressed in the form Q'(t) = kcp (T2 − T1 )r¯ 1 u'1 (t − t); k and the time delay t depend on the frequency. Over a range of measured frequencies, k varied between 0·5 and 1·0. There are, therefore, also circumstances under which the Case II form of the heat input in expression (2.11) is achievable in practice. The aim of this section is not to obtain detailed predictions for a particular geometry of flame or heated grid, but rather to demonstrate that the form of the unsteady heat input can have a significant effect on the frequency of a thermoacoustic

562

. . 

oscillation. This leads on to an investigation of how this effect is accounted for (if at all) in some of the methods of solution in the literature. For the rate of heat input in Case I, described by equation (2.10), substitution for p'(x, t) from equations (2.5) and (2.6) into the jump conditions (2.7) leads to an equation for v. After some algebra the frequency of oscillation is found to be given implicitly by tan (a) tan (b) = r¯ 1 c¯1 /r¯ 2 c¯2

for Case I,

(2.12)

where a = vb/c¯1 and b = v(l − b)/c¯2 . It is a straightforward matter to solve equation (2.12) numerically. The roots are real and the lowest frequency solution is shown in Figure 2 for various values of T2 /T1 and for the particular geometry b = l/2. When the rate of heat input in Case II described by equation (2.11) is used, a different equation for the frequency is obtained. It is evident from equation (2.5) that u'1 (t) = A eivt(e−ivb/c¯1 − eivb/c¯1 )/r¯ 1 c¯1 .

(2.13)

When this is used in equation (2.11), substitution into equation (2.7) leads to tan (a) tan (b) = c¯1 /c¯2

for Case II.

(2.14)

Again the roots are real, and the lowest frequency solution is plotted in Figure 2 for comparison with the solution of equation (2.12). When the temperature is uniform along the duct (T2 = T1 ), the predicted frequencies just correspond to that for which the duct length is a quarter-wavelength. For other temperature ratios, the relationship between the unsteady heat input and the flow crucially affects the frequency of oscillation. Indeed by a temperature ratio of six, there is nearly a 60% difference between the frequency predicted for no unsteady heat input per unit mass and no unsteady heat input per unit volume. Of course, it is well known that the unsteady heat input affects the frequency. Rayleigh [1] commented that ‘‘the pitch is raised if heat be communicated to the air a quarter period before the phase of greatest condensation: and the pitch is lowered if the heat be communicated a quarter period after the phase of greatest condensation’’. In the geometry considered here, the velocity fluctuation u'1 lags the pressure perturbation and so the rate of heat input in equation (2.11) lags the greatest condensation and leads to a lower frequency than that for Q' = 0. The numerical calculations show that this shift in frequency

Figure 2. The lowest frequency of oscillation as a function of temperature ratio T  2 /T  1 for b = l/2. ——, Case I, no unsteady rate of heat input (roots of equation (2.12)); – – – – , Case II, no unsteady heat input per unit mass (roots of equation (2.14)).

 

563

can be significant, and it is appropriate to discuss how this effect is accounted for in some of the methods of solution in the literature. 2.2.      Many authors just assume that thermoacoustic oscillations occur at the ‘‘acoustic frequency’’ of a duct. By this, they invariably mean the frequency at which a solution of the homogeneous wave equation

0 1

1 1 1 2p' − r¯ 9· 9p' = 0 r¯ c¯ 2 1t 2

(2.15)

satisfies the boundary conditions. This frequency takes no account of the relationship between the unsteady heat release rate and the flow. For example, for the particular geometry in Figure 1 it leads to tan (a) tan (b) = r¯ 1 c¯1 /r¯ 2 c¯2 ,

(2.16)

for any form of the unsteady heat release rate. A comparison of equations (2.4) and (2.15) shows why this is so. Equation (2.4), which was derived from the equations of motion, reduces to the homogeneous form (2.15) only when there is no unsteady rate of heat input per unit volume. When the rate of heat input is unsteady, the pressure perturbations satisfy the inhomogeneous wave equation (2.4) and the form of the heat input must be considered in a calculation of the frequency of oscillations in the duct. 2.3.     Hedge et al. [9] used a Green function technique to determine the one-dimensional pressure perturbation generated by unsteady combustion within a duct. They essentially solved the inhomogeneous wave equation (2.4) for q'(x, t) = qˆ (x) eivt, by introducing a Green function G(x=x0 ) which satisfies r¯ c¯ 2

0 1

d 1 dG + v 2G = −cp (g − 1)d(x − x0 ) dx r¯ dx

(2.17)

and the boundary conditions. The solution of equation (2.4) is then given by p'(x, t) = pˆ (x) eivt, where pˆ (x) =

iv cp

g

G(x=x0 )qˆ (x0 ) dx0 .

(2.18)

In particular, when the heat input is concentrated on the plane x = b, qˆ (x0 ) = Q d(x0 − b), and pˆ (0) = ivQ G(0=b)/cp .

(2.19)

A general form for the Green function was given by Hedge et al. [9, equation 14]. For the geometry illustrated in Figure 1, their expression simplifies to sin (b) c¯1 c¯2 . v T1 c¯2 sin (a) sin (b) − T2 c¯1 cos (a) cos (b)

G(0=b) = −

(2.20)

Hedge et al. found the resonance frequencies of the duct by minimizing the magnitude of the denominator of the Green function. For the particular form (2.20), the denominator vanishes at frequencies that satisfy tan (a) tan (b) = T2 c¯1 /T1 c¯2 = r¯ 1 c¯1 /r¯ 2 c¯2 ,

(2.21)

564

. . 

a condition identical to that in equation (2.12). If the unsteady heat input is a specified broadband source, equation (2.19) shows that large pressure perturbations are generated at this frequency. However, it is misleading to interpret this as the frequency of a thermoacoustic instability, as suggested in reference [6]. A thermoacoustic instability involves coupling between the heat input and the flow. The unsteady heat input occurs as a response to fluctuations in flow. This can be expressed in the form Q = Z(v)pˆ (0), (2.22) for some function Z(v). Substitution for Q in equation (2.19) leads to pˆ (0)(1 − ivZ(v)G(0=b)/cp ) = 0.

(2.23)

The frequency of the thermoacoustic oscillation is the frequency at which self-sustaining oscillations can occur. It is clear from equation (2.23) that pˆ (0) can be non-zero only if {1/G(0=b)} − ivZ(v)/cp = 0.

(2.24)

This reduces to 1/G(0=b) = 0 only if Z(v) = 0. This explains why the zeros of 1/G(0=b) were found in equation (2.21) to be identical to frequencies calculated for the particular case of no unsteady rate of heat input. When the heat input responds to the flow, Z(v) is non-zero and the frequency is shifted. As an example, consider Case II with no unsteady heat input per unit mass. The mode shape in x Q b, as described by equations (2.5) and (2.13), shows that u'1 (t) = −i sin (a)pˆ (0) eivt/r¯ 1 c¯1 .

(2.25)

Hence the unsteady heat input described by equation (2.11) can be written in the form Q = −icp (T2 − T1 ) sin (a)pˆ (0)/c¯1 , leading to Z(v) = −icp (T2 − T1 ) sin (a)/c¯1 .

(2.26)

Substitution of the Green function and Z(v), from equations (2.20) and (2.26), respectively, into equation (2.24) leads, after some algebra, to the same equation for v as that given in equation (2.14). In summary, a Green function can be a useful tool in the investigation of thermoacoustic oscillations. However, it is not sufficient simply to inspect the Green function alone. The form of the coupling between the heat input and the flow must also be considered. In the notation of this section, this coupling is described by the function Z(v). The frequency of the thermoacoustic oscillation is then a zero of {1/G(0=b)} − ivZ(v)/cp . The form of the relationship between the unsteady heat input and the flow affects the frequency. 2.4.    Culick and his colleagues used a Galerkin method to solve the inhomogeneous wave equation (2.4). This method acknowledges that the form of the coupling between the unsteady heat input and the flow affects the frequency, but linearizes the change in frequency and mode shape. Reference [3] clearly describes the method and gives a review of earlier work. For the one-dimensional problem illustrated in Figure 1, the method involves expanding the pressure perturbation as a Galerkin series: a

p'(x, t) = s hm (t)cm (x). m=1

(2.27)

  The functions cm (x) are eigensolutions of the homogeneous wave equation,

0 1

vm2 d 1 dcm c + r¯ = 0, c¯ 2 m dx r¯ dx

m = 1, 2, . . . ,

565

(2.28)

that satisfy the boundary conditions dcm /dx = 0 at x = 0 and cm = 0 at x = l. (2.29) It is a straightforward matter to show from this definition that the functions cm are orthogonal,

g

l

cm (x)cn (x) dx = 0

for m $ n.

(2.30)

0

For the geometry in Figure 1, the mean temperature is uniform in the two regions x Q b and x q b. The homogeneous wave equation in (2.28) then simplifies to d2cm vm2 d2cm vm2 in 0 E x Q b, + 2 cm = 0 in b Q x E l, (2.31a, b) 2 cm = 0 2 + dx c¯1 dx 2 c¯2 together with two jump conditions across x = b: +

b [cm ]xx = = b− = 0

and

$ % 1 dcm r¯ dx

x = b+

= 0.

(2.32)

x = b−

It is a matter of straightforward algebra to check that cm (x) =

6

cos (vm x/c¯1 ) in x Q b (cos am /sin bm ) sin (vm (l − x)/c¯2 ) in x q b

7

(2.33)

in a solution of equations (2.31) and (2.32) that satisfies the boundary conditions (2.29). vm is the mth root of tan am tan bm = (r¯ 1 c¯1 )/(r¯ 2 c¯2 ) and am = vm b/c¯1 and bm = vm (l − b)/c¯2 . Substitution for p'(x, t) from equation (2.27) into equation (2.4) leads to

0

a

s m=1

1

d2hm 1q' + vm2 hm cm (x) = (g − 1) . dt 2 1t

(2.34)

After multiplication by cn (x) and integration with respect to x, the orthogonality condition (2.30) shows that equation (2.34) simplifies to (g − 1) d2hn + vn2hn = En dt 2

g

l

0

1q' c (x) dx, 1t n

(2.35)

where En =

g

l

cn2 dx = 12 (b + (l − b) cos2 am /sin2bm ),

(2.36)

0

after substitution for cn from equation (2.33). No approximation has yet been made, but in the next stage in the analysis it is assumed that 1q'/1t is small in magnitude. It is then argued that since 1q'/1t is small, it need only be evaluated approximately. The acoustic approximations hn (t)cn (x) and (h˚ n (t)/r¯ vn2 ) dcn / dx for the pressure and velocity perturbation are then used when calculating 1q'/1t. If second derivatives of the amplitudes arise, they are replaced by the zeroth order approximation h¨ n (t) 2 −vn2hn (t). The errors induced by these approximations can be checked by applying the method to Case II to find the lowest frequency of combustion oscillation when there is no unsteady heat input per unit mass.

. . 

566

Then the rate of heat input has the form given by equation (2.11) and q'(x, t) = cp (T2 − T1 )r¯ 1 u'1 d(x − b). Rewriting the velocity perturbation in the acoustic form (h¨ 1 (t)/r¯ v12) dc1 /dx leads to q'(x, t) = cp (T2 − T1 )

h˚ 1 (t) dc1 − (b ). v12 dx

(2.37)

The right side of equation (2.35) therefore involves the second derivative, h¨ 1 (t). This is replaced by −v12h1 (t), following Culick’s rules, to give h¨ 1 + v12h1 = −h1

cp dc (g − 1)(T2 − T1 ) 1 (b−)c1 (b) E1 dx

= h1 cp v1 (g − 1)(T2 − T1 ) sin a1 cos a1 /E1 c¯1 ,

(2.38)

after substitution for c1 from equation (2.33). Equation (2.38) is a straightforward second order differential equation for h1 (t). The frequency of oscillation of h1 (t), v, can be found by substituting h1 (t) = C eivt into equation (2.38) to give v 2 = v12 − cp v1 (g − 1)(T2 − T1 ) sin a1 cos a1 /E1 c¯1 .

(2.39)

As in this method it is assumed that the difference v − v1 is small, Culick [3] therefore recommended that the square root in equations such as (2.39) be evaluated approximately by using the binomial theorem. Equation (2.39) then simplifies to v = v1 − cp (g − 1)(T2 − T1 ) sin a1 cos a1 /2E1 c¯1 .

(2.40)

This approximation to the lowest frequency of combustion oscillation is plotted in Figure 3 for comparison with the exact value. The linearized Galerkin method calculates the change in frequency due to the unsteady heat input, but treats the shift as small. The form of heat input in equation (2.11), i.e., no unsteady heat input per unit mass, vanishes identically for T2 = T1 and is infinitesimally small for small T2 − T1 . Hence, the frequency predicted from the linearized Galerkin method and the true frequency agree exactly for T2 = T1 . The two curves of frequency against temperature ratio also have the same slope at this point, but they diverge significantly for larger temperature ratios, as shown in Figure 3. By a temperature ratio of six, the Galerkin method predicts a frequency which is less than half the exact value.

Figure 3. The lowest frequency of oscillation for Case II with no unsteady heat input per unit mass for b = l/2. ——, Correct solution; ····, acoustic resonance frequency and the value predicted by investigating zeros of 1/G(0=b); – – – – , frequency predicted by linearized Galerkin theory; · · · · , frequency predicted by ray theory.

 

567

2.5.   The phase change of a high frequency sound wave, as it travels through a region of gradually varying sound speed c¯ (x), can be described by ray theory. According to ray theory, negligible sound energy is reflected as the ray propagates through such a region, and a sound wave of frequency v, propagating in the x-direction undergoes a phase change v f0l dx/c¯ (x) over a distance l [10]. Ray theory is appropriate for short wavelength waves travelling through a region in which the sound speed varies gradually and there is no unsteady heat input per unit volume. It is therefore entirely the wrong physical limit for low frequency thermoacoustic oscillations, where the wavelength is longer than, or comparable with, the length of the region of heat input and the oscillations are driven by the unsteady heating. Nevertheless, ray theory is used in industry to estimate the frequency of thermoacoustic oscillations and so it is worthwhile to compare its predictions with the exact frequency in some model calculations. According to the ray theory, the phase of a ray is changed by v(b/c¯1 + (l − b)/c¯2 ), after travelling along the duct in Figure 1. The pressure wave undergoes phase changes at 0 and p on reflection at the closed and open ends of the duct respectively. An estimate of the resonance frequency then follows from a requirement that the total phase change after propagation up and down the duct be an integer multiple of 2p: i.e., 2v(b/c1 + (l − b)/c¯2 ) + p = 2pn,

n integer.

(2.41)

This clearly gives frequencies which are independent of the form of the rate of heat release. The lowest frequency solution of equation (2.41) is plotted in Figure 3. It agrees with the true frequency of oscillation only for the trivial case with no mean temperature gradient, T2 = T1 . When the mean temperature varies along the duct over a length scale which is much shorter than the wavelength, ray theory is inadequate. This is because it does not take reflections from the region of heat input into account and so is unable to describe the flow perturbations correctly. 2.6.    Thermoacoustic oscillations have been investigated for the geometry illustrated in Figure 1. Results for the frequency of oscillation for two particular forms of the unsteady heat release rate are summarized in Figure 2. They show that the form of the relationship between the unsteady heat input and the flow crucially affects the frequency of oscillation. Many authors have neglected this effect or treated it as small. It has been said [3] that good approximations to the frequency of a thermoacoustic oscillation can often be obtained by calculating the acoustic frequencies of the geometry. That is equivalent to solving the case of no unsteady heat release rate per unit volume. The results in Figure 3 serve as a warning that this procedure is not universally appropriate. Hedge et al. [6] recommended the use of a Green function G to investigate one-dimensional perturbations in a duct. They then interpreted the roots of 1/G = 0 as the frequency of the thermoacoustic oscillation. In Figure 3 it is shown that this leads to results which can be significantly in error for the particular case of no unsteady heat input per unit mass. The analysis in section 2.3 explains why it is not sufficient simply to inspect the Green function alone. The frequency of a thermoacoustic oscillation is a zero of 1/G only if the rate of heat input per unit volume does not vary significantly in response to the flow perturbations. If it does, the form of the coupling between the heat input and the flow affects the frequency and must be considered. Culick and his colleagues have developed a Galerkin method to investigate thermoacoustic oscillations. This method acknowledges that the form of the coupling between the

568

. . 

unsteady heat input and the flow affects the frequency, but treats the shift in frequency, from that for no unsteady heat release rate per unit volume, as linear. The method therefore works exactly when there is no unsteady heat release rate per unit volume. In section 2.4 the method was applied to the case of no unsteady heat input per unit mass and once again the results are plotted in Figure 3. The predicted frequency and true frequencies agree exactly for T2 = T1 . The two curves of frequency against temperature ratio also have the same slope at this point, but they diverge significantly for larger temperature ratios. By a temperature ratio of six, the Galerkin method predicts a frequency which is less than half the exact value. In section 2.5 the frequency ‘‘predicted’’ by ray theory was compared with the true frequency. There is no theoretical basis for believing that ray theory can accurately describe low frequency thermoacoustic oscillations. Indeed, the plots in Figure 3 show that ray theory is unable to describe correctly the variation of the frequency of oscillation with the temperature ratio T2 /T1 .

3. THE EFFECT OF A MEAN FLOW

Most practical thermoacoustic oscillations involve a mean flow. However, often the Mach number of the oncoming flow is so small that it is tempting to neglect the mean velocity entirely. The errors introduced by such an approximation are investigated in this section. A mean flow has two main consequences. Trivially, it affects the speed of the propagation of the acoustic waves, which then travel downstream with speed c¯ + u¯ and upstream at c¯ − u¯ . In addition, a mean flow admits a second type of linear waves. These are convected with the mean velocity and involve convected entropy or vorticity. As a consequence, the equations of conservation of mass, momentum and energy become coupled. These effects may be illustrated by considering one-dimensional thermoacoustic oscillations in the geometry in Figure 4. The inlet is choked so that the mass flow rate is constant there: (r'/r¯ 1 ) + (u'/u¯1 ) = 0

at x = 0.

(3.1)

The boundary condition of a downstream open end can again [10] be taken as p' = 0

at x = l.

(3.2)

The heat input is concentrated at the fixed plane x = b, the rate of heat input per unit cross-sectional area being denoted by Q'(t). Consider disturbances with time dependence eivt. Upstream of the zone of heat input, there are acoustic waves propagating in both directions and the flow is isentropic. The

Figure 4. One-dimensional disturbances in a duct with a mean flow and mean rate of heat release.

 

569

perturbation in flow quantities can be written in the form p'(x, t) = eivt(A e−ivx/c¯1 (1 + M 1 ) + B eivx/c¯1 (1 − M 1 ) ),

(3.3a)

u'(x, t) = eivt (A e−ivx/c¯1 (1 + M 1 ) − B eivx/c¯1 (1 − M 1 ) )/r¯ 1 c¯1 ,

(3.3b)

r'(x, t) = p'(x, t)/c¯ , 2 1

cp T'(x, t) = p'(x, t)/r¯ 1 ,

(3.3c, d)

for 0 E x Q b, where A and B denote the strengths of the acoustic plane waves and M  1 = u¯1 /c¯1 . Downstream of the region of heat input, there might be a convected hot spot in addition to plane sound waves, and so p'(x, t) = eivt (C e−ivx/c¯2 (1 + M 2 ) + D eivx/c¯2 (1 − M 2 ) ), ivt

u'(x, t) = e (C e

−ivx/c¯2 (1 + M 2)

r'(x, t) = cp T'(x, t) =

−D e

ivx/c¯2 (1 − M 2)

(3.4a)

)/r¯ 2 c¯2 ,

(3.4b)

p'(x, t) Sr¯ 2 iv(t − x/u¯2 ) − e , c¯22 cp

(3.4c)

p'(x, t) Sc¯22 + eiv(t − x/u¯2 ), r¯ 2 (g − 1)cp

(3.4d)

for b Q x E l. C and D are the amplitudes of the acoustic waves, S is that of the entropy wave or convected hot spot and there are no vorticity waves in this one-dimensional example. M  2 = u¯2 /c¯2 . The constants A, B, C, D and S are related by conservation conditions across the zone of heat input. If the heat input is at a constant position x = b, then conservation of mass across this discontinuity shows that r¯ 1 u'1 + r'1 u¯1 = r¯ 2 u'2 + r'2 u¯2 ,

(3.5) −

+

where the suffices 1 and 2 denote the flow quantities at x = b and b Similarly, it follows from conservation of momentum across x = b that

respectively.

p'1 + r'1 u¯12 + 2r¯ 1 u¯1 u'1 = p'2 + r'2 u¯22 + 2r¯ 2 u¯2 u'2 .

(3.6)

The perturbation in the energy flux out of a control volume around x = b exceeds the energy flux into it by Q'(t): i.e., cp T01 (r¯ 1 u'1 + r'1 u¯1 ) + r¯ 1 u¯1 (cp T'1 + u¯1 u'1 ) + Q' =cp T02 (r¯ 2 u'2 + r'2 u¯2 ) + r¯ 2 u¯2 (cp T'2 + u¯2 u'2 ),

(3.7)

where T0 is the mean stagnation temperature T + u¯ /cp . Care needs to be taken to recover the jump conditions for zero mean flow from equations (3.5)–(3.7). The strength S of the entropy wave enters these equations only in the combination u¯2 S. In the limit u¯2 :0, S tends to infinity, while the product u¯2 S remains finite. It is convenient to begin an investigation of low mean flow Mach numbers by noting that equation (3.5) may be used to recast the energy equation (3.7) in the form: 1 2

2

r¯ 2 u¯2 (cp T'2 + u¯2 u'2 ) = Q' − cp (T02 − T01 )(r¯ 1 u'1 + r'1 u¯1 )+r¯ 1 u¯1 (cp T'1 + u¯1 u'1 ).

(3.8)

After using equation (3.4d) to expand cp T'2 and taking the limit u¯ :0, this simplifies to {r¯ 2 u¯2 c¯22/(g − 1)cp }S eiv(t − b/u¯2 ) = Q' − cp (T2 − T1 )r¯ 1 u'1 .

(3.9)

It then follows from a combination of equations (3.4c) and (3.9) that, for low mean flow Mach numbers, (g − 1)Q' (g − 1)cp (T2 − T1 )r¯ 1 u'1 (g − 1)Q' + =− + (r¯ 1 − r¯ 2 )u'1 , c¯22 c¯22 c¯22

u¯2 r'2 = −

(3.10)

. . 

570

since (g − 1)cp (T2 − T1 )/c¯ = (T2 − T1 )/T2 = 1 − r¯ 2 /r¯ 1 . Finally, substitution for u¯2 r¯ 2' in the equation of mass conservation (3.5) leads to 2 2

r¯ 2 u'1 = r¯ 2 u'2 − (g − 1)Q'/c¯22.

(3.11)

The no flow jump condition (2.8) then follows directly from equation (3.11). The pressure continuity condition (2.7a) can be deduced in a straightforward way from the low Mach number limit of equation (3.6). When there is a mean flow, the strengths of the acoustic and entropy waves are coupled by the continuity conditions (3.5)–(3.7) and need to be solved for together. For an open end the downstream boundary condition has a particularly simple form (see equation (3.2)), that just relates the strengths of the upstream and downstream propagating acoustic waves. For more complex geometries when, for example, the duct is terminated by a nozzle, the entropy and acoustic waves may also be coupled by the exit boundary condition. A description of how the unsteady flow perturbs the rate of heat input is needed before the frequency of the thermoacoustic oscillation can be determined. As in section 2, two particular forms for the relationship between the rate of heat input and the flow will be considered in detail. The first case is that of no unsteady rate of heat input; Case I,

Q'(t) = 0,

no unsteady rate of heat input.

(3.12)

The second case is that of no unsteady heat per unit mass, the constant heat input per mass being that required to raise the mean stagnation temperature from T01 to T02 : Case II, Q'(t) = cp (T02 − T01 )(r¯ 1 u'1 + r'1 u¯1 ),

no unsteady heat input per unit mass. (3.13)

Substitution for the unsteady flow from equations (3.3) and (3.4) into the boundary conditions (equations (3.1) and (3.2)), the jump conditions (equations (3.5)–(3.7)) and the heat release equation (either equation (3.12) or (3.13)) leads to six homogeneous equations for six unknowns. These equations can be written in the form A F B G G C XG D G 2 −ivb/u¯ G S(r2 c¯2/cp ) e f Q /c¯1

2

J G G G = 0, G G j

(3.14)

where X is a 6 × 6 matrix. For Case I, 1) −(1 − M 0 F 1 + M 1 G (1 + M 1 )e1 −(1 − M 1 )e2 −(1 + M 2 )e3 c¯1 /c¯2 G (1 + M 1 )2e1 (1 − M 1 )2e2  2 )2e3 −(1 + M X =G G(1 + M 1 )a1 e1 (1 − M 1 )a2 e2 −(1 + M 2 )a3 e3 c¯2 /c¯1 0 0 e5 G f 0 0 0

0

0

(1 − M  2 )e4 c¯1 /c¯2

M  2 c¯1 /c¯2

−(1 − M  2 )2e4

M  22

−(1 − M  2 )a4 e4 c¯2 /c¯1

e6 0

1 2

M  32c¯2 /c¯1

0 0

0J 0G

G G 1G 0G 1j 0

(3.15)

  

571

with e1 = e−ivb/c¯1 (1 + M 1 ),

e2 = eivb/c¯1 (1 − M 1 ),

e3 = e−ivb/c¯2 (1 + M 2 ),

e4 = eivb/c¯2 (1 − M 2 ),

e5 = e−ivl/c2 (1 + M 2 ),

e6 = eivl/c¯2 (1 − M 2 ),

a1 = M  1 + 12 M  21 + (g − 1)−1, a3 = M  2 + 12 M  22 + (g − 1)−1

a2 = M  1 − 12 M  21 − (g − 1)−1, and

a4 = M  2 − 12 M  22 − (g − 1)−1.

For Case II, with no unsteady heat input per unit mass, the last row in the matrix defining X is to be replaced by −cp (T02 − T01 )(1 + M  1 )e1 /c¯12

cp (T02 − T01 )(1 − M  1 )e2 /c¯12

0

0

0

1.

(3.16)

Equation (3.14) has a unique trivial solution A = B = C = D = S = Q = 0 unless the determinant of X vanishes. Therefore only disturbances for which det X = 0 can propagate. For given mean flow properties, this is an equation for v and, in general, it has complex roots. Since the time dependence is of the form eivt, the sign of Imaginary v determines whether disturbances grow or decay. Real v gives the frequency of the mode. It is a straightforward matter to determine the zeros of det X numerically, and some results are shown in Figure 5. Inspection of the form of the matrix X in equation (3.15) shows that the leading order effect of the mean flow is O(M ). However, it is apparent from the numerical results that the importance of the the mean flow depends not only on the inlet Mach number and stagnation temperature ratio (which affect the value of M  2 ), but also on the form of the coupling between the heat input and the unsteady flow. It is striking how significant Mach number effects can be. For example, one might be tempted to think that an inlet Mach number of 0·15 is sufficiently small for the mean flow to be neglected. But at this Mach number, the frequency of the thermoacoustic oscillation for Case II, T02 /T01 = 6, is reduced to half its no flow value. Care needs to be taken before simply concluding on the basis of a low inlet Mach number that the mean flow may be neglected.

Figure 5. The lowest frequency of oscillation as a function of inlet Mach number for b = l/2. T  01 = 288 K, p¯2 = 1 bar and T  02 /T  01 = 6. ——, Case I (no unsteady rate of heat input); – – – –, Case II (no unsteady heat input per unit mass); W, Case I with no mean flow (root of equation (2.12)); Q, Case II with no mean flow (root of equation (2.14)).

. . 

572

4. THE EFFECTS OF DRAG

Most practical ways of adding heat also exert a drag force on the flow. In this section, the effect of this drag on the frequency of thermoacoustic oscillations is investigated. If the instantaneous drag can be modelled quasi-steadily, the drag force per unit duct area may be written in the form 12 CD r1 u12(t), where CD is a drag coefflcient. In particular, when a bluff body or grid is used for heat transfer or to stabilize a flame, its blockage reduces the duct area available for the flow. If the expansion from this reduced area is abrupt, an elementary one-dimensional quasi-steady analysis can be used to relate the drag coefficient to the blockage ratio r. This gives [12] CD = (1 − (1 − r)−1 )2.

(4.1)

For linear perturbations the momentum equation (3.6) becomes p'1 + r'1 u¯12 + 2r¯ 1 u¯1 u'1 = p'2 + r'2 u¯22 + 2r¯ 2 u¯2 u'2 + 12 CD r'1 u¯12 + CD r¯ 1 u¯1 u'1 ,

(4.2)

and the third row of the matrix X in equation (3.15) is replaced by (1 + 2M 1s + M  21s)e1

(1 − 2M 1s + M  21s)e2

−(1 + M  2 )2e3

−(1 − M  2 ) 2e 4

M  22

0, (4.3)

where s = 1 − 12 CD = 12 + (1 − r)−1 − 12(1 − r)−2,

(4.4)

after substitution for CD from equation (4.1). Results for various inlet Mach numbers and blockage ratios are shown in Figure 6. The drag force exerted by a grid or flame holder with a blockage ratio of 25% or less is found to have a negligible effect on the frequency of thermoacoustic oscillations for the range of inlet Mach numbers considered. A blockage ratio of 50% alters the frequency for Mach numbers greater than 0·1. 5. THE EFFECTS OF DISTRIBUTED HEAT INPUT

So far, the heat input has been considered as concentrated at a single plane, but in many practical applications it is distributed over a significant length of the resonator. For

Figure 6. The lowest frequency of oscillation as a function of inlet Mach number for b = l/2. T01 = 288 K, p¯2 = 1 bar and T02 /T01 = 6 and various values of the blockage ratio. ——, No drag; – – – – , 25% blockage ratio; ····, 50% blockage ratio.

 

573

example, when a flame is stabilized in a duct in the wake of a bluff body, the combustion is initiated as the flow passes the flame holder, but the fluid continues to burn throughout the downstream portion of the duct [4]. The conditions under which a distributed heat input can be considered as concentrated are discussed in this section. It is sometimes said that heat input with axial extent d can be treated as lumped at a single plane provided d is small in comparison with the acoustic wavelength (or more precisely provided d W c¯ /v). However, this is an oversimplification when there is a mean flow. Then there is also an entropy wave, with the much shorter length scale u¯ /v (see equation (3.4)). In a region of distributed heat input’ the strengths of the acoustic and entropy waves change continuously with position. It is then easier to determine the flow by working directly from the equations of motion. The one-dimensional equations of mass, momentum and energy conservation for the mean flow are respectively: d (r¯ u¯ ) = 0, dx

d (p¯ + r¯ u¯ 2) = 0, dx

d q¯ (x) (c T + 12 u¯ 2) = , dx p r¯ u¯

(5.1–5.3)

where q(x, t) is the rate of heat input per unit volume. When q¯ (x) is specified as a function of x, and the flow at one axial position is known, equations (5.1)–(5.3) can be integrated with respect to x in a straightforward way to determine the mean flow throughout the duct. Linear perturbations can be calculated in a similar way. For oscillations proportional to eivt, the unsteady one-dimensional equations of mass, momentum and energy may be written in the form (d/dx)(r¯ u' + r'u¯ ) = −ivr',

(5.4)

(d/dx)(p' + r'u¯ + 2r¯ u¯u') = −iv(r¯ u' + r'u¯ ), 2

(5.5)

(d/dx)[(r¯ u' + r'u¯ )(cp T + u¯ ) + r¯ u¯ (cp T' + u¯u')]=q' − iv[r'(cv T + u¯ ) + r¯ (cv T' + u¯u')]. (5.6) 1 2

2

1 2

2

The upstream boundary conditions determine the inlet flow perturbations. Then, once the mean flow, the frequency v and the relationship between q' and the unsteady flow has been specified, equations (5.4)–(5.6) can be integrated along the duct, thus determining the flow perturbations at all positions. The details are given in reference [5]. At a general value of v, the exit boundary condition is not satisfied. It is therefore necessary to iterate in v, at each state calculating the flow in the duct, until the complex values of v, for which the exit boundary condition is met, are determined. These are the frequencies of the thermoacoustic oscillations. Only disturbances with these particular frequencies satisfy all the boundary conditions and can exist as free modes of the heated duct. As an illustrative sample, consider the duct geometry of section 3, where the inlet flow is isentropic and choked, and the exit is open. The inlet boundary condition is described by equation (3.1), together with the isentropic relationships r' = p'/c¯12 and cp T' = p'/r¯ 1 . Equation (3.2) shows the exit boundary condition to be p'(l, t) = 0. The mean heat input is uniformly distributed over a length d, centred on x = b: q¯ (x) =

6

cp r¯ 1 u¯1 (T02 − T01 )/d 0

7

for =x − b= Q 12 d . for =x − b= q 12 d.

(5.7)

The lowest frequency of oscillation has been calculated for this distributed heat input by integrating equations (5.1)–(5.6), for Case I: q'(x, t) = 0,

no unsteady rate of heat input per unit volume.

(5.8)

. . 

574

Figure 7. The lowest frequency of oscillation as a function of stagnation temperature ratio for a uniform distribution of mean heat input over a distance d. Case I, no unsteady rate of heat input, q' = 0, with M 1 = 0·1, b = l/2, T01 = 288 K, p¯2 = 1 bar. ——, Concentrated mean heat input; – – – –, d=0·05l; ····, d = 0·25l.

The frequencies predicted for two values of d/l are shown in Figure 7. Also plotted for comparison are results for the concentrated mean heat input q¯ (x) = cp r¯ 1 u¯1 (T02 − T01 )d(x − b). These are the roots of det X = 0, where X is the 6 × 6 matrix in equation (3.15). It is apparent from these plots that there are significant differences between the predicted frequencies for concentrated and distributed mean heat input. For example, distributing the mean heat over a length as short as 5% of the duct can lead to a 25% change in the frequency of the oscillation at the higher temperature ratios, even though vd/c¯1 is only 0·1. The detailed form of the axial distribution of mean heat input is not so important. In Figure 8, results are given for a triangular distribution in the rate of mean input:

8

4cp r¯ 1 u¯1 (T02 − T01 )(x + 12 d − b)/d q¯ (x) = 4cp r¯ 1 u¯1 (T02 − T01 )(b + 12 d − x)/d 0

9

for b − 12 d Q x E b for b E x Q b + 12 d . for =x − b= q 12 d.

(5.9)

The plots are virtually indistinguishable from those in Figure 7, where the mean input is distributed uniformly over the same axial extent d and with the same centroid b.

Figure 8. As Figure 7, but with a triangular distribution of mean heat input.

 

575

Hence, even if the mean heat has only a modest spatial distribution over a length d E 0·1c¯ /v, it can still lead to appreciably different frequencies of oscillation from that when the heat input is concentrated at a single axial plane. These numerical results can be interpreted and explained analytically. In the Appendix, the conservation equations (5.4)–(5.6) are integrated across a region of heat input. There it is shown (see equation (A13)) that the pressure and velocity fluctuations on either side of a region of heat input are unaffected by its spatial distribution over a length d provided vd W1 c¯1

v cp

and

g

2

r¯ (u¯ − u¯2 )s' dx W p'.

(5.10)

1

Here s' denotes the perturbation in entropy and 1 and 2 are the positions x = b − 12 d and b + 12 d on either side of the region of heat input. Entropy fluctuations, with their shorter length scale u¯ /v are more strongly influenced by a distribution in heat input. If the entropy waves propagating downstream of a region of heat input of axial length d are to have the same strength as those produced by concentrated heat addition, it is evident from equation (A15) that d must be small in comparison with u¯ /v. For a low Mach number mean flow this is a strong constraint. In the examples considered in this paper, the pipe exit is open and the boundary condition of zero pressure perturbation does not involve the entropy waves. Hence, the two conditions (5.10) are sufficient to ensure that the frequency of oscillation is unaffected by a spatial distribution of heat input. The entropy fluctuations in the region of heat input must be determined before the integral in conditions (5.10) can be estimated. For a mean flow and linear perturbations, the entropy equation rT Ds/Dt = q(x, t) may be written in the form

0

1

1s' q¯R q' u' p' 1s' − − . + u¯ = q¯ u¯ p¯ 1t 1x p¯

(5.11)

Consider Case I, in which q'(x, t) = 0. For disturbances of frequency v and low mean flow Mach numbers, equation (5.11) simplifies to ivs' + u¯ 1s'/1x = −q¯Ru'/p¯u¯ .

(5.12)

When vd/u¯ is very small in comparison with unity, u¯ 1s'/1x is the largest term on the left side of equation (5.12). After substitution for q¯ (x) from equation (5.7), integration leads to

00

u¯s'(x) 0 O cp

1

T2 x − b + 12 d −1 u'1 d T1

1

for

vd W 1, u¯

(5.13)

where the initial condition s'1 = 0 has been used. Note that the estimate of u¯2 s'2 obtained by putting x = b + 12 d in expression (5.13) agrees with the value calculated in equation (3.9) for d = 0. Expression (5.13) may be used to estimate the integral in conditions (5.10) and gives v cp

g

2

1

r¯ (u¯ − u¯2 )s' dx 0 O

0 0

1 1

vd T2 − 1 p' c¯1 T1 2

for

vd W 1. u¯

(5.14)

Hence the two conditions (vd/c¯1 )(T2 /T1 − 1)2 W 1 and vd/u¯ W 1 are sufficient to ensure that the frequency of oscillation is not altered by the distribution of heat input over a small length d.

. . 

576

When the heat input is caused by combustion, the mean flow velocities are often so low that the limit vd/u¯ w 1 is more appropriate. Then ivs' is the largest term on the left side of equation (5.12) and substitution for q¯ (x) from equation (5.7) leads to the estimate ivs'2 0 O

00

11

cp T2 − 1 u'1 , d T1

(5.15)

and hence v cp

g

0 0

2

r¯ (u¯ − u¯2 )s' dx 0 O r¯ 1 u¯1

1

1 1

T2 − 1 u'1 . T1 2

(5.16)

It then follows from the inequalities in (5.10) that the conditions for the frequency of oscillation to be unaffected by a distribution in heat input are vd W1 c¯1

and

M 1

0

T2 −1 T1

1

2

W 1,

with

vd w 1. u¯

(5.17)

It is now possible to explain the numerical results in Figure 7. For d = 0·25l, vd/c¯1 is not small and the frequency of oscillation is substantially changed by the axial distribution in heat input. For d = 0·05l, vd/c¯1 is much smaller than unity but vd/u¯1 is not. It is then clear from expression (5.17) that for temperature ratios near unity the effects of the distribution should not be significant but that it becomes important at higher temperature ratios. These predictions are confirmed by the numerical results shown in Figure 7. The second of conditions (5.17) implies that the sensitivity of the frequency of oscillation to the axial extent of the heat input is reduced at low Mach numbers. That is confirmed by the results in Figure 9 for M  1 = 0·01. A comparison with Figure 7 shows that the distribution of heat input does indeed have less effect on the frequency of oscillation at very low Mach numbers. It is evident from equation (5.11) that the particular form of unsteady heat input for which q'/q = (u'/u¯ ) + (p'/p¯ ) (5.18) leads to s'(x, t) 0 0. When there are no entropy waves the conditions (5.10) reduce to the simple requirement that d be acoustically compact. The frequencies of thermoacoustic oscillations for this form of unsteady heat input are shown in Figure 10. For distributed heat input, the frequency was calculated by integrating

Figure 9. As Figure 7, but with M 1 = 0·01.

 

577

Figure 10. As Figure 7, but for q'/q¯ = (u'/u¯ ) + (p'/p¯ ).

equations (5.1)–(5.6), with the unsteady heat input (5.18). Results for the concentrated mean heat input were found by replacing the last row in the matrix X in equation (3.15) by −cp (T02 − T01 )(1 + gM  1 )e1 /c¯12

cp (T02 − T01 )(1 − gM  1 )e2 /c¯12

0

0

0

1,

(5.19)

and finding the roots of det X = 0. As predicted, for the form of unsteady heat input (5.18), the frequency of oscillation is less affected by the distribution of heat input than that for Case I with q' = 0 shown in Figure 7. The results for Case II, in which the rate of release is proportional to the mass flow rate, show similar trends. For a concentrated heat input, the frequency of oscillation is a root of det X = 0, where X is described by equation (3.16). When the heat input is distributed, the unsteady heat input generalises to Case II,

q' u' r' = + , q¯ u¯ r¯

no unsteady heat input per unit mass.

(5.20)

This is close in form to equation (5.18) and only weak entropy waves are generated. Once again, the frequency of oscillation is altered little by the distribution in heat input. In summary, entropy waves have a short length scale, u¯ /v, and are more affected by an axial distribution in the heat input than the acoustic waves with their longer length scale

Figure 11. As Figure 7, but for Case II with q'/q¯ = (u'/u¯ ) + (r'/r¯ ).

578

. . 

c¯ /v. Acoustic and entropy waves are coupled through a region of heat input. The two conditions vd/c¯1 W 1 and v f r¯ (u¯ − u¯2 )s' dx/cp W p' must be satisfied, if the acoustic waves are to be unaffected by the distribution of heat input over an axial distance d. Results in Figure 7, for M  1 = 0·1 and Case I with no unsteady rate of heat input per unit volume, show that distributing the heat input over a length as short as 5% of the duct can lead to a 25% change in the frequency of oscillation at the higher temperature ratios. To guarantee in general that heat input distributed over a length d can be treated as concentrated, one requires both vd/u¯1 W 1 and vd/c¯1 W 1. This is very restrictive. However, flow conditions that reduce the volume term f r¯ (u¯ − u¯2 )s' dx, that accounts for the coupling between the acoustic and the entropy waves, reduce the sensitivity of the thermoacoustic oscillations to the axial extent of the heat input. The results in Figures 7–11 show that this can be achieved in several ways: (i) A reduction in the mean temperature ratio. The results for small (T02 − T01 )/T01 in Figure 7 show that then the frequency of oscillation is only slightly affected by a modest distribution in heat input. (ii) A reduction in mean inlet Mach number. A comparison of Figures 7 and 9 shows that this reduces the effect of the axial extent of the heat input. (iii) Particular forms for the unsteady heat input. The unsteady heat input in equation (5.18) generates no entropy waves, while only weak waves are produced when the heat input per unit mass is constant. The results in Figures 10 and 11 show that then the frequencies of the thermoacoustic oscillations are insensitive to a spatial distribution in heat input, when the downstream boundary condition involves only the acoustic waves. 6. CONCLUSIONS

Model problems with very simple geometries have been considered to investigate the influence of various flow effects on the frequency of thermoacoustic oscillations. The form of the coupling between the heat input and the unsteady flow has been demonstrated to have a crucial effect on the frequency of oscillation. Indeed, for the particular case of heating at a single axial plane in a duct with one open and one closed end, there is nearly a 60% difference between the frequency predicted for no unsteady heat input per unit mass and no unsteady rate of heat input. A number of calculation methods recommended in the literature have been tested by applying them to model problems. This has shown that they do not account fully for this effect. The presence of a mean flow significantly complicates the analysis of a thermoacoustic oscillation, since an entropy wave becomes coupled to the acoustic field. It is, therefore, tempting to neglect the mean flow, whenever the inlet Mach number is low. However, great caution must be exercised before making such an assumption. Mean flow effects are found to be surprisingly significant. For example, with a mean stagnation temperature ratio of six, the frequency of a thermoacoustic oscillation for an inlet Mach number of 0·15 can be reduced to half its no flow value. Many practical ways of adding heat also present a blockage to the flow and exert a drag force on the fluid. A drag force exerted by a grid or flame holder with a blockage ratio of 25% or less is found to have a negligible effect on the frequency of thermoacoustic oscillations for inlet Mach numbers in the range 0 E M  1 E 0·15. A blockage ratio of 50% alters the frequency for inlet Mach numbers greater than 0·1. Even a modest distribution of the heat input over an axial distance d can lead to a significantly different frequency of oscillation from that when the heat input is concentrated. For example, even an axial extent of heat input as short as 5% of the duct length can lead to a 25% change in the frequency of oscillation. Entropy waves have a short

 

579

length scale u¯ /v and are more affected by an axial distribution in the heat input than the acoustic waves with their longer length scale c¯ /v. Hence, to guarantee, in general, the heat input distributed over a length d can be treated as concentrated, one requires both vd/u¯ W 1 and vd/c¯ W 1. However, these conditions can be relaxed when the downstream boundary condition involves only acoustic waves, and the acoustic waves (and hence the frequency of oscillation) are not affected by the entropy waves. This is the case at low mean flow Mach numbers or for mean temperature ratios near unity. It also occurs for certain forms of the unsteady heat input, which lead to no or only very weak entropy waves. In all these cases heat input, distributed over a distance d, can be considered as concentrated provided that d is compact. REFERENCES 1. L R 1896 The Theory of Sound, London: Macmillan. 2. S. M. C and T. J. P 1988 Proceedings of the Institute of Acoustics 10, 103–153. Interactions between acoustics and combustion. 3. F. E. C. C 1988 AGARD-CP-450. Combustion instabilities in liquid-fuelled propulsion systems—an overview. 4. P. J. L 1988 Journal of Fluid Mechanics 193, 417–443. Reheat buzz—an acoustically coupled combustion instability, part I: experiment. 5. G. J. B, A. P. D and P. J. L 1988 Journal of Fluid Mechanics 193, 445–473. Reheat buzz: an acoustically coupled combustion instability, part 2: theory. 6. U. G. H, D. R, B. T. Z and B. R. D 1987 AIAA-87-0216. Fluid mechanically coupled combustion-instabilities in ramjet combustors. 7. A. P. D 1988 AGARD-CP-450. Reheat buzz—an acoustically coupled combustion instability. 8. F. A. W 1965 Combustion Theory. Reading, Massachusetts: Addison-Wesley. 9. U. G. H, D. R and B. T. Z 1988 American Institute of Aeronautics and Astronautics Journal 26, 532–537. Sound generation by ducted flames. 10. A. P. D and J. E. F W 1983 Sound and Sources of Sound. Chichester: Ellis Horwood. 11. A. M. C 1982 Journal of Fluid Mechanics 121, 59–105. Low frequency sound radiation and generation due to the interaction of unsteady flow with a jet pipe. 12. A. P. D and G. J. B 1984 AIAA-84-2321. Reheat buzz—an acoustically driven combustion instability. APPENDIX

The aim of this Appendix is to determine conditions under which equations (5.4)–(5.6), describing changes through a distributed region of heat input, integrate to give the jump conditions (3.5)–(3.7) for concentrated heat input. Integration of equations (5.4)–(5.6) across the region of heat input, from position 1, x = b = 12 d, to 2, x = b + 12 d, leads to

g g g

2

[r¯ u' + r'u¯ ]21 = −iv

r' dx,

(A1)

(r¯ u' + r'u¯ ) dx,

(A2)

1

2

[p' + r'u¯ 2 + 2r¯ u¯u']21 = −iv

1

[(r¯ u' + r'u¯ )cp T0 + r¯ u¯ (cp T' + u¯u')]21 − Q'=−iv

2

(r'(cv T + 12 u¯ 2) + r¯ (cv T' + u¯u')) dx,

1

(A3) where Q' = f q'(x) dx. It is evident that equations (A1)–(A3) reduce to the concentrated heat input jump conditions when their right sides are negligible. 2 1

. . 

580

In a perfect gas, the entropy s = cv ln p − cp ln r + constant. For linear perturbations, this may be rearranged to show that r' = (p'/c¯ 2) − (r¯ /cp )s'. (A4) Similarly, cp T' = (p'/r¯ ) + Ts'. (A5) The inlet boundary conditions ensure that there are no incoming entropy waves and so s'1 = 0. Equations (A1)–(A3) will now be rearranged so that they relate p'2 , u'2 and s'2 to the inlet flow with additional integrals over the region of heat input. First note that r'2 can be eliminated by subtracting the product of equation (A1) and u¯2 from equation (A2) to give p'2 + r¯ 2 u¯2 u'2 = p'1 + r'1 u¯1 (u¯1 − u¯2 ) + r¯ 1 (2u¯1 − u¯2 )u'1 −iv

g

2

(r'(u¯ − u¯2 ) + r¯ u') dx.

(A6)

1

A second equation for p'2 and u'2 is derived by noting that the energy equation, [ru(cp T + 12 u 2)]21 = Q −

g

1 (r(cv T + 12 u 2)) dx, 1t

2

1

(A7)

may be rewritten in the form

$

g pu + 12 ru 3 g−1

%

2

g 0 2

=Q− 1

1

since cp rT = gp/(g − 1). Equation (A8) is equivalent to

$

g pu g−1

%

1

1 p + 1 ru 2 dx, 1t g − 1 2

2

+ r1 u1 [12 u 2]21 + [ru]21 12 u22 = Q − 1

g 0

(A8)

1

1 p + 1 ru 2 dx. 1t g − 1 2

(A9)

For linear perturbations of frequency v this leads to g (u' p¯ + u¯2 p'2 ) + r¯ 1 u¯1 u¯2 u'2 g−1 2 2 g =Q' + (u' p¯ + u¯1 p'1 ) + r¯ 1 u¯12u'1 + (r¯ 1 u'1 + r'1 u¯1 )12 (u¯12 − u¯22) g−1 1 1

g0 2

−iv

1

1

p' + r¯ u¯u' + r'12 (u¯ 2 − u¯22 ) dx, g−1

(A10)

after substitution for [r¯ u' + r'u¯ ]21 from equation (A1). When r' is expanded in terms of pressure and entropy from equation (A4), equations (A6) and (A10) become, respectively, p'2 + r¯ 2 u¯2 u'2 = p'1 + r'1 u¯1 (u¯1 − u¯2 ) + r¯ 1 (2u¯1 − u¯2 )u'1

g 00 2

−iv

1

1

1

p' r¯ s' − (u¯ − u¯2 ) + r¯ u' dx, c¯ 2 cp

(A11)

and g (u' p¯ + u¯2 p'2 ) + r¯ 1 u¯1 u¯2 u'2 g−1 2 2 g =Q' + (u' p¯ + u¯1 p'1 ) + r¯ 1 u¯12u'1 + (r¯ 1 u'1 + r'1 u¯1 )12 (u¯12 − u¯22) g−1 1 1

g0 2

−iv

1

0

1

1

p' p' r¯ + r¯ u¯u' + 2 − s' 12 (u¯ 2 − u¯22) dx. g−1 c¯ cp

(A12)

 

581

Equations (A11) and (A12) describe the changes in p' and u' across a region of heat distributed over a length d. It is evident that these are equal to the changes for a concentrated heat input provided that the contributions from the integrals are negligible: i.e., provided that vd W1 c¯1

and

v cp

g

2

r¯ (u¯ − u¯2 )s' dx W p',

(A13)

1

and the flow is subsonic. These conditions are sufficient to ensure that the acoustic waves on either side of a region of heat input are unaffected by its distribution over a finite length. To determine s'2 , the strength of the entropy wave, note that equations (A1) and (A3) may be combined to give a modified energy equation: r¯ 1 u¯1 (cp T'2 + u¯2 u'2 − cp T'1 − u¯1 u'1 ) + (r¯ 1 u'1 + r'1 u¯1 )cp (T02 − T01 ) − Q'

g0

1

2

=iv

r'(cp T02 − cv T − 12 u¯ 2) − r¯ (cv T' + u¯u') dx.

1

(A14)

Substitution for r' and T' from equations (A4) and (A5) into equation (A14) leads to

0

r¯ 1 u¯1 T2 s'2 +

1

g0 2

+iv

1

0

1

p' p'2 + u¯2 u'2 = Q' + r¯ 1 u¯1 1 + u¯1 u'1 − (r¯ 1 u'1 + r'1 u¯1 )cp (T02 − T01 ) r¯ 2 r¯ 1

−s'

1

r¯ p' (c T − 1 u¯ 2) + 2 cp (T02 − T0 ) − r¯ u¯u' dx. cp p 02 2 c¯

(A.15)

Estimation of the integral in equation (A15) shows that vd/u¯ W 1

(A16)

is required if the strength of an entropy wave produced by heat addition is to be unaffected when the heat input is distributed over a length d. The implications of the conditions (A13) and (A16) are discussed in section 5.