Early detection of combustion instabilities in liquid rocket engines: Assessment of statistical,recurrence,and fractal analyses for sub- and supercritical pressure conditions

Abstract

Safe shutdown of a liquid rocket engine in ground testing prior to the onset of damaging combustion instabilities through reliable detection of instability precursors would translate to time and cost savings in engine development programmes. Methods derived from statistical, recurrence, and fractal analysis have been successfully applied in the literature to detect precursors in unsteady pressure signals from canonical combustion experiments, gas-turbine combustion experiments, and sub-scale rocket combustion experiments operated at low pressures. In the present work, several such methods were applied to data from two cryogenic oxygen-natural gas rocket experiments operated at higher pressures than previously reported; both sub- and supercritical with respect to oxygen. The goal was to identify methods that can discern limit-cycle instabilities from intermittently unstable operation and are sufficiently responsive to be applied as emergency shut-down criteria in engine tests. Among the methods applied were the standard deviation, variance of the auto-correlation, the second spectral moment, the ratio between determinism and recurrence rate, the Hurst-exponent, and the multifractal range. The second spectral moment, the Hurst-exponent, and a measure derived from the multifractal spectrum all have short detection delays for instability onset and short-lived could be discerned from self-sustaining instabilities with an appropriate choice of threshold value. They also have moderate computation cost which makes them of interest for potential real-time implementation. The Hurst-exponent has the additional advantage of a common threshold value for all test cases addressed, demonstrating its potential for broader application independent of combustion device or operating conditions.

Keywords

Rocket engine subcritical combustion supercritical combustion liquid oxygen/natural gas single-injector multi-injector Hurst-exponent early-warning-signals

Introduction

The occurrence of high-frequency (HF) combustion instabilities poses a high risk to the safe operation of liquid-propellant rocket engines.¹ Generally stable combustion features pressure fluctuations below 5% of the mean chamber pressure. If the amplitude of the pressure oscillations exceeds this threshold but in a chaotic manner this is depicted as rough combustion. Unstable combustion displays organized oscillations and can be easily recognized against the broadband noise of the highly turbulent combustion.²

The Rayleigh-Criterion provides a theoretical foundation for thermoacoustic driving, which can lead to combustion instabilities.³ The mathematical description of this criterion is given in equation 1.

\int_{0}^{T} \int_{V} p^{'} (\vec{x}, t) {\dot{q}}^{'} (\vec{x}, t) d V d t > \int_{0}^{T} \int_{V} \sum_{i} L_{i} (\vec{x}, t) d V d t

(1)

