Dynamical systems approach to study thermoacoustic transitions in a liquid rocket combustor Free

Liquid rocket propulsion has become indispensable in aiding mankind for space exploration. In the late 1950s, when the space race was at its peak between the then superpowers: United States and the Soviet Union, NASA was desperately pushing to place a man on the moon to respond to the resounding success of the Soviet Union’s Sputnik satellite. The Apollo program, as it was called, took a serious setback when the F-1 liquid rocket exhibited instabilities containing oscillations of perilous magnitudes to render the mission impossible. A special program known as Project First was established and after six more years of approximately 2000 full-scale tests, these instabilities were mitigated to make the history we know today.¹ Similarly, these combustion instabilities have come to haunt several other missions, notably the Soviet RD-0110 engines and the Ariane space program in Europe.² These instabilities came to be widely known as combustion instabilities or thermoacoustic instabilities and have been observed in liquid rocket engines, solid rocket motors, tactical and strategic missiles, aeroderivative gas turbine engines, power-producing gas turbines, industrial boilers, etc. Thermoacoustic instability is characterized by large harmful oscillations in pressure and heat release rate arising due to the positive feedback between the acoustic pressure oscillations in a confinement and the heat release rate oscillations in the flame.³ The occurrence of thermoacoustic instabilities can potentially lead to partial or total mission failure of rockets, rising costs, increased delays, and developmental setbacks. Specifically in rockets, they give rise to large amplitude thrust oscillations compromising controllability and structural integrity. As a result, thermoacoustic instability has been studied widely by academia and various national and private institutions.

Rayleigh criterion³ was one of the earliest measures to assess the stability of combustion systems.⁴ This criterion judges the competition between the acoustic driving and damping in the system. Earlier approaches based on linear stability analysis^5–9 assessed the stability of liquid rocket engines and solid rocket motors for different configurations. Crocco and Cheng¹⁰ built the $n−τ$ model accounting for various time lags pertaining to different processes to assess the stability of the system. This theory made a great step toward characterizing thermoacoustic instability in rocket engines of various configurations.¹¹ Crocco et al.¹² showed that the time lags of both acoustic pressure and gas velocity need to be used in the $n−τ$ model for transverse oscillations. A similar approach was followed by Zinn¹³ to model continuous transverse modes in three dimensions, which are more prominent in rockets. However, linear stability analysis could not conclusively explain the various dynamics observed in rocket combustors. Crocco and co-workers have presciently predicted the importance of nonlinear dynamics.^10,14,15 Zinn and Powell¹⁶ used the $n−τ$ model with Galerkin method to analyze nonlinear combustion instability in liquid rocket engines. Estimates for limit cycle amplitude and triggering thresholds were obtained analytically by Mitchell et al.¹⁷ for a longitudinal mode rocket motor with shock waves. We refer interested readers to an excellent recent review by Sirignano¹¹ for an elaborate discussion on these techniques. Apart from these brilliant analytical attempts, focus was largely skewed toward identifying the physical processes behind these combustion instabilities and attenuating these harmful oscillations with the help of baffles, Helmholtz resonators, and spray and impingement alterations.² Another section of research in rockets defined stability boundaries for the system for different configurations^2,9,18 and changes in the controllable operating parameters such as injector spacing inside the combustion chamber, propellant temperature, etc.

Rocket combustion is a highly nonlinear and dynamic process. The nonlinearities^9,18,19 may arise out of gas dynamic processes, flame interactions, boundary interactions, high thermal energy density²⁰ (⁠$O∼30$ GW/m$3$⁠), and the turbulent base flow. Furthermore, extreme rates of heat addition in rockets is a major source of nonlinearities in rockets. The magnitudes of the oscillations of the system variables approach the order of magnitude of the mean variables. As a consequence, the nonlinearities in the system become significant and promote the transfer of energy across higher modes. The turbulent base flow induces wrinkles along the flame boundaries which are smoothed out at different rates depending on their length scales. The presence of flow separation at sharp edges, rapid flow expansions, and interaction of the acoustic oscillations with the coherent structures^21,22 in the reactive flow-field add upon the nonlinearities in the system. Furthermore, the wave steepening mechanism causes acoustic waves to turn into shockwaves.^23–25 The usage of nonlinear theory correctly predicts the sawtooth wave profiles in pressure for cases containing shock discontinuities, while linear theory predicts smooth sinusoidal waveforms.^11,23,26

Chester,²⁷ in a classic study, investigated the conditions required for the occurrence of flow discontinuities in the form of shock waves²⁸ in a closed duct using a rigorous gas dynamics framework. Inspired by the experimental work of Saenger and Hudson,²⁸ he showed that the shock waves arise as the natural solution for frequencies close to the resonant frequency since the nonlinear terms from acoustics, viscosity, and heat conduction contribute significantly. Furthermore, he proved that the number of shock waves in the duct increase linearly as the oscillations occur at harmonic frequencies. He further showed the effects of bulk viscosity of the fluid and boundary layer on the oscillations at near resonant frequencies.

On top of all these events, there exist several interactions across various subsystems such as injector hydrodynamics and flame dynamics, rendering the system complex.^{7,9,18,22,29,30} Several processes occurring in rocket engines are artifacts of the nonlinearities in the system.^9,25 Limit cycle oscillations could arise due to the balance between the acoustic driving and damping mechanisms in the system along with other limiting mechanisms like propellant flow.

A stable combustor can be excited with a finite amplitude disturbance to trigger self-sustained oscillations of considerable amplitudes. This phenomenon is known as triggering instabilities in rockets. During triggering, the system dynamics transitions to high amplitude state of oscillations through a finite amplitude perturbation above a threshold amplitude, called triggering amplitude. When the amplitude of the initial condition is less than the triggering amplitude, the system behavior decays asymptotically to a stable state. The phenomenon of triggering is observed when the system is operating in the bistable zone. In rockets, the time history of acoustic pressure oscillations is usually accompanied by a rise in the mean pressure levels. This phenomenon, known as DC shift²⁵ in the rocket propulsion literature, exposes the rocket to dangerous pressure levels. As a result of these nonlinear behaviors, it is vital to understand the dynamics exhibited by a rocket combustor from the perspective of nonlinear dynamics.

