Abstract
A Rijke tube with a distributed heat source is investigated. Driven by the widely existing thermoacoustic instability in lean premixed gas turbine combustors, this work aims to explore the physicochemical underpinning and assist in the elucidation and analysis of this problem. The heat release model consists of a row of distributed heat sources with individual heat release rates. The integrated heat release rate is then coupled with the acoustic perturbation for thermoacoustic analysis. A continuation approach is employed to conduct the bifurcation analysis and capture the nonlinear behaviour inherent in the system. Unlike the conventional approach by the Galerkin method, the acoustic equations are originally discretized using the Method of Lines (MOL) to build up a dynamic system. The model is first validated and shown to yield good predictions with available experimental data. Influences of multiple heat sources, time delay, and heat release distribution are then studied to reveal the extensive nonlinear characteristics involved in the case of a distributed heat source. It is found that distributed heat source plays an important role in determining the stability of a thermoacoustic system.
1. Introduction
Thermoacoustic instability has been a serious impediment to develop NOx tolerant combustion systems for both aircraft propulsion and power generation gas turbines including rocket motors and industrial burners [1–6]. It results from the interaction between the heat release and acoustic pressure or velocity oscillations within the combustion system. Rijke tube, a typical time-delayed thermoacoustic system, has been a classical tool employed for the study of thermoacoustic instability. It usually consists of an open-end tube and heat source inside it. When the heat source is placed in certain positions along the tube, sound would emit from the tube. The sound is generated as a result of the transfer from unsteady heat release to acoustic energy. In spite of its simple structure, it could demonstrate abundant nonlinear phenomena, such as fixed point, bifurcation, limit cycle, quasiperiodicity, and chaos. Comparing with the expensive and complicated experiments, its simplicity and operability have made it an excellent example to investigate thermoacoustic instability in practical problems [7–10].
Rijke tube has been broadly studied to understand the intrinsic nonlinear behaviors of the thermoacoustic instability during the past decades. Several review works have been conducted by [11–13] regarding the theoretical, experimental, and numerical research on Rijke tube perturbations. The theoretical work by Culick [14, 15] focused on the nonlinear behavior of acoustic waves within a combustion chamber and was extensively referred to by a number of subsequent works on the study of Rijke tube. Heckl [16] developed an empirical heat release model based on experimental results and predicted limit cycle amplitudes and nonlinear effects in the tube. Hantschk and Vortmeyer [17] investigated self-excited thermoacoustic instabilities in the Rijke tube using a commercial CFD code and nonlinearity in the heat flux from the heating source to the flow was found to determine the limit cycle amplitudes. Matveev [18, 19] introduced a nonlinear heat transfer function for the nonlinear stability analysis of a horizontal Rijke tube. Hysteresis phenomenon was reported in the stability boundary and limit cycles were predicted as observed in experiments. Ananthkrishnan et al. [20] captured the dynamics of the nonlinear acoustic waves with reduced-order models via truncating the modal expansions and determined the number of modes required for accurate predictions. Balasubramanian and Sujith [21] studied the role of nonnormality and nonlinearity in thermoacoustic system in a Rijke tube and found out that the nonnormality could result in transient growth of oscillations which can trigger nonlinearities in the system. Subramanian et al. [22] conducted bifurcation analysis of the dynamic behaviors of a horizontal Rijke tube as a function of different system parameters. Nonlinear stability boundaries and other nonlinear phenomena were also observed in the analysis of the thermoacoustic system. Juniper [23] employed flame describing function, adjoint-based approach, and matrix-free continuation to predict the limit cycles and bifurcations in the Rijke tube.
Most of the previous research has focused on the classical Rijke tube with a single compact heat source, either a flame or a hot-wire gauze, allowing for the fact that the heat source is small enough compared with the acoustic wavelength. However, in practical combustion systems, the flame scale is usually considerable under most circumstances, especially with larger fuel flow or output power load, smaller excess air ratio, smaller nozzle spray angle, and so forth [24–26]. In this case, it is not quite proper to deal with the heat source as a single point source.
There have been various previous studies on the effect of distributed heat source which have extended the study of thermoacoustic instability to a wider scope. For example, Dowling’s theoretical work [2] proves that distributed heat input over an axial distance can lead to a significantly different frequency of oscillation from that when the heat input is concentrated. Kim and his coworkers [3] investigated the influence of the spatial heat release distribution using local flame transfer function (FTF) measurements. Heckl [27] employed Green’s function to model the behaviour of a Rijke tube with a distributed heat source and concluded that the heat source distribution has a first order influence on the stability of the thermoacoustic system and could not be ignored.
However, most of them are centred on the frequency-domain analysis. This paper focuses specifically on the time domain and conducts the bifurcation analysis regarding the thermoacoustic system, which has not been attempted previously. In this study, a horizontal Rijke tube with distributed heat source is employed for stability study. The distributed heat source consists of a row of single heat source. A modified form of Heckl’s model is utilized to couple the heat release with acoustics. Method of Lines (MOL) is adopted to discretize the governing equations and a linear multistep method (LMS-method) and Newton iteration method [28, 29] are used for linear stability study. A numerical continuation method is then employed to obtain bifurcation diagrams for investigating nonlinear behavior including the Hopf bifurcation, fold bifurcation, and limit cycles. This model is shown to be accurate and could contribute to understanding the physical nature of the thermoacoustic instability and provide guidance concerning the design and operation of a practical thermoacoustic system.
2. Physical Model
Rijke tube is often oriented vertically in practice, in which the base flow is driven by natural convection. With the view of neglecting this complicated convection, in this paper, a horizontal Rijke tube is of major concern for the instability study of the thermoacoustic system.
Figure 1 shows a schematic of a horizontal Rijke tube with a distributed heat source. The base flow driven by an external fan passes through the tube and is heated up by hot wire gauzes (as discrete heat sources) placed along the tube at locations , where . Naturally, the tube displays an infinite number of acoustic modes. According to Rayleigh criterion, the thermal energy could be transferred to acoustic energy as long as they are in phase [30]. The influence of the mean flow and mean temperature gradients is excluded.
The one-dimensional governing equations for the acoustic momentum and energy are where and are streamwise location and time, respectively; , , and are the fluid density, velocity, and pressure, respectively; the subscript 0 represents the unperturbed variables and the tilde denotes the perturbed variables; , , , and are the heat capacity ratio, generic acoustic damping coefficient, speed of sound, and length of the tube, respectively; is the integrated rate of heat release perturbation per unit volume.
The heat release model is relatively less complex than real thermoacoustic system. It consists of a row of distributed heat sources with individual heat release rate. Here, the heat release rate depends on acoustic fluctuations in the tube, which is a function of the acoustic velocity at the source, delayed by a time delay . It is noted that each heat source has its own heat release rate, time delay, and position along the tube. Thus, by appropriately and reasonably choosing proper values, it could flexibly model different conditions. The interactions between the heat sources are quite complicated and are neglected in this model. One can select different heat release rates and time delays to represent the interaction, though it is factitious. For each discrete heat source, a modified form of King’s law [16] is utilized to model the heat release ratewhere , , and represent the length, diameter, and temperature of the hot wire gauze, respectively; and are the fluid thermal conductivity and heat capacity at constant volume, respectively; is the temperature of the unperturbed base flow; is the cross-sectional area of the tube; is the time delay; is the Dirac Delta function to narrow the heat-release region specifically at the wire position .
Therefore, the integrated unsteady heat release can be given asand then coupled with the acoustic perturbation for thermoacoustic analysis.
The generic damping coefficient [31] is defined aswhere is the kinematic viscosity and is the thermal diffusivity. It comprises the acoustic dissipations resulting from radiation at both ends and the acoustic boundary layer along the tube.
The nondimensional variables are defined as where for the ideal gas, is the universal gas constant, , and is the Mach number (). With these definitions, one can substitute and obtain the governing equations (1) in nondimensional form aswhere is the heat release coefficient defined as It represents the strength of the discrete heat source, incorporating all the details of the fluid, the hot wire gauze, and the tube.
The boundary conditions are treated as and at , 1, which enable the concerned problem to be a well-posed initial value problem within time domain. The acoustic energy per unit volume is selected as the indicator of the acoustic amplitude inside the tube. It is one of the most convenient ways to measure the size of acoustic perturbation and also has physical interpretation [32]. The dimensionless form is defined as Particularly, the minimum acoustic energy on periodic solutions is chosen for the bifurcation analysis as used by [23, 33].
The geometrical parameters and physical properties used in this paper correspond to a typical laboratory Rijke tube as defined in the previous works [23, 33]. Specifically, m, m, m2, m, m/s, , kg/m3, m/s, W/(m·K), m2/s, and m2/s. Correspondingly, the heat release coefficient is set in the range from 0 to 2. The time delay is selected from 0.01 to 0.04, which is slightly smaller than the value 0.05 found in [16]. However, this has minor effect on the results since the time delay is very much smaller compared with the period over which the transient growth occurs (which has order of 2–20) [23]. In the paper, to depict a realistic fame shape, a parabolic distribution is defined by specifying different in terms of different heat sources, which can be practically implemented by varying the power input of hot wire gauzes as well as the time delay by changing the speed of base flow [34].
3. Numerical Method
To investigate the variety of nonlinear behaviors displayed by the Rijke tube, the partial differential equations (6) should be reduced into a series of ordinary differential equations (ODEs) to form a dynamic system for bifurcation study. Among the past studies, Galerkin method [35] has been most commonly used, for which two sets of presumed basis functions for and are employed to convert (6) to ODEs for each acoustic mode. Nevertheless, in this paper, the Method of Lines (MOL) approach was utilized for the various advantages it provides as exemplified in [36]. The basic idea of the MOL is to directly discretize the spatial independent variable in the PDEs but keep the temporal variable continuous. Hence, one obtains a series of ODEs naturally without introducing any approximations. It is worth mentioning that this approach has seldom been used in the past for such applications due to the complexity of the discretization and expensive computational costs.
Here, a standard finite difference scheme was applied for the discretization of governing equations (6). The domain along the length of the tube was divided into points indexed as , which split the direction into parts with an identical interval of . Second-order central difference scheme was used to discretize the spatial derivatives and . For the treatment of boundary conditions, second-order upwind difference scheme was employed to deal with the Neumann boundary conditions for at and , while specific values were given for the Dirichlet boundary conditions of . Consequently, we could set up the dynamic system written symbolically as where is the argument vector of size , is the vector function of size , and is the bifurcation parameter vector .
A numerical continuation method [28, 29] was utilized to investigate the bifurcation analysis of the time-delayed thermoacoustic system. First, the system was linearized around the preresolved equilibrium (steady state) solution and the corresponding eigenvalues of the linearized system, namely, the roots of the characteristic equation, were calculated subsequently. These roots were first approximated by a linear multistep method (LMS-method) and then corrected using a Newton iteration method [28]. Mathematically, the rightmost characteristic root, or the characteristic root with the maximal real part, conclusively indicates the linear stability of the system. With the variation of one or more bifurcation parameters, bifurcation may occur as the rightmost characteristic root crosses zero, as shown in Figure 2(a). Moreover, all the branches emitting from the bifurcation points were captured by a prediction-correction approach [37, 38]. Specifically, the stability of each branch was determined by the Floquet-multiplier scheme, for which the solution is stable provided the moduli of all the multipliers are less than unity, as demonstrated in Figure 2(b).
(a)
(b)
In the above section, the governing equations have been discretized into a set of ODEs using the finite difference method. As is well known, the discretization resolution is of significant importance on the convergence and accuracy of the numerical solutions. Before conducting further analysis, a grid independence check was primarily done to assure the accuracy. From Figure 3, it could be seen that there is tiny change for the bifurcation diagrams when the number of discretization points increases to and upwards. To compromise between the accuracy and speed, the domain was thus discretized into 120 points for the following bifurcation analysis.
4. Results and Discussions
In this section, effects of multiple heat sources, time delay, and heat release distribution are investigated in detail to reveal the extensive nonlinear characteristics involved in the case of a distributed heat source.
The numerical approach stated above was first validated through comparison with the experimental data [18, 39] and numerical calculation [22] in the case of a single heat source. This curve was obtained by varying the heat release coefficient and time delay simultaneously. Each point on the curve stands for a Hopf bifurcation point. When the heat release coefficient is larger than the values of these points, an infinitesimal perturbation would result in a large-amplitude limit cycle, defined globally to be unstable. Therefore, the Hopf bifurcation curve denotes the linear stability boundary of the system. Any point within the inner region of the curve is unstable, which should be avoided in practical operation. From Figure 4, it can be concluded that the current approach yields very good prediction of the stability boundary and is more accurate than the results provided in [22] using the Galerkin method.
The effect of varying the heat release parameter on the dynamic behavior of the system was presented by the bifurcation diagram shown in Figure 5.
It is a typical subcritical Hopf bifurcation and similar to the bifurcation diagrams reported in [21, 23, 33]. Hopf bifurcation and fold bifurcation can both be observed in this diagram, with the heat release parameter being and , respectively. The stability of each periodic solution was determined by Floquet multipliers, as discussed in above section. The minimum acoustic energy on periodic solutions was selected to measure the size of acoustic perturbations. The equilibrium solution was stable for . When increased to , Hopf bifurcation arose. The system lost stability and small amplitude periodic solutions (limit cycle) emerged from the point. These limit cycles were unstable until reached , which stabilized the system. Hence, for , the steady state solution was stable for perturbations of any magnitude. For , an infinitesimal perturbation would result in a large-amplitude limit cycle. In the region of , linearly stable equilibrium solutions, small-amplitude unstable limit cycles, and large-amplitude stable limit cycles coexist to form a bi-stable region. Hence, the fold bifurcation and the Hopf bifurcation can indicate the linear and nonlinear stable region of a system.
It is well known that the generation of sound in a Rijke tube is a result of the standing acoustic wave being set up in the tube [40]. The heat source plays the role of exciting the duct mode acoustic waves as well as reinforcing and sustaining the excited waves. According to the Rayleigh criterion [30], “when the pressure oscillation is in phase with the unsteady heat release, thermal energy could be transferred to acoustic energy and the acoustic oscillations can be further strengthened.”
With the result from a single heat source, the influence of multiple heat sources was investigated as shown in Figure 6. It can be seen that there are striking differences regarding and between distributed heat sources and concentrated heat sources (refer to the different positions along the tube). In the presence of two distributed heat sources in Figure 6(a), both and have increased, which means the system has become more stable. This is due to the fact that the acoustic perturbation excited by the foregoing heat source could interact with the subsequent unsteady heat release in the case of distributed heat sources. As long as they are (partly) in phase, the heat would be added during compression and removed during expansion. Consequently, the acoustic perturbation could be further strengthened. Otherwise, it would be damped. In other words, as the integrated heat release acts as a source of acoustic energy, the thermoacoustic system becomes less stable, and vice versa. Furthermore, the enhancement of and would be different for different distributed positions. Similar results have also been found for three and four heat sources in Figures 6(b) and 6(c). The results presented could adequately illustrate the importance of multiple heat sources in the system stability. The presence of distributed heat source would influence the stability of system vis-a-vis the compact heat source. Besides, the system stability also depends on the distribution length. It has different influence on the nonlinear dynamics via locating the distributed heat sources at different positions along the tube. By choosing proper values for the heat release rates, this model could be employed for multiple heat sources distributed anywhere along the tube.
(a) Two heat sources
(b) Three heat sources
(c) Four heat sources
The influence of time delay was examined as displayed in Figure 7. Subcritical Hopf bifurcations were also observed with different time delays. When the time delay remains the same at one heat source, the increase at the other heat source would weaken the system stability (see , 0.02 and , 0.04), whereas the decrease would stabilize the system significantly (see , 0.02 and , 0.02). Moreover, the proportional increase or decrease of both time delays would make the system less or more stable, respectively (see , 0.02 and , 0.03). Time delay is the result of thermal inertia between the heat transfer and the velocity [40]. When the velocity varies, the corresponding heat release cannot change instantaneously. A smaller time delay means that it takes shorter time for the variation of velocity to get reflected by heat release. Thus, the perturbation would be diminished more quickly and the system becomes more stable.
Furthermore, to depict a realistic flame shape, a qualitatively suggestive parabolic-type distribution is considered (see Figure 8). The heat sources were placed at , , , , and separately. Different heat release coefficient describing parabolic curve was used to represent different intensity/source strength of discrete heat sources.
Figure 9 reveals that the Hopf bifurcation at was linearly decreasing as increased. It was due to the fact that the parabolic distribution reinforced the coupling between integrated heat release and the acoustic field and made the system less stable. In this case, the system will be more susceptible to triggering.
5. Conclusions
In this study, the dynamic behavior of a horizontal Rijke tube with a distribute heat source is studied in detail. The distributed heat source comprises a row of discrete heat sources and a modified form of Heckl’s model is adopted to model the integrated heat release rate. Unlike the frequently used Galerkin method, Method of Lines (MOL) technique is employed to discretize the governing equations. A numerical continuation method is used for the study of thermoacoustic instability within the context of bifurcation and nonlinear analysis. It is found that distributed heat source plays an important role in determining the stability of a thermoacoustic system. The results reveal if the heat source is of a distributed nature as opposed to being compact, it would influence the stability of system. By locating the distributed heat sources at different positions along the tube, the nonlinear characteristics of system would also be affected differently. It could either strengthen or weaken the acoustic energy based on specified distribution. Besides, the decreased (increased) time delay is proven adequately to have a stabilizing (destabilizing) effect on the system. Moreover, the distribution of heat release parameter also has an effect on the system stability. A parabolic-type distribution could diminish the stable region of the system remarkably.
Due to the lack of experimental/numerical data in the open literature regarding the presence of a distributed heat source, the model could not be verified quantitatively. However, the results are shown to be physically reasonable and sound and could provide the basis for comparable data involving future experiments/predictions. In the follow-on work, a two-dimensional flame model will be implemented and proper heat release function will be identified to depict the combustion in practical systems.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
X. Yang acknowledges gratefully the support of President’s Doctoral Scholar Award from the University of Manchester. S. Lei would like to gratefully acknowledge the financial support by Industrial Development Agency (IDA) Ireland.