Formation mechanism of high-frequency combustion oscillations in a model rocket engine combustor

Combustion oscillations are manifested as organized combustion oscillations with large-amplitude and high-frequency and arise from the strong mutual interaction between unsteady heat release rate and acoustic pressure fluctuations in a combustor. The onset of combustion oscillations causes unacceptable structural damage and shortens the life span of combustors, obstructing the development of combustors involving land-based gas-turbine power plants, aircraft engines, and rocket engines. Many physical mechanisms responsible for the excitation and sustainment of combustion oscillations in various types of turbulent combustor have been examined in detail and are summarized in several books^1,2 and comprehensive review articles.^3,4

Advanced analytical methods inspired by dynamical systems theory and complex networks have recently shed light on the study of combustion dynamics in thermoacoustic systems^5–14 and have proactively been employed with the two main aims of (i) gaining a physical understanding of the transition process from the stable state to combustion oscillations and (ii) developing potential detectors for capturing the precursor of combustion oscillations. Murugesan and Sujith⁵ studied the dynamic behavior of combustion noise and combustion oscillations in a bluff-body-type turbulent combustor using the visibility graph.¹⁵ They showed that a scale-free structure related to fractality appears in the degree distribution in the graph during combustion noise, while the scale-free structure collapses during combustion oscillations. Okuno et al.⁶ reported the possible presence of high-dimensionality and small-world-like nature in the pressure fluctuations during combustion oscillations in a swirl-stabilized turbulent combustor using cycle networks¹⁶ and phase space networks.¹⁷ Murayama et al.⁹ clarified the network structure in the flow velocity field during combustion oscillations in the same combustor⁶ using the turbulent network.¹⁸ The possible presence of a primary hub in the hydrodynamic shear layer between the inner vortex breakdown region and the outer recirculation region, which is an important factor in the formation of combustion oscillations, can be extracted from the spatial distribution of the vertex strength and the community structure. Mondal et al.⁸ studied the synchronized state between the pressure and heat release rate fluctuations in a bluff-body-type combustor by estimating the order parameter¹⁹ in terms of collective synchronization. Combustion noise is a state of phase-asynchronous aperiodic oscillations, while combustion oscillations form a state of phase-synchronous periodic oscillations. They also showed that the mean and variance of the order parameter significantly increase as the combustion state undergoes a major transition from combustion noise to combustion oscillations. Guan et al.¹² studied the synchronized state between the pressure and heat release rate fluctuations in a liquid-fueled diffusion-flame combustor by estimating the probability of recurrence.²⁰ Both fluctuations are desynchronized during combustion noise, while they undergo a weak phase synchronization during combustion oscillations with a weak limit-cycle state. The directional coupling between the pressure and heat release rate fluctuations in a bluff-body-type combustor has recently been studied by Godavarthi et al.¹¹ during a transition and the subsequent combustion oscillations. They showed that the strong directional dependence of the heat release rate fluctuations on the pressure fluctuations during an amplification of combustion oscillations can be explained by using a methodology combining recurrence plots²⁰ and a recurrence network.²¹ 

Nguyen and co-workers^22–24 have recently carried out numerical studies on the triggering mechanism of combustion oscillations in a model liquid rocket engine combustor with a single injector. We have more recently studied the dynamic behavior of high-frequency combustion oscillations in a cylindrical combustor with an off-center-installed coaxial injector.²⁵ The appearance and disappearance of the scale-free structure are found to be formed in the turbulence network during combustion oscillations. Both the order parameter and the phase parameter quantifying the degree of phase synchronization between the pressure and heat release rate fluctuations in the combustor are useful for detecting a precursor of combustion oscillations. The switching of the directional dependence between the pressure and heat release rate fluctuations plays an important role in the transition to combustion oscillations, which is clearly identified by the transfer entropy²⁶ in the framework of information theory. Our previous study²⁵ clearly revealed the applicability of the turbulence network, the order parameter, and the transfer entropy in dealing with the dynamic behavior of combustion oscillations. However, an in-depth elucidation of the physical mechanism behind the formation and sustainment of combustion oscillations has not been provided in our previous study, focusing particularly on the relevance of the pressure and flow velocity fluctuations in fuel/oxidizer injectors to the heat release rate fluctuations in combustors. The previous numerical studies^22–24 have not clarified the directional coupling during combustion oscillations in a model liquid rocket engine combustor on the basis of dynamical systems and information theories. Note that Kasthuri et al.²⁷ have experimentally studied the dynamical transitions from a stable combustion to combustion oscillations in a liquid rocket combustor using the recurrence quantification analysis based on dynamical systems theory.