Furthermore, much of the research focused on the combustion stability assessment of rockets are based on building accurate models and computational fluid dynamics (CFD) simulations, owing to the high costs involved in testing even subscale hardware.^{7,9,22,31–33} Most of the studies have focused solely on understanding the dynamics during thermoacoustic instability. However, to characterize the various dynamics present in the system and build models that capture the relevant features, we need to characterize the dynamics during the transition from stable operation to thermoacoustic instability. Many studies exist in the rocket literature^34–38 pointing out to the exponential growth of amplitudes of acoustic pressure oscillations during the transition from stable operation to unstable operation. Recently, Selvakumaran et al.³⁹ detected the signature of intermittent oscillations in heat release rate fluctuations of a composite solid propellant. Adopting tools from dynamical systems theory, Guan et al.⁴⁰ showed switching between period-2 and period-3 oscillations during the state of thermoacoustic instability in a full-scale solid rocket motor.

Another challenge faced by the rocket propulsion community is to identify the dynamical transitions from stable operation to thermoacoustic instability. Conventional measures^41,42 such as root mean square of the oscillations, maximum amplitude from the amplitude spectrum through Fourier transform, etc., cannot be applied universally to all rocket combustors without a priori knowledge of the oscillation amplitudes. Different rocket combustors vary in their mean operating pressure, choice of propellants, geometry, etc. Given the widely different operating conditions and amplitudes observed for different rocket combustors, a measure which is bounded within a certain range of values for different rocket combustors would be a better candidate than the conventional measures to track the transition. Recently, Orth et al.⁴³ used the maximum of cross correlation, which is bounded between $−1$ and 1 to distinguish the stable operation from thermoacoustic instability. However, this measure cannot isolate intermittency from thermoacoustic instability. Furthermore, it requires careful selection of the two signals to be cross correlated. In spite of several advances, the rocket propulsion community is still in need of robust measures from dynamical systems theory to characterize and also to detect the transition from stable operation to thermoacoustic instability.

Adopting the framework of dynamical systems and complex systems theory, recent studies in the gas turbine literature have shown immense progress toward understanding several dynamical states of combustor operation such as chaos, period-$n$ limit cycle, and quasiperiodicity.^44–46 Several measures such as Hurst exponent,⁴⁷ recurrence quantification measures,^48–52 and measures from complex network analysis^53–55 have been deployed to detect the proximity to the onset of thermoacoustic instability. Furthermore, synchronization theory has been exploited to study the coupling between the acoustic pressure and the heat release rate oscillations.^56–58 In light of these advancements, it would be interesting to analyze the dynamics of a liquid rocket combustor using tools from complex systems and dynamical systems theory.

In this paper, by adopting various tools from nonlinear time series analysis, we detect the different dynamical states and also characterize the dynamical transitions observed in acoustic pressure oscillations of a liquid rocket combustor. We observe that the transition from stable operation to thermoacoustic instability occurs via intermittency, a state consisting of alternate occurrence of bursts of periodic oscillations among epochs of low amplitude aperiodic oscillations. Through the use of first return map, we unravel intricate features of thermoacoustic instability, where the periodic dynamics alternates between period-3 and period-4 oscillations. We show that the measures based on recurrence quantification analysis and multifractal analysis can aid in detecting the dynamical transitions in the acoustic pressure oscillations, which is not possible through conventional measures.

The rest of the paper is outlined as follows. In Sec. II, we briefly discuss the various tools employed in this study. In Sec. III, we concisely describe the model liquid rocket combustor used in this study. In Sec. IV, we perform time series analysis of pressure oscillations for the different dynamical states. Then, we perform recurrence quantification analysis and multifractal analysis. We conclude by detailing several measures that can be used to detect the dynamical transitions across different dynamical states in the system.

In this section, we briefly describe the methodology used to perform the nonlinear time series analysis throughout the rest of the paper.

In practical applications, such as thermoacoustic instability in the combustion chamber of a rocket, it is difficult to obtain data of all the independent variables that govern the dynamics of the system. In such situations, usually only a handful of system variables (in the limiting case, at least one) are available to be acquired by an experimentalist. The dynamics of a liquid rocket combustor in the higher dimensional phase space can be reconstructed from a state variable (for example, acoustic pressure: $p′$⁠) by Takens’ delay embedding theorem.⁵⁹ Such a reconstruction involves converting the univariate time series data into a set of delayed vectors from the appropriate choices of time delay (⁠$τ$⁠) and embedding dimension (⁠$d$⁠). We construct the vectors $x′(d)=[(p′(t),p′(t+τ),p′(t+2τ),…,p′(t+(d−1)τ))]$ from the measured pressure signal, $p′(t)$⁠. Here, $t$ is varied from 1 to $n−(d−1)τ$⁠, where $n$ is the total number of data points in the signal. Each delay vector corresponds to a state point in the phase space and the combination of all these vectors constitutes a phase space trajectory. To perform an appropriate phase space reconstruction for a particular state of the system, we need to obtain the optimum time delay (⁠$τ$⁠) and the minimum embedding dimension (⁠$d$⁠) for the given signal. Here, $τ$ can be estimated using average mutual information⁶⁰ or autocorrelation function.⁶¹ The minimum embedding dimension (⁠$d$⁠) can be obtained using false nearest neighbor method⁶¹ or alternately Cao’s method,⁶² which we use in this study.

Recurrence of state points in the phase space is a fundamental property of deterministic dynamical systems. Recurrence plots (RPs) are used to visually identify the time instants at which the phase space trajectory of the system revisits roughly the same area in the phase space.⁶³ The patterns present in a recurrence plot allow us to characterize the features of the signal embedded in the $d$-dimensional phase space. The construction of the recurrence plot requires a prior knowledge of the optimum time delay (⁠$τ$⁠) and minimum embedding dimension (⁠$d$⁠).

