Numerical Study on the Route of Flame-Induced Thermoacoustic Instability in a Rijke Burner

The self-excited thermoacoustic instability in a two-dimensional Rijke-type burner with a center-stabilized premixed methane–air flame is numerically studied. The simulation considers the reacting flow, flame dynamics, and radiation model to investigate the important physical processes. A finite volume-based approach is used to simulate reacting flows under both laminar and turbulent flow conditions. Chemical reaction modeling is conducted via the finite-rate/eddy dissipation model with one-step reaction mechanisms, and the radiation heat flux and turbulent flow characteristics are determined by using the P-1 model and the standard k-ε model, respectively. The steady-state reacting flow is first simulated for model verification. Then, the dynamic pressure, velocity, and reaction heat evolutions are determined to show the onset and growth rate of self-excited instability in the burner. Using the fast Fourier transform (FFT) method, the frequency of the limit cycle oscillation is obtained, which agrees well with the theoretical prediction. The dynamic pressure and velocity along the tube axis provide the acoustic oscillation mode and amplitude, also agreeing well with the prediction. Finally, the unsteady flow field at different times in a limit cycle shows that flame-induced vortices occur inside the combustor, and the temperature distribution indicates that the back-and-forth velocity changes in the tube vary the distance between the flame and honeycomb in turn, forming a forward feedback loop in the tube. The results reveal the route of flame-induced thermoacoustic instability in the Rijke-type burner and indicate periodical vortex formation and breakdown in the Rijke burner, which should be considered turbulent flow under thermoacoustic instability.


Introduction
In 1859, Rijke discovered strong acoustic oscillations when heat was added to the lower half of a vertical tube opened at both ends [1]. Several similar devices were proposed to study the heat-driven acoustic phenomenon in past decades, and these devices were grouped under the general heading "Rijke tube" [2]. Among these devices, one type in which the flow is combustible gas and the heat is produced by the flame is called a Rijke burner or a Rijke tube burner.
Experimental research on the Rijke burner is intuitive and can yield highly reliable information about thermoacoustic instability, such as temperature, pressure, and velocity. Abundant nonlinear behaviors, such as limit cycle oscillation, quasiperiodic, frequencylocked, and chaotic behaviors, were discovered [3][4][5]. Based on experimental results, some classical flame models have been proposed [6][7][8][9][10]. However, experimental methods can collect data only at finite points. The parameters of the entire field cannot be obtained by experimental methods, and a comprehensive understanding of the flow field of thermoacoustic instability is still lacking.
Analytical studies have been conducted to explain the mechanisms of thermoacoustic instability since the 1960s [11,12]. Culick established the linear theory of combustion instabilities and first derived the nonlinear acoustic model in combustion chambers based on experimental research [13][14][15][16]. Balasubramanian and Sujith investigated the role of nonnormality and nonlinearity in flame-acoustic interactions and found that nonnormality and nonlinearity lead to triggering and redistribution of energy between eigenmodes [17]. Yoon et al. predicted the bootstrapping instability mechanism observed in rocket motors by approximated modal analysis [18,19]. Li et al. observed the stability switch phenomenon with a one-dimensional (1-D) Rijke-type thermoacoustic system [20]. Through the analytical method, we have obtained a thoughtful understanding of the mechanism of thermoacoustic instability, especially for predicting the frequency, amplitude, and stability of self-excited thermoacoustic instability.
With the development of computational technology, computational fluid dynamics (CFD) simulations have increasingly helped to describe the mechanisms of thermoacoustic instability and to provide guidance for experimental research, making up for the shortage of experimental work.
Many numerical studies have been conducted on combustion instabilities. Because of the nonlinearity in the combustion system, the interactions between reaction mechanisms, hydrodynamics, and acoustics reveal different spatial and temporal scales and are sensitive to various parameters. Chatterjee investigated the occurrence of combustion instabilities in a Rijke tube-type combustor using a two-dimensional (2-D) finite volume method, captured the instability growth to a limit cycle of the pressure oscillation, and predicted the frequency and magnitude of the thermoacoustic instability [21]. Chatterjee et al. simulated a three-dimensional (3-D) premixed chamber with ANSYS FLUENT software and studied the dependence of a homogeneous mixture on various parameters such as the premixing length, airflow rate, and equivalence ratio to optimize the combustor [22]. Wang also successfully simulated the lean premixed combustion instabilities in a 2-D Rijke-type burner by CFD simulation and found that the radiation and reaction mechanism play critical roles in predicting the occurrence of thermoacoustic instability [23]. Blanchard et al. simulated a bluff-body flame holder, and the proper orthogonal decomposition modes of the simulated results and experiments showed quantitative agreement [24]. Ghulam et al. used COMSOL software to investigate the driving mechanism of combustion instabilities in a small rectangular combustion duct and obtained the pressure waves in the longitudinal and transverse directions [25]. On the research of the interaction of turbulence and combustion, numerical simulation has also shown its extraordinary ability [26]. These CFD simulations help in understanding the detailed flow parameters, such as pollution control, efficiency enhancement, and mass transport processes, for further analysis and can thus serve as an aid to engineers and designers for building combustion systems with a decreased possibility of thermoacoustic instabilities.
The present study concentrates on the numerical simulation of a 2-D Rijke burner to understand the dynamics of self-excited thermoacoustic instability. A finite volumebased approach is used to simulate reacting flows under both laminar and turbulent flow conditions. The onset of growth of the pressure, velocity, and reaction heat are captured. The oscillation frequency and mode are calculated and compared with the theoretical prediction, and they show good agreement. Finally, the unsteady flow fields at different times show that flame-induced vortices occur inside the combustor, and the temperature distributions show that the velocity perturbation in the tube changes the distance between the flame and honeycomb in turn, forming a forward feedback loop in the tube. The results of this work show that there are vortices in the Rijke-type burner when thermoacoustic instability occurs, and it should be considered a turbulent reacting flow field.

