Numerical study of the linear and non-linear damping in an acoustically forced cold-flow test rig with coupled cavities

Abstract

Quantifying the oscillation amplitude of thermoacoustic instabilities remains a critical and challenging issue, as it is a complex balance between driving and damping processes. The New Pressurized Coupled Cavities (NPCCs) setup designed for the study of acoustic damping is analyzed in this work. It is a cold-flow test rig mimicking the geometry of a liquid rocket engine and equipped with an acoustic forcing device. The chamber 1T mode triggers a strong non-linear harmonic response, while the 1T1L and 1T2L exhibit weak non-linearities. Disturbance energy budgets are used in large-eddy simulations to characterize the damping phenomena with the 1T2L and 1T1L forcing. The correct global damping of the system is retrieved, and local damping contributions are extracted. Then, a non-linear term representing the energy transfer to the harmonics is derived from non-linear acoustics theory. Combined with a linear model, this model correctly retrieves the limit-cycle of the 1T mode.

Keywords

Liquid rocket engines high-frequency combustion instabilities large eddy simulation damping disturbance energy balance

Introduction

Thermoacoustic instabilities remain a fundamental issue when designing a liquid rocket engine (LRE). The extremely high energy delivered by the combustion can lead to destruction of the engine when it is coupled with acoustic modes of the chamber. Quantifying the oscillation amplitude remains a critical and challenging issue, as it is a complex balance between driving and damping processes. This is true in the small-amplitude linear phase of the process, associated with the triggering of the self-sustained oscillation, and in the non-linear one, leading to the so-called limit-cycle related to various saturation and energy transfer processes. To make decisions in the design phase, different low-order approaches have been proposed in recent years.^1–3 In this framework, combustion dynamics models have been developed to represent the unsteady heat release under acoustic solicitation,^1,4 but few studies focus on the damping mechanisms that play a significant role both in the early stage of the instability when the pressure perturbation remains low and in the final stage when the pressure oscillations reach a limit-cycle potentially leading to the system’s destruction. This is the purpose of the present work.

At the beginning of the instability, the pressure oscillations remain low compared to the mean chamber pressure. Assuming small perturbations, the disturbance energy budget⁵ seems an adapted framework to study and quantify the transfer of energy responsible for the growth or damping of an instability. A description of the budget’s implementation in a numerical framework can be found in the PhD thesis of Giauque.⁶ A notable application on the DLR (Deutsches Zentrum für Luft- und Raumfahrt) research combustor model ‘D’ (BKD)^7,8 was presented by Urbano and Selle.⁹ This work revealed the importance of taking into account hydrodynamic phenomena such as vortex shedding in a complex reactive case. However, the closure of the budget on complex cases remains very challenging.

When the instability grows, the small perturbations hypothesis holds no more, preventing the use of the disturbance energy budget. The damping is then mostly driven by non-linear energy transfer towards harmonics of the excited frequency. Models have been proposed to represent this phenomenon, as the two-mode approximation presented by Culick.¹⁰

In this article, previous experimental and numerical results on the New Pressurized Coupled Cavities (NPCCs) rig are first recalled. The NPCC rig⁴ is a cold-flow test rig that mimics the geometry of a rocket engine combustion system with a dome, three injectors and a chamber. It can be acoustically forced via two nozzles at the end of the chamber excited by the Very High Amplitude Modulator (VHAM),¹ a perforated wheel mechanism that was used in various experiments to successfully modulate the system’s eigenmodes. Then, a specific focus is made on the linear damping mechanisms, present in the early stage of an instability. Disturbance energy budgets⁵ are applied to large eddy simulations (LESs) of the NPCC rig, leading to a local quantification of the linear damping mechanism as well as a physical understanding of the phenomena involved. The next section focuses on non-linear energy transfer towards harmonics, a key damping phenomenon dominant in the last stages of the instabilities. A damping model is derived from non-linear acoustics theory and applied with success to a simple one-dimensional (1D) numerical case. Finally, the last section uses both linear and non-linear components of the damping to give an accurate model of the limit-cycle found in the most complex NPCC experimental case.

Limit-cycle observed during forcing of the NPCC test rig

The NPCC test rig is a non-reactive coupled-cavity configuration running with pure air and operating at $p_{0} =$ 3.5 bar. A modulation device (VHAM) permits the excitation of selected eigenmodes. In this section, the main feature of the simulation and results are recalled (the reader is referred to Marchal et al.¹¹ for more details). The leading phenomena are summarized to drive the rest of the analysis.

Forcing of the NPCC test rig

The numerical domain is presented in Figure 1. The mass flow rate is constant for all the following analyses. The pressurized air goes through the three injectors with a velocity of around 10 m/s, corresponding to $R e \approx 14, 000$ . In the following, the pressure fluctuations are evaluated at the probe HFc1, located in a pressure antinode for all the considered eigenmodes.