The recurrence plot of any time series is constructed by computing the pairwise distances between the state points of the reconstructed phase space. For a time series of length $n$⁠, the recurrence matrix is given by the following equation:

where $Θ$ is the Heaviside step function and $ϵ$ is a threshold to define the neighborhood of a state point in the reconstructed phase space. $‖xi′−xj′‖$ is the Euclidean distance between any two state points, $i$ and $j$⁠, on the reconstructed phase space. Whenever a state point in the phase space recurs in the predefined threshold, it is marked as a black point. Nonrecurring points are marked as white points in the recurrence plot. $Rij$ is one for a black point and zero for a white point. Thus, a recurrence plot is a two-dimensional arrangement of black and white points that exhibits different patterns characterizing different dynamics of the signal.

Several statistical measures can be derived from the organization of such black and white points in the recurrence plots. Such an analysis is known as the recurrence quantification analysis of a measured signal. Measures such as determinism (DET), recurrence rate (RR), trapping time, entropy, laminarity, and average diagonal length can be used to study the recurrence behavior of the phase space trajectory.^64,65 These measures could further be used to distinguish between the various dynamical states exhibited by the system. Here, we discuss the usage of determinism (⁠$DET$⁠), recurrence rate (⁠$RR$⁠), and the ratio between these quantities (⁠$RATIO$⁠) in the analysis of acoustic pressure data obtained experimentally from the model liquid rocket combustor.

Recurrence rate measures the density of black points in a recurrence plot and can be obtained as

where $N$ = $n−(d−1)τ$ is the number of state vectors in the reconstructed phase space.

Determinism measures the percentage of black points in a recurrence matrix, which form diagonal lines of minimum length $lmin$⁠,

where $P(l)$ is the probability distribution of diagonal lines having length $l$ and $lmin$ = 2.

The ratio of determinism and recurrence rate (⁠$RATIO$ = $DET/RR$⁠) has been introduced by Webber and Zbilut⁶⁶ to discover transitions in physiological systems.

Classical Euclidean geometry deals with smooth objects which have an integer dimension. However, many things in nature contain wrinkles when observed at different levels of magnification. Such objects or signals are classified as fractals, and they exhibit self-similar features at various observational scales.⁶⁷ Measures such as length, area, and volume for such objects are dependent on the scale at which the measurements are performed. The logarithmic plot of the measure of the object vs the scale at which the object is measured would give a straight line with an inverse power law.⁶⁷ The absolute value of the slope of this line is known as fractal dimension (⁠$D$⁠). The framework of fractal theory can be used to describe a fractal time series which exhibits self-similarity at various time scales.⁶⁸ For a fractal time series, $H$ quantifies the amount of correlation in the signal and is related to the fractal dimension⁶⁹ of the time series as $D=2−H$⁠. If $p(t)$ is a fractal time signal whose Hurst exponent is $H$⁠, then $p(ct)=p(t)/cH$ is another fractal signal preserving the same statistics⁴⁷.

Certain complex signals cannot be described using a single fractal dimension. These signals can be described with a range of fractal dimensions and such signals are classified as multifractals. In this study, we use multifractal detrended fluctuation analysis (MFDFA)⁷⁰ to study the multifractal characteristics of the time series of acoustic pressure oscillations. To estimate the Hurst exponent, the time series [$p(t)$] is mean [$⟨p(t)⟩$] adjusted to get a cumulative deviate series $yi$ as

The deviate series is then separated into an integer number $nw$ nonoverlapping segments of equal span $w$⁠. To look for trends in the segments, a local polynomial fit (⁠$y¯i$⁠) is made to the deviate series $yi$ and the fluctuations about the trend are obtained by subtracting the polynomial fit from the deviate series. Next, a quantity known as structure function (⁠$Fwq$⁠) of order $q$ and span $w$ can then be obtained from the fluctuations for $q≠0$ as

For $q=0$⁠, we have

The generalized Hurst exponents [$H(q)$] is then obtained from the slope⁷¹ of the linear regime in a log-log plot of $Fwq$⁠, for a range of span sizes, $w$⁠. In this study, we obtain this linear regime for 2–10 cycles⁷² of the acoustic oscillations observed with a frequency of 2650 Hz during thermoacoustic instability. Henceforth, the generalized Hurst exponents can be represented as a spectrum of singularities, $f(α)$⁠, via a Legendre transform⁷³ 

This spectrum, represented as a plot of $f(α)$ against $α$⁠, is known as the multifractal spectrum. The multifractal spectrum provides information on the fractal characteristics of the data. Further details regarding MFDFA can be found in the work of Kandelhart et al.⁷⁰ and Ihlen.⁷¹ 

In the literature, the generalized Hurst exponent $H(q)$ for $q=2$ is popularly known as the Hurst exponent (⁠$H$⁠). For $q=2$⁠, $H$ becomes the scaling of the root mean square of the standard deviation of the fluctuations with the window size. Since its introduction, $H$ has been used for various applications.^74,75 In thermoacoustics, Nair and Sujith⁴⁷ have used $H$ to capture the transition from stable operation to thermoacoustic instability via intermittency in a laboratory-scale turbulent combustor. Also, Unni and Sujith⁴⁹ have used $H$ as a precursor to detect blowout in a turbulent combustor.