Our main purpose in this study is to reveal the physical mechanism behind the formation and sustainment of high-frequency combustion oscillations in a cylindrical combustor with an off-center-installed coaxial injector using sophisticated analytical methods inspired by symbolic dynamics and complex networks. In particular, we attempt to clarify the directional feedback process between the flow velocity fluctuations in the fuel/oxidizer injectors, the pressure fluctuations, and the heat release rate fluctuations in the combustor during a transition and the subsequent combustion oscillations. The symbolic transfer entropy²⁸ as an extended version of the transfer entropy enables the quantification of the dominant information flow between two variables. In this study, we estimate the symbolic transfer entropy of the flow velocity fluctuations in the fuel/oxygen injectors and the heat release rate fluctuations in the combustor. The symbolic recurrence plots^20,29 enable us to examine the similarity of the dynamic behaviors of two variables. In addition to the symbolic transfer entropy, we adopt symbolic recurrence plots to clarify the mutual synchronization between the pressure and heat release rate fluctuations during the transition to combustion oscillations. We finally obtain a deeper understanding of the physical mechanism behind the formation of combustion oscillations by introducing an unweighted and directed spatial network consisting of the pressure and heat release rate fluctuations in the combustor.

This paper is organized into four sections. A brief description of the numerical computation and the framework of time series analysis is given in Sec. II. We present the results and discussion in Sec. III. A summary is given in Sec. IV.

In this study, we employ the analytical methods based on symbolic dynamics and complex networks for the spatiotemporal data³⁰ during a transition and the subsequent high-frequency combustion oscillations previously obtained by high-resolution large-eddy simulation (see Ref. 25 for details). This large-eddy simulation produces the first tangential (1T)-mode oscillations exceeding ±3% with respect to the average pressure in the combustor [Fig. 1(a)], with the flame dynamics switching back and forth between the attached and lifted flame edge from the injector rim.³⁰ Large-scale vortex rings^25,30 are produced from the injector rim during combustion oscillations [Fig. 1(b)].

The transfer entropy²⁶ is a useful measure for investigating the directional information flow between two variables. A point of care in the estimation of the transfer entropy is the construction of the distribution in the joint and conditional probabilities of the variables. Various methods of approximating the probability distribution have been proposed,^31,32 but include the problem of setting many parameters used as thresholds. Bandt and Pompe³³ proposed the permutation entropy based on a symbolic dynamics approach, which considers the probability distribution of rank-order patterns in a time series. The symbolic transfer entropy²⁸ based on Bandt and Pompe's concept does not need the parameters to be set. Note that the importance of the rank-order patterns has been shown by Amigó.³⁴ In this study, we estimate the directionality index $ΔSS$⁠, defined as the difference between the symbolic transfer entropies, in three cases: (i) spatial average cross-sectional fuel flow velocity fluctuations in the injector $w¯I,F$ and heat release rate fluctuations $q′$⁠, (ii) spatial average cross-sectional oxidizer flow velocity fluctuations in the injector $w¯I,O$ and $q′$⁠, and (iii) pressure fluctuations in the combustor $p′c$ and $q′$ during a transition and the subsequent combustion oscillations,