Figure 1.

Numerical domain for the LES simulation of the NPCC rig. The acoustic modulation is provided by the VHAM.¹ The position of the HFc1 pressure sensor is displayed. LES: low eddy simulation; NPCC: New Pressurized Coupled Cavity; VHAM:Very High Amplitude Modulator.

The modulation of the rig is studied for three eigenmodes: 1T (chamber-only mode), 1T1L (in the chamber, coupled with the 1T mode of the dome) and 1T2L (chamber-only mode). The forcing frequencies are shown in Table 1.

Table 1.

Modulation frequencies of the NPCC rig.

Eigenmode	1T	1T1L	1T2L
Frequency (Hz)	1226	1468	2036
Type	ChO	CC	ChO

NPCC: New Pressurized Coupled Cavity; ChO: chamber-only; CC: coupled-cavity.

It was demonstrated in Marchal et al.¹¹ that the forcing of the VHAM could be modelled within the LES using only geometrical assumptions, closing the system. The forcing amplitude $Δ p^{'} = 0.45 bar$ depends only on the mean pressure in the chamber and the distance between the VHAM and the nozzles.

LES results

The modulated cases are simulated using LES with the AVBP solver.^12,13 The resulting pressure perturbation $p^{'}$ plotted against the experimental curve are recalled in Figure 2.¹¹ The simulations are in excellent agreement with the experiments. Notably, the spectral content of the limit-cycle seems well-reproduced for the 1T mode which presents a strong harmonic content. The associated spectrogram (Figure 3) shows the emergence of the harmonics one after the other. The limit-cycle is clearly the result of a competition between forcing and damping, the former being imposed following a geometrical reasoning without underlying assumptions as described in Marchal et al.¹¹ It follows that the LES gives a good representation of the damping phenomena occurring in the NPCC rig.

Figure 2.

Acoustic pressure fluctuations in the chamber (HFc1 probe), large eddy simulation (LES); experiment; top: 1T, mid: 1T1L, bot: 1T2L. Both signals are high-pass filtered with a cut-off frequency of 400 Hz.

Figure 3.

Spectrogram of the pressure response in the chamber, 1T mode, large eddy simulation (LES) results, corresponding to the signal displayed Figure 2 top.

Damping in the NPCC rig

It is our understanding that the acoustic waves in the NPCC test rig are damped by three main mechanisms:

The interaction between the acoustic waves and the flow at different locations of the rig: injector exit, dome, and nozzles. This is extensively studied using disturbance energy budgets in the next section and will be shown to be responsible for a small fraction of the overall damping in experimental cases.

The viscous and thermal dissipations at the walls, studied in Searby et al.¹⁴ To capture these phenomena in the LES, the mesh should precisely represent the viscous and thermal boundary layer. It is not in the scope of this work, as it would imply a huge central processing unit cost due to the heavy mesh and the long averaging time required for good statistics. This contribution is left for future studies.

The non-linear energy transfer towards harmonics. This is the subject of the last section of this paper and will be shown to be the main contributor to the overall damping. As the acoustic pressure’s spectral content seems well captured by the LES (Figure 2), it is reasonable to believe that the simulation reproduces this cascade phenomenon with accuracy.

Linear damping characterized using disturbance energy budgets

The disturbance energy budget⁵ is well-adapted to study the interaction between the acoustics and the flow. However, it is only valid in the assumption of small perturbations, which is not the case of the experiments and simulations presented in the previous section ( $(Δ p / p_{0}) \approx 13 %$ ). Thus, additional LESs were performed on the NPCC geometry reducing the forcing amplitude to 10% of the experimental amplitude, that is, $(Δ p / p_{0}) \approx 1.3 %$ (1T mode: 1T_lo, 1T1L mode: 1T1L_lo and 1T2L mode: 1T2L_lo, see Table 2).

Table 2.

LES computations of the NPCC rig.

Computation	Forcing amplitude	Presence of
	(bar)	harmonics
1T	0.45	Yes
1T_lo	0.045	Yes
1T1L	0.45	yes
1T1L_lo	0.045	No
1T2L	0.45	Yes
1T2L_lo	0.045	No

LES: large eddy simulation; NPCC: New Pressurized Coupled Cavities.

For the 1T1L_lo and 1T2L_lo cases, the harmonic content disappears as illustrated for case 1T1L_lo on Figure 4. Therefore, it is expected that these computations represent fully linear cases, well-suited for the application of the disturbance energy budget. In addition, the computation presenting the non-modulated situation will also be used in this part and referred to as NM.

Figure 4.

Spectrograms for the numerical cases 1T1L (top) and 1T1L_lo (bottom).

Note: Remarkably, despite the reduced modulation, the 1T_lo case still contains a non-negligible contribution of the harmonics. This is a very interesting feature that is left for future investigations.

Disturbance energy budgets in the LES code

The disturbance energy budget was implemented in the LES code as presented in Giauque.⁶ Making some assumptions adapted to the case of the NPCC rig (adiabatic case, single species, and no chemical reaction), it can be formulated and integrated on a volume $V$ with boundary $δ S$ as:

\int_{V} \frac{\partial E_{2}}{\partial t} d V = \int_{V} D_{2} d V + \int_{δ S} W_{2} \cdot n d S

(1)

with

{\begin{matrix} E_{2} = \frac{p^{' 2}}{2 \bar{ρ} {\bar{c}}^{2}} + \frac{1}{2} \bar{ρ} u^{' 2} + ρ^{'} \bar{u} \cdot u^{'} \\ W_{2} = (p^{'} + \bar{ρ} \bar{u} \cdot u^{'}) (u^{'} + \frac{ρ^{'}}{\bar{ρ}} \bar{u}) + \bar{(ρ u)^{'} H^{'}} - T \bar{(ρ u)^{'} s^{'}} \\ D_{2} = D_{ζ} + D_{ψ} \end{matrix}

(2)

and

\begin{aligned} {\begin{matrix} D_{ζ} & = - (ρ u)^{'} \cdot ζ^{'} - \bar{(ρ u)^{'} \cdot ζ} \\ D_{ψ} & = (ρ u)^{'} \cdot ψ^{'} + \bar{(ρ u)^{'} \cdot ψ^{'}} \end{matrix} \end{aligned}

where

() = \bar{()} + ()^{'}

p

is the pressure,

u

the velocity,

ρ

the density,

H

the enthalpy,

T

the temperature,

ζ = (\nabla \times u) \times u

and

ψ = (1 / ρ) (\partial τ_{i j} / \partial x_{i})

with

τ_{i j}

the shear-stress tensor.

In the formulation of the disturbance energy $E_{2}$ presented in equation (2), one recognizes the classical acoustic energy density (the term featuring entropy disturbances being neglected here). In the flux term $W_{2}$ , the first two brackets encompass contributions from acoustic waves and their interaction with the mean flow. The last two terms represent the interaction between acoustic and thermal energies. The source/damping term $D_{2}$ features the contributions of the vorticity in $D_{ζ}$ and the dissipation by shear stress in $D_{ψ}$ .

Choice of a reference solution

When computing disturbance energy budgets, the choice of the reference solution is crucial. It is used to get averaged values $\bar{()}$ and in the computation of perturbed values $()^{'}$ . Good statistics are required to get a clean reference solution. For the NM case, the reference is taken as the average over 200 ms starting after a steady state is achieved. For the modulated cases (1T1L_lo and 1T2L_lo), the reference is the average over 200 periods of oscillations. The perturbations are then computed with respect to this reference. It was not possible to get rid of all the ultra-low frequency effects, any remaining convergence issues will be mentioned when encountered.

Closure of the budget on linear NPCC simulations

A mapping of the rig enables the local analysis (1, dome; 2, Injectors; 3, jets; 4, chamber; and 5, nozzles) of the disturbance energy budget. In the limit-cycle regime, the disturbance brought by the modulation is compensated (over each period) by the damping and the change in flux distribution:

\int_{V} \bar{\frac{\partial E_{2}}{\partial t}} d V = 0

(3)

leading to

\int_{V} \bar{D_{ψ}} d V + \int_{V} \bar{D_{ζ}} d V + \int_{δ S} \bar{W_{2} \cdot n} d S = 0

(4)

The closure of the budget is verified in each case by checking the sum of these three terms. Figure 5 displays the averaged budget for cases 1T2L_lo and NM for each area of the mapping and for the whole domain (a similar behaviour is observed for 1T1L_lo compared to 1T2L_lo). The closure error (i.e. the sum of the terms that should be zero for a perfectly closed budget) is roughly 10% of

D_{ζ}

when looking at the full domain which is extremely satisfactory taking into account the relative complexity of this 3D case. The vorticity (

D_{ζ}

, black rectangles) and diffusion (

D_{ψ}

, dotted rectangles) terms are dominant on the averaged budget and nearly compensate themselves in each subdomain. The averaged fluxes are one order of magnitude below. Thus, the budget closes reasonably well on the whole system and in each zone of the mapping.

Figure 5.

Averaged disturbance energy budget for all subdomains and the full domain in the permanent regime. 1T2L_lo and NM case are displayed. Each histogram represents the three terms for the corresponding area of the mapping: the dome, the three injectors taken separately, the three jets taken separately and all together, the right end side of the chamber, the full chamber, the two nozzles taken separately and the full domain.

Local damping contribution

To isolate the effect of the modulation, it is useful to observe the difference between modulated and non-modulated cases. This is possible thanks to the linear nature of the phenomena encountered with the lower amplitude modulation. Figure 6 presents the difference 1T2L_lo – NM. A relatively similar behaviour is observed for 1T1L_lo – NM (the differences between the two modulated cases are developed later on). The acoustic modulation increases both $| D_{ζ} |$ (especially in the jets) and $| D_{ψ} |$ (especially in the dome). This proves that the jet-acoustic damping mechanism is linked to the production of vorticity ( $D_{ζ}$ , which is always positive) and its subsequent dissipation ( $D_{ψ}$ , always negative).

Figure 6.

Difference of the budget’s terms 1T2L_lo – NM. Each histogram represents the three terms for the corresponding area of the mapping: the dome, the three injectors taken separately, the three jets taken separately and all together, the right end side of the chamber, the full chamber, the two nozzles taken separately and the full domain.

To highlight the difference between chamber-only (1T2L_lo) and coupled-cavity (1T1L_lo) modulated cases, the budget difference 1T2L_lo – 1T1L_lo is displayed in Figure 7. Doing the difference as such implies that a positive value for the resulting term indicates a greater value of this term for 1T2L_lo (chamber-only case). Beware that the dissipation term $D_{ψ}$ being negative, a higher value of $D_{ψ}$ means less dissipation. Comparing the two modulated cases is trickier than it may seem. Indeed, the modulation in the 1T2L_lo case proved to be slightly more effective than the coupled-cavity case (1T1L_lo ), as seen on the overall acoustic flux (small square pattern) entering the domain which is greater (‘Full domain’ histogram). We could associate this ‘efficiency’ effect with a difference in the nozzle flow for each mode. This question is still under investigation and beyond the scope of the present paper.

Figure 7.

Difference of the budget’s terms 1T2L_lo – 1T1L_lo. Each histogram represents the three terms for the corresponding area of the mapping: the dome, the three injectors taken separately, the three jets taken separately and all together, the right end side of the chamber, the full chamber, the two nozzles taken separately and the full domain.

Several features can be observed when analysing Figure 7:

For the 1T2L_lo (chamber-only) case, the acoustic flux is oriented preferentially towards the central jet.

For the 1T1L_lo (coupled-cavity) case, on the contrary, and even though the mode pattern in this area of the chamber is similar for both cases, the lateral jets get more acoustic flux than the central one, tending to indicate an interaction with the feeding part of the domain for that case.

$D_{ψ}$ in the injectors is lower in the chamber-only case, indicating that more acoustic damping takes place in the injectors in the coupled-cavity one, probably for the same reason;

There is a discrepancy between the two lateral injectors suggesting a small convergence issue, that should be carefully looked at, but we strongly believe that this does not alter our main conclusions.

Quantifying the damping phenomena means, in a linear framework, finding a damping coefficient

σ

. Figure 8 illustrates the decay of the envelope of the disturbance energy when the modulation is suddenly stopped. The damping coefficient

σ

appears in the expression modelling the decaying part represented in red¹ :

E (t) = E_{NM} + (E_{modulated} - E_{NM}) e^{σ t}

(5)

Figure 8.

Assumption on the behaviour of disturbance energy when cutting the modulation.

Its derivative is

\frac{\partial E (t)}{\partial t} = σ (E_{modulated} - E_{NM}) e^{σ t} = σ (E (t) - E_{NM})

(6)

It follows

σ = \frac{1}{E (t) - E_{NM}} \frac{\partial E (t)}{\partial t}

(7)

In the following, it is convenient to split the flux term for each subdomain into (a) the part coming from the forcing

W_{forcing}

and (b) the part coming from the flow

W_{flow}

. Applied to our specific case, the decaying phase (without modulation,

W_{forcing} = 0

) would translate to:

\begin{aligned} \int_{V} \frac{\partial E}{\partial t} d V |_{decay} & = \int_{V} D_{ζ} d V |_{decay} + \int_{V} D_{ψ} d V |_{decay} \\ + \int_{δ S} W_{flow} \cdot n d S |_{decay} \end{aligned}

(8)

On the other hand, the situation in limit-cycle (equation (4)) gives:

\begin{aligned} \int_{V} \bar{D_{ψ}} d V |_{limit - cycle} + \int_{V} \bar{D_{ζ}} d V |_{limit - cycle} \\ + \int_{δ S} \bar{W_{flow} \cdot n} d S |_{limit - cycle} = \int_{δ S} \bar{W_{forcing} \cdot n} d S |_{limit - cycle} \end{aligned}

(9)

In the following, it is assumed that

σ

is constant in the decaying phase. It will be evaluated at the instant the modulation is stopped. Assuming equality of the different terms at the instant the modulation is switched off, a reasonable hypothesis, one gets:

\int_{V} \frac{\partial E}{\partial t} d V |_{decay} = \int_{δ S} \bar{W_{forcing} \cdot n} d S |_{limit - cycle}

(10)

The denominator of equation (11) can be expressed in a control volume

V

at the instant the modulation is stopped:

\int_{V} E (t) - E_{NM} d V = \int_{V} E_{modulated} - E_{NM} d V

and in our case evaluated simply by integrating the disturbance energy over

V

in modulated and non-modulated cases. In steady states, it is constant in both cases. It comes finally:

σ = \frac{\int_{δ S} \bar{W_{forcing} \cdot n} d S |_{limit - cycle}}{\int_{V} E_{modulated} - E_{NM} d V}

(11)

Linking

σ

in a control volume

V

with the incoming fluxes

W_{forcing}

and the surplus of acoustic energy coming from the modulation

\int_{V} E_{modulated} - E_{NM} d V

. According to Searby et al.,¹⁴ the acoustic pressure damping

σ_{p}

is one half of

σ

, leading to:

σ_{p} = \frac{\int_{δ S} \bar{W_{forcing} \cdot n} d S |_{limit - cycle}}{2 (\int_{V} E_{modulated} - E_{NM} d V)}

(12)

The challenge then lies in isolating

W_{forcing}

. In the NPCC rig, in a first approximation, the flux coming from the forcing can be associated with the acoustic flux

\int_{δ S} \bar{p^{'} u^{'}} \cdot n d S

. As displayed in Figure 9, it forms the majority of the flux coming from the nozzles into the domain (red patch) in the modulated cases.

Figure 9.

Fluxes entering and leaving the top nozzle in case 1T2L_lo. A positive flux means energy incoming into the nozzle subdomain, a negative flux translates to energy leaving the subdomain. The energy transmitted to the New Pressurized Coupled Cavities (NPCC) rig is then the negative flux on the left side of the nozzle.

Isolating $W_{forcing}$ could be done by subtracting $\int_{δ S} \bar{p^{'} u^{'}} \cdot n d S$ before and after stopping the modulation. In our case, it is unnecessary as the acoustic flux of the non-modulated case is very low compared to one of the modulated cases (see Table 3).² This is true over the full domain as well as locally. Thus, the approximation:

\int_{δ S} \bar{W_{forcing} \cdot n} d S |_{limit - cycle} \approx \int_{δ S} \bar{p^{'} u^{'}} \cdot n d S |_{limit - cycle}

(13)

is made in the following, locally and for the whole domain.

Table 3.

Comparison of the total flux term and the incoming acoustic flux in the full domain for cases NM, 1T1L_lo and 1T2L_lo. The difference in the total flux betwen 1T1L_lo and 1T2L_lo is traced to the flow discrepancies inside the nozzles between these two cases, with more outgoing flux in the 1T1L_lo case.

Case	$\int_{δ S} \bar{W_{2} \cdot n} d S$	$\int_{δ S} \bar{p^{'} u^{'} \cdot n} d S$
	(mW)	(mW)
NM	3	0.7
1T2L_lo	7	29
1T1L_lo	1	29

Applied to the mapping of the rig, equation (12) then gives local damping contributions for each subsection (Figure 10). A large part of the damping in the coupled-cavity case (1T1L_lo ) comes from the dome, whereas this part is reduced for the chamber-only case (1T2L_lo). In the coupled-cavity case, the section of the chamber featuring the jets shows strong damping (except for the top jet, but we believe it is a convergence issue), whereas the end of the chamber shows $σ_{p} > 0$ with a similar amplitude, making the chamber overall neutral with regard to the acoustic damping. The central injector is a source of acoustic power whereas the lateral injectors feature a slight damping.

Figure 10.

Local damping contributions found with the fluxes of the disturbance energy budgets; top: 1T1L_lo, bot: 1T2L_lo. Red stands for $σ > 0$ , blue for $σ < 0$ . Absolute amplitude is represented by opacity, transparent being 0.

For the chamber-only case, the lateral jets are damping acoustics, whereas the central jet acts as a source of acoustic power. Overall the chamber is a place of damping, complementing the dominant damping effect of the dome. The injectors are all slightly acting as sources of acoustic power.

The global damping found using this methodology is close to the global damping extracted from the LES with a simple linear model as displayed in Table 4, especially for the coupled-cavity case.

Table 4.

Global acoustic damping $σ_{p}$ found using disturbance energy budget and linear fit on large eddy simulation (LES) results.

Case	1T1L_lo	1T2L_lo
$σ_{p}$ [ $s^{- 1}$ ] using disturbance energy budget	$-$ 3.1	$-$ 3.0
$σ_{p}$ [ $s^{- 1}$ ] using a linear fit on LES results	$-$ 3.0	$-$ 6.0

Applying these damping coefficients on a simple linear model permits to properly recover the limit-cycle obtained numerically in the low amplitude (linear) cases (1T1L_lo, 1T2L_lo).¹¹ However, when applied to the actual experimental cases (i.e. with larger modulation amplitudes), the model predicts a larger amplitude at the limit-cycle, especially for case 1T (see Figure 11), suggesting that a physical phenomenon is not accounted for. This aspect is discussed in the next section.

Figure 11.

Comparison between the large eddy simulation (LES) signal at probe HFc1 (opacity = 70%, bandpass filtered around the 1T frequency) and fitted linear model. Excitation amplitude: 0.45 bar. Case: 1T. Reproduced from Marchal et al.¹¹

Derivation of a non-linear model for the energy transfer to harmonics

The pressure field for the 1T case when modulated at the experimental amplitude reveals in the LES the formation of a pressure wave bouncing on the walls of the chamber at the modulation frequency (Figure 12). This explains the peculiar shape of the signal observed at the probe location in the chamber (Figure 2 top). In this section, the transition from a standing mode to a shock-like behaviour is detailed and a quantitative model is derived for the energy transfer from the fundamental to the first harmonic frequency.

Figure 12.

Large eddy simulation of the 1T case with a modulation amplitude extracted from the experiment. The pressure field is observed at three different phases taking the modulation frequency as reference.

Weak shock formation in a high-amplitude travelling wave

The following development is inspired by Garrett’s book.¹⁵ The logic was developed for travelling waves, however, we believe that it can be used also for standing waves, which are the object of the present work. Consider a wave of high-amplitude travelling in a medium, for example, air, down a duct of constant cross-section. If the amplitude is not sufficiently damped by thermoviscous effects, it will influence the local sound speed and induce a steepening of the wave. This disturbance to the sound speed is characterized with the help of the local particle velocity $ν$ , which depends on the amplitude of the wave. The local amplitude-dependent sound speed is then defined by:

c (ν) = c_{0} + (1 + \frac{\partial c}{\partial y} \frac{\partial y}{\partial ν}) \equiv c_{0} + Γ ν

(14)

where

c_{0}

is the equilibrium sound speed,

y

the amplitude and

Γ

is called the Grüneisen parameter. It is useful to introduce the acoustic Mach number associated with the local velocity perturbation

ν

M a_{ac} = \frac{ν}{c_{0}} \approx \frac{p^{'}}{γ p_{0}}

(15)

using the travelling waves’ property

p^{'} = ρ_{0} c_{0} ν

(extended in a first approximation to standing waves).

After some distance, the travelling wave will eventually form a shock. The shock inception distance $D_{S}$ is expressed with regard to the wavelength $λ$ of the wave:

D_{S} = \frac{λ}{2 π Γ M a_{ac}}

(16)

A Fourier decomposition¹⁶ of the wave progressively deformed after a travelling distance

x

using the scaling

λ = \frac{x}{D_{S}}

gives the following harmonics’ amplitudes

C_{m}

| C_{m} | = \frac{2}{λ m} J_{m} (m λ) for λ < 1 and m = 1, 2, 3, \dots

(17)

where

J_{m}

is the Bessel function of the first kind.

In the first approach, only the energy transfer from the fundamental to the first harmonic frequency is considered. Indeed, the spectrograms (Figure 2) show that in our case the fundamental and the first harmonic frequencies retain most of the energy. The quantity of interest is then $\partial | C_{1} | / \partial t$ .

Considering equation (17) and the fact that $λ < 1$ , the first two terms of the development of $| C_{1} |$ give a good approximation of its value³ :

| C_{1} | \approx \frac{2}{λ} (\frac{λ}{2} - \frac{λ^{3}}{16}) = 1 - \frac{λ^{2}}{8}

(18)

Leading to:

\frac{\partial | C_{1} |}{\partial t} = - \frac{1}{4} \frac{\partial λ}{\partial t} = - \frac{1}{4 D_{S}} \frac{\partial x}{\partial t}

(19)

Using the definition of

D_{S}

given in equation (16) and the fact that the wave travels at the speed of sound

c_{0}

\frac{\partial | C_{1} |}{\partial t} = - \frac{1}{4} Γ ω M_{ac}

(20)

with

ω = 2 π f

the pulsation.

The variation of the pressure perturbation amplitude $A_{p}^{1}$ of the fundamental is then expressed as:

\frac{\partial A_{p}^{1}}{\partial t} = ρ_{0} c_{0} ν \frac{\partial | C_{1} |}{\partial t} \approx - p^{' 2} \frac{Γ ω}{4 ρ_{0} c_{0}^{2}}

(21)

showing a quadratic link between the variation of the perturbation and its amplitude, with an analytic coefficient. In the following, let

b_{1} \equiv Γ ω / 4 ρ_{0} c_{0}^{2}

Application to a purely numerical 1D case

The model is first tested on a simple 1D canonical configuration. The setup consists of a 1 m tube acoustically forced at one end (with a totally reflective boundary condition – acoustic waves enter the system but cannot exit) with the other end being a wall. The domain is filled with pure oxygen at 1 bar and 1050 K, leading to $c_{0} \approx 600 m \cdot s^{- 1}$ and $ρ_{0} \approx 0.37 kg \cdot m^{- 3}$ . Only Euler equations are solved and the system is stabilized using artificial viscosity when the shock is formed. The excitation frequency is set to 300 Hz, so that the first longitudinal mode is excited. After 0.3 s of simulation, the forcing is stopped so that only non-linear damping is present in the system.

This case exhibits roughly the same spectral behaviour as the NPCC rig in the growing phase (Figure 13).

Figure 13.

Spectrogram obtained from the signal of a probe located at the inlet in the LES of the simplified 1D case with excitation $u^{'} = 0, 5 {m.s}^{- 1}$ . LES: large eddy simulation; 1D: one-dimensional.

Thus, it is interesting to test the non-linear model in this case which gets rid of most of the complex phenomena studied before (no viscosity, no jets, and no coupled cavities). The modelling of this case also incorporates a linear forcing term to account for the energy gained during the modulation.

The quadratic non-linear model predicts the limit-cycle of the mode with only a slight underestimation and retrieves well the slope of the amplitude decrease when the modulation is stopped (Figure 14).

Figure 14.

Temporal evolution of the amplitude of the 1T mode’s pressure perturbation in the simplified one-dimensional (1D) case with $u^{'} = 0.5 m \cdot s^{- 1}$ . When $t < 0.3$ s, the system is modulated. At $t = 0.3$ s, the modulation is stopped. Thick line: 1D simulation; thin line: corresponding model including non-linear effects.

Application to the NPCC test rig

It was shown in Marchal et al.¹¹ that a linear model for the global damping was able to retrieve the evolution of the pressure at the beginning of the VHAM modulation with satisfactory accuracy. When harmonic contents appear, the numerical results differ from the linear model and converge towards a limit-cycle below the one predicted by the linear model. Extending the model with the non-linear contribution analytically derived above, the modelled amplitude is expected to give a closer approximation of the limit-cycle observed experimentally.

The reduced-order model developed by Marchal et al.¹¹ and extended with this non-linear contribution can be applied to the NPCC case. The equation for the amplitude of the fundamental mode now writes:

{\ddot{η}}_{1} (t) + 2 (α_{lin} + b_{1} * | η_{1} (t) |) \dot{η_{1}} (t) + ω_{0}^{2} η_{1} (t) = S_{exc}^{1} (t)

(22)

with

α_{lin}

the linear component of the damping,

b *_{1} \approx \frac{b_{1}}{0, 64}

, with the factor 0.64 arising from the switch to amplitude formulation.⁴ For the excitation term

S_{exc}^{1}

, the parameters are the same as in Marchal et al.¹¹ The linear damping coefficient

α_{lin}

is taken from Table 2. For the non-linear part,

Γ_{air} = 6 / 5

ρ_{0} = 4.14 kg \cdot m^{- 3}

and

c_{0} = 340 m \cdot s^{- 1}

are used for the computation of

b_{1}

Tested against LES data, the model produces very satisfying results for the 1T mode (Figure 15). In contrast with the purely linear model (Figure 11), the predicted amplitude of $p^{'}$ stays close to the LES data from the start of the modulation to the limit-cycle regime. The model underpredicts the steady-state regime amplitude by about 13%. This slight discrepancy could indicate that (a) the model induces a stronger damping compared to the result of the LES and/or (b) the modelled excitation source term $S_{exc}^{1}$ is a bit too weak. An explanation for the stronger damping may come from the non-ideal geometry of the rig (the shock wave is bouncing on the wall and interacting with the jets) compared to the travelling wave’s theory from which the model was derived.

Figure 15.

Comparison between large eddy simulation (LES) signal in the chamber and reduced-order model including linear and non-linear components, 1T mode. Excitation amplitude: 0.45 bar. Top: from the start of the modulation. Bottom: zoom on limit-cycle.

For the 1T1L and 1T2L cases, the results are not as good (Figure 16). The model underpredicts the amplitude of the pressure perturbation by a large value (around 50%). This is probably due to the fact that the modes’ shapes include a longitudinal component, which prevents their non-linear development from being described by the 1D analogy presented above. On the contrary, the 1T mode could be successfully seen as a 1D wave bouncing on the walls. This also explains why the 1T modulation produces non-linearities so easily while the other two modulation cases remain mostly linear.

Figure 16.

Comparison between large eddy simulation (LES) signal in the chamber and reduced-order model including linear and non-linear components, case 1T1L (top), case 1T2L (bot). Excitation amplitude: 0.45 bar.

Thus, the non-linear model developed here should be restricted to the first approach to transverse oscillations. This is still an interesting result considering the tendency of LREs’ combustion chambers to exhibit 1T oscillations.¹⁰

Conclusion

This work explored the question of acoustic damping in LREs using a fundamental approach. By getting rid of combustion in the NPCC rig, the damping phenomena that would otherwise be hard to detect crushed by the flame–acoustics coupling were identified and quantified.

In linear cases, or at the very start of the instabilities where the amplitude remains low, the disturbance energy budget seems a promising tool. The application presented here shows the possibility of conducting this analysis on a 3D LES, with all the relevant terms and with a closed budget, over a long physical time simulated. By extracting the fluxes from the budget, a good approximation of the damping was found. The next step is to extend this methodology to a reactive case.

To account for the wave-steepening phenomenon that leads to the formation of a weak shock, observed in the NPCC experiment and the corresponding LES for the 1T mode, a non-linear model was developed. It seems very promising to help in the determination of the limit-cycle amplitude in LREs, where very strong fluctuations lead to the saturation of the fundamental frequency. The extension of this model to multi-dimensional modes is in progress.

It is the author’s understanding that the present numerical simulations are unable to capture the damping associated with the viscous and thermal boundary layers, due to a mesh that is too coarse at the boundaries. This does not hinder the reproduction of the experimental cases where the non-linear phenomena are preponderant. However, it would be interesting to simulate the rig with a very fine mesh to capture the damping due to the boundary layers in purely linear cases. Towards the same goal, it could be feasible to derive a wall function that takes into consideration the acoustic–boundary layer interaction in the same idea as presented in Berggren.¹⁷ This could reduce the computational cost of such simulations significantly.

Finally, it could be useful to go back to the bench and find a way to excite the NPCC rig with a low amplitude, in order to get experimental results in the purely linear regime, as the VHAM proved to be too powerful in that regard. Implementing speakers on the bench could be a way, and this is considered for future work.

Footnotes

Acknowledgements

This work was granted access to the HPC resources of IDRIS made available by Grand Equipement National de Calcul Intensif (GENCI) under allocation A0042B06176. A part of this work was performed using High Performance Computing resources from the Mesocentre Computing Center of CentraleSupelec and Ecole Normale Superieure Paris-Saclay supported by CNRS and Region Ile-de-France. We would like to thank CERFACS for the access to the software AVSP and AVBP, as well as IFPEN for the latter.

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

Thomas Schmitt

Sébastien Ducruix

Notes

References

Mery

et al. Injection coupling with high amplitude transverse modes: Experimentation and simulation. C R Mécan. 2009; 337: 426–437.

Laurent

et al. A novel modal expansion method for low-order modeling of thermoacoustic instabilities in complex geometries. Combust Flame 2019; 206: 334–348.

Fougnie

Schmitt

Ducruix

. Reduced order modeling for the prediction of combustion instabilities in liquid rocket engines: coaxial injection response and acoustic damping modeling. In: Aerospace Europe Conference 2023 – 10^{T
H} EUCASS – 9^{T
H} CEAS, July 2023, Lausanne, Switzerland.

Gonzalez-Flesca

et al. Reduced order modeling approach to combustion instabilities of liquid rocket engines. AIAA J 2018; 56: 4845–4857.

Brear

et al. Disturbance energy transport and sound production in gaseous combustion. J Fluid Mech 2012; 707: 53–73.

Giauque

. Fonctions de transfert de flamme et énergies des perturbations dans les écoulements réactifs. PhD Thesis, Institut National Polytechnique de Toulouse, 2007.

Gröning

, et al. Influence of hydrogen temperature on the stability of a rocket engine combustor operated with hydrogen and oxygen. CEAS Space J 2017; 9: 59–76.

Gröning

et al. Injector-driven combustion instabilities in a hydrogen/oxygen rocket combustor. J Propul Power 2016; 32: 1–14.

Urbano

Selle

. Driving and damping mechanisms for transverse combustion instabilities in liquid rocket engines. J Fluid Mech 2017; 820: 1–12.

10.

Culick

FEC

. Unsteady motions in combustion chambers for propulsion systems. Technical report. NATO, 2006.

11.

Marchal

Schmitt

Ducruix

. Numerical characterization of the damping in a cold-flow experiment with coupled cavities. In: Space Propulsion 2020+1 100% Virtual Conference, 17 to 19 March, 2021.

12.

Schonfeld

Rudgyard

. Steady and unsteady flows simulations using the hybrid flow solver avbp. AIAA J 1999; 37: 1378–1385.

13.

Moureau et al

. High-order methods for DNS and LES of compressible multicomponent reactive flows on fixed and moving grids. J Comput Phys 2005; 202: 710–736.

14.

Searby et al

. Prediction of the efficiency of acoustic damping cavities. J Propul Power 2008; 24: 516–523.

15.

Garrett

. Understanding acoustics. 2nd ed. Springer, 2020.

16.

Hargrove

. Fourier series for the finite amplitude sound waveform in a dissipationless medium. J Acoust Soc Am 1960; 32: 511–512.

17.

Berggren

Bernland

Noreland

. Acoustic boundary layer as boundary conditions. J Comput Phys 2018; 371: 633–650.