A schematic diagram of the multielement model liquid rocket combustor is presented in Figs. 1(a) and 1(b). A detailed description of the flow conditions and the experimental hardware can be found in Orth et al.⁴³ and Harvazinski et al.⁷⁶ The oxidizer is supplied by an oxidizer rich preburner that uses hydrogen as the fuel and is located upstream of an oxidizer manifold. The preburner provides oxygen with 4%–5% mass fraction of water vapor to the experiment at a mean chamber pressure of 6.55 MPa and a mean temperature of 635 K. The oxidizer manifold is sized to minimize dynamic pressure losses and provide uniform flow to each of the injection elements downstream of it. Each injector has a choke plate upstream of it to decouple any feed system dynamics from the experiment and vice versa. Methane is injected through shear coaxial injector elements at the downstream end of the oxidizer posts through a manifold with a choked inlet. The mean Mach number in the oxidizer posts is 0.25 at nominal operating conditions. A centerline-centerline injector spacing of 25.7 mm is used in the test case presented in this study. The injector exit diameter is 15.8 mm. The exit nozzle has the same aspect ratio as the combustion chamber and is designed to obtain a mean chamber pressure of approximately 1140 kPa during the test. The chamber width is chosen to drive self-excited transverse mode dynamics at a frequency of 2.65 kHz at nominal test conditions. The length of the chamber is designed for a fundamental longitudinal (1L) mode of 3475 Hz so that the transverse mode harmonic frequency does not coincide with the 1L mode or its harmonics.

The combustion chamber is instrumented densely with high frequency pressure transducers and K-Type thermocouples. Piezoresistive Kulite WCT-312M sensors (observational error: $±6.895$ kPa) sampled at 250 kHz are used to measure the pressure in the propellant manifolds, oxidizer post, and the combustion chamber. The sensors are mounted in a recess cavity to avoid thermal saturation effects. The cavity is designed as a Helmholtz resonator with a resonance frequency of 22.4 kHz. The location of the pressure transducers used for the analysis are labeled in Fig. 1(b).

A representative time series of acoustic pressure oscillations obtained from the pressure transducer located at the right side wall of the combustor is shown in Fig. 1(c). The time interval in region I corresponds to starting of the preburner and the ignition of the main chamber. The first jump in the pressure signal close to 1 s corresponds to the start of the preburner and the second jump around 2.5 s corresponds to the ignition of the main chamber. Region III pertains to the shutdown of the engine. The acoustic pressure oscillations in region II are of prime interest in this study, as this interval of the signal represents the actual dynamical transitions from stable operation to thermoacoustic instability in the liquid rocket combustor.

In this section, we characterize the temporal behavior of acoustic pressure oscillations observed during the onset of thermoacoustic instability in the liquid rocket combustor. Toward this purpose, we examine the time series of the chamber acoustic pressure oscillations, as shown in Fig. 2, acquired for the same operating conditions (working fluids, flow rates, upstream pressures, and temperatures) and the injector configurations. However, we notice that although the operating conditions are the same during experiments, the dynamics arising out of the combustor is different during each trial. The datasets chosen for the analysis along with the dynamical transitions observed are summarized in Table I.

For test A [Fig. 2(a)], we observe that the time series is entirely composed of stable operation, exhibiting low amplitude aperiodic oscillations. For test B [Fig. 2(b)], we observe small epochs of marginally large amplitude periodic oscillations interspersed within the aperiodic oscillations of the signal. We refer to this dynamical state as intermittency. In general, intermittency refers to a dynamical state composed of high amplitude bursts of periodic oscillations amidst epochs of low amplitude aperiodic oscillations in an apparently random manner.⁴⁸ Next, we obtain a transition from stable operation to thermoacoustic instability via intermittency for test C [Fig. 2(c)]. Here, thermoacoustic instability is comprised of large amplitude periodic oscillations. For test D [Fig. 2(d)] and test E [Fig. 2(e)], we detect only two dynamical states: intermittency followed by thermoacoustic instability without the occurrence of a stable combustor operation. However, the time spent in the periodic epoch of intermittency is higher during test E than that for test D. The reasons behind such a difference in the dynamics of the combustor behavior for the same operating conditions remain unanswered.

A careful observation of the dynamics of the liquid rocket combustor shows the existence of three primary dynamical states in the acoustic pressure oscillations. These states are stable operation (low amplitude aperiodicity), intermittency (epochs of periodicity interspersed between epochs of aperiodicity in an apparently random manner), and thermoacoustic instability (epochs of sustained periodicity). During the periodic epochs of intermittency and thermoacoustic instability, we observe that the periodic waveform nearly takes the shape of a sawtooth wave profile. Furthermore, we notice that the state of intermittency always precedes the onset of thermoacoustic instability. Such an observation is different from previous descriptions of the onset of thermoacoustic instability in rocket combustors where the transition from small amplitudes to large amplitudes is reported to occur through an exponential growth.^34–38 Recently, Orth et al.⁴³ bandpass filtered the time series of acoustic pressure oscillations in the same model multielement combustor, used in the present study. When the frequencies pertaining to the fundamental mode are bandpassed, they observed the presence of an exponential growth rate in the amplitude of oscillations. They also observed a similar exponential growth rate when the harmonic frequencies are bandpassed. However, in the present study, we analyze the time series with its entire frequency content preserved. In this study, we characterize the dynamical features of the representative portions of the time series pertaining to these three dynamical states observed during different trials of experiments. We choose stable operation of test A, intermittency from test E, and thermoacoustic instability from test E. Next, we will look into the frequency content present in these three dynamical states.

The amplitude spectrum with a frequency resolution of 12 Hz generated out of the fast Fourier transform (FFT) algorithm is plotted in Fig. 3. For stable operation [Fig. 3(a)], we observe that the amplitude spectrum is broadband, containing a wide range of frequencies at smaller amplitudes. During intermittency [Fig. 3(b)], we observe a dominant peak emerging around 2500 Hz amidst the neighboring band of frequencies. During thermoacoustic instability [Fig. 3(c)], we notice a sharp peak at $f1=2650$ Hz along with several of its harmonics (⁠$nf1$⁠) of considerable amplitudes. We have marked only the first ten harmonics (⁠$f2=2f1$ to $f10=10f1$⁠) for conciseness. The presence of several harmonics of considerable amplitudes during thermoacoustic instability is due to the spiky nature of the signal caused by the steepening of the compression wave front into a shock wave.^25,26 The shift in the dominant frequency in time is attributed to the increase in mean temperature during the transition.