Here, I denotes the injector, F denotes the fuel, O denotes the oxidizer, and $SS,w¯I,F→q′$ is the symbolic transfer entropy for the direction of $w¯I,F$ to $q′$⁠. $πw¯I,FD(tm)$ is the rank-order pattern of $w¯I,F$ and $πq′D(tm)$ is the rank-order pattern of $q′$ at time t_m, where D is the embedding dimension. $P(πq′D(tm+1+τa),πq′D(tm+τa),πw¯I,FD(tm))$ is the joint probability of $πq′D(tm+1+τa),πq′D(tm+τa)$⁠, and $πw¯I,FD(tm)$⁠, and $P(πq′D(tm+1+τa)|πq′D(tm+τa))$ is the conditional probability of $πq′D(tm+1+τa)$ given $πq′D(tm+τa)$⁠. Note that $ΔSS,w¯I,F→q′$ takes a positive value when $w¯I,F$ drives $q′$⁠. In this study, we consider the advection time of the vortex rings from the z position of the spatial average cross-sectional flow velocity fluctuations in the injector to each grid of the flame region as the delay time τ_a (= $(Li+Lc)/va$⁠, where v_a is the advection velocity of the vortex rings from inside the injector). v_a is estimated to be 108 and 168 m/s at the fuel and oxidizer injectors, respectively. L_c is the distance from the injector rim to each grid of the flame region in the combustor and L_i is the distance from the z position of the spatial average cross-sectional flow velocity fluctuations to z = 0 mm, where z = 0 mm corresponds to the injector rim. Note that τ_a is set to zero for the estimation of $ΔSS,p′c→q′$ in this study to examine the synchronized state of pressure and heat release rate fluctuations simultaneously.

The symbolic recurrence plots^20,29 incorporating the rank-order patterns in a time series can quantify the similarity of the dynamic behaviors of two variables and have been adopted for examining the synchronized state between two variables.^29,35 In this study, we adopt the symbolic recurrence plots for $p′c$ and $q′$ during a transition and the subsequent combustion oscillations. They are represented by the matrix consisting of $RS,ij$ as:

Here, $πp′cD(ti)$ is the rank-order pattern of $p′c, πq′D(tj)$ is the rank-order pattern of $q′$⁠, and D is the embedding dimension. A long diagonal structure appears in $RS,ij$ when the dynamic behavior of a time series coincides with that of another time series. $p′c$ and $q′$ have a phase difference when the diagonal structure is formed parallel to the main diagonal line of i = j. τ corresponds to the phase difference from the main diagonal line. We estimate the determinism in the symbolic recurrence plots D_R defined as

Here, N is the total number of data points of $p′c$ and $q′$⁠. l is the diagonal length, l_min is the minimum allowable length of the diagonal line, and $Pτ$ is the frequency distribution of length l in the diagonal line. D_R represents the ratio of the number of dots in the symbolic recurrence plots forming the diagonal line structures to the total number of dots. Note that D_R takes a high value as two time series are synchronized with each other.

The product of the pressure and heat release rate fluctuations in the combustor, denoted as $p′cq′$⁠, is an important physical quantity for estimating the driving and damping regions of combustion oscillations. The region with $p′cq′>$ 0 corresponds to the acoustic power source that drives combustion oscillations. Krishnan et al.³⁶ have recently proposed a weighted and undirected spatial network consisting of $p′cq′$ and extracted an acoustic power source of combustion oscillations in a bluff-body-type turbulent combustor. Following a similar methodology,³⁶ we adopt an unweighted and directed spatial network, namely, a thermoacoustic power network, in this study. Each grid in the numerical domain corresponds to a node in the network. If node j where $p′c,jq′j>0$ exists within radius r_i, the node i where $p′c,iq′i>0$ is connected to node j. Here, r_i of node i is the distance that the sound wave at node i propagates during the time resolution (⁠$1×10−5 s$⁠). The sound velocity at node i is defined as $ci(=γRTi$⁠, where γ is the specific heat ratio, R is the gas constant, T_i is the temperature at each grid, $γ=1.4$⁠, and R = 287 $J/(kg K$⁠). Note that small changes in γ and R do not significantly affect the construction of the network under the conditions in this study: thus, we set $γ=1.4$ and R = 287 $J/(kg K$⁠) for the construction of the network. We estimate two important network properties: the degree and the cluster coefficient in the network. The degree k_i is obtained using the adjacency matrix A_ij,

The cluster coefficient C_i is estimated as the average transitivity over all nodes,

Here, E_i is the number of triangles including node i.

Figure 2 shows the temporal evolutions of the fuel pressure fluctuations $p¯I,F$⁠, fuel flow velocity fluctuations $w¯I,F$⁠, oxidizer pressure fluctuations $p¯I,O$⁠, and oxidizer flow velocity fluctuations $w¯I,O$ in the injectors. Here, $p¯I,F, w¯I,F, p¯I,O$⁠, and $w¯I,O$ are spatially averaged in the ranges of x = 0 mm, 67.7 mm $≤y≤$ 69 mm, and 69.55 mm $≤y≤$ 74.45 mm for the fuel and oxygen injector cross sections, respectively. $p¯I,F$ and $w¯I,F$ start to oscillate at $t≥$ 10 ms, and high-amplitude and strong periodic fluctuations are formed as t exceeds approximately 20 ms. A similar trend is observed for $p¯I,O$ during a transition and the subsequent combustion oscillations. In contrast, $w¯I,O$ remains almost unchanged with small-amplitude fluctuations regardless of the z location inside the oxygen injector. Figure 3 shows the variations in the power spectral densities of $p¯O,F, w¯O,F, p¯I,O$⁠, and $w¯I,O$⁠. Similarly to that in Fig. 2, $p¯I,F, w¯I,F, p¯I,O$⁠, and $w¯I,O$ are spatially averaged in the ranges of x = 0 mm, 67.7 mm $≤y≤$ 69 mm, and 69.55 mm $≤y≤$ 74.45 mm for fuel and oxygen injector cross sections, respectively. Here, the power spectral density is calculated for 20 ms $≤t≤$ 40 ms, which corresponds to the formation region of combustion oscillations. As reported in a previous study,³⁰ the dominant 1T mode is observed for pressure fluctuations $p¯c$ in the combustor. A predominant peak corresponding to the 1T mode clearly appears in the power spectral densities of $p¯I,F$ and $w¯I,F$ at around 1.1 kHz. The dominant peak is also formed in $p¯I,O$⁠, but the magnitude of the power spectral density is sufficiently smaller than that of $p¯I,F$⁠. The peak value of the 1T mode for $w¯I,O$ is small, indicating that $w¯I,O$ is negligibly affected by the 1T mode. Matsuyama et al.³⁰ reported that the pressure losses of the fuel and oxidizer injectors are approximately 1 and 20 kPa, respectively. The pressure loss in the fuel injector is significantly smaller than that in the oxidizer injector since both the density and flow velocity of the fuel are lower than those of the oxidizer. As a result, the pressure fluctuations in the fuel injector due to the propagation of pressure fluctuations in the combustor strongly affect the flow velocity fluctuations in the fuel injector, compared with those in the oxidizer injector. This causes the significant difference in power spectral density in the flow velocity fluctuations between the fuel and oxidizer injectors. The pressure fluctuations in the combustor strongly affect the flow velocity fluctuations in the fuel injector but weakly affect those in the oxidizer injector. These results show that the pressure fluctuations in the combustor propagate into the fuel injector, resulting in the significant flow velocity fluctuations in the fuel injector.

Figure 4 shows the temporal evolution of the directionality index $ΔSS,w¯I,F→q′$ between $w¯I,F$ and $q′$⁠. Here, $ΔSS,w¯I,F→q′$ and $ΔSS,w¯I,O→q′$ are spatially averaged in the ranges of x = 0 mm, 65 mm $≤y≤$ 79 mm, and 30 mm $≤z≤$ 50 mm. $ΔSS,w¯I,F→q′$ always takes a positive value and remains nearly unchanged as a function of time. The flow velocity fluctuations in the fuel injector significantly drive the heat release rate fluctuations during a transition and the subsequent combustion oscillations, accompanied by a time delay between the two fluctuations. In contrast, $ΔSS,w¯I,O→q′$ always takes a negative value, indicating that the flow velocity fluctuations in the oxidizer injector do not drive the heat release rate fluctuations. Matsuyama et al.³⁰ reported that the periodic ignition of the unburnt $H2/O2$ mixture occurs after the unburnt $H2/O2$ mixture moving in the downstream direction comes into contact with high-temperature products of the $H2/air$ flame. The onset of the periodic ignition is associated with the fuel flow velocity fluctuations at the inlet of the combustor. They also mentioned that the periodic fluctuations of the ignition position act as an important factor that drives combustion oscillations. On the basis of the previous study³⁰ and the results shown in Figs. 2–4, the flow velocity fluctuations in the fuel injector generated by the pressure fluctuations in the combustor give rise to the periodic ignition of the unburned $H2/O2$ mixture, resulting in a significant change in heat release rate fluctuations in the combustor. Urbano et al.³⁷ examined the dynamics of combustion oscillations in a turbulent combustor with multiple coaxial injectors. They reported two important points: (i) the pressure fluctuations in the combustor produce the flow velocity fluctuations of the hydrogen jet at the injector outlet and (ii) the fuel flow fluctuations strongly affect the heat release rate fluctuations. Although the physical mechanism behind the onset and sustainment of combustion oscillations observed in this study differs from that in their study, large flow velocity fluctuations of the fuel flow play an important role in the significant change in heat release rate fluctuations.

Figure 5 shows the variation in the determinism D_R of the symbolic recurrence plots as a function of time. Here, D_R is spatially averaged in the ranges of x = 0 mm, 65 mm $≤y≤$ 79 mm, and 30 mm $≤z≤$ 50 mm. D_R starts to increase as t exceeds 14 ms. It monotonically increases with time and finally reaches approximately 0.72. This indicates that $p′c$ and $q′$ gradually synchronize with each other at 14 ms $≤t≤$ 25 ms. D_R can capture the initiation of the synchronized state during a transition to combustion oscillations. A highly synchronized state arises at $t≥$ 25 ms owing to the formation of limit cycle oscillations with strong periodicity. These results show that the onset of combustion oscillations is associated with the synchronization of the pressure and heat release rate fluctuations in the combustor. Figure 6 shows the temporal evolution of $ΔSS,p′c→q′$ during a transition and the subsequent combustion oscillations. Here, $ΔSS,p′c→q′$ is spatially averaged in the ranges of x = 0 mm, 54 mm $≤y≤$ 90 mm, and 0 mm $<z≤$ 50 mm. $ΔSS,p′c→q′$ is approximately −0.06 at $t≤$ 9.5 ms, which indicates that the heat release rate fluctuations drive the pressure fluctuations in the combustor. $ΔSS,p′c→q′$ gradually increases as t exceeds 9.5 ms. $ΔSS,p′c→q′$ becomes nearly zero at t = 17 ms, indicating that the degree of mutual coupling between the pressure and heat release rate fluctuations in the combustor remains almost the same during combustion oscillations. When t exceeds approximately 18 ms, $ΔSS,p′c→q′$ becomes positive and the pressure fluctuations in the combustor start to affect the heat release rate fluctuations during the transition to combustion oscillations. The change in the sign of $ΔSS,p′c→q′$ almost corresponds to that obtained by Hashimoto et al.²⁵ When the combustion state leads to combustion oscillations, a highly synchronized state is formed between the pressure and heat release rate fluctuations in the combustor.

Figure 7 shows the extracted $p′c$ and the degree distribution of the thermoacoustic power network during a transition and combustion oscillations. Here, the network is constructed from $pc,j′$ and $q′j$ in the ranges of x = 0 mm, 65 mm $≤y≤$ 80 mm, and 0 mm $<z≤$ 50 mm. As the amplitude of $p′c$ increases with t during the transition state, high-degree nodes are widely distributed along the hydrodynamic shear layer between the inner $O2$ and outer $H2$ jets [see Figs. 7(A) and 7(B)]. The presence of the high-degree nodes indicates the formation of the primary hub in the network. The primary hub suddenly disappears at t = 16.70 ms but re-emerges in the upstream region at t = 16.80 ms. The sudden collapse and re-emergence of widely distributed thermoacoustic power sources in the shear layer play an important role in the driving of combustion oscillations during the transition state. In contrast, the primary hub appears extensively along the shear layer region when $p′c$ takes a local minimum value at t = 28.85 ms during combustion oscillations. The formation region of thermoacoustic power sources is larger than that of the transition state. The formation and collapse of thermoacoustic power sources periodically occur in the shear layer region during combustion oscillations. The Rayleigh index is a well-known indicator of the driving region of combustion oscillations.² The formation location of thermoacoustic power sources in the shear layer region corresponds reasonably to the region where the local Rayleigh index is positive.³⁰ The important point to note here is that the degree distribution of the thermoacoustic power network can reveal the interaction region between thermoacoustic power sources using only the positive and negative values of the product of the pressure and heat release rate fluctuations. One cycle length of pressure and heat release rate fluctuations is required at least for the estimation of the Rayleigh index. In contrast, the degree distribution in the thermoacoustic power network can be estimated without the cycle length, allowing us to obtain the instantaneous driving region. For the average degree distribution obtained from one cycle of $pc,j′$ and $q′j$ during combustion oscillations, we observe a high degree distribution in the residence region of high-temperature products of the $H2/air$ flame associated with the attachment and detachment of the flame base to the injector rim. The thermoacoustic power network enables us to find this important region that is difficult to identify by only the Rayleigh index. These results show that, in addition to the estimation of the instantaneous driving region of combustion oscillations, the degree distribution of the thermoacoustic power network is a useful index in discussing the physical mechanism behind a transition and the subsequent combustion oscillations. Figure 8 shows the temporal evolutions of the spatially averaged cluster coefficient $C¯$ and $p′c$ during combustion oscillations, together with the spatial distribution of C_i. Here, C_i is spatially averaged in the ranges of x = 0 mm, 65 mm $≤y≤$ 79 mm, and 0 mm $<z≤$ 50 mm. When t = 29.40 ms, C_i takes high values near the injector rim. This means that the clusters of thermoacoustic power sources are produced near the injector rim when $p′c$ takes a local maximum. In contrast, the thermoacoustic power clusters are widely formed throughout the shear layer and the downstream region when $p′c$ becomes minimum at t = 29.75 ms. $C¯$ significantly changes in antiphase with $p′c$⁠, indicating the periodic formation and collapse of the thermoacoustic power clusters. The periodicity corresponds to the 1T mode in the combustor. These results show that the cluster coefficient of the thermoacoustic power network is also useful in discussing the driving mechanism of combustion oscillations.

From these results, the formation mechanism of high-frequency combustion oscillations can be summarized as follows. The propagation of the pressure fluctuations in the combustor due to the sudden onset of large heat release rate fluctuations induces the flow velocity fluctuations in the fuel injector. The periodic contact of the unburnt $H2/O2$ mixture with high-temperature products of the $H2/air$ flame owing to fuel flow velocity fluctuations at the injector gives rise to significant fluctuations in the ignition location, resulting in the heat release rate fluctuations. As the heat release rate fluctuations affect the pressure fluctuations in the combustor, both fluctuations synchronize with each other. The repetition of the formation and collapse of thermoacoustic power source clusters drives combustion oscillations throughout the synchronization. When the combustion state leads to combustion oscillations, the formation and collapse of the clusters of thermoacoustic power sources occur periodically, accompanied by a highly synchronized state. These physical processes play an important role in a transition and the subsequent combustion oscillations.

We have clarified the formation mechanism of high-frequency combustion oscillations in a model rocket combustor using sophisticated analytical methods based on symbolic dynamics and complex networks. The flow velocity fluctuations in the fuel injector generated by the pressure fluctuations in the combustor give rise to the periodic ignition³⁰ of the unburnt fuel/oxidizer mixture, resulting in a significant change in heat release rate fluctuations in the combustor. The heat release rate fluctuations drive the pressure fluctuations in the combustor before a transition to combustion oscillations, while the pressure fluctuations in the combustor gradually begin to affect the heat release rate fluctuations during a transition to combustion oscillations. When the combustion state leads to combustion oscillations, a highly synchronized state is formed between the pressure and heat release rate fluctuations in the combustor. The directional feedback process during the transition and subsequent combustion oscillations is identified by the directionality index in the symbolic transfer entropy. The mutual synchronization between the pressure and heat release rate fluctuations in the combustor is clearly shown by the determinism in the symbolic recurrence plots. The degree distribution of the thermoacoustic power network can reveal the interaction region between the thermoacoustic power sources in the combustor. The repetition of the formation and collapse of thermoacoustic source clusters in the hydrodynamic shear layer region between the inner oxidizer and outer fuel jets plays an important role in the driving of combustion oscillations. The degree distribution and cluster coefficient of the thermoacoustic power network provide a better understanding of the physical mechanism behind the formation of combustion oscillations.

The spatiotemporal data obtained by large eddy simulation³⁰ that support the findings of this study are available from JAXA. Restrictions apply to the availability of these data, which were used under license for this study. Data are available from the authors upon reasonable request and with permission from JAXA. Other data that support the findings of this study, which are obtained by time series analysis based on symbolic dynamics and complex networks, are available from the corresponding author upon reasonable request.