Rijke Burner Configuration
The Rijke burner used in this study consists of a 1572 mm-long steel tube with an inner diameter of 73.7 mm, as shown in Figure 1. A premixed methane-air mixture is fed into the inlet end of the tube, and a honeycomb flame stabilizer with a thickness of 25.4 mm is positioned in the middle of the tube. A flat flame is anchored at the top of the honeycomb flame stabilizer, and releases unsteady heat into the flow, inducing pressure oscillation in the tube. A 2-D axisymmetric flow field is used to approximate the 3-D flow field to simplify the computational model as shown in Figure 2.  The model of the Rijke burner is generated, and boundary conditions are set, including the mass flow rate inlet, pressure outlet, coupled wall between the honeycomb and flow, coupled wall between the inside tube wall and flow, outer tube walls, and axisymmetric centerline. The input parameters used in setting up the Rijke burner mode are as follows. The incoming flow is set to 1.27 × 10 −4 kg/s with a fuel/air equivalence ratio of 0.75. The equivalence ratio is defined as the value of the actual fuel/air mass ratio relative to stoichiometric fuel/air ratio. If it is less than 1, the combustion is lean with excess air. The convective heat transfer coefficient for the outside tube wall is estimated to be 20 W/(m 2 ·K), and the reference condition is set as normal temperature and pressure. The boundary conditions and initial conditions are summarized in Tables A1 and A2 of Appendix A, respectively.

Mathematical Model
The governing equations of this model include conservation equations of mass, kinetics, species, and energy, combined with turbulence equations and radiation equations. The conservation equations are expressed by where , , and are the density, time, and velocity, respectively; represents the dependent transport variables [27]; represents the diffusion coefficient; and is the source term of the variable .
The governing equations are solved with the finite volume method. Although at the steady state of the combustion field the average velocity is approximately 0.026 m/s, the corresponding Reynolds number is approximately 120, and the Mach number is 1.5 × 10 −4 in this model, indicating that this is a laminar flow for the non-reacting condition. However, influenced by the highly nonlinear chemical reaction, it is inappropriate to ignore the turbulence effect in this model. The k-ε model [28] is the most common model used in CFD to simulate mean flow characteristics for turbulent flow conditions. Here, the k-ε model with enhanced wall treatment is used to simulate the reacting flow.
The species is determined by the species transport model, while the finite-rate/eddy dissipation model is applied to compute the interaction between turbulence and chemistry. Besides, the material properties used are listed in Table A3 of Appendix A.
The methane/air reaction is complicated and involves dozens of reactions and hundreds of elementary reactions. Several reaction mechanisms have been proposed to accurately describe the combustion process. The detailed reaction mechanisms are highly precise [29], while the computational costs are increased. In a numerical comparison between the different methane/air combustion mechanisms, Li [31] pointed out that the one-step mechanism agrees with experimental results quantitatively for fuel-lean conditions. Acampora et al. [31] concluded that moving from a one-step global mechanism to a twostep mechanism does not significantly increase the accuracy of predictions. According to the mentioned points, the one-step reaction mechanism is used for the reaction of methane and air, and is expressed as The rate of fuel consumption is given by the global Arrhenius equation: where the parameters A, ⁄ , m and n are given in Table A4 of Appendix A [32,33].
Due to the high temperatures of the flame, radiative heat transfer becomes one of the dominant mechanisms of heat transfer in addition to conduction and convection heat transfer [34,35]. Hence, proper simulation of the radiative heat transfer is of great importance to this problem. The P-1 model is based on the expansion of the radiation intensity I into an orthogonal series of spherical harmonics and has high accuracy and the ability to consider the radiation exchange in the flame environment [36][37][38][39] For this advantage, the P-1 model is widely used in the simulation of combustion [40,41]. In this work, the P-1 radiation model was chosen to model the high-temperature flow problem.