To probe the hidden features of the dynamics during each state, we reconstruct the phase space traced by the acoustic pressure oscillations. Toward this purpose, we need to evaluate the optimum time delay and minimum embedding dimension for each state. Furthermore, to estimate the optimum time delay, we plot the average mutual information (AMI) for different time lags⁶⁰ as shown in the first column of Fig. 4. AMI measures the mutual dependence of the signal and its delayed version at two different time instants. The first minima of the AMI can be used as the optimum time delay for the construction of the phase space. However, we observe that the optimum time delay cannot be unambiguously determined using AMI [Figs. 4(a)–4(c)], due to the difficulty in clearly identifying the first local minima, especially in Figs. 4(b) and 4(c). Hence, we turn to the autocorrelation function (ACF) to estimate the optimum time delay.⁶¹ 

Autocorrelation function (ACF) calculates the linear correlation between a time series and its delayed copy of the same time series. The value of ACF ranges between $−1$ and 1. The optimum time delays obtained from ACF corresponds to the first zero crossing in the plot, which are denoted by dashed lines in Figs. 4(d)–4(f). The corresponding optimum time delays for stable operation, intermittency, and thermoacoustic instability are 0.04 ms, 0.136 ms, and 0.084 ms, respectively.

Furthermore, we need to estimate the minimum embedding dimension required for phase space reconstruction. We rely on Cao’s method⁶² to identify the minimum embedding dimension. The two parameters, $E1$ and $E2$⁠, are evaluated for a range of embedding dimensions from 1 to 20. $E1$ measures the ratio of mean distances between two points in the phase space in two successive embedding dimensions. When a sufficient embedding dimension is attained, $E1$ attains a value close to 1 and remains constant for further increments in embedding dimension. $E2$ is a quantity, which can distinguish between deterministic and stochastic signals. For a completely random signal, $E2$ remains nearly unity for any embedding dimension.⁶² For deterministic signals, $E2$ varies for lower embedding dimensions and saturates beyond a certain embedding dimension.

The optimum embedding dimension is the dimension, denoted by dashed lines in Figs. 5(a)–5(c), for which $E1$ and $E2$ starts to become invariant with further increase in dimension (⁠$d$⁠). In addition, we observe that $E2$ is not unity for some embedding dimensions, denoting that the dynamics during stable operation are not completely stochastic. The minimum embedding dimension chosen is 13 for stable operation [Fig. 5(a)] and 10 for both intermittency [Fig. 5(b)] and thermoacoustic instability [Fig. 5(c)].

With the optimum time delay obtained for each state, we plot the three-dimensional phase portraits for stable operation, intermittency, and thermoacoustic instability in Figs. 6(a)–6(c), respectively. We observe that the phase portraits during stable operation in Fig. 6(a) is cluttered and has no distinct repeating pattern corresponding to the low amplitude aperiodic oscillations. However, during thermoacoustic instability in Fig. 6(c), we obtain a pattern (marked 1–7 in order), which repeats at equal intervals of time. The phase portrait of this state shows a stretched trefoil-knot like structure, similar to that observed in gas phase detonations.⁷⁷ This structure is radically different from the phase portrait of thermoacoustic instability observed for gas turbine combustors, which mostly trace out a ring or elliptical orbit.^44,56 During thermoacoustic instability in this rocket combustor, due to an increase in the speed of sound because of rising temperature and convective effects in the compression phase, the waveform tends to catch up with the expansion front.^23,24 This leads to the steepening of the compression wave front into a shock wave. As a result, the pressure wave front has a faster growth in the amplitude during the compression phase compared to the slow decay of the oscillation in the expansion phase. This characteristic behavior is captured faithfully in the corresponding phase portrait wherein the phase space trajectory spends relatively shorter times during the compression phase [points 1–2 in Fig. 6(c)] compared to the expansion phase [points 2–7 in Fig. 6(c)] of the signal. During intermittency in Fig. 6(b), we obtain a phase portrait bearing some resemblance to the phase portrait during thermoacoustic instability. The presence of amplitude modulation during periodic oscillations and the aperiodic oscillations corrugates the phase portrait of intermittency.

A Poincaré map or the first return map preserves many properties of periodic, quasiperiodic, and chaotic orbits.⁶¹ Hence, we use a return map, tracking the successive local maxima of the signal, to probe the dynamics. In Fig. 7, the first return map tracking the local maxima of the acoustic pressure oscillations during stable operation, intermittency: aperiodic and periodic epochs, and thermoacoustic instability are plotted.

The trajectory traced by the return map helps us in identifying the precise dynamical state which is sometimes not apparent from the visual inspection of the three-dimensional phase portrait. In a first return map, a point is observed for limit cycle oscillations with period-1, a ring is observed for the quasiperiodic oscillations, and a clutter of points is observed for a chaotic signal.⁷⁸ Also, if the consecutive dots traced in the return map of period-$n$ oscillations are joined, it results in the trajectory of a $n$-sided polygon.

The aperiodic oscillations [see Figs. 7(a) and 7(b)] during stable operation and intermittency show a clutter of trajectories without exhibiting any specific pattern. However, for periodic oscillations [Figs. 7(c) and 7(d)] during both intermittency and thermoacoustic instability, we observe the random occurrence of period-3 and period-4 oscillations as shown by triangles (I-II-III) and quadrilaterals (1-2-3-4), respectively, in their first return maps. This further suggests that the state of thermoacoustic instability is nontrivial and is not the same as the period-1 limit cycle oscillations, which is usually observed for gas turbine engines. It is particularly interesting to note that a similar switching between period-2 and period-3 limit cycle dynamics have been reported recently for a full-scale solid rocket motor.⁴⁰ At this juncture, we must note that caution must be exercised while applying tools designed to detect conventional period-1 limit cycle oscillations as they might fail for such complex period-3 and period-4 oscillations.