If pressure oscillations (

p^{'}

) are in phase with heat release oscillations (

{\dot{q}}^{'}

) and exceed the losses (

\sum_{i} L_{i}

) the necessary condition for growing amplitudes is met. The stability of a rocket engine can be affected by various factors such as unexpected resonances originating from the turbopumps or feed-lines, propellant two-phase effects, and injector characteristics. Usually, the dynamic stability of the thermoacoustic system can only be reliably determined in costly full-scale testing.

While extensive health monitoring with emergency shut-down logic is standard practice in ground testing, detecting and avoiding combustion instabilities before an engine sustains damage remains a challenge. Reliable detection of precursors of oscillatory combustion in an engine with unknown stability characteristics has enormous cost-saving potential for rocket engine development and operation.

Two approaches are prominent in predicting the onset of combustion instabilities, the first being model-driven. Sujith et al.⁴ applied measures from recurrence and multifractal analysis to a multi-injector rocket combustor operating at chamber pressures below 30 bar. Within this study, a measure derived from recurrence quantification analysis (RQA) was identified as capable of distinguishing between the different operating states of the combustor. From fractal analysis, the Hurst-exponent and the multi-fractal spectrum width showed the same capabilities. Banerjee et al.⁵ additionally investigated measures derived from autocorrelation and power-spectral-densities in a laboratory-scale combustor. It was found that, again, the Hurst-exponent, but also the variation of the auto-correlation (AC) and second spectral moment could be applicable as early-warning-signals (EWSs).

The second approach originates from machine learning. Sengupta et al.⁶ combined nonlinear autoregressive models with exogenous variables and Bayesian neural network ensembles for a rocket thrust chamber to estimate the pressure amplitude with an uncertainty quantification. Waxenegger-Wilfing et al.⁷ combined measures from recurrence analysis with a support vector machine to classify different stability regimes and applied the class prediction as an EWS.

The current work focuses on model-driven approaches from statistical, recurrence, and fractal analysis and aims at assessing the applicability of these EWSs in cryogenic oxygen-natural gas combustion under conditions representative of launcher engines. Data from two different sub-scale combustion devices were used as test cases. The first device is a single-injector, optically accessible combustor and the second is a multi-injector thrust chamber. Both devices used the propellant combination of liquid oxygen and compressed or liquefied natural gas (LOX/CNG/LNG). The inclusion of two different combustors with different numbers of injector elements and different power levels is motivated by the aim to have an EWS not limited to a specific combustion device.

In these experiments, self-excited thermoacoustic instabilities occurred, with both phases of short-lived (SL) bursts of oscillatory combustion as well as sustained, limit-cycle instabilities. In contrast to stable combustion which is dominated by low-amplitude broadband noise of the highly turbulent combustion, thermoacoustic instabilities are accompanied by a tremendous increase of mechanical and thermal loads which the structure containing the combustion is exposed to and thus the potential to destroy the rocket engine in short period of time. While the intermittent bursts of oscillatory combustion are not at all a desired combustion behaviour their destructive impact is limited due to the short period of heightened thermal loads. A desired EWS for testing at the test bench would therefore be able to distinguish between SL and self-sustaining instabilities. In the latter case, a test abort needs to be triggered, whereas SL instabilities should not lead to resource-intensive false test aborts.

This paper is structured as follows. First, the test cases are described, comprising two experimental setups and five tests in total. Then, the precursor detection methods assessed in this work are introduced. The methods applied for statistical analysis were the standard deviation, variance of the AC, and the second spectral moment. For recurrence analysis, the ratio between determinism and recurrence rate was utilized. For fractal analysis, the Hurst exponent and the multifractal range were employed. Next, all the methods are applied to time series unsteady pressure data from a single test with the single-injector combustor and their utility for precursor detection is assessed. Finally, three of the most promising methods are applied to the remaining sets of test data covering a wider set of operating conditions and including the multi-injector configuration. Advantages and disadvantages of these selected methods are discussed and recommended threshold values for instability detection are offered.

Experimental method

The hot-fire tests presented in this work were conducted at the European Research and Technology Test Facility P8 for cryogenic rocket engines with the DLR research combustion chamber BKN and research thrust chamber BKD. While BKN is an extensively instrumented single-element combustor limited in overall power density compared to flight engines, the multi-element BKD features power levels and densities comparable to small upper stage rocket engines but renders itself less accessible for diagnostics.

Combustion chamber N

For these tests, BKN consisted of a single-element injector head, an optical chamber segment with 50-mm diameter, up to three additional chamber segments and a convergent divergent nozzle.

In the presented tests three different combustion chamber lengths ( $l_{c c}$ ) were used: 359, 539 and 609 mm. The nozzle throat diameter ( $d_{t h r o a t}$ ) is 14.5 mm for all tests. A schematic representation of the combustor also indicating the diagnostic capabilities is illustrated in Figure 1(a). Within this work, the acoustic measurements of the unsteady pressure sensor ( $p_{c c, 1}^{'}$ ) closest to the injection plane, distanced 34.5 mm further downstream, are analyzed. The signal was sampled at a rate of 100 kHz with a 30 kHz second-order anti-aliasing filter and have a measurement range of $\pm 30$ bar. The shear coaxial injection element features a tapered and recessed LOX post. The schematic representation of the injection element with its dimensions are given in Figure 1(b). A detailed description of the experiment can be found in previous publications reporting on data gained with BKN and LNG⁸ as fuel.

Figure 1.

Combustor and injection element schematics in BKN. (a) Experimental combustor model ‘N’; (b) Injection element.

Thrust chamber D

BKD consists of a multi-element injector head, a measurement ring, at least one cylindrical chamber segment and a convergent divergent nozzle.

For the two presented tests with BKD it was configured with two cylindrical chamber segments resulting in an overall cylindrical length of 400 mm. Again, a schematic representation of the combustor also indicating the diagnostic capabilities is illustrated in Figure 2(a). The measurement ring mounted between the injector head and the first chamber segment is instrumented with four unsteady pressure sensors at the circumferential angular positions of 0^∘, 80^∘, 180^∘, and 300^∘. The acoustic measurements of the first unsteady pressure sensor were analyzed within this work. The sample rate of the signal is also 100 kHz with a 30 kHz second-order anti-aliasing filter. The measurement range was set to only $\pm 10$ bar. This resulted in cropped signals during the intervals of increased thermoacoustic activity. LOX and LNG are injected through 42 shear coaxial injection elements with a tapered LOX post, as illustrated in Figure 2(b) and an inner (exit) diameter of the LOX posts of $d = 3.7$ mm, a tip width of 0.15 mm and a size of the annular gap of 0.5 mm. The analyzed tests featured a recess of 2 mm, equivalent to 0.54 $d$ . A detailed description of the experiment can be found in a previous publication reporting on injection-coupled (IC) and limit-cycle combustion instabilites with BKD and LNG⁹ as fuel.

Figure 2.

Combustor and injection element schematics in BKD. (a) Experimental thrust chamber (schematic)⁹ model ‘D’; (b) Injection element.

Operating conditions

Acoustic measurements of five test runs are analyzed in this work. Within these tests the chamber pressure ( $P_{c c}$ ) ranges between 40 and 60 bar and ratio of oxidizer to fuel mass flow rate (ROF) at the main injector between 1.5 and 3.5. A description of the operating conditions is given in Table 1 for the tests with BKN and in Table 2 for the tests with BKD. The operating conditions are described with the chamber pressure, ROF, the temperatures in the manifolds of oxidizer ( $T_{o}$ ) and fuel ( $T_{f}$ ). Additionally, the chamber length is provided since it is varied for the tests with BKN. It is further indicated if the occurred HF phenomena were IC, SL or sustained in a limit cycle (LC).

Table 1.

Injection conditions (envelope), chamber length and HF characteristics for combustor ‘N’.

	$P_{c c}$	ROF	$T_{o}$	$T_{f}$	$l_{c c}$	HF
Test ID	(bar)	–	(K)	(K)	(mm)	–
N.1	40–52	3.0–3.4	$\approx$ 105	200–250	429	SL
N.2	40–50	3.2–3.6	$\approx$ 105	190–250	539	SL + LC
N.3	40–55	2.0–3.2	$\approx$ 105	250	609	LC

HF: high-frequency; ROF: ratio of oxidizer to fuel mass flow rate; SL: short lived; LC: limit cycle.

Table 2.

Injection conditions (envelope), chamber configuration and HF characteristics for combustor ‘D’.

	$P_{c c}$	ROF	$T_{o}$	$T_{f}$	$l_{c c}$	HF
Test ID	(bar)	–	(K)	(K)	(mm)	–
D.1	40–60	1.5–2.0	110–120	200–210	426	IC
D.2	50	1.7–2.7	110–120	200–215	426	SL + LC

HF: high-frequency; ROF: ratio of oxidizer to fuel mass flow rate; IC: injection coupled; SL: short lived; LC: limit cycle.

The test sequences of the tests with BKN are shown in Figure 3 and the tests with BKD in Figure 4. The plots show a spectrogram (top) of the unsteady pressure signal (middle), along with traces of the static chamber pressure, LOX and CNG/LNG injection temperature and ROF (bottom). These traces of parameters describe the performed sequence of operating conditions, while the raw unsteady pressure signal (middle) and the spectrogram (top) describe the acoustic character.

Figure 3.

Test sequence spectrograms of Tests N.1-N.3 showing solely short-lived bursts of oscillatory combustion in case of N.1, additional self-sustaining instabilities in N.2 and solely limit-cycle instabilities in N.2.

Figure 4.

Test sequence spectrograms of Test D.1 and D.2 showing high-amplitude combustion instabilities for distinct operation conditions in D.1 and short-lived in combination with self-sustaining instabilities in D.2.

The test sequence spectrogram for the BKN test N.1 shows solely SL bursts of oscillatory combustion. These bursts resemble phenomena depicted as coherence resonance¹⁰ or intermittency.¹⁰ While both phenomena share common signatures in the pressure data they differ in their underlying mechanisms. In the case of test N.1 it is likely that coherence resonance is observed in the early phase of the test. Perhaps in the later phase, a transition to intermittency takes place with the recurrence rate of these bursts increasing as a system parameter is, most likely, approaching its critical value. However, due to the limited data a clear distinction between coherence resonance and intermittency is not feasible and does also not hinder the main goal of this work. The test N.2 also features those SL bursts but ultimately one of these bursts transitions into limit-cycle thermoacoustic instabilities. An abrupt change from stable to unstable combustion appears in Test N.3.

The test sequence in the BKD test D.1 shows slightly heightened amplitudes of the chamber’s first tangential mode which transition to high amplitude combustion instabilities due to a change of the operating conditions. BKD test D.2 shows similar characteristics as the BKN test N.2 even though BKD features a quite different combustion hardware.

In combination, all five tests cover a wide range of operating conditions, combustion stability dynamics and hardware. It should also be noted that combustion was initiated with a torch igniter. However, the analyzed pressure readings are extracted seconds after the completion of the ignition sequence. The residence time of propellants in these tests is in the order of miliseconds assuming no closed recirculation zones. Therefore it can be ensured that the ignition is not affecting the evaluation. The mass flows determining chamber pressure and ROF are controlled via regulation valves and flow metres.

Methodology

In the following the methods applied to analyze the data are presented. Analogue to a potential real-time application all measures were calculated with measurement values collected prior to the corresponding evaluation time. Further, no filters were applied to the data being evaluated by the following measures. However, in some of the following figures the instantaneous amplitude is presented as an indication of the level of thermoacoustic activity next to the corresponding measures. Here the Hilbert transform requires a monofrequent signal and thus a bandpass filter, including the resonance frequency of the combustion device, was applied to the pressure data for this visualization purpose.

Since flight engines are often equipped with only one unsteady pressure sensor an important requirement for all considered measures is the limitation to data from one measurement source.

Statistical measures

A popular approach for monitoring the combustion and detecting unfavourable behaviour is the calculation of statistical measures indicating the variance ( $σ^{2}$ ) or its square root, the standard deviation ( $σ$ ) of a given number of data points.

σ (p_{N}^{'}) = \sqrt{\frac{1}{N} \sum_{i = 1}^{N} (p_{i}^{'} - \bar{p^{'}})^{2}}

(2)

Here

p_{i}^{'}

denotes the

i

th and

\bar{p^{'}}

the mean value of

N

pressure samples from the acoustic pressure measurements used to calculate

σ

. A disadvantage of this approach is that the standard deviation cannot distinguish between pressure oscillations originating from combustion noise or (potentially more harming) thermoacoustic oscillations. Further, if the pressure measurement values fed to the calculation of

σ

are sampled at a certain rate, the standard deviation also highly depends on the probe rate.

Banerjee et al.⁵ proposed the variance of the $A C$ and the second spectral moment ( $μ_{2}$ ) as more reliable EWS of impending combustion instabilities. The AC for a time series $p (t)$ and time-lag $τ$ is defined as

A C (τ) = \frac{\sum_{t_{i} = 1}^{N - k} (p^{'} (t_{i}) - \bar{p^{'}}) (p^{'} (t_{i} + τ) - \bar{p^{'}})}{\sum_{t_{i} = 1}^{N} (p^{'} (t_{i}) - \bar{p^{'}})^{2}}

(3)

with

t_{i}

as the

i

th time index of

p (t)

and

k

as any positive integer smaller than the amount of samples

N

in the time series. The variance of the AC is then given as

σ^{2} (A C) = \frac{1}{N - 1} \sum_{τ = 1}^{k} (A C (τ) - \bar{A C})^{2} .

(4)

For perfectly periodic signals

σ^{2} (A C)

is 0.5.

While the previous measures can be interpreted as the second central moments about their mean, the $n$ th order spectral moment based upon a one-sided power spectral density (PSD; $P_{x x} (f)$ ) is defined as

μ_{n} = m_{n} (2 π)^{n} = (2 π)^{n} \sum_{i = 1}^{N} f_{i}^{n} P_{x x} (f_{i}) Δ f_{i}

(5)

with

f_{i}

as the frequency of the

i

th bin with the frequency spacing

Δ f_{i}

between each frequency bin.¹¹ Thermoacoustic instabilities in rocket engines exhibit oscillations with high amplitude at high frequencies compared to the broadband combustion noise. Therefore

μ_{2}

should increase drastically with the onset of combustion instabilities.

Recurrence measures

In deterministic dynamical systems the recurrence of state points ( $x^{'}$ ) in the phase space is a fundamental property.⁴ The base for each recurrence measure is the recurrence matrix. The entry of the binary matrix is calculated via

R_{i j} = Θ (ϵ - ‖ x_{i}^{'} - x_{j}^{'} ‖)

(6)

with

i, j = 1, \dots, N - (d - 1) τ_{o p t}

d

corresponding to the minimum embedding dimension and

τ_{o p t}

to the optimum time delay. Further is

Θ

the Heaviside step function,

ϵ

defines the threshold for the neighbourhood’s radius of a state point, here the acoustic pressure, in the phase space and

‖ x_{i}^{'} - x_{j}^{'} ‖

denoting the Euclidian distance between a pair of state points. Therefore

R_{i j}

equals one if a state point recurs within the set neighbourhood and zero otherwise.

τ_{o p t}

can be estimated as the first zero crossing in the AC of the analyzed time series¹² and the minimum embedding dimension can be approximated by Cao’s method.¹³ For sufficient embedding dimensions, the ratio of the mean distances between two state points in the phase space is close to one for two successive dimensions. Several rules of thumb exist for the choice of the threshold as collected by Schinkel et al.¹⁴ Here the threshold value was chosen based on the standard deviation of the analyzed time series. Since this approach suffers from low signal-to-noise ratios the value of the standard deviation was additionally scaled by the ratio of the standard deviation of the high-pass filtered series (as estimation of the signal) and low-pass filtered series (as estimation of the noise).

From the recurrence matrix, a recurrence plot can be generated. Here every value in the matrix equalling one is plotted as a black dot and white otherwise. Measures can be derived from the organization of black and white dots in a recurrence plot. Sujith et al.⁴ identified the ratio of the determinism and the recurrence rate as potential EWS. The determinism describes the percentage of ones in a recurrence matrix forming diagonal lines of a minimum length $l_{m i n} = 2$

DET = \frac{\sum_{l = l_{m i n}}^{N_{s v}} l P (l)}{\sum_{l = 1}^{N_{s v}} l P (l)}

(7)

where

P (l)

is the probability distribution of diagonal lines with length

l

and

N_{s v} = N - (d - 1) τ_{o p t}

the number of state vectors in the phase space. The recurrence rate indicating the density of ones in a recurrence matrix can be calculated as

RR = \frac{1}{N_{s v}} \sum_{i, j = 1}^{N_{s v}} R_{i j}

(8)

Fractal measures

The Hurst-exponent originates from the fractal theory, which can be used to describe a fractal time series¹⁵. It quantifies the amount of correlation in the signal and relates to the fractal dimension.¹⁶ H can be estimated in the following way⁴:

First the times series, here $p (t)$ representing the signal of the acoustic pressure sensor $p_{c c, 1}^{'}$ , is mean adjusted

y_{i} = \sum_{t = 1}^{i} (p (t) - ⟨ p (t) ⟩)) i = 1, 2, \dots, n

(9)

with

⟨ p (t) ⟩ = \frac{\sum_{t = 1}^{n} p (t)}{n} .

(10)

Next the deviate series is separated into a number

n_{w}

of non-overlapping segments with equal span

w

. In order to detrend the segments a local, here second order, polynomial (

{\bar{y}}_{i}

) is fitted to the deviate series

y_{i}

and subtracted. The structure function of order

q \neq 0

and span w for the signal is obtained by

F_{w}^{q} = {[\frac{1}{n_{w}} \sum_{i = 1}^{n_{w}} {(\sqrt{\frac{1}{w} \sum_{t = 1}^{w} (y_{i} (t) - {\bar{y}}_{i})^{2}})}^{q}]}^{1 / q} .

(11)

For

q = 0

the structure function is

F_{w}^{q} = \exp [\frac{1}{2 n_{w}} \sum_{i = 1}^{n_{w}} \log (\sqrt{\frac{1}{w} \sum_{t = 1}^{w} (y_{i} (t) - {\bar{y}}_{i})^{2}})] .

(12)

Finally, the slope¹⁷ of the linear regime in a log-log plot of

F_{w}^{q}

for a range of span sizes

w

provides the generalized Hurst exponents

H^{q}

. The generalized Hurst-Exponent of second order (

H^{2}

) is commonly denoted as Hurst-exponent. Following Sujith et al.⁴ the span size for calculating

H

varied between two and ten acoustic cycles. Generally for a non-stationary as well as unbounded signal

H

> 1

. Further reference points are

H ≃ \frac{3}{2}

for brownian noise,

H ≃ 1

for pink noise and

H ≃ \frac{1}{2}

for gaussian white noise. In case of periodic signals

H

is approaching zero. Evaluation of the generalized Hurst exponents for various orders enables to obtain the multifractal spectrum via a Legendre transform. The mass exponents can be derived from

H^{q}

with

τ_{q} = q H^{q} - 1,

(13)

and the singularities with

α = \frac{δ τ_{q}}{δ q} .

(14)

The spectrum of singularities

f (α)

is received by

f (α) = q α - τ_{q} .

(15)

f (α)

plotted against

α

then visualizes the multifractal spectrum.

Results and discussion

Assessment of early warning signals

First the applicability of the measures presented in the previous section on combustion devices operating above 40 bar mean chamber pressure is assessed. All measures are evaluated each 4 ms, which is in the range of the minimal evaluation time of computer hardware used at corresponding test benches. The amount of samples (evaluation window) is varied between the methods adjusted to the necessary amount of samples with respect to each method.

The assessment will be conducted on the basis of the BKN test N.2, since BKN is more comparable in power (density) than BKD to the experiments, which were used in literature to assess the aforementioned measures as EWSs. The pressure trace of this test is characterized by stable combustion, bursts of oscillatory combustion, and high-amplitude limit-cycle instabilities, as illustrated in Figure 5.

Figure 5.

Unsteady pressure signal ( $p_{c c, 1}^{'}$ ) and instantaneous amplitude of bandpass (BP) filtered analytic unsteady pressure signal ( $\hat{H i l b e r t} (p_{c c, 1, B P}^{'})$ ) for the 1L frequency in test N.2. The red line is indicating the time when 90% of the maximum acoustic amplitude during limit-cycle instabilities is reached.

The instantaneous amplitude based upon the band-pass (BP) filtered signal containing the frequency of the first longitudinal resonance mode (1L mode) is also visualized. Further three intervals of stable, intermittent and unstable combustion are highlighted in green, orange and red.

Figure 6 shows the phase portrait of the acoustic pressure in these three intervals in test N.2.

Figure 6.

Phase portraits of the acoustic pressure for the stable, intermittent and unstable interval (Figure 5) in Test N.2.

As expected in the phase portrait for stable combustion no specific pattern can be recognized. This is due to the lack of an organized structure in the pressure trace which is dominated by aperiodic oscillation and broadband combustion noise. The phase portrait for the unstable interval draws a completely different picture. Here analogue to the data analyzed by Sujith et al.⁴ a three-dimensional deformed triangle is observable. This structure is known from gas phase detonations.¹⁸ The intermittent phase portrait shows periodic triangle-like structures similar to the unstable interval, but superimposed by a strong drift of the triangle’s centre. This can be explained by the behaviour of the chamber pressure. While for the sustained limit-cycle instabilities the mean chamber pressure stays constant, the chamber pressure rises with increasing HF amplitude and decreases afterwards.

Another presentation of the combustion behaviour are the Poincare sections in Figure 7.

Figure 7.

Poincare sections of the acoustic pressure for the stable, intermittent and unstable interval (Figure 5) in Test N.2.

Here the successive local maxima of the acoustic pressure are traced. In the case of stable combustion the poincare section is a nearly straight line. For the unstable interval, a quadrilateral structure can be recognized, which implies that next to the dominant unstable frequency (1L) also subharmonics play a role.¹⁹ The straight line with a high point density again shows that during the phase of decaying acoustic pressure, the signal is also superimposed by noise. Finally, the poincare section for the intermittent case shows a structure mixed from the stable and unstable case. Here the aforementioned pressure rise and decline during the phase of thermoacoustic oscillations is visible by the quadrilateral structure with its centre moving along the straight line.

Figure 8 visualizes the temporal development of the standard deviation, the variance of the AC and the second spectral moment in blue and the instantaneous amplitude of the 1L mode in black. It can be noticed that $σ$ lacks the ability to distinguish between short bursts of oscillatory combustion and limit-cycle instabilities.

Figure 8.

Statistical measures ( $σ$ , $σ^{2} (A C)$ , $μ_{2}$ ) compared to the instantaneous amplitude of BP filtered analytic unsteady pressure signal ( $\hat{H i l b e r t} (p_{c c, 1, B P}^{'})$ ) for the 1L frequency in test N.2. The red lines are indicating the time when a fictive red-line would have been triggered. AC: auto-correlation; BP: bandpass.

The highest value of $σ$ arises during the phase of stable combustion interrupted by SL thermoacoustic instabilities. If, for example, the critical $σ$ value would have been set to 8 (as an initial guess), this would have aborted the test during the third burst of oscillatory combustion. However, it is not possible to select and threshold, which keeps the combustor running during oscillatory combustion and aborts the test throughout sustained limit-cycle instabilities with the potential to harm the combustion device.

$σ^{2} (A C)$ on the other hand is capable to detect the sustained periodic oscillations during instabilities, since its value increases after the onset of instabilities and is moving closer towards 0.5. As previously mentioned a perfect periodic signal would result in $σ^{2} (A C) = 0.5$ . Analogue to Banerjee et al.⁵ more than 200 acoustic cycles were needed for $σ^{2} (A C)$ to converge. This results in an evaluation window of 250 ms. Therefore the detection of the limit-cycle instabilities is delayed. The delayed detection in combination with the high necessary calculation time diminishes its suitability as EWS.

The second spectral moment is capable of distinguishing between stable combustion, SL instabilities and unstable combustion since the spectral moment is slightly increased during the sustained instabilities compared to the intermittent instabilities. A threshold value of $μ_{2} = 9.5 \times^{11}$ would not be violated during SL, but throughout self-sustaining limit-cycle instabilities. However, it would be difficult to define a meaningful threshold for specific hardware prior to the first test without knowledge of the unstable frequencies and amplitudes. The evaluation window for this measure was set to 6 ms enabling to collect enough samples for a frequency resolution of the PSD to result in different frequency bins for the low-frequency content and the HF instabilities.

The recurrence plots in Figure 9 visualize the recurrence matrix. Ideally for every interval which is evaluated the minimal embedding dimension and the optimal time lag should be calculated. This would result in time inefficient calculation process of the derived measures. Therefore it was decided to maintain $τ_{o p t} = 25$ for all interval which corresponds to the minimal $τ_{o p t}$ during thermoacoustic oscillations and $d = 9$ which corresponds to the maximum $d$ during the phase of stable combustion interrupted by SL thermoacoustic instabilities. Zbilut et al.²⁰ concluded that a “correct” delay or embedding is inappropriate. It is more important to have a sufficiently large embedding “containing” the relevant dynamics. Choosing a time delay $τ < τ_{o p t}$ leads to oversampling of the data including less significant features at smaller scales.

Figure 9.

Recurrence plots of the acoustic pressure for the stable, intermittent and unstable interval (Figure 5) in Test N.2.

Apart from the centreline mainly isolated short diagonal lines can be seen for stable combustion. This is characteristic for aperiodic fluctuations. Perhaps the slight periodicity in the diagonal lines could indicate unstable periodic orbits between multiple attractors.

During predominant stable combustion with intermittent bursts of instabilities diagonal lines parallel to the centreline are visible indicating at least for a short time a periodic signal. After the decline of the thermoacoustic oscillations still some periodic structures are recognizable revealing some cyclicities in the process.

For unstable combustion, the whole recurrence plot is filled with parallel diagonal lines characterizing a stationary periodic signal with an identical evolution of state at different times and thus high recurrence rate and determinism.

Figure 10 illustrates the RQA measures discussed in the previous section.

Figure 10.

RQA measures ( $R R$ , $D E T$ , $R D R$ ) compared to the instantaneous amplitude of BP filtered analytic unsteady pressure signal ( $\hat{H i l b e r t} (p_{c c, 1, B P}^{'})$ ) for the 1L frequency in test N.2. The red line is indicating the time when a fictive red-line would have been triggered.RQA: recurrence quantification analysis; RR: recurrence rate; DET: determinism; RDR: ratio of determinism to recurrence rate; BP: bandpass.

During bursts of oscillatory combustion there is only a minor rise in the recurrence rate (RR), its value approaches one during sustained instabilities. The determinism (DET) in the signal is quite high throughout the whole data but still increased during oscillatory combustion. Both values combined in the ratio of determinism to recurrence rate (RDR) show decreased values during SL instabilities compared to stable combustion and a drastically reduced value during sustained instabilities. If the threshold for the $R D R$ value would have been set to 4 this approach would be able to distinguish between all three combustion states. Analogue to $σ^{2} (A C)$ a relatively large window of sampled data with 150 ms is necessary due to the minimum embedding dimensions. Again, this leads to a delayed detection of the sustained limit-cycle HF combustion instabilities.

Next, the multifractal spectrum is plotted in Figure 11.

Figure 11.

Multifractal spectrum of the acoustic pressure for the stable, intermittent and unstable interval (Figure 5) in Test N.2.

Here the different multifractal characteristics for stable intermittent and unstable combustion can be extracted. Stable combustion features a high mean value for the singularities $\bar{α}$ and moderate range $\hat{α} = α_{m a x} - α_{m i n}$ for different orders $q$ . Intermittent combustion shows a reduced value for $\bar{α}$ but a relatively large range of singularities indicating strong nonlinearities and oscillations at different timescales. Finally during combustion instabilities $\bar{α}$ as well as $\hat{α}$ are approaching 0 characterizing a stationary state with a non-varying dominant frequency.

The temporal development of the Hurst-Exponent ( $H^{2}$ ), the multifractal measures $\bar{α}$ as well as $\hat{α}$ and a combined parameter of the latter are visualized in Figure 12. The evaluation window was set to 12 ms enabling at least 10 acoustic cycles to be captured. Sujith et al.⁴ recommended a variation of the span $w$ between 2 and 10 acoustic cycles. The same range of $w$ was also applied in this work.

Figure 12.

Multifractal measures ( $H^{2}$ , $\bar{α}$ , $\hat{α}$ , $\hat{\bar{α}}$ ) compared to the instantaneous amplitude of bandpass (BP) filtered analytic unsteady pressure signal ( $\hat{H i l b e r t} (p_{c c, 1, B P}^{'})$ ) for the 1L frequency in test N.2. The red lines are indicating the time when a fictive red-line would have been triggered.

$H^{2}$ is quite noisy during stable combustion but with a clear tendency of approaching zero at the onset of oscillatory combustion. While $H^{2}$ remains at an intermediate state during short oscillatory bursts, its values are closer to zero during longer bursts comparable to its value during sustained instabilities. If a threshold for $H^{2} = 0.1$ is chosen this would lead to a test abortion during one of these prolonged bursts. While no damage was sustained by the current test setups from these SL instability events, this level of detection sensitivity may be useful in testing of engines that may be harmed by this type of instability. The same applies to $\bar{α}$ since this corresponds to $H^{0}$ in most cases. As observable in Figure 11 the range of the singularities $α_{m a x} - α_{m i n}$ shows a broader spectrum. Therefore the combined parameter $\hat{\bar{α}} = \bar{α} + \hat{α}$ may help to distinguish more clearly between SL and self-sustaining instabilities. A threshold of $0.05$ for this measure would result in test abortion after the onset of the limit-cycle instabilities.

In summary, the Hurst-exponent ( $H^{2}$ ), the combined multifractal measure ( $\hat{\bar{α}}$ ) and the second order spectral moment ( $μ_{2}$ ) showed potential in identifying a harmful combustion state and even enabling the distinction between stable combustion, SL instabilities and unstable combustion with a reasonable amount of computation costs, making them interesting for potential real-time implementation.

Application of early warning signals

The three measures identified as most promising in the previous section are now applied to further tests. In Figure 13 the Hurst-exponent, the combined multifractal measure and the second spectral moment are illustrated for the BKN test N.1 which features only intermittent bursts of oscillatory combustion but no sustained thermoacoustic instabilities.

Figure 13.

Hurst-exponent ( $H^{2}$ ), multifractal parameter ( $\hat{\bar{α}}$ ) and second spectral moment ( $μ_{2}$ ) compared to the instantaneous amplitude of bandpass (BP) filtered analytic unsteady pressure signal ( $\hat{H i l b e r t} (p_{c c, 1, B P}^{'})$ ) for the 1L frequency in test N.1. The red lines are indicating the time when a fictive red-line would have been triggered.

Both fractal measures show the desired behaviour of maintaining values above their respective thresholds ( $H^{2} > 0.1$ and $\hat{\bar{α}} > 0.05$ ), while during sustained HF instabilities the thresholds were consistently crossed. The second spectral moment would abort the test at 13.488 s with the threshold identified previously for N.2.

Test N.3 (Figure 14) features a direct onset of limit-cycle instabilities. All three measures exceed their critical thresholds during the transition to HF instabilities.

Figure 14.

Hurst-exponent ( $H^{2}$ ), multifractal parameter ( $\hat{\bar{α}}$ ) and second spectral moment ( $μ_{2}$ ) compared to the instantaneous amplitude of bandpass (BP) filtered analytic unsteady pressure signal ( $\hat{H i l b e r t} (p_{c c, 1, B P}^{'})$ ) for the 1L frequency in test N.3. The red lines are indicating the time when a fictive red-line would have been triggered.

Next the behaviour of $H^{2}$ , $\hat{\bar{α}}$ and $μ_{2}$ will be assessed for tests D.1 and D.2 featuring the BKD multi-injector thrust chamber. First test D.1 in Figure 15 is analyzed.

Figure 15.

Hurst-exponent ( $H^{2}$ ), multifractal parameter ( $\hat{\bar{α}}$ ) and second spectral moment ( $μ_{2}$ ) compared to the instantaneous amplitude of bandpass (BP) filtered analytic unsteady pressure signal ( $\hat{H i l b e r t} (p_{c c, 1, B P}^{'})$ ) for the 1T frequency in test D.1. The red lines are indicating the time when a fictive red-line would have been triggered.

In $H^{2}$ as well as $\hat{\bar{α}}$ a gradual decrease of the values towards the interval of high amplitude instabilities can be seen. This probably indicates a gradual approach of the system towards it thermoacoustical unstable attractor.

Shortly prior to the onset of HF instabilities, between 27 and 29 s, a rise in $H^{2}$ , $\hat{\bar{α}}$ can be detected. The reason for this behaviour is that the spans $w$ for evaluation of these parameters were optimized to detect the first transversal resonance (1T) mode of the chamber. The change of the operating conditions preceding unstable combustion is accompanied by a brief excitation of the 1L mode leading to additional frequency content in the signal disturbing the 1T periodicity and broadening the multifractal range. This temporarily shifts $\bar{α}$ and thus $H^{2}$ as well as $\hat{α}$ to higher values.

Finally, these values decrease below the critical threshold during the high-amplitude interval. It should be noted here that the threshold for the multifractral measure was doubled to $\hat{\bar{α}} > 0.1$ due to the generally increased multifractal range in BKD as a result of higher power density in this thrust chamber. Since no common threshold is applicable to all cases this presents an additional parameter to be estimated prior to the test. In the case of the second spectral moment a gradual change can be seen as well, reaching its highest value the high-amplitude interval. Again, a suitable threshold $μ_{2} = 5 \times^{12}$ has to be chosen, which can could not have been estimated in advance.

The behaviour of the combustion in D.2 (Figure 16) is similar to N.2. A phase of stable combustion precedes a phase with intermittent instabilities followed by sustained limit-cycle high-amplitude LC instabilities.

Figure 16.

Hurst-exponent ( $H^{2}$ ), multifractal parameter ( $\hat{\bar{α}}$ ) and second spectral moment ( $μ_{2}$ ) compared to the instantaneous amplitude of bandpass (BP) filtered analytic unsteady pressure signal ( $\hat{H i l b e r t} (p_{c c, 1, B P}^{'})$ ) for the 1L frequency in test D.2. The red lines are indicating the time when a fictive red-line would have been triggered.

Analogue to the behaviour in N.2, the fractal measures decrease momentarily during intermittent bursts of oscillatory combustion and finally fall below the critical threshold of 0.1. The second spectral moment has a steep rise at the onset of instabilities. A suitable threshold would be $μ_{2} = 3.6 \times^{12}$ .

A summary of the EWS statistics is given in Table 3.

Table 3.

Thresholds for EWS, time of threshold violation with respect to each test and evaluation settings.

		$0.9 \hat{p_{L C}^{'}}$	$σ (p^{'})$	$σ^{2} (A C)$	$μ_{2}$	$R D R$	$H^{2}$	$\hat{\bar{α}}$
N.1	Threshold				$9.5 \times 10^{11}$		0.1	0.055
	Time				13.588		–	–
N.2	Threshold	4.6	8	0.25	$9.5 \times 10^{11}$	4	0.1	0.055
	Time	7.764	6.336	8.012	7.767	7.848	6.908	7.796
N.3	Threshold	4.5			$9.5 \times 10^{11}$		0.1	0.055
	Time	7.298			7.400		7.344	7.388
D.1	Threshold				$3.6 \times 10^{12}$		0.1	0.1
	Time				26.012		26.016	26.012
D.2	Threshold				$5.0 \times 10^{12}$		0.1	0.1
	Time				28.340		29.556	29.344
Rate			4 ms	4 ms	4 ms	4 ms	4 ms	4 ms
Window			4 ms	250 ms	6 ms	150 ms	12 ms	12 ms

EWS: early-warning-signal; LC: limit cycle; AC: auto-correlation; RDR: ratio of determinism to recurrence rate.

The table incorporates the threshold and time when 90% of the limit-cycle amplitude ( $0.9 \hat{p_{L C}^{'}}$ ) was reached with respect to each test. For the test N.1 no limit-cycle oscillations occurred and for D.1 and D.2 no final limit-cycle amplitude could be determined due to saturated signals. Since the preliminary assessment of the measures was conducted for N.2, these rows feature no blank spots. For the remaining tests solely the thresholds and violation-time are given for the preselected measures $μ_{2}$ , $H^{2}$ and $\hat{\bar{α}}$ . Finally, the evaluation rate and the necessary evaluation window in order to sample sufficient data for each method are indicated.

In summary, the preselected measures seem to be, at least to some extent, applicable for a wider set of operating conditions, including various chamber pressures and ROFs, and different combustion devices. While for the Hurst-exponent a common threshold applies for both devices, the threshold for the multifractal parameter has to be adjusted for each combustion device and the second spectral moment even for each test.

Conclusion

The suitability of methods for detecting combustion instability precursors in cryogenic oxygen-natural gas rocket combustion at sub- and supercritical measures were assessed. Different methods from statistical, fractal and recurrence analysis were applied to five different tests from two different combustion devices featuring stable combustion, short-lived bursts of oscillatory combustion, and sustained, high-amplitude combustion instabilities.

The measure based on a standard deviation ( $σ$ ) within the samples of the acoustic pressure sensor was not able to distinguish between intermittent and sustained combustion instabilities.

The measures based on the variance of the AC ( $σ^{2} (A C)$ ) and RQA ( $R D R$ ) were able to distinguish between the states of stable, intermittent and unstable combustion, but were identified as not practicable since the large evaluation window necessary to get meaningful results resulted in delayed detection and higher computational cost.

The fractal measures ( $H^{2}$ and $\hat{\bar{α}}$ ) as well as the second spectral moment ( $μ_{2}$ ) were also able to distinguish the different states of combustion at moderate computational expense. The suitability of the fractal measures is slightly superior compared to the second spectral moment since a common detection threshold for instability onset could be identified independent of the test characteristics for $\hat{\bar{α}}$ and additionally independent on the combustion device for $H^{2}$ .

Finally, it should be noted that the Hurst-exponent $H^{2}$ also fell below the critical threshold during longer bursts of oscillatory combustion. This characteristic may be desirable in certain situations, since even short bursts high-amplitude instability may endanger some combustion devices. Furthermore, the Hurst-exponent relies on the evaluation of the generalized Hurst-exponent for only one order and thus features reduced computational costs.

Footnotes

Acknowledgements

The work is associated with the Franco-German Rocket Engine Stability iniTiative (REST). The authors would also like to thank the crew of the P8 test bench. Special thanks to Alex Grebe for his assistance in preparing and performing the experiments.

Declaration of conflicting interests

The authors declare no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The authors received no financial support for this article’s research, authorship, and/or publication.

ORCID iDs

Jan Martin

Wolfgang Armbruster

Justin S Hardi

References

Harrje

Reardon

. (eds.) Liquid Propellant Rocket Combustion Instability. Washington, DC: NASA SP-194, 1972.

Sutton

Biblarz

. Rocket Propulsion Elements. 8 ed. New York: John Wiley & Sons, 2010.

Rayleigh

. The explanation of certain acoustical phenomena. Nature 1878; 18: 319–321.

Sujith

Pawar

Krishnan

et al. Application of dynamical systems theory and complex systems theory to combustion instability in liquid rocket engines. Technical report, 2020.

Banerjee

Pavithran

Sujith

. Early warnings of tipping in a non-autonomous turbulent reactive flow system: efficacy, reliability, and warning times, 2023. DOI: https://doi.org/10.48550/ARXIV.2306.03486.

Sengupta

Waxenegger-Wilfing

Martin

et al. Forecasting thermoacoustic instabilities in liquid propellant rocket engines using multimodal bayesian deep learning. Int J Spray Combust Dyn 2022; 14: 218–228.

Waxenegger-Wilfing

Sengupta

Martin

et al. Early detection of thermoacoustic instabilities in a cryogenic rocket thrust chamber using combustion noise features and machine learning. Chaos: Interdisciplinary J Nonlinear Sci 2021; 31. DOI: https://doi.org/10.1063/5.0038817.

Martin

Armbruster

Stützer

et al. Flame dynamics of an injection element operated with LOX/h2, LOX/CNG and LOX/LNG in a sub- and supercritical rocket combustor with large optical access. Int J Spray Combust Dyn 2023; 15: 147–165. DOI: https://doi.org/10.1177/17568277231172302.

Martin

Armbruster

Hardi

et al. Experimental investigation of self-excited combustion instabilities in a LOX/LNG rocket combustor. J Propul Power 2021; 37: 944–951. DOI: https://doi.org/10.2514/1.b38289.

10.

Kabiraj

Steinert

Saurabh

et al. Coherence resonance in a thermoacoustic system. Phys Rev E 2015; 92: 042909. DOI: https://doi.org/10.1103/physreve.92.042909.

11.

Sweitzer

Bishop

Genberg

. Efficient computation of spectral moments for determination of random response statistics. In Proceedings of ISMA.

12.

Nayfeh

Balachandran

. Applied Nonlinear Dynamics. Weinheim: Wiley, 1995. DOI: https://doi.org/10.1002/9783527617548.

13.

Cao

. Practical method for determining the minimum embedding dimension of a scalar time series. Physica D: Nonlinear Phenomena 1997; 110: 43–50. DOI: https://doi.org/10.1016/s0167-2789(97)00118-8.

14.

Schinkel

Dimigen

Marwan

. Selection of recurrence threshold for signal detection. Eur Phys J Special Topics 2008; 164: 45–53. DOI: https://doi.org/10.1140/epjst/e2008-00833-5.

15.

Mandelbrot

. How long is the coast of britain? statistical self-similarity and fractional dimension. Science 1967; 156: 636–638. DOI: https://doi.org/10.1126/science.156.3775.636.

16.

Mandelbrot

. The Fractal Geometry of Nature. vol. 1. New York: W.H. Freeman & Co., 1982.

17.

Ihlen

EAF

. Introduction to multifractal detrended fluctuation analysis in matlab. Front Physiol 2012; 3. DOI: https://doi.org/10.3389/fphys.2012.00141.

18.

Abderrahmane

Paquet

. Applying nonlinear dynamic theory to one-dimensional pulsating detonations. Combust Theory Modell 2011; 15: 205–225. DOI: https://doi.org/10.1080/13647830.2010.535566.

19.

Kabiraj

Saurabh

Wahi

et al. Experimental study of thermoacoustic instability in ducted premixed flames: periodic, quasi-periodic and chaotic oscillations. In n3l-Int’l Summer School and Workshop on Non-Normal and Nonlinear Effects In Aero-and Thermoacoustics. p.12.

20.

Zbilut

Zaldivar-Comenges

Strozzi

. Recurrence quantification based liapunov exponents for monitoring divergence in experimental data. Phys Lett A 2002; 297: 173–181.