Grid Independence and Model Verification
To verify the computational model and grid independence, the numerical results of the axial temperature distribution at the centerline under steady-state combustion are compared. Three different numbers of nodes, namely, 26,312, 67,056, and 195,233 nodes, are compared with the previous numerical modeling of Wang [23], and the results are presented in Figure 3. The convergence limit is assumed when the residual values are less than 10 −6 .
As shown in Figure 3, the temperature distributions before the flame stabilizer almost converged to the same value, while in the flame stabilizer area and the flame zone, the temperature increased sharply because of the heat released by the flame. The result of the 67,056-node scheme is closer to the result presented by Wang [23] than the other two grid schemes. The peak temperatures of the computational grids with 26,312, 67,056, and 195,233 nodes were approximately 1804 K, 1847 K, and 1908 K, respectively, and the outlet temperatures were 408 K, 331 K, and 354 K, while the peak temperature and outlet temperature in Wang's simulation were 1927 K and 320 K, respectively. In addition, for the outlet temperature, the result of the 67,056-node scheme is much closer to the result of Wang. Consequently, the mesh with 67,056 nodes was selected as the best computation grid for computational modeling.

Steady-State Solution
The modeling of thermoacoustic instability proceeded in two stages-the first stage was the simulation of the steady-state reacting flow, and the second stage was the simulation of the transient reacting flow. For the steady-steady reaction simulation, the flow field was initialized with the inlet parameters. The solution of the steady-state reacting flow included the accurate capture of the flame anchoring process shown in Figures 4 and 5, the consumption of reactants and product formation shown in Figure 6, and the steadystate temperature distribution along the combustor shown in Figure 3. Each process/distribution is important for verification of the computational results against the theoretical and experimental results. The reaction rate and temperature contour plots downstream of the flame stabilizer are shown in Figures 4 and 5, in which the location of the highest reaction rate and the temperature indicate the flame location. As seen in the figures, the flame sits just on the outlet of the flame stabilizer at a distance of 1-2 mm away from it. The temperature contours show that the incoming flow was preheated to approximately 800 K at the outlet of the honeycomb and then heated to a peak temperature of 1828 K by the reaction. Then, a steep temperature gradient, which was dominated by radiation heat transfer, occurred downstream of the peak temperature point. With decreasing temperature, the influence of radiation weakens, and the temperature gradient becomes increasingly gentle until the flow reaches the outlet.
The mass fractions of reactants and products along the centerline in the flame stabilizer area and the flame zone are shown in Figure 6. Due to the preheating effect of the flame stabilizer and flame radiation, the reaction started before the flame stabilizer outlet and rapidly increased to a peak downstream of it. At the flame location, the reaction rate reached its maximum value; the CH4 fuel was completely consumed, and its mass fraction became zero. The steady-state species distribution and kinetic reaction rate show the exact location of the flame and demonstrate the effectiveness of the numerical simulation in this work.