The phase portraits of high-dimensional attractors are usually visualized by projecting them into the lower dimensions. However, a lot of information will be lost when the phase space is condensed into lower dimensions. Eckmann et al.⁶³ proposed a visual representation tool, known as recurrence plot that enables us to investigate the behavior of $n$-dimensional phase space trajectory through a two-dimensional representation of its recurrences. The recurrence plot contains unique patterns for each kind of oscillation. For example, periodic oscillations are represented by continuous diagonal lines because the trajectory of such signals revisits roughly the same region of phase space in equal intervals of time. For random signals, we obtain a grainy structure in the recurrence plot. For chaotic signals, unlike random signal, one would obtain isolated short lines parallel to the main diagonal line.⁷⁹ For a detailed description on recurrence plots, we encourage the reader to see Marwan et al.⁸⁰ 

Recurrence plots (RPs) for the acoustic pressure oscillations during the stable operation, intermittency (both aperiodic and periodic epochs), and thermoacoustic instability are shown in Fig. 8. The black patches during the occurrence of aperiodic oscillations in Figs. 8(a) and 8(b) correspond to the trajectory trapped within a small region in the phase space. The short (or broken) lengths of diagonal lines in RP [see the zoomed inset in Figs. 8(a) and 8(b)] during both stable operation and aperiodic region of intermittency imply deterministic behavior and could possibly suggest chaotic dynamics for the aperiodic oscillations. However, dedicated tests have to be performed before confirming chaotic dynamics. The recurrence plots during periodic oscillations of intermittency [Fig. 8(c)] and that of thermoacoustic instability [Fig. 8(d)] show continuous diagonal lines, indicating strong deterministic characteristics in the dynamics. However, during the periodic portion of intermittency, the diagonal lines are relatively broken due to the gradual decrease in the amplitude of oscillations in the signal.

Many complex signals exhibiting aperiodic oscillations contain certain structural characteristics, which are difficult to be captured by various tools discussed so far. Fractal theory can be used to describe such complex signals that are composed of multiple time scales. By applying fractal analysis to thermoacoustic systems, Nair and Sujith⁴⁷ showed that the stable operation (i.e., a state of combustion noise) in a turbulent combustor has multifractal features and these multifractal signatures vanish at the onset of thermoacoustic instability. By following their approach, we study the multifractal behavior of acoustic pressure oscillations observed in the model liquid rocket combustor.

In Fig. 9(a), we plot the variation of generalized Hurst exponents with the variation in the order-$q$ for different dynamical states observed during the onset of thermoacoustic instability. We notice that, during stable operation and intermittency, the large-scale fluctuations and small scale fluctuations scale differently as the variation of $H(q)$ shows a different trend for both the states. Contrary to this, $H(q)$ shows a negligible change with variation in $q$ during thermoacoustic instability, indicating the existence of a single scale during this state.

Furthermore, we observe a nonlinear variation of the mass exponents, $τ(q)$⁠, with scaling order $q$ in Fig. 9(b) for all the states except thermoacoustic instability. Generally, a linear and nonlinear variation of $τ(q)$ represents monofractal and multifractal behavior of the signal, respectively.⁷¹ This indicates that the states of stable operation and intermittency exhibit multifractal behavior that reduces to a monofractal-like behavior during thermoacoustic instability. Also, the resulting multifractal spectra shown in Fig. 9(c) for stable operation and intermittency exhibit a wide spectrum spanning several values of singularity exponents (⁠$α$⁠). Thus, the variation of generalized Hurst exponents, mass exponents, and the multifractal spectrum strongly point out to the presence of multifractal nature in these oscillations.

During thermoacoustic instability, this multifractality is lost. This loss of multifractality is evident from the invariant nature of $H(q)$⁠, the linear variation of $τ(q)$ with $q$⁠, and the collapse of the multifractal spectrum to a shorter arc centered around a nonzero $α$⁠. This nonzero value of $α$ and the noninteger value of the $H(q)$ further confirm the monofractal-like behavior of acoustic pressure signals during thermoacoustic instability. Such a monofractal behavior for periodic signals has been shown previously for plasma signals.⁸¹ 

Additionally, the multifractal spectra during stable operation and intermittency display a right skewed behavior [Fig. 9(c)]. This right skewness suggests that the multifractal dynamics of the pressure oscillations is determined predominantly by the small-scale fluctuations. It is also reflected in the reduction in the slope of generalized Hurst exponents for positive order $q$⁠, indicating that the $q$th-order root mean square values are insensitive to the local fluctuations with large magnitudes.⁷¹ Having studied the dynamical features of acoustic pressure oscillations during the onset of thermoacoustic instability, we now proceed to characterize the dynamical transitions observed in the system dynamics of liquid rocket combustor quantitatively.

We have shown that a thermoacoustic system can exhibit different dynamical states such as stable operation, intermittency, and thermoacoustic instability. A measure which can distinguish between these different dynamical states would be an ideal tool for engineers and simulators to help in assessing the stability of a rocket combustor.

In Fig. 10, we show several measures that exhibit a quantitative change during the transition from stable operation to thermoacoustic instability. In Fig. 10(a), we plot the time series of acoustic pressure without removing the mean pressure, during test C containing the transition from stable operation to thermoacoustic instability via intermittency, for which the measures are evaluated. The variation of conventional measures employed to detect the transition to thermoacoustic instability such as root mean square value [Fig. 10(b)], the variance of the oscillations [Fig. 10(c)], and magnitude of the dominant frequency from the amplitude spectrum [Fig. 10(d)] is plotted. The entire time series is split into 100 segments of 6 ms interval each for plotting Figs. 10(b) and 10(c). Due to the compromise in the frequency resolution with shorter window size, we use a relatively larger window interval of 55.6 ms, which resulted in 8 segments of the actual time series, for plotting Fig. 10(d).

The variation of both root mean square and variance of the acoustic pressure oscillations increases progressively as the system dynamics approaches thermoacoustic instability. The nonmonotonic trend in the variation of these measures prior to thermoacoustic instability is due to the presence of intermittency. The magnitude of the dominant frequency in the amplitude spectrum calculated with a frequency resolution of 18 Hz exhibits a gradual variation from stable operation to thermoacoustic instability. However, to determine the onset of thermoacoustic instability from these measures, an a priori knowledge of the expected amplitude levels out of the combustor is required. Armed with the knowledge of the amplitude levels during the onset of thermoacoustic instability in a combustor, one can determine whether thermoacoustic instability is attained or not. However, in most scenarios, the amplitude levels in a combustor are difficult to predict as they depend highly on the operating conditions, working fluids, etc. Even if this is overlooked, using these measures, we cannot robustly distinguish the transition between the states of stable operation, intermittency, and thermoacoustic instability.

In an attempt to overcome the shortcomings of these conventional measures, Orth et al.⁴³ introduced the maximum of cross correlation (⁠$CCmax$⁠) as a measure to distinguish between stable operation and thermoacoustic instability. $CCmax$⁠, bounded between $−1$ and 1, captures the highest similarity between two time series. In Fig. 10(e), we show the variation of the maximum value of the cross correlation (⁠$CCmax$⁠) between the acoustic pressure signals acquired at two different locations in the combustor [labeled as “Fuel Manifold Pressure” and “Right Wall Pressure” in Fig. 1(b)]. In this study, we find that $CCmax$ is unable to distinguish between intermittency and thermoacoustic instability as the values of $CCmax$ are nearly the same during intermittency and thermoacoustic instability. Next, we show the variation in the recurrence based measure: the ratio of determinism to recurrence rate (⁠$RATIO$⁠) in Fig. 10(f). We note that the value of $RATIO$ starts decreasing with the onset of intermittency and decays to almost zero during thermoacoustic instability. $CCmax$ and $RATIO$ are plotted for a window size of 7.5 ms corresponding to 20 cycles of oscillations. The robustness of $RATIO$ in distinguishing the different dynamical states for a range of recurrence thresholds is discussed in  Appendix A.

Finally, the variation of fractal measures, Hurst exponent (⁠$H$⁠) in Fig. 10(g) and multifractal spectrum width (⁠$α2$ − $α1$⁠) in Fig. 10(h) are plotted to distinguish the dynamical transitions across stable operation, intermittency, and thermoacoustic instability. Here, $α2$ and $α1$ are the extreme values of the singularity exponents covered by the multifractal spectrum. The multifractal spectrum width (⁠$α2−α1$⁠) is calculated by measuring the range of singularity exponents covered by the spectrum. For the multifractal measures, a window size of 37.6 ms corresponding to 100 cycles of oscillations with an overlap of 90 cycles is used. The multifractal spectrum width drops from near 0.4 to lower than 0.02 during the onset of thermoacoustic instability. However, the presence of intermittency cannot be detected by this measure. The value of Hurst exponent (⁠$H$⁠) varies from around 0.5 during stable operation to less than 0.1 during the onset of thermoacoustic instability. During intermittency, if the value of $H$ drops below 0.1, this model rocket combustor can be considered to be in the proximity of an impending thermoacoustic instability. However, the critical Hurst exponent below which thermoacoustic instability is imminent may vary from system to system. Hence, $RATIO$⁠, Hurst exponent, and multifractal spectrum width collectively can be used to distinguish the combustor operation across all three states for a rocket combustor, as they possess fixed values for a particular type of dynamics, unlike traditional measures such as $rms$ value, amplitude of frequency peaks, and variance of the oscillations. The statistical significance and robustness of the multifractal measures for different parameters are described in  Appendix B.

Next, in Fig. 11, we show that the same measure $RATIO$ can also be used to detect the transitions from aperiodic to periodic oscillations, and vice versa, in a signal [see Fig. 11(a)]. We compare the efficacy of $RATIO$ as compared to $CCmax$ in detecting such transitions. We also show the variation of $DET$ and $RR$ in Figs. 11(c) and 11(d), respectively. We observe that through a windowed variation of $CCmax$ [Fig. 11(b)] and $RATIO$ [Fig. 11(e)], we can detect the switching between periodic and aperiodic behavior during intermittency. Here, $CCmax$ is obtained by cross correlating the same two pressure signals used to calculate $CCmax$ plotted in Fig. 10. Zoomed views of the normalized pressure time series of the two signals (⁠$pn,fuel′$ and $pn′$⁠) are plotted for an aperiodic epoch of intermittency, a periodic epoch of intermittency, and thermoacoustic instability in Figs. 11(i)–11(iii). A window size of 2.3 ms corresponding to two hundred slices of the actual time series is used to calculate all measures in Fig. 11. A smaller window size is necessary to detect the aperiodic-periodic transitions. $DET$⁠, $RR$⁠, and subsequently, $RATIO$ are obtained by calculating the recurrences of the phase trajectories within a threshold of 20% of the maximum size of the corresponding attractor. The time delay and embedding dimension are calculated for the entire time series and are found to be 0.196 ms and 10, respectively.

We observe an uncharacteristically higher value of $DET$ for the aperiodic oscillations compared to other combustors⁸². The value of $DET$ for both aperiodic and periodic dynamics in this data [see Fig. 11(c)] remains nearly the same. The value of $DET∼1$ suggests the possibility of high deterministic features⁸⁰ in the aperiodic oscillations of the rocket combustor dynamics. This high determinism value could be a result of the dynamics of the flame front, arising from the globally unstable hydrodynamic field.⁸³ 

On the other hand, the value of $RR$ exhibits a significant increase during the transition from aperiodic to periodic oscillations [see Fig. 11(d)]. Hence, $RATIO$ exhibits a higher value for aperiodic oscillations and a lower value for periodic oscillations. On the other hand, for $CCmax$⁠, we expect a value close to 0 for aperiodic oscillations with low similarity and a higher value close to 1 for periodic oscillations with large similarity. The blue and green shaded regions in Figs. 11(a)–11(e) represent an aperiodic epoch and a periodic epoch, respectively, during intermittency. During the aperiodic epoch, we observe that $CCmax$ shows lower values, while $RATIO$ exhibits larger values. We observe the opposite behavior in both $RATIO$ and $CCmax$ during the periodic epoch of intermittency. During thermoacoustic instability (see red shaded region in Figs. 11(a)–11(e)], the values of both these measures are largely invariant, denoting sustained periodic behavior in the system. For this state, we observe that the values of both $CCmax$ and $RATIO$ are low. The lower value of $CCmax$ is unexpected during thermoacoustic instability as the dynamics during this state is periodic.

The reason behind the lower value of $CCmax$ for both periodic and aperiodic oscillations can be understood from the overlapped plot of the two pressure signals used for the calculation of $CCmax$ [see Figs. 11(i)–11(iii)]. To aid us in detecting the similarity, the two time series (⁠$pn,fuel′$ and $pn′$⁠) are normalized. For the aperiodic epoch of intermittency, we do not observe any similarity between the two signals [Fig. 11(i)]. During the periodic epoch of intermittency [Fig. 11(ii)], we observe that the two signals follow a nearly similar trend at a finite nonzero time lag, leading to higher values in $CCmax$⁠. On the contrary, during the state of thermoacoustic instability [Fig. 11(iii)], we notice that the time series of $pn,fuel′$ contains significantly higher frequencies, whereas that of $pn′$ contains lower frequency corresponding to fundamental mode of the combustor (2650 Hz). This difference in the oscillations of acoustic pressure at different locations contribute to lower the value of $CCmax$⁠. Unlike $CCmax$⁠, we observe that the lower values of $RATIO$ correctly captures the periodic oscillations during thermoacoustic instability as well as during intermittency. This suggests that using $RATIO$ is better than $CCmax$ to unambiguously determine the periodic-aperiodic-periodic transitions in the acoustic pressure signal observed during the onset of thermoacoustic instability. We also remark that $RR$ can be a good candidate to distinguish the aperiodic-periodic transitions if there is a significant variation in $RR$ during the aperiodic-periodic transitions.

The dynamics of acoustics pressure oscillations during the transition from stable operation to thermoacoustic instability in a model multielement liquid rocket combustor is analyzed. We observe that the transition from small amplitude stable operation to large amplitude thermoacoustic instability occurs through intermittency. Intermittency is a dynamical state wherein bursts of high amplitude periodic oscillations appear amidst epochs of low amplitude aperiodic oscillations, distributed in a seemingly random manner.

The waveform during thermoacoustic instability is highly nonlinear, consisting of typically steepened pressure wavefronts leading to the formation of shock waves, and is significantly different from the sinusoidal limit cycle oscillations typically seen in gas turbine combustors. As a result, we obtain a characteristic trefoil knotlike shape of the phase space attractor during thermoacoustic instability. Furthermore, we detect the dynamical switching between possibly period-3 and period-4 oscillations in an apparently random manner during thermoacoustic instability and the periodic epochs of intermittency. Such complex limit cycle dynamics are seldom seen in gas turbine combustors.

Through a suitable multifractal analysis, we detect the collapse of multifractality during the onset of thermoacoustic instability. We present a recurrence based measure (⁠$RATIO$⁠) and two fractal based measures (multifractal spectrum width and the Hurst exponent) that can be used to distinguish between different states of combustor operation. We found that these measures are more robust than the existing measures such as root mean square of the oscillations, spectral amplitude, maximum of cross correlation, etc., in distinguishing the dynamical state of a rocket engine. The measures illustrated in this study can be used to validate the CFD multifidelity simulations used for optimizing the stability and performance metrics of the rocket combustor. Such an approach can reduce the developmental time scales of a rocket engine. To summarize, the signals pertaining to rocket combustors are different from their gas turbine counterparts and other derived laboratory combustors due to the significant contribution of nonlinearities in the rocket combustor. Hence, extreme care must be exercised while extending the results obtained for gas turbine combustors to rocket combustors.

The authors gratefully acknowledge the funding provided by Air Force Office of Scientific Research under Award No. FA2386-18-1-4116 (Grant Program Manager: Lt. Col. Sheena Winder). R.I.S. thanks Dr. Venkateswaran Sankaran (AFRL) for his numerous productive discussions. We thank Michael R. Orth for the experimental data. P.K. and I.P. are grateful to MHRD-India and Indian Institute of Technology Madras for providing research assistantship.

The temporal variation of the recurrence measure, $RATIO$⁠, for different recurrence thresholds (⁠$ϵ$⁠) is shown in Fig. 12. The recurrence thresholds are selected as a proportion of the maximum size (⁠$s$⁠) of the phase space attractor of acoustic pressure oscillations reconstructed using Takens’ delay embedding theorem for test C. Owing to the dependence of the magnitudes of the determinism (⁠$DET$⁠) and recurrence rate (⁠$RR$⁠), the absolute values of $RATIO$ vary during the transition. However, the underlying trend required to distinguish the different dynamical states remains intact for the all the recurrence thresholds shown.

Prior to performing the multifractal analysis, we need to estimate the range of scales necessary to capture the multifractal characteristics of the acoustic pressure oscillations. As explained in Sec. II B, we plot the structure function (see Fig, 13) against the range of binarized scales necessary to capture the small-scale and large-scale fluctuations. From the plot, we examine the trends for the acoustic pressure fluctuations during the dynamical states of stable operation, intermittency, and thermoacoustic instability. We observe a linear regime for the range of scales from 2 to 10 cycles of the dominant instability frequency (2650 Hz) for all the dynamical states. Hence, the multifractal measures such as Hurst exponent (⁠$H$⁠) and multifractal spectrum width (⁠$α2−α1$⁠) are computed with this range of scales.

In Fig. 14, we show the temporal variation of Hurst exponent (⁠$H$⁠) estimated for test C. The error bars are estimated with a confidence of 90% based on the goodness of the fit to measure the slope in the plot of structure function. We observe that $H$ can be used to robustly demarcate the onset of thermoacoustic instability from the states of stable operation and intermittency.

The sensitivity of $q$-range in the computation of the width of the multifractal spectrum (⁠$α2−α1$⁠) is plotted for test C in Fig. 15. Here, we observe that $α2−α1$ is fairly robust in exhibiting similar trends in the variation from stable operation to thermoacoustic instability via intermittency.