CC BY
71Electronic Journal «Technical acoustics» http://webcenter.ru/~eeaa/ejta/
2003, 12
K. I. Matveev and F. E. C. Culick
California Institute of Technology, Pasadena, CA, 91125, USA, [email protected]
Limit-Cycle Properties of a Rijke Tube
Received 24.05.2003, published 12.06.2003
Thermoacoustic instability appears when unsteady heat release is favourably coupled with acoustic pressure perturbations. The important technical applications involving thermoacoustics are combustion instability in rocket motors and low-pollutant lean flames; noisy industrial burners; pulsed combustors; and thermoacoustic engines. The simplest device for studying thermoacoustic instability is a Rijke tube. In this work, a series of experiments is carried out to determine the nonlinear behavior of the transition to instability and the excited regimes for an electrically driven Rijke tube. A hysteresis effect in the stability boundary is observed. A mathematical theory involving heat transfer, acoustics, and thermoacoustic interactions is developed to predict the transition to instability and limit-cycle properties.
INTRODUCTION
Thermoacoustic instability is defined as excitation of acoustic modes in chambers with heat sources due to the coupling between acoustic perturbations and unsteady heat addition. When the heat release has a component in phase with the pressure fluctuation, then, according to Rayleigh’s criterion [1], acoustic oscillations are encouraged. If the power transformed from the heat released to the acoustical motions exceeds the losses, then instabilities grow until nonlinear factors limit this development. Thermoacoustic instabilities appearing in thermal systems with a mean flow, such as various combustors and heat exchangers, can be conveniently studied by employing a Rijke tube, first built by Rijke [2] in the 19th century. This device consists of a duct with a permeable heater inside. In the case of horizontal orientation, airflow through the tube is provided a fan. For some combinations of the main system parameters (airflow rate, heater location, and power supplied), a high-intensity sound is generated in the tube.
The qualitative reasons for Rijke oscillations are well understood at the present time [3]. Unsteady heat transfer from the heater to the airflow depends on the flow oscillations. In unsteady flow, the instantaneous heat transfer rate lags the current flow velocity due to the thermal inertia of the heat transfer process. At certain positions of the heater, an unsteady heat addition component appears that is in phase with pressure fluctuation, and hence, thermoacoustic instability is encouraged in accordance with Rayleigh’s criterion [1]. Despite the long history of experiments with Rijke tubes, little attention has been previously paid to the assessment of the role of experimental procedures and the main parameter fluctuations that may significantly affect the operating point defined by a given set of main parameters. The authors have been conducting a project, aimed at determining the transition to instability
in a Rijke tube with known uncertainties [4, 5]. The primary goal of the present work is to investigate, both experimentally and mathematically, hysteresis and the limit-cycle properties of the electrically driven Rijke tube built at Jet Propulsion Center, Caltech.
1. EXPERIMENTAL SYSTEM
The scheme of the Rijke tube studied in this work is shown in Figure 1a. Details of the experimental setup are available in [5]; here, only a brief description of the system is given. The main part of the apparatus is a one-meter aluminum tube having square cross section area
0.009 m . A fan provides a mean flow. The flow direction is from right to left in the figure. A damping chamber is installed to prevent interaction between the blower and tube acoustics. The tube is placed horizontally to avoid the influence of natural convection on the mass flow rate. Therefore, the airflow rate is an independent system parameter for this configuration and can be controlled using a valve in front of the blower.
The heat source is electrically powered 40-mesh nichrome gauze. A sketch of the heating element is shown in Figure 1b. The gauze is mounted on a frame made of mecor to avoid thermal and electric contacts with the tube walls. Copper rods are used to feed the gauze with electric current delivered from a power supply. The constructed heating element possesses properties required for both experimental and theoretical study: the heater can repeatedly withstand high temperatures for long time intervals; it is able to release large amounts of thermal power (over 1 kW in our system) without geometrical changes; the air flow is heated more or less uniformly over the cross-section; and the longitudinal size of the heater is much smaller than characteristic dimensions in the system, so it can be considered infinitely thin in acoustic modeling.
INSULATION FRAME
COPPER STRIPS
DAMPING CHAMBER
RIJKE TUBE
7
HEATER
BLOWER
POWER SUPPLY
GAUZE
(a)
(b)
COPPER RODS (TO POWER SUPPLY)
+
Figure 1. (a) Experimental setup; (b) heater configuration
Heater location, airflow rate and supplied power are controlled parameters. For certain values of these parameters, the system becomes unstable and a high-intensity sound is generated. The goal of the experiment is to quantify the transition to instability and limit-cycle properties as functions of main system parameters. A computerized data acquisition system acquires time resolved measurements of pressure fluctuation, airflow, and power with
a scan rate up to 8000 Hz; and a thermocouple array provides the spatially resolved temperature field along the tube. A special experimental procedure to obtain robust measurements is described in [5].
2. LINEAR MODELING
A one-dimensional linear wave equation for the pressure fluctuation p in the medium, having non-uniform temperature field and containing volumetric sources of intensity Q and heat addition per unit of time and volume Q', is applied
d2p 2 d2p a2 dp0 dp' ^.dQ' 2 dQ
-2- - a +-----------— — = (y- 1) — + p0a ----. (1)
dt dx p0 dx dx dt dt
The sound velocity a and mean gas density po are time-averaged values, but they are
functions of temperature, which in turn is dependent on a horizontal coordinate. The Mach number of the mean flow in our system is less than 10-3, and its effect on acoustics is much smaller than the influence of the forcing terms in Eq. (1). Solution to the wave equation is sought in the form
p(x, t) = ea/(x), (2)
where a is a complex number; and /(x) is the complex-valued mode shape. By applying boundary conditions at the tube ends and matching conditions at the heating element, the problem is closed. Each nontrivial solution, found numerically, corresponds to a certain eigen
value a = 5 + io; its real part 5 is the growth rate and imaginary part O is the frequency of
an acoustic mode.
The dominant component of the unsteady heat addition in the system is the heat released by the grid. We make use of matching conditions at the grid developed by Merk [6] for the upstream and downstream velocity perturbations in the case of the low Mach number
f f m \ \
u'_, (3)
. V
u + = —
cp+
cpT-1
vcpX
Tr
1 ’lin
1+
V
where Trlin is the linear transfer function between a reduced velocity and heat release
fluctuations; cp is the heat capacity. Subscripts + and - stand for downstream and upstream
sides of the heater. A transfer function, dependent on the mean flow velocity, frequency, heater geometry and fluid properties, is defined by the formula
Q' u' _
or = TrM — ■ (4)
Q 0 u0-
The complex-valued transfer function Trlin = \Triin\e~lip relates the unsteady component of the
heat transfer rate Q' to the fluctuating velocity u'_ upstream of the heater. Phase delay (p of the heat transfer with respect to the velocity oscillation is a consequence of thermal inertia in the heat transfer process.
The heating element in the system studied is the gauze consisting of thin wires. For this configuration and characteristic flow regime, the most relevant transfer function can be extracted from the data reported by Kwon and Lee [7]. They studied numerically the unsteady heat transfer from an isothermal cylinder in the range of Reynolds number (from 1 to 10) similar to our case. In this regime, Oseen’s flow approximation is no longer valid and a boundary layer is not yet formed; therefore, there are no analytical expressions available for the transfer function. The magnitude and phase delay of the transfer function TrKL, calculated from data given in [7], are shown in Figure 2. The transfer function depends on the nondimensional mean flow velocity u0 = u0 / -Jox and the non-dimensional radius of the wire
r *= r^Jo/ x ; X being the thermal diffusivity. However, this is not a linear transfer function, since it was computed at a finite amplitude of the reduced velocity fluctuation u'/ylox (which was always equal to one in the referenced study). In order to find the transition to linear instability accurately, this function is transformed in our study to the true linear transfer function using nonlinear theory developed in the next section.
<D "Ö =5 -i—■
£=
CD
cu
Figure 2. Transfer function: (a) absolute value; (b) phase delay
70
60
S’ 50 p.
■ = 40
30
20
10
. w
\ \
\x r =2.0 ■
' r*=0.5 \ \ \ \ '
(b)
3. NONLINEAR THEORY
Exact nonlinear modeling of the Rijke tube would require a prohibitively detailed flow simulation due to the system complexity. Our goal is to stay within the approach of reduced-order modeling, retaining the factors of the most importance. The complex processes involved are approximated introducing reasonable assumptions and hypotheses.
It was shown in [8] that neither nonlinear gas dynamics nor nonlinear boundary conditions at the tube ends are of importance for determining the limit cycle properties in the system considered. The fact that the oscillating velocity magnitude estimated in the vicinity of the heater tends to be stabilized near the mean flow velocity (slightly exceeding it) in the unstable regimes, and grows slowly with the increase of the supplied power, suggest that flow reversal at the heater location is critical for nonlinear modeling. When the flow is reversed, i.e. the total velocity is negative; application of formula (4) may result in convective heat transfer
directed from the cooler flow to the hotter grid, which is impossible. Therefore, Eq. (4) is not valid for the nonlinear case, and must be modified.
Let us first consider a simpler situation, namely heat transfer from a heated isothermal object to a flow having mean and acoustic components. When the characteristic size of the object is much less than an acoustic wavelength and the Mach number is small, the flow may be treated as incompressible. The commonly used form for the convective heat transfer rate in a steady flow, valid over a wide range of the Reynolds number, is
Q 0 = kV^, (5)
where k is a coefficient proportional to the difference between the temperatures of the flow and the object. The heat transfer rate from the heat source in our system obeys this law [4].
When flow is non-steady, consisting of superimposed steady and harmonically fluctuating flow components, Eq. (5) remains valid in the low-frequency and low-amplitude limits of oscillations. If the amplitude of velocity perturbations exceeds a mean flow velocity, then some pockets in the airflow cross the heater plane more than once. Since at the flow reversal the incident air is already pre-heated, the effectiveness of the time-averaged heat transfer decreases. In that case, the quasi-steady heat transfer rate can be approximated by a simple rule:
Qqs (t) = kJ| u(t) |, when u(t) > 0,
(6)
= 0, when u(t) < 0,
where index qs means quasi-steady; and u(t) = u0 + u'(t) is the total flow velocity in the
heater vicinity. Equation (6) is related to the hypothesis proposed by Lehmann [9]. He assumed that the cold air experiences a constant increase in temperature in passing through the heating element the first time, and during flow reversal no heat transfer occurs. If the velocity perturbation is small with respect to the mean flow velocity, then Eq. (6) can be transformed into form of Eq. (4), and the corresponding transfer function will be Trq = 1/2 .
In the nonlinear problem, a sinusoidal variation of the flow velocity with frequency o inevitably produces a heat response with other frequency components. We define the nonlinear transfer function Tr in the following way:
= Tr , (7)
uo
o
where (f )O designates a component of the function f that oscillates with frequency O. The natural choice for the parameter measuring nonlinear behavior is the ratio of the oscillating velocity amplitude and the mean velocity a = ux' / u 0, estimated in the vicinity of the heater.
Now we propose the main hypothesis in our analysis: the frequency dependence of the nonlinear transfer function is assumed to be the same as for the linear transfer function; and the nonlinearity is taken into account by a coefficient G(a), which is a function of the parameter a,
Tr («a) = G(a)Trlm (ffl). (8)
The mean heat transfer rate in the nonlinear case is assumed to be equal to the quasi-steady heat transfer rate averaged over a period of oscillation:
Q Un (a) = Q Is (a). (9)
Symbol un is introduced to differentiate the heat transfer in the real unsteady flow from the steady flow and the quasi-steady flow model.
Applying hypothesis (8) to the quasi-steady process, we obtain the formula
Trqs (o,a) = G(a)TrZ («). (10)
Taking into account that Trqsn = 1/2, the nonlinear component of the transfer function is
expressed as follows:
G(a) = 2Trqs (a). (11)
Therefore, with the assumptions used here, a quasi-steady transfer function and a linear transfer function completely describe a nonlinear transfer function.
Let us calculate Trqs using Eq. (6) and a sinusoidal function for the oscillating component in the total velocity u'(t) = u 0asin(«t). From Eq. (6) it follows that the mean quasi-steady
heat transfer rate Q 0qs is
2n
Qqs = JV1 + asin<p|d<p, (12)
2n 0
where Q0 is defined by expression (5). The unsteady component of the heat transfer rate is
Q'qs (t) = Qo
1 + a sin (ot\
1 2n
)— J^|l + a sin <p|dq
2n 0
(13)
The quasi-steady nonlinear transfer function defined by Eq. (7) can be computed from the following equation:
(t) = Trqs a sin(ft)t). (14)
Q T
vQ 0 )«
Substituting Eq. (12) and (13) into expression (14) and taking a Fourier component corresponding to frequency (O, we find the quasi-steady nonlinear transfer function:
1 2 T f A. \qs
Trqs (a) = --aT
QL Q o
(r)sin(oK)dT, (15)
0 ( Q0 )
Nonlinear modifications calculated for the mean heat transfer rate and the transfer function are shown in Figure 3. From Eq. (15) it follows that Trqs (a) ^ 1/2 when a ^ 0,
i.e. G(0) = 1. The nonlinear mean heat transfer rate decreases when the nonlinear parameter increases from zero to one; this is caused by the square root law (Eq. 6). After a certain point,
the heat exchange starts to increase, since the heat transfer rate is enhanced by the growing magnitude of acoustic velocity. The behavior of the transfer function is opposite: it increases at the beginning, achieving maximum near unity of the nonlinear parameter, and then decreases, since in the limit of high oscillating velocity, or, equivalently, low mean flow velocity, the unsteady heat transfer does not possess a component with the frequency of flow oscillation.
Non-monotonic dependence of the transfer function Tr on the fluctuating velocity amplitude results in interesting implications for the excited regimes of operation of a Rijke tube. From Eq. (7) it follows that a magnitude of the unsteady heat release is proportional to the transfer function. The same holds for thermoacoustic energy conversion, if the heat release has a component in phase with pressure perturbation. Therefore, for oscillating velocity amplitudes smaller than the mean flow velocity (compared at the heater location), the instability amplification increases with acoustic amplitude. That means that nonlinear thermoacoustic instability is possible when the system is linearly stable, so multiple system states can exist for the same values of the system parameters. The choice of a particular state is determined by history of parameter variation and other factors, e.g. a level of noise disturbances. For sufficiently high amplitudes of the oscillating velocity, a magnitude of the transfer function decreases; hence, further development of the instability is discouraged. This factor determines the amplitudes of the limit cycles.
Figure 3. Nonlinear corrections for mean heat transfer rate Q0un /Qst, solid line; and for
transfer function Tr / Tr lin, dashed line
Established transformation of the linear to nonlinear transfer function (Eq. 8) allows us to relate the available transfer function TrKL (Figure 2) to the truly linear transfer function Trlin using the following expression:
Tr^n («) = TrG «/y. (16)
G(1/ U0)
The inverse value of the reduced mean flow velocity plays the role of the nonlinear parameter, because the amplitude of the reduced oscillating velocity was always equal to unity in calculating TrKL [7].
4. RESULTS AND DISCUSSION
Experimental investigation of the transition to instability and the limit cycles has been carried out in a way consistent with modeling approach. In individual runs, the heater position and the mass flow rate are fixed, and power is treated as a variable parameter. The power is varied in both directions, increasing and decreasing, in order to capture history-dependent properties of the system. Two heater locations are chosen, 1/4 and 5/8 of the tube length from the upstream end. The first position is close to the optimal one for excitation of the first acoustic mode of the tube. When the heater is in the downstream section of the pipe, it is known that the first mode is not excitable because of the unfavorable coupling between pressure and heat release fluctuations [5]; however, the second mode can become unstable. A set of mass flow rates is selected to cover the range where the system can be unstable within limitations set by the available power.
4.1. Stability boundaries
During the experimental runs aimed at capturing the transition to instability, the variable parameter, power input, was varied slowly enough to provide sufficient time for the temperature field settlement and to avoid unsteadiness in the system conditions, which might result in nonlinear triggering of instabilities [5]. The stability boundary is shown in Figure 4 for both heater locations. Crosses, whose sizes are equal to the experimental errors, present test data. Bold and light crosses designate the boundaries when power is increased and decreased respectively. At the heater location x / L = 1/ 4, the stability curve is split into two branches at high mass flow rates, each corresponding to one of the directions of power variation. Existence of two branches is a manifestation of a hysteresis effect: two different states (stable and excited) can exist at the same values of the system parameters. The presence of a particular state is dependent on the method of varying the parameters, i.e. on the history of the system. For the heater location x / L = 5 / 8, the system is even more history-dependent. When power is increased further after transition to instability, sound disappears at some point, and does not appear when the variation of power is reversed. That is why experimentally obtained stability boundary is shown in Figure 4b for increasing power only. The characteristic value of power needed for system excitation is higher for x / L = 5/8 than that for x / L = 1/4; hence the system is less prone to instability when the heater is in the downstream section of the tube.
The calculated stability boundaries are shown by curves in Figure 4. The details of the computational algorithm are described in [5]. For the heater location x / L = 1/4 and the range of moderate and high flow rates, the model results agree well with the experiment. Deviations noticeable at low flow rates may be due to increased roles of natural convection and possible air leakage in the system. For the heater location x/L = 5/8, accuracy is compromised at high mass flow rates. Since the second mode is responsible for instability, the characteristic frequency is about two times higher than that for the upstream heater location; and the assumption in the acoustic model that the tube and the damping chamber are narrow may be marginal due to the high frequency of the sound waves in this case.
1500
1000
o
Cl
500
1 1 1 1400 I 1 I I
(a) (b)
J 1300 • -
I 1200 } -
EXCITED 4- 1
■ - 1100 ■ -
+ / £ 1000 . _
l / + <5
\ tV i T -+- \ f -l" o 900 Cl ■ \ + -
\ f 800 \
700 ■ +/ -
~"H"~ STABLE 600 - -
STABLE
1 1 1 500 1 1 1 1
1 2 3
Mass Flow Rate [gm/sec]
1.5 2 2.5 3 3.E
Mass Flow Rate [gm/sec]
Figure 4. Stability boundary for heater positions: (a) x / L = 1/4; (b) x / L = 5/8. Experimental data are given with error bars; bold crosses correspond to power increase, light crosses to power decrease. Calculated stability boundary is shown by a curve
4.2. Hysteresis effect
The computational procedure for determining the equilibrium-excited states is similar to calculation of linear stability, except in two respects. First, the time-averaged heat transfer is modified as shown by Eq. (12). Second, the transfer function is dependent on the nonlinear parameter a (Eq. 8); this factor is dominant in determining the limit-cycle amplitude. The results obtained are presented for pressure amplitude at location 3/20 of the tube length from the upstream end of the tube, the location of one of the pressure transducers during the tests.
For the transition to instability, we consider the heater location x / L = 1/4, where test records are available for both directions of power variation. At high flow rates, a significant gap appears between the critical powers, corresponding to the transitions between stable and excited regimes (Figure 5c). Thus, two states can be observed for the same set of system parameters: one state is associated with a limit cycle having finite amplitude; in the other state the amplitude is zero. This effect is identified as hysteresis at the stability boundary. The choice of a certain state is determined by the direction of the power variation. The results obtained by modeling show a similar behavior (Figure 5d). Bold lines correspond to the stable states, with zero and finite amplitudes; and the dashed line designates an unstable equilibrium. The hysteresis pattern is explained by a non-monotonic dependence of the magnitude of the transfer function on the amplitude of oscillating velocity, as discussed in Section 4. A size of the hysteresis loop predicted by the model is significantly larger than that measured.
At low mass flow rates, the test does not exhibit a significant hysteresis. In fact, there is a narrow range of power between stable and unstable regimes where a steady state is not well defined. A tone can appear and disappear, and beating in the sound is sometimes observed.
Pressure magnitudes in the range of the power 165-170 W (Figure 5a) correspond to some arbitrary recorded states. When the power exceeds 170 W, the sound amplitude is more or less stable. The model still produces hysteresis, although with a smaller gap between the critical powers, corresponding to the transitions between stable and excited states (Figure 5b).
Thus, the model developed here explains the mystery of hysteresis, although it significantly overpredicts the effect. The most likely reason for this discrepancy is a noisy environment in the vicinity of the heater inside the tube. Variations of the flow rate even in stable regimes are around 10%, corresponding to the size of the error bars shown in Figure 4. Those oscillations may be caused by the geometry of the system, e.g., the frame that supports the heated grid. Vortices can shed from this frame and disturb the flow. Such perturbations may initiate early transitions between stable and excited states, effectively reducing the gap between critical power levels. That explains the overprediction of the hysteresis effect at high mass flow rates. At the lower flow rates, two additional factors appear that force the hysteresis gap to contract even more: first, the relative error in mass flow rate increases when the mean velocities are lower (Figure 4); second, the predicted gap is smaller. That can lead to total disappearance of hysteresis, similar to the experimental observations (Figure 5a).
Figure 5. Transition to instability: (a) and (c) experiment (A, power increase; V, power decrease); (b) and (d) modeling (bold line, stable limit cycle; dashed line, unstable). Mass flow rate: (a) and (b) 0.89 gm/sec; (c) and (d) 2.75 gm/sec. Heater location x / L = 1/4
4.3. Limit-cycles far from stability boundary
Now consider the limit-cycle properties deep inside the excited regimes for both heater positions. In order to carry out the experiment in a reasonable time, we had to decrease the duration between power increments to cover a wide range of parameters. That inevitably
produces some delay in the temperature field settling due to the system thermal inertia. The accuracy of the results, obtained both experimentally and theoretically in this work for the excited regimes of operations of the Rijke tube far from stability boundaries, is accepted to be lower than the accuracy of the stability boundary. With regard to the limit cycles deep in the excited regimes, we are mostly interested in estimating their properties and approximate modeling.
The calculated and experimental results for the amplitude and frequency of the self-excited (first) acoustic mode in the unstable regimes of operation of the Rijke tube are shown in Figure 6 for the heater location x / L = 1/4 and two mass flow rates, which are fixed in individual runs. The amplitude and frequency of the excited mode are given versus supplied power. Experimental trends in amplitudes and frequencies are predicted correctly. The difference in experimental data for opposite directions of power variation is caused mainly by the thermal inertia of the system; the average temperature was higher when the power was being decreased; that caused higher frequencies.
Figure 6. Properties of the excited acoustic eigen mode at heater location x / L = 1/ 4. Mass flow rate: (a) and (b) 1.63 gm/sec; (c) and (d) 2.40 gm/sec. Measured data: A, power increase;
V, power decrease. Lines correspond to the model results
The calculated and experimental results for the limit-cycle properties of the self-excited mode at the heater location x / L = 5/8 are shown in Figure 7 for two mass flow rates. At these system conditions, the excitable mode is the second acoustic eigen mode of the tube, and its frequency is about twice as big as that of the first mode. The errors are much larger than those for the heater location x / L = 1/ 4, probably due to inaccuracy in acoustic modeling at high frequencies. Although the model predicts marginally the maximum level of pressure magnitude in the excited state, it does not capture the amplitude dependence on the power input as seen in the test data. The tendency in the frequency behavior is modeled well, but the numerical value is smaller by a few percent than the measured frequency.
Figure 7. Properties of the excited acoustic eigen mode at heater location x / L = 5 / 8 . Mass flow rate: (a) and (b) 2.05 gm/sec; (c) and (d) 2.20 gm/sec. Measured data: A, power increase.
Lines correspond to the model results
CONCLUSIONS
An experimental study of the transition to instability and the nonlinear regimes of operation of an electric horizontally oriented Rijke tube has been carried out. A rather simple theory is proposed in this paper for the nonlinear modeling of a complex thermoacoustic phenomenon in the presence of a mean flow. The results confirm the speculation that nonlinear behavior of Rijke tubes of the sort investigated in this work is dominated by nonlinear characteristics of heat transfer from the source of the instabilities.
A hypothesis for the general form of the heater transfer function that accounts for finite flow disturbances is proposed. A quasi-steady model for nonlinear heat transfer, which is valid in the low-frequency limit of oscillations, is utilized for determining the nonlinear correction to the transfer function.
The theory constructed demonstrates satisfactory agreement with experimentally obtained stability boundaries for two locations of the heating element, corresponding to different excitable modes of the Rijke tube. The hysteresis effect in the transition between stable and excited state is quantified experimentally and explained theoretically in terms of nonlinear heat transfer.
Providing the theory of heat transfer is valid, the model developed here can be recommended for approximate analysis and preliminary design of thermal devices where thermoacoustic instability is a concern. For more accurate knowledge of the process in a particular real-world system, sophisticated studies using Computational Fluid Dynamics or detailed experiments must be carried out.
REFERENCES
1. J. W. S. Rayliegh. The Theory of Sound. New York: Dover Publications, 1945 re-issue.
2. P. L. Rijke. Notiz über eine neie art, die luft in einer an beiden enden offenen Röhre in Schwingungen zu versetzen. Annalen der Physik, 1859, vol. 107, pp. 339-343.
3. R. L. Raun, M. W. Beckstead, J. C. Finlinson, and K. P. Brooks. A review of Rijke tubes, Rijke burners and related devices. Progress in Energy and Combustion Science, 1993, vol. 19, pp. 313-364.
4. K. I. Matveev and F. E. C. Culick. Experimental and mathematical modeling of thermoacoustic instabilities in a Rijke tube. 40-th Aerospace Sciences Meeting and Exhibit, Reno, NV, USA, AIAA Paper 2002-1013, 2002.
5. K. I. Matveev and F. E. C. Culick. A study of the transition to instability in a Rijke tube with axial temperature gradient. Journal of Sound and Vibration, 2003, vol. 264, pp. 689706.
6. H. J. Merk. Analysis of heat-driven oscillations of gas flows. Applied Scientific Research, 1956, A6, pp. 402-420.
7. Y.-P. Kwon and B.-H. Lee. Stability of the Rijke thermoacoustic oscillation. Journal of the Acoustical Society of America, 1985, vol. 78, pp. 1414-1420.
8. K. I. Matveev and F. E. C. Culick. Modeling of unstable regimes in a Rijke tube. 5-th International Symposium on Fluid-Structure Interactions, New Orleans, LA, USA, ASME Paper IMECE 2002-33369, 2002.
9. K. O. Lehmann. Uber die theory der netztone. Annalen der Physik, 1937, vol. 29, pp.527-555.