Self-Excited Pressure Oscillation-The Exponential Growth of Thermoacoustic Instability
To capture the dynamic characteristics of thermoacoustic insatiability, a transient simulation was conducted. Here, the steady-state reacting flow field was set as the initial condition of the dynamic flow. The time step was set to between 1 × 10 s and 1 × 10 s to obtain enough information in the time domain while meeting the convergent requirements in previous research [21,23,42]; here, it was set to 1 × 10 s. The transient solution includes the growth of the self-excited instability, the subsequent limit cycle behavior of pressure and velocity, the oscillating frequencies and modes, and the periodic behavior of the flow field.
The instability growth mechanism is characterized by increasing oscillations of pressure inside the combustor. Figure 7a shows the time evolutions of average pressure and average velocity 5 cm downstream of the flame stabilizer. Figure 7b shows the growth of total reaction heat released by the flame. The steady-state value of pressure, velocity, and heat release were 0 Pa, 0.1226 m/s, and 63.68 W, respectively. Approximately 0.05 s after the time integration was initialized, the pressure, velocity, and reaction heat started to increase exponentially. The oscillation grew until approximately 0.86 s and reached maximum amplitudes. After several periods of amplitude adjustment from 0.86 s to 1.1 s, the system reached a "saturated" state, which is referred to as the limit cycle. The pressure, velocity, and reaction heat maintained at amplitudes of approximately 208 Pa, 0.39 m/s, and 10.79 W, respectively.  Figure 7c illustrated the phase relationship between the average pressure and velocity 5 cm downstream of the flame stabilizer at the limit cycle oscillation state. It can be seen that the velocity had a 88.11° phase advance relative to the pressure. Figure 7d illustrated the phase relationship between the average pressure 5 cm downstream of the flame stabilizer and total heat release of the flame at the limit cycle oscillation state. The reaction heat release had a 24.3° phase delay relative to the pressure, which fulfills the Rayleigh criterion that oscillations are encouraged when heat fluctuates in phase with pressure perturbation [43]. It can be calculated that the reaction heat had a 114.34° phase delay relative to the velocity. This time delay between reaction heat and velocity has been proved to play a critical role in the stability of the Rijke-type thermoacoustic system [20].
Two important parameters to evaluate the growth of thermoacoustic instability are the growth rate of the acoustic disturbances α and G. The definition and detail derivation of α and G can be found in the reference [44]. The time evolution of the pressure and heat release of the flame from 0 s to 0.4 s are illustrated in Figures 8 and 9. It can be seen that both the pressure and heat release were growing exponentially after an adjustment of about 0.05 s. Two exponential functions are constructed to best fit the series of local maximum and minimum points in Figures 8 and 9, and they are written as: where is the series of local maximum and minimum pressure, is the series of local maximum and minimum reaction heat, A, B, C, D, , α , and are parameters to be solved.
The Equations (4) and (5) are solved by iteration, and the calculated growth rate is 13.77 s , is 13.85 s . The corresponding G is 1.22 × 10 rad/s , which is within the range of the maximum possible growth rate for this type of combustion-driven acoustic disturbance, =9.2 × 10 rad/s .  From the analysis of Yoon et al. [18,19], the analytic pressure modes and velocity models of the acoustic perturbation for the acoustically closed-open Rijke-type burner can be expressed by Equations (6) and (7): = cos ((2 − 1) π /2 L) = 1,2,3 ⋯ ∞ = sin ((2 − 1) π /2 L) = 1,2,3 ⋯ ∞ and the corresponding acoustic frequencies can be calculated by: where is the n-th acoustic frequency, and a is the sound speed. According to the Rayleigh criterion, oscillations are encouraged when heat fluctuates in phase with pressure perturbation [43]. It can be determined from Equation (6) that the second acoustic mode may be excited when the unsteady heat is located at the center of the tube length because the flame is near the antinode of the second acoustic mode, which has a three-quarter wave mode shape. The frequency of this mode is 168.8 Hz with the inlet parameter by solving Equation (8).
The frequency component of the pressure oscillation from 1.4 s to 1.6 s is calculated and displayed in Figure 10. The excited frequency is 179 Hz, which is 6% higher than the predicted acoustic frequency of 168.8 Hz for the non-reacting flow. The small discrepancy is due to the flow in the tube being heated by the flame, which induced an increase in the sound speed. Furthermore, the increased sound speed resulted in an increase in the oscillation frequency. Thus, researchers should take the frequency increase effect into consideration when predicting the thermoacoustic instability frequency with the inlet parameters. To compare the simulated oscillating mode with the prediction, the dynamic pressure along the centerline at time 1.400445 s was normalized by the maximum pressure amplitude, meanwhile the velocity along the center line at time 1.40585 s was normalized by the maximum velocity amplitude. In Figure 11, the normalized pressure and axial velocity of this work are compared with the analytical second acoustic pressure mode and velocity mode, and the normalized pressure and axial velocity in reference [23]. The simulation pressure distribution along the centerline is in good agreement with the analytical results and the results of Wang. While for the velocity distribution, two areas appear to have significant differences between the numerical results and the analytic results. The first area is in the flame stabilizer, where the simulated velocity is two times the analytical velocity because the flame stabilizer occupied half of the total volume of the entire flame stabilizer zone. The other inconsistent area is on the outlet face, because it was set as a pressure outlet, and the mass flow in the tube must be conserved. To further study this problem, it will be better to impose constraints on the velocity derivative on the outlet boundary. Aside from these two areas, the simulated velocity distribution agrees with the prediction results. Figure 11. Comparison of the simulated pressure and velocity along the centerline with the analytical prediction results. Figures 12 and 13 show the instantaneous axial pressure distribution and axial velocity along the centerline of the Rijke burner at 8 different times of a limit cycle: namely, at times T, T/8, T/4, 3T/8, T/2, 5T/8, 3T/4, and 7T/8. In Figure 12, it can be seen that the pressure node point was at 1/3 length, and the pressure antinode point occurred downstream of the flame stabilizer at approximately 2/3 length, which matches the theoretical pressure maximum location for a three-quarter wave shown in Figure 11. From Figure 13, one can see that the velocity node point occurred slightly before the 2/3 tube length, which is also the pressure summit point of the pressure distribution.   Figure 14 illustrates the axial velocity distribution and streamline at different times within a period of one limit cycle, where the yellow regions represent velocities above zero, the blue regions represent velocities below zero, and the gray regions represent the solid zone. Once the system is excited into the periodical oscillating state, the flow in the tube does not remain steady and uniform; rather, it begins oscillating back and forth. At the beginning of a period, the velocity was positive in the left 2/3 of the tube length and negative in the right 1/3 of the tube length, as shown in Figure 14a. At approximately 2/3 of the tube length, which is a node point of the velocity mode in the 1-D analysis, the axis velocity was zero. This result is in good agreement with the predicted velocity mode, and the velocity along the axis is in agreement with the velocity at time T in Figure 13. In the flame zone, due to the influence of the high-temperature flame, the flow was compressed and distorted, inducing the formation of a clockwise vortex after the flame area, as shown in Figure 14c. With the growth of this clockwise vortex, a counterclockwise vortex was formed at the inner wall of the left half of the tube in Figure 14d. These vortices were enlarged and changed the direction of flow in Figure 14e, making the flow negative on the left 2/3 of the tube length and positive on the right of the 1/3 tube length at time T+T/2, as shown in Figure 14f.

Unsteady Flow Field in the Limit Cycle State
After time T+T/2, the flow was compressed near the inlet face and was expanded in the flame zone. When the negative flow was heated by the flame, followed by the formation of new vortices, and the subsequent development and collapse of these new vortices as shown in Figure 14g-j, the direction of flow changed in the opposite direction. A new period began.
The flow field variation in a limit cycle shown in Figure 14 demonstrates that vortices emerged and collapsed when thermoacoustic instability occurred. Thus, the computation of thermoacoustic instability should consider a turbulent flow, even in low Reynolds number cases. Thus, it can be concluded that the distance between the flame and flame stabilizer varies with the change in the axial velocity in the limit cycle state; the flow between the flame stabilizer and the 2/3 tube length is compressed when the axial velocity is positive, while it is expanded when the axial velocity is negative. During this periodic compression and expansion, the streamline is stretched and distorted, leading to the periodic formation and disappearance of vortices in the tube. This periodic behavior influences both the shape and distance of the flame via feedback; as a result, the forward feedback between the heat release and flow perturbation results in thermoacoustic instability.

Conclusions
A two-dimensional CFD simulation of a Rijke burner was conducted in this work. The self-excited thermoacoustic instability was numerically captured and analyzed, and the following conclusions are drawn.
1. The onset and growth of the thermoacoustic instability were successfully captured.
The simulated growth of the instability G was 1.22 × 10 rad/s , which is less than the maximum possible growth rate. 2. The simulated oscillation frequency and mode were in good agreement with the 1-D analytical results. The simulated pressure oscillation frequency was 6% higher than the analytical results based on the inlet parameters, since the flow temperature and the sound speed were increased by the reaction heat. The normalized pressure and normalized velocity along the center line were in agreement with the prediction results of the second acoustic mode. 3. As the system was excited into a limit cycle oscillation state, the amplitude of the velocity perturbation was larger than the average steady-state flow; thus, in some local areas, the flow velocity may be below zero, resulting in local backward flow in the burner. 4. In the limit cycle state, the distance between the flame and honeycomb outlet varied with the periodic change in the axial velocity. The flow between the flame stabilizer and the 2/3 tube length was periodically compressed and expanded, while the streamline was periodically stretched and distorted, inducing the periodic formation and disappearance of vortices in the tube. This is a typical behavior of flame-induced vortices, which leads to the flow field transition from laminar flow to turbulent flow. 5. The periodic behavior of the flow field induced by the unsteady heat release influenced both the shape and distance of the flame via feedback; thus, the forward interaction between the heat release and flow perturbation resulted in thermoacoustic instability. Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare that they have no conflict of interest.