the heat-release linear response to acoustic modulation at the fuel inlet.

This analysis has the potential to guide future development of theoretical
models for LRE flame response, which are expected to account for the 3 major
contributions identified. Additionally, the FTF data is available, and
expected to be embedded in acoustic solvers for the prediction of
thermoacoustic instabilities in multi-injectors LREs.

###  Chapter 7 Nonlinear dynamics in a doubly-transcritical LO2/LCH4 coaxial
jet flame

The previous chapter investigated the linear response of a LO2/LCH4 LRE flame,
subjected to small amplitude acoustic modulation imposed at the fuel inflow.
This perturbation regime is typical of the onset of a thermoacoustic
instability, but significantly differs from an established large-amplitude
limit-cycle, which may trigger nonlinear phenomena in the flame response. The
present chapter therefore extends the analysis initiated in the previous one,
by assessing the nonlinear dynamics of the same doubly-transcritical flame in
the Mascotte test-rig, subjected to strong fuel inflow acoustic modulation,
over a wide frequency range spanning from $\mathcal{O}(1\ \textrm{kHz})$ to
$\mathcal{O}(20\ \textrm{kHz})$. It starts with the identification of two key
nonlinear features affecting the heat-release, namely the response intensity
saturation and the generation of higher harmonics. A local Flame Describing
Function is computed and commented. The vortex nonlinear dynamics, which are
strongly correlated to the flame response, are then discussed. Finally, the
physical mechanisms contributing to the heat-release response saturation are
identified, and their respective nonlinear behaviors are compared.

####  7.1 Introductory remarks

The numerical setup is identical to the one of Chapter 6. The geometry is that
of the single-injector Mascotte test rig (see in Fig. 5.6 and Fig. 6.1),
operating in the doubly-transcritical regime. Both propellants (LO2 and LCH4)
are injected as dense cryogenic fluids into the combustion chamber at $P=75$
bar (see Tab. 5.1). A series of LES are performed to impose harmonic acoustic
modulation of the fuel inflow, at 16 distinct frequencies $f$ ranging from
$\mathcal{O}(1\ \textrm{kHz})$ to $\mathcal{O}(20\ \textrm{kHz})$. This
frequency range is representative of thermoacoustic instabilities usually
observed in LREs. The target fluctuating inlet velocity, imposed thanks to the
methodology of [239], writes:

$\displaystyle u_{CH4}^{\prime}=Au_{CH4}^{0}\sin(2\pi ft)$ (7.1)

Nondimensionalized forcing frequency and axial location are defined based on
the flame length and the fuel inlet bulk velocity:
${f}^{*}=fL_{f}/u_{CH4}^{0}$, ${x}^{*}=x/L_{f}$. Several other
nondimensionalized quantities are introduced later and detailed when
necessary. Simulations are run for at least 20 oscillations periods to ensure
the dissipation of transients and the convergence of Fourier statistics. In
addition, every time signal is systematically detrended before computing
Fourier coefficients, in order to mitigate the effect of intrinsic low-
frequency oscillations detailed in Sec. 6.2.3. The notation $\hat{v}$ refers
to the Fourier coefficient of the quantity $v$ at the forcing frequency
${f}^{*}$, unless specified otherwise.

At the onset of a thermoacoustic instability, the perturbations are small in
comparison to the mean flow quantities, and they therefore induce a flame
response that scales linearly with their amplitude. This linear regime was
emulated in Chapter 6 by imposing small-amplitude acoustic modulations
($A=0.025$). During a combustion instability, the onset stage is often
followed by a continuous growth of the fluctuations intensity, until a stable
limit-cycle is reached. In this latter phase, the perturbations strength is
significantly larger: for instance in [207], the bulk fuel inflow velocity was
observed to oscillate from 10% to 15% of its mean value. These strong
disturbances significantly alter the flame response, that then entails
nonlinear features. Those can induce elaborate phenomena, such as chaotic
oscillations, harmonics and subharmonics generation, or even flame blow-off
and flash-back. In the present work, two common nonlinear flame response
mechanisms are investigated:

  * •

The flame response saturation, or in other words a relative decrease in the
strength of the heat-release fluctuations. This saturation is a key ingredient
in the establishment of a stable limit-cycle, as it slows down and ultimately
stops the growth of the thermoacoustic oscillations when an equilibrium state
is reached. The heat-release saturation was evidenced in 6.14, where the local
FTF gain is shown to decrease with $A$.

  * •

The generation of higher harmonics, that is a known route to chaos in laminar
combustion systems, through the period-doubling process that it induces [69].
Subharmonics generation may also occur, but are not considered here due to the
considerable associated computation cost.

In the light of these preliminary remarks, imposing a 10% amplitude modulation
($A=0.1$) appears as a reasonable choice to trigger nonlinear flame response
mechanisms representative of a LRE thermoacoustic limit-cycle. The following
Sections provide an analysis of the LES results based on comparisons with the
linear response computed in Chapter 6 (where the forcing amplitude was 4 times
smaller than that used here).

####  7.2 Nonlinear heat-release response

The two-dimensional maps of the azimuthal vorticity and heat-release Fourier
coefficients shown in Fig. 7.1 reveal a flow response spatial pattern similar
to the one observed with a low-amplitude forcing (see Fig. 6.12 for
comparison).

Figure 7.1: (a) Maps of the Fourier coefficients real parts for a high-
frequency forcing ($\bm{f^{*}=23.6}$). Dark lines: contour of
$\bm{\hat{Q}}$-criterion, colored by the normalized azimuthal vorticity
$\bm{\hat{\omega}_{\theta}^{*}=(\hat{\omega}_{\theta}W_{F})/(Au_{CH4}^{0})}$.
Blue-red: Fourier coefficient of the heat-release rate
$\bm{\hat{\omega}_{T}}$. (b) Same, for a low-frequency forcing
($\bm{f^{*}=2.9}$).

As they leave the fuel injector, acoustic waves generate a variable shear that
in turn produces pairs of vortex-rings which are convected along the mean flow
streamlines. At high frequency, these vortices are generated directly at the
injector exit, while at low frequency they are spawned downstream in the
second half of the flame. As they travel in the vicinity of the reaction
front, these vortices perturb it, which produces in-phase heat-release
fluctuations. The absolute local heat-release oscillations are naturally more
intense for this large amplitude forcing than in the linear case
($\hat{\omega}_{T}(A=0.1)\sim 5\times 10^{11}$ vs.
$\hat{\omega}_{T}(A=0.025)\sim 2\times 10^{11}$ at ${f}^{*}=23.6$, and
$\hat{\omega}_{T}(A=0.1)\sim 2\times 10^{11}$ vs.
$\hat{\omega}_{T}(A=0.025)\sim 1\times 10^{11}$ at ${f}^{*}=2.9$). The
relative vortex response is however roughly diminished by half
($\hat{\omega}_{\theta}^{*}(A=0.1)\sim 7$ vs.
$\hat{\omega}_{\theta}^{*}(A=0.025)\sim 14$ at ${f}^{*}=23.6$, and
$\hat{\omega}_{\theta}^{*}(A=0.1)\sim 2$ vs.
$\hat{\omega}_{\theta}^{*}(A=0.025)\sim 4$ at ${f}^{*}=2.9$).

#####  7.2.1 Local Flame Describing Function

The nonlinear dependence of the heat-release fluctuations $\hat{\omega}_{T}$
with respect to the forcing amplitude $A$ suggested in Fig. 7.1 can be more
precisely quantified through the computation of a frequency-domain Flame
Describing Function (FDF), which is a generalization of the FTF concept to
nonlinear systems [64, 250]. Since the modulation frequencies are high, and
since the flame is longer than in a usual gaseous case, it is not considered
compact. Consequently, a local FDF depending on both the forcing frequency
${f}^{*}$ and the axial location ${x}^{*}$ is computed. This type of non-
compact FDF was for instance computed by Li et al. [251] for a 1 meter long
turbulent premixed gaseous flame. The local FDF gain and phase are defined in
a similar fashion as in Sec. 6.3.1:

$\displaystyle
n({x}^{*},{f}^{*})e^{j\varphi({x}^{*},{f}^{*})}=\dfrac{L_{f}\hat{q}(x^{*})}{Q_{0}\left(\hat{u}/u_{0}\right)}$
(7.2)

where $\hat{q}(x^{*})$ is the Fourier coefficient of the fluctuating heat-
release rate par unit length, $Q_{0}$ is the mean flame power, and $u_{0}$ is
the mean bulk velocity in the fuel injection line. The associated cumulative
FDF gain $N_{c}({x}^{*},{f}^{*})$, which quantifies the overall heat-release
oscillations, is still defined as:

$\displaystyle
N_{c}({x}^{*},{f}^{*})=\left|\dfrac{L_{f}}{Q_{0}\left(\hat{u}/u_{0}\right)}\int_{0}^{{x}^{*}}\hat{q}(x^{*})d{x^{\prime}}^{*}\right|$
(7.3)

It is reminded that $N_{c}({x}^{*}=1,{f}^{*})$ is the global gain commonly
used for compact flames. Note that the computation of a complete FDF would
require to perform additional sets of simulations at different amplitude
levels. In the present work, only a small portion of this FDF at a fixed
amplitude $A=0.1$ is considered.

The local FDF gain, phase, and cumulative gain are represented over the
frequency range of interest and the flame length in Fig. 7.2.

Figure 7.2: Maps of the local FDF for a forcing at $\bm{A=0.1}$. (a) Local FDF
gain $\bm{n({x}^{*},{f}^{*})}$. (b) Cumulative FDF gain
$\bm{N_{c}({x}^{*},{f}^{*})}$. (c) Local FDF phase
$\bm{\varphi({x}^{*},{f}^{*})}$. The dark lines are points directly computed
from the 16 forced frequencies. For visualization purposes the surfaces were
obtained by cubic spline interpolation along the $\bm{{f}^{*}}$ axis.

The local FDF gain (Fig. 7.2-(a)) presents a preferential response region
whose topology is similar to that observed in the linear case (see Fig. 6.13
for comparison). At high-frequency, the heat-release response concentrates in
a small region at the injector exit, while for higher ${f}^{*}$ this region
widens and is shifted downstream towards the flame tip. Differences are
however observed on the magnitude of this preferential response region. In the
linear regime, the FTF gain reaches maximal values of the order of $12$ near
the half of the flame, for intermediate forcing frequencies ($8<{f}^{*}<11$).
In contrast, large amplitude modulations yield a FDF gain peaking at roughly
$5$ in the near-injector region for high ${f}^{*}$. This decrease of the
relative heat-release oscillations intensity clearly evidences a nonlinear
saturation of the flame response, which can be more quantitatively assessed
through the direct comparisons between the local FDF and FTF shown in Fig.
7.3.

Figure 7.3: Comparison between the local FDF (case $\bm{A=0.1}$) and the local
FTF (case $\bm{A=0.025}$), over the entire flame length, for a few select
forcing frequencies $\bm{{f}^{*}}$. (a) to (d): local FDF/FTF gains
$\bm{n({x}^{*},{f}^{*})}$. (e) to (h): local FDF/FTF phases
$\bm{\varphi({x}^{*},{f}^{*})}$.

At low frequency, the FDF saturation is rather marked, as the local gain in
the preferential region is roughly divided by a factor of 3 between the low-
amplitude and the large-amplitude modulations (Fig. 7.3-(a)). It is even more
pronounced at intermediate frequencies, where $n(A=0.1)\sim n(A=0.025)/5$ near
the half of the flame (Fig. 7.3-(b,c)). On the contrary, the local heat-
release response does not saturate for high-frequency forcing in the
preferential response region located near the injector (Fig. 7.3-(d)). In this
latter case, the only noticeable difference between the linear and the
nonlinear regime is a slight displacement of the response peak upstream,
towards the injector exit.

The FDF phase (Fig. 7.2-(c) and Fig. 7.3-(e) to (h)) is also similar to that
obtained in the low-amplitude case: it evolves linearly with ${x}^{*}$ from
the exit of the injector to a downstream location where it reaches a plateau
and stagnates. This pattern indicates a first region of constant propagation
speed $u_{\omega_{T}}$ of the heat-release disturbances along the axial
direction, followed by a second zone presenting bulk uniform heat-release
oscillations. Note that the plateau does not exist at low ${f}^{*}$,
suggesting a constant propagation speed $u_{\omega_{T}}$ over the entire flame
length. Direct comparisons between the propagation speeds $u_{\omega_{T}}$ in
the linear and the nonlinear regimes are provided in Fig. 7.4.

Figure 7.4: Comparison between the propagation speeds
$\bm{u_{\omega_{T}}^{*}=u_{\omega_{T}}/u_{CH4}^{0}}$ of the heat-release
disturbances along the flame axis, as a function of the forcing frequency
$\bm{{f}^{*}}$, for the large-amplitude and the small-amplitude modulations
($\bm{A=0.1}$ and $\bm{A=0.025}$, respectively). The equivalent mean flow
velocity $\bm{\langle\overline{u}^{*}_{L_{i}}\rangle}$ along the inner
streamline $\bm{L_{i}}$ (defined in Eq. (6.13)) is also indicated for
reference.

Comparable values are obtained in both cases for frequencies ${f}^{*}>5$, thus
indicating that the propagation of the heat-release disturbances along the
flame axis is not noticeably affected by nonlinear response mechanisms.
Remarkable differences are however observed for ${f}^{*}<5$, where the heat-
release perturbations propagate twice as fast in the nonlinear regime as in
the linear one. Interestingly, in the former case $u_{\omega_{T}}^{*}$ appears
more strongly correlated to the equivalent mean flow velocity
$\langle\overline{u}^{*}_{L_{i}}\rangle$ in the inner streamline $L_{i}$
(where vortex rings are convected), whereas in the linear regime
$u_{\omega_{T}}^{*}$ is significantly lower than the equivalent mean flow
velocity.

These dissimilarities in the propagation speeds $u_{\omega_{T}}^{*}$ have
first-order effects on the saturation of the cumulative FDF gain $N_{c}$ (Fig.
7.2-(b)). Because of the faster propagation at ${f}^{*}<5$, the heat-release
fluctuates more in-phase over the flame length in the nonlinear regime than in
the linear one. This spatio-temporal coherence of the heat-release
fluctuations directly results in a global flame response that is only slightly
weakened in comparison to that obtained with small-amplitude modulations
($N_{c}({x}^{*}=1,A=0.1)\sim 0.4$ vs. $N_{c}({x}^{*}=1,A=0.025)\sim 0.6$ at
${f}^{*}=2.9$). On the contrary, even though the local saturation at high
${f}^{*}$ is not as pronounced as that at lower ${f}^{*}$ (see Fig. 7.3-(a) to
(d)), the global response saturation is more marked in the former case. This
is due to a propagation speed $u_{\omega_{T}}$ nearly identical to the one in
the linear regime (Fig. 7.3-(f) to (h)), which produces out of phase heat-
release oscillations. This weak spatio-temporal coherence yields global FDF
gains for ${f}^{*}>5$ significantly lower than that of the FTF.

To conclude, the computation and analysis of the local FDF leads to three
major observations:

  * •

An significant saturation of the local heat-release response in the
preferential response region is evidenced for low to intermediate forcing
frequencies. Such local saturation does not exist for high-frequency
modulation.

  * •

The FDF and the FTF phases are nearly identical, at the exception of a faster
propagation of the heat-release perturbations along the flame axis in the
nonlinear regime for ${f}^{*}<5$.

  * •

Counterintuitively, the global flame response saturation differs from the
local one. The global gain saturation is more pronounced for high forcing
frequencies, while it is only slightly marked for low to intermediate
${f}^{*}$, which directly results from the faster propagation speed at
${f}^{*}<5$.

#####  7.2.2 Higher harmonics generation

Another nonlinear response feature commonly encountered in flame subjected to
strong acoustic perturbations is the generation of higher harmonics, at
frequencies multiple of the fundamental forcing frequency ${f}^{*}$. Because
of obvious limitations related to simulations cost and data storage, only the
first harmonic at $2{f}^{*}$ is considered in the present work. The possible
occurrence of such harmonic response in the modulated LES is assessed by
performing Fourier analyses at $2{f}^{*}$. In this section (and only this
one), the notation $\hat{v}$ therefore refers to the Fourier coefficient of
the variable $v$ computed at the frequency $2{f}^{*}$. Figure 7.5 shows maps
of the vorticity and heat-release harmonic responses for low-frequency and
high-frequency modulations.

Figure 7.5: (a) Maps of the Fourier coefficients real parts, computed at the
1st-harmonic $\bm{2{f}^{*}}$, for a high-frequency forcing
($\bm{f^{*}=23.6}$). Dark lines: contour of $\bm{\hat{Q}}$-criterion, colored
by the normalized azimuthal vorticity
$\bm{\hat{\omega}_{\theta}^{*}=(\hat{\omega}_{\theta}W_{F})/(Au_{CH4}^{0})}$.
Blue-red: Fourier coefficient of the heat-release rate
$\bm{\hat{\omega}_{T}}$. (b) Same, for a low-frequency forcing
($\bm{f^{*}=2.9}$).

In comparison to the fundamental response at ${f}^{*}$ (see in Fig. 7.1), the
1st-harmonic flow response displays a similar pattern comprising pairs of
vortex rings in the fuel stream that perturb the reactive layer and therefore
result in heat-release fluctuations at $2{f}^{*}$. Identical frequency-
dependent response regions are revealed, with a near-injector preferential
response for a high-frequency modulation and a wider downstream preferential
response for low-frequency excitation. Expectedly, in both cases the 1st-
harmonic flow structure presents shorter wavelengths that are halved in
comparison to the fundamental, indicating the occurrence of a secondary
vortex-shedding mode in the fuel stream. Note that the harmonic response
relative amplitudes are lower than the fundamental ones: for instance at
${f}^{*}=2.9$ the 1st-harmonic vorticity response is roughly halved in
comparison to the fundamental, while at ${f}^{*}=23.6$ it is approximately
only 50% weaker.

More quantitatively, harmonics generation in the flame nonlinear response can
be evaluated by a generalized frequency-domain transfer function. In this
matter, Haeringer et al. [252] introduced the Extended Flame Describing
Function (xFDF), which relies on the observation that if the flame response
produces higher harmonics, the forcing velocity is then likely to also contain
those. Thus, the heat-release 1st-harmonic $\hat{q}_{2}$ can originate from
two distinct contributions: (1) the 1st-harmonic FDF that transforms the
fundamental of the velocity signal $\hat{u}_{1}$ into a heat-release
oscillation at $2{f}^{*}$, and (2) the fundamental FDF that converts the
velocity harmonic $\hat{u}_{2}$ oscillating at $2{f}^{*}$ into a heat-release
fluctuation also at $2{f}^{*}$. Formally, this yields the definition of the
1st-harmonic FDF:

$\displaystyle
n_{2}({x}^{*},{f}^{*},A)e^{j\varphi_{2}({x}^{*},{f}^{*},A)}=\dfrac{L_{f}\hat{q}_{2}({x}^{*})}{Q_{0}(\hat{u}_{1}/u_{0})}-n({x}^{*},{f}^{*},A)e^{j\varphi({x}^{*},{f}^{*},A)}\left(\dfrac{\hat{u}_{2}}{\hat{u}_{1}}\right)$
(7.4)

where $n$ and $\varphi$ refer to the fundamental FDF defined in Eq. (7.4), and
$A$ is still the magnitude of the fundamental velocity $\hat{u}_{1}$
normalized by $u_{CH4}^{0}$. However in the present case, Fourier analyzes of
velocity time-series in the fuel injection line did not reveal the appearance
of any harmonic $\hat{u}_{2}$. The 1st-harmonic FDF can therefore be computed
by only retaining the first term in the right-hand side of Eq. (7.4). In Fig.
7.6, the strength of the harmonic flame response is assessed in the linear and
nonlinear regimes by comparing the 1st-harmonic gain $n_{2}({x}^{*},{f}^{*})$
to the fundamental gain $n$ in the preferential response region (computed as
$\max_{{x}^{*}}(n)$).

Figure 7.6: Comparison of the 1st-harmonic gain $\bm{n_{2}}$ with respect to
the fundamental gain $\bm{n}$, between the nonlinear ($\bm{A=0.1}$) and the
linear ($\bm{A=0.025}$) regimes, at a few select forcing frequencies
$\bm{{f}^{*}}$.

Unsurprisingly, large-amplitude modulations yield significantly stronger
harmonic responses than the small-amplitude ones, over most of the frequency
range. Note that this point is not true at ${f}^{*}=10.3$, where the harmonic
response magnitude of both cases are comparable (Fig. 7.6-(b)). This is
however not due to the generation of the harmonic by nonlinear mechanisms in
the flame response, but rather to the presence of the self-sustained vortex-
shedding mode at ${f}^{*}_{i}=21.4\approx 2{f}^{*}$, as discussed in Sec.
6.2.3 and Sec. 6.5. In the nonlinear regime, the 1st-harmonic response is the
strongest at low forcing frequencies near the middle of the flame (Fig.
7.6-(a)), where $n_{2}$ reaches about 40% of the fundamental response $n$. The
1st-harmonic response strength is lower at higher ${f}^{*}$, but still of
significant importance as it peaks at 20% of the fundamental gain. The 1st-
harmonic propagation speed along the flame axis (computed from the FDF phase
$\varphi_{2}({x}^{*},{f}^{*})$) does not noticeably differ from the
fundamental one, with values of $1.2u_{CH4}^{0}$ at ${f}^{*}=23.6$ and
$0.7u_{CH4}^{0}$ at ${f}^{*}=1.5$.

####  7.3 Nonlinear vortex dynamics

Section 6.4 showed that in the linear regime, the heat-release response is
tightly related to the dynamics of vortex-rings pairs that are generated and
convected in the CH4 annular jet. Resemblances in the spatial response
patterns of the flame and of the vorticity observed in Fig. 7.1 suggest a
similar mechanism in the nonlinear regime. This section therefore aims at
assessing nonlinear vortical waves dynamics in the fuel stream, with an
emphasis on the relation between the respective saturations of the vorticity
and flame responses, as well as on the convection speed of the vortex-rings.
Note that the generation of vorticity higher harmonics is likely, as evidenced
in Fig. 7.5, but is not discussed in this work. The convection trajectories of
the vortices pairs are observed to be identical to those in the linear regime:
vortex-rings formed at the outer corner of the CH4 injector travel along a
straight line $L_{o}$, while those engendered at the inner corner move along
the path $L_{i}$, as illustrated in Fig. 6.18. Nonlinear vortical waves
dynamics are here quantified thanks to two frequency-domain Vorticity
Describing Function (VDF), which are defined analogous to a FDF, as the
generalization of the Vorticity Transfer Function of Eq. (6.15):

$\displaystyle|\hat{\omega}_{\theta}^{*}|({x}^{*},{f}^{*})\
e^{j\varphi_{\omega_{\theta}}({x}^{*},{f}^{*})}=\dfrac{R({x}^{*})}{Au_{CH4}^{0}}\hat{\omega}_{\theta}({x}^{*},{f}^{*})$
(7.5)

where the normalization by the typical vortex radius $R({x}^{*})$ allows the
direct comparison between the strengths associated to vortices of different
sizes. A VDF is computed on each trajectory $L_{i}$ and $L_{o}$.

#####  7.3.1 Vorticity preferential response regions

The characteristic shape of the heat-release response region, observed in both
linear and nonlinear regimes (Fig. 7.2-(a)), was shown in Sec.6.4.2 to
correspond to a vorticity preferential response region. This latter results
from a loading and firing mechanism that governs the vortex-shedding dynamics,
yielding a formation area located at ${x}^{*}\sim u_{CH4}^{0}/{f}^{*}$, for
small-amplitude modulations. Nonlinear phenomena altering the vortex
generation are here estimated in Fig. 7.7 by comparing the locations of the
heat-release and vorticity response regions between the linear and the
nonlinear regimes. Those are computed as the locales of the peaks in the FDF
gain and VDF gain, respectively.

Figure 7.7: Comparison of the axial locations of the heat-release and
vorticity preferential response regions in the nonlinear and the linear
regimes, as a function of the modulation frequency $\bm{{f}^{*}}$. The
vorticity preferential region is the one of the inner layer $\bm{L_{i}}$.

In both linear and nonlinear cases, the vorticity and the heat-release
response regions appear strongly correlated, which confirms the direct
influence of the vortex-rings generation on the flame response, even at large
modulation amplitudes. Similar vortex generation locations roughly scaling as
${x}^{*}\sim u_{CH4}^{0}/{f}^{*}$ are observed, with the difference that the
response regions in the nonlinear regime are slightly shifted upstream towards
the injector. This deviation was already visible on the local FDF in Fig. 7.3.
Note that the shifting length $\Delta x$ between the large and the low
amplitude cases is the same for the vorticity and the heat-release response
regions, and that it seems independent of the forcing frequency, with an
approximate value $\Delta x\approx 0.1L_{f}$. Such translation is not
negligible, since it approaches one-third of the heat-release perturbation
wavelength at ${f}^{*}=2.9$, and more than two heat-release wavelengths at
${f}^{*}=23.6$.

#####  7.3.2 Vortex convection

The nonlinear convection of the vortex-rings along their respective
trajectories is estimated thanks to the equivalent propagation speed $\langle
v_{\omega_{\theta}}^{*}\rangle({f}^{*})$ in the direction of the path $L_{i}$
or $L_{o}$. It is computed with the gradient of the VDF phase
$\varphi_{\omega_{\theta}}({x}^{*},{f}^{*})$ based on the definition of Eq.
(6.16) and on the $|\hat{\omega}_{\theta}^{*}|$-weighted spatial average
detailed in Eq. (6.17). The coefficient $\alpha\approx\langle
v_{\omega_{\theta}}\rangle/\langle\overline{v}\rangle$, defined in Eq. (6.14),
is of particular interest to quantify the difference between the true vortex
convection speed and the local mean flow velocity. The local VDF gains
$|\hat{\omega}_{\theta}^{*}|$ are used to measure the vortices strength.

Figure 7.8-(a,c) compares the linear and nonlinear vortex-rings propagation
speeds along $L_{i}$ and $L_{o}$.

Figure 7.8: Comparison of the vortex responses in the nonlinear and the linear
regimes, on both trajectories $\bm{L_{i}}$ (left column) and $\bm{L_{o}}$
(right column), as a function of the modulation frequency $\bm{{f}^{*}}$. (a)
and (c): Comparison of the equivalent vortex convection speeds, defined by Eq.
(6.17). (b) and (d): Comparison of the maximal VDF and VTF peaks magnitudes
for each forcing frequency.

The inner vortex propagation speed $\langle v_{\omega_{\theta}}\rangle$
presents a comparable growth in the linear and the nonlinear cases, as it
increases by approximately $0.5u_{CH4}^{0}$ between ${f}^{*}=1.5$ and
${f}^{*}=23.6$ (Fig. 7.8-(a)). It is however lower for large-amplitude than
for small-amplitude modulations, with a deviation approaching $0.1u_{CH4}^{0}$
over most of the frequency range. As a result, the nonlinear convection speed
never exceeds the bulk injection velocity, whereas it reaches
$1.15u_{CH4}^{0}$ at high frequencies in the linear case. This slower vortex
propagation in the nonlinear regime is likely to be responsible for the
upstream translation of the preferential response regions towards the injector
evidenced in Fig. 7.7. The coefficient $\alpha$ was found to stay close to
unity in the linear regime (varying from $0.9$ at low ${f}^{*}$ to $1.1$ at
high ${f}^{*}$ in Fig. 6.20). This sensitivity of the convection speed with
respect to the local mean flow was attributed to the vortex inertia that
quickly decays after its shedding, under the influence of the intense pseudo-
boiling due to the proximity of the flame. In contrast, at larger modulation
amplitude the variation of $\alpha$ is more pronounced, since it grows from
$0.75$ at ${f}^{*}=1.5$ to $1.25$ at ${f}^{*}=23.6$. This shows that in the
nonlinear regime the vortex-rings motion in the inner layer $L_{i}$ is less
sensitive to the local mean velocity. Such difference is imputed to the true
convection speed that is affected not only by the local mean velocity, but
also by the local velocity fluctuations induced by the vortex pairs themselves
(see the term $u^{\prime}(x,t)$ in Eq. (6.14)). In contrast, the outer
vortices convection speed $\langle v_{\omega_{\theta}}\rangle$ is only
slightly influenced by the larger amplitude forcing (Fig. 7.8-(c)). It indeed
remains close to $0.9u_{CH4}^{0}$ over most of the frequency range,
independently of the region where the vortices travel. Similarly to the linear
regime, this low sensitivity to the local mean flow is due the large inertia
of the outer vortex-rings, as they trap dense pockets of cryogenic fuel. This
invariance of $\langle v_{\omega_{\theta}}\rangle$ combined with the
translation of the preferential response region towards the injector results
in an increase of $\alpha$ at high ${f}^{*}$ even more dramatic than in the
linear case. Indeed, the outer vortex-rings are in this case forced through
the PIRZ (Fig. 6.6) where $\langle\overline{v}\rangle$ is low, which in turn
amplifies $\alpha$.

The saturation of the vorticity response strength is revealed in Fig.
7.8-(b,d), where the magnitude of the peaks in the VDF and VTF gains are
compared. Both inner and outer vortex-rings responses are considerably weaker
at $A=0.1$. Remarkably, the resonant peaks that were observed around
${f}^{*}=10$ in the linear regime have vanished, and the vorticity response
strength only presents a slow, continuous growth with ${f}^{*}$. In Sec. 6.5
the apparition of these peaks was imputed to equal vortex convection speeds in
$L_{i}$ and $L_{o}$ around ${f}^{*}=10$, resulting in a synchronization of
both vortices, that in turn produces stronger vorticity fluctuations. On the
contrary, in the nonlinear regime this synchronization does not occur at
${f}^{*}\approx 10$ ($\langle
v_{\omega_{\theta}}\rangle_{{f}^{*}=10}=0.78u_{CH4}^{0}$ in $L_{i}$ vs.
$\langle v_{\omega_{\theta}}\rangle_{{f}^{*}=10}=0.9u_{CH4}^{0}$ in $L_{o}$),
but rather near ${f}^{*}\gtrsim 18$ (Fig. 7.8-(a,c)), which corroborates the
shifting of the maximal VDF gain towards higher frequencies. Most importantly,
this profound modification of the vorticity response intensity is directly
related to the heat-release response saturation (Fig. 7.3-(a) to (d)). At low
to intermediate frequencies, the heat-release and inner vortex-rings
saturations are similar. At ${f}^{*}=3$, the FDF and VDF gains are roughly
halved in comparison to the linear regime. The maximum saturation occurs
around ${f}^{*}=10$, where both FDF and VDF gains are divided by a factor 5.
Some dissimilitudes are however visible at high frequency, since there is no
noticeable saturation of the heat-release, whereas the vorticity response on
$L_{i}$ displays a saturation factor of 2 to 3, and that on $L_{o}$ a
saturation factor of 1.5.

To conclude, the analysis of the Vorticity Describing Functions computed on
the trajectories $L_{i}$ and $L_{o}$ leads to three major observations:

  * •

The propagation of vortex-rings along $L_{i}$ is less sensitive to the local
mean flow velocity in the nonlinear regime than in the linear one. The true
vortex convection speed is believed to be reduced under the effect of the
local velocity fluctuations induced by the vortices themselves. On the
contrary, the propagation of the outer vortex-rings is not altered.

  * •

A direct consequence of the previous point is a translation of the heat-
release and vorticity preferential response regions towards the injector.

  * •

The saturation of the vorticity response is sufficient to explain that of the
heat-release for low to intermediate modulation frequencies. Because of the
inner-vortex rings lower convection speed, the minimal vorticity saturation
occurring when both inner and outer vortices synchronize, is now observed at
high frequencies. It is however still noticeable, which cannot explain the
absence of saturation of the flame response at high ${f}^{*}$.

####  7.4 Contributions to HR nonlinear dynamics

In this section, the nonlinear dynamics of the distinct contributions to the
heat-release response are analyzed to identify potential saturation
mechanisms. The reasoning of Sec. 6.3.2 based on the flame sheet assumption is
repeated. The heat-release fluctuations are decomposed into five distinct
contributions, which are recalled here for convenience:

$\displaystyle\dfrac{\hat{q}}{q_{0}}=\left[\dfrac{\hat{\rho}}{\rho_{0}}+\dfrac{\hat{\mathcal{D}}}{\mathcal{D}_{0}}+\dfrac{\hat{\Theta}_{Z}}{\Theta_{Z0}}+\dfrac{\hat{\Psi}_{F}}{\Psi_{F0}}+\dfrac{\hat{\xi}}{\xi_{0}}\right]_{Z=Z_{st}}$
(7.6)

where $\mathcal{D}$ is the species diffusion coefficient, $\Theta_{Z}$
represents the mixture fraction gradient, $\Psi_{F}$ is related to the
internal flame structure, and $\xi$ is the flame surface area per unit length
(i.e. proportional to the local mean flame radius). Once again, Eq. (7.6) is
used to define 5 Contribution Describing Functions (CDFs), as amplitude-
dependent generalizations of the Contribution Transfer Functions (CTFs)
detailed in Eq. (6.10). For each contribution $v$ of Eq. (7.6), the associated
CDF writes:

$\displaystyle
n_{v}({x}^{*},{f}^{*})e^{j\varphi_{v}({x}^{*},{f}^{*})}=\dfrac{L_{f}q_{0}\hat{v}(x^{*})}{Q_{0}v_{0}\left(\hat{u}/u_{0}\right)}$
(7.7)

The reader is reminded that this equation is obtained through a linearization,
followed by a Fourier transform at the modulation frequency ${f}^{*}$. This
approach may be justified for small-amplitude fluctuations, but not in a
nonlinear regime where higher-order terms $\widehat{v_{i}v_{j}}$ may be of
primary importance. Note however that such terms are likely to produce higher-
harmonic oscillations, and since the present section only aims at assessing
the saturation mechanisms, Eq. (7.7) is still a reasonable approximation.

Figure 7.9 compares the 5 CDF gains along the flame axis for a few forcing
frequencies.

Figure 7.9: CDF gains $\bm{n_{v}}$ along the flame axis, at 4 forcing
frequencies.

The relative weights of the contributions are highly similar to those obtained
in the linear regime (see Fig. 6.15 for comparison). The diffusion coefficient
fluctuations are the overall major contributor to the heat-release response.
The second most dominant contribution stems from the mixture fraction gradient
fluctuations in the near-injector region at high frequency, or from the flame
surface area in the second half of the flame at lower frequency. The internal
flame structure and the density do not vary sufficiently to significantly
contribute to the flame response. The clear similitudes observed between the
linear and the nonlinear regimes suggest that the physical mechanisms
responsible for these relative contributions are identical to those already
described in Sec. 6.3.2.

In order to quantify the saturation of each one of the three dominant
contributions ($\mathcal{D}$, $\Theta_{Z}$, and $\xi$), their respective CDF
gains are compared to the CTF gains in Fig. 7.10.

Figure 7.10: Comparison between the CDF and the CTF gains along the flame axis
for a few select modulation frequencies $\bm{{f}^{*}}$. First row: species
diffusion coefficient contribution. Second row: mixture fraction gradient
contribution. Third row: flame surface area contribution.

At low frequency (Fig. 7.10-(a,e,i)), the contributions’ responses are
approximately half as strong in the nonlinear case than in the linear one,
which also corresponds to the saturation observed on the FDF gain (Fig.
7.3-(a)). However, the flame saturation is not in this case due to the
contributions’ saturation, since the vorticity (Fig. 7.8-(a)), considered as
the origin of the flow disturbances affecting the reaction front, is also
twice smaller for $A=0.1$ than for $A=0.025$. The same remark is true at
intermediate frequency (second and third columns in Fig. 7.10), where the
heat-release, its contributions, and the vorticity (Fig. 7.8-(b,d)) display a
similar saturation factor comprised between 4 and 5. In contrast, the
contributions responses are remarkably different at high frequency (Fig.
7.10-(d,h,l)): the diffusion and the mixture-fraction fluctuations are
respectively 50% and 20% more intense in the nonlinear than in the linear
regime. This trend goes against that visible on the vorticity response, which
still saturates at high ${f}^{*}$. The enhanced contributions of $\mathcal{D}$
and $\Theta_{Z}$ to the heat-release fluctuations may counterbalance the
vortex-rings saturation, and would therefore explain the only slight flame
saturation observed at high frequencies in Fig. 7.3-(d).

To conclude, the examination of the Contribution Describing Functions leads to
three major observations:

  * •

The contributions to the heat-release fluctuations, their relative weights,
and their governing mechanisms appear to be similar in the nonlinear and
linear regimes.

  * •

The dominant contributors are not responsible for the flame saturation for low
to intermediate modulation frequencies, since their respective responses are
in line with that of the vorticity.

  * •

Most remarkably, at high frequencies the species diffusion and mixture-
fraction gradient variations are more intense in the nonlinear than in the
linear regime. They therefore offset the vorticity response saturation, and
result in the absence of heat-release saturation.

####  7.5 Conclusion

A series of LES was performed to assess the nonlinear dynamical response of a
doubly-transcritical LO2/LCH4 cryogenic flame to harmonic acoustic modulation
imposed at the fuel injector. Particular emphasis was put on extending the
linear analysis of Chapter 6 to the case of large-amplitude forcing. Two
nonlinear mechanisms were especially targeted, namely the flame response
saturation, and the generation of higher-harmonics. These computations
necessitated considerable computational resources, including 15 million CPU
hours on a machine equipped with Intel Xeon 8168 processors (2.7GHz).

A portion of a non-compact Flame Describing Function was computed at a fixed
modulation amplitude, equal to 10% of the bulk injection velocity, and over a
wide frequency range spanning from $\mathcal{O}(1\ \mathrm{kHz})$ to
$\mathcal{O}(20\ \mathrm{kHz})$. The examination of this FDF evidenced a more
pronounced heat-release response saturation for low to intermediate forcing
frequencies, in the downstream half of the flame. Conversely, at high
frequencies the heat-release preferential response region, located near the
injector exit, does not noticeably saturate, but is shifted upstream towards
the injector. An analysis of the nonlinear vortex dynamics and of the
contributions to the heat-release fluctuations led to a physical
interpretation of these results. More precisely, the flame saturation at low
to moderate frequencies was attributed to that of the vortex-rings pairs
traveling at the edges of the CH4 annular jet. The pronounced vorticity
response saturation at moderate frequencies was in turn imputed to the lower
propagation speed of the inner vortices in comparison to the linear regime,
thus resulting in a desynchronization of the vortex-rings pair and in a weaker
response. The absence of heat-release saturation at high frequency is believed
to originate from a moderate vorticity saturation, counterbalanced by enhanced
fluctuations of the species diffusion coefficient and of the mixture-fraction
gradient.

The generation of a first harmonic at twice the modulation frequency was
observed in the heat-release and vorticity fields. A portion of a non-compact
Extended Flame Describing Function was computed to estimate the strength of
this harmonic in comparison to the fundamental. It was shown that the harmonic
generation occurs mostly at low to intermediate frequencies near the middle of
the flame, where the 1st-harmonic magnitude approaches 40% of the fundamental.

The proposed analyses are expected to be useful to guide future theoretical
studies interested in the nonlinear response of LRE coaxial jet flames to
acoustic perturbations. In particular, the strong influence of the complex
vortex dynamics on the flame response was unveiled, which emphasizes the
necessity to develop advanced shedding and convection models for vortex-ring
pairs in cryogenic variable-density flows. Finally, the computed FDF and xFDF
are made available and can be embedded in thermoacoustic low-order models of
multi-injectors LREs.

###  Chapter 8 Flame-wall interaction effects on the flame root stabilization
mechanisms of a doubly-transcritical LO2/LCH4 coaxial jet-flame

In the previous chapters, as in the majority of existing studies related to
LRE combustion, the boundary condition at the injector lip is supposed
adiabatic. This assumption is highly disputable, considering the
submillimetric lip thickness and the proximity of both cold cryogenic
reactants and hot combustion products. The present chapter therefore aims at
thoroughly assessing the effect of Flame-Wall Interaction (FWI) on the
anchoring mechanism of a doubly-transcritical LO2/LCH4 jet-flame in the
geometry of the Mascotte test rig. It uses high-fidelity quasi-DNS with
complex chemistry to study the flame root stabilization mechanisms, and the
FWI effects are accounted for by considering the unsteady conjugate heat-
transfer problem. It starts by introducing the multi-physics numerical
framework employed to resolve such a strongly coupled problem. A simulation
with adiabatic boundary conditions is also performed for comparison.
Simulation results are then presented. An analysis of the flame root structure
provides a detailed insight into its stabilization mechanisms and shows that
the large wall heat losses at the lips of the coaxial injector are of primary
importance. It is found that adiabatic walls simulations lead to enhanced
cryogenic reactants vaporization and mixing, and to a quasi-steady flame,
which anchors within the oxidizer stream. Then, the analysis of the non-
adiabatic flame dynamics evidences self-sustained oscillations of both lip
temperature and flame root location at similar frequencies. The flame root
moves from the CH4 to the O2 streams at approximately $450$Hz, affecting the
whole flame structure. This distinctive anchoring mechanism is expected to
have first-order effects on the global flame topology and dynamics, including
its response to acoustic perturbations. To the knowledge of the author, this
is one of the first attempts to simulate a doubly-transcritical shear-coaxial
flame including FWI and complex CH4/O2 chemistry.

####  8.1 The conjugate heat-transfer problem

Flame-Wall Interaction (FWI) is known to strongly affect flame anchoring, for
instance by inducing flame quenching and lift-off. These effects are in turn
of prevalent importance on the flame global topology, on its intrinsic
dynamics, as well as on its response to acoustic perturbations [61].
Conversely, any modification of the flame features may result in a
modification of the heat flux that the flame provides to solid boundaries,
which may perturb the temperature field within the wall. This two-way coupling
between reactive flow phenomena and heat conduction in the solid constitutes a
conjugate heat transfer problem, as illustrated in Fig. 8.1.

Figure 8.1: Schematic of the conjugate heat transfer problem, or in other
words the two-way coupling due to Flame-Wall Interaction.

For most combustion applications, this coupling is rather loose: indeed, the
heat diffusion time-scale in the solid can be as large as $\mathcal{O}(1\
\textrm{s})$, while the flow characteristic time-scale is in $\mathcal{O}(1\
\textrm{ms})$. The temperature within the wall is then insensitive to fast
oscillations in the fluid (but still depends on the mean flow and on slower
oscillations), such that only the average solid temperature field is involved
in the conjugate heat transfer problem. However, as briefly mentioned in
[191], the dynamic flame-wall coupling is expected to be stronger in the case
of a LRE coaxial injector, due to the combination of very high heat fluxes and
thin lips. This effect is expected to be even more pronounced in the case of
doubly transcritical combustion, since both reactants are injected in a
cryogenic state, which will be confirmed later. Note that the numerical
modeling of FWI usually requires complex chemistry, due to the relatively cold
temperature in the vicinity of the solid boundary. Thus, most previous
numerical studies of LRE ccombustion overlooked the FWI problem at the coaxial
injector lips. They instead relied on two crude simplifications: (1) adiabatic
boundary conditions were assumed at the coaxial injector lip [192, 191], and
(2) simplified chemistry, tabulated [228, 199] or based on two-step reduced
schemes [198], was employed. In the light of these observations, the present
work assesses the FWI influence on the anchoring mechanisms of the doubly-
transcritical LO2/LCH4 flame studied in the previous chapters. Most
importantly, it considers an unsteady, synchronized conjugate heat transfer
problem, where the flame and the lip temperature field may be strongly
coupled.

#####  8.1.1 Numerical setup

As the present study focuses on flame root stabilization, computations are
performed on a two-dimensional slice of the Mascotte geometry (see Fig. 5.6).
The computation domain, shown in Fig. 8.2 is also truncated in the axial
direction.

Figure 8.2: The two-dimensional computational domain. Thick dark lines
represent no-slip adiabatic walls for the fluid domain. Dashed lines are non-
adiabatic no-slip walls were fluid and solid domains are coupled. The basis of
the lip (hashed line) is an adiabatic boundary condition for the solid domain.

The area of interest, of dimensions $2d_{O2}\times 1d_{O2}$, is located near
the injector lip of thickness $\delta_{l}$, and consists of cells of dimension
0.03$\delta_{l}$ (i.e. a few $\mu\rm m$). The diffusion flame reactive
thickness (of order $100\ \mu$m for the strain rate levels considered here)
implies a high resolution of the front with at least 10-16 cells across the
reactive thickness. Injection conditions (recalled in Tab. 8.1 for
convenience) are the same as those employed in the previous chapters: they are
similar to the operation point T1 in [194]. The only difference is that the
chamber pressure has been increased to 75 bar (instead of 54 bar in the
experimental study), giving higher reduced pressures.

| Injection conditions |  Thermal properties
---|---|---
|  $\dot{m}_{inj}/\dot{m}_{inj}^{O2}$ $T_{R}^{inj}$ $P_{R}^{inj}$ |  $e_{f}$ (SI) $\kappa_{f}$
O2 |  $1.0$ $0.55$ $1.49$ |  $260$ $75.5$
CH4 |  $3.59$ $0.62$ $1.64$ |  $2200$ $8.9$
HG | - |  $80$ $245$
Table 8.1: Injection conditions (mass flow rate, reduced pressures and
temperatures), and fluid thermal properties (effusivity $\bm{e}$ (SI), and
effusivity ratio $\bm{\kappa_{f}}$) for reactants and hot gases (HG). The
global equivalence ratio is $\bm{\phi_{g}=14.3}$, the momentum flux ratio is
$\bm{J=(\rho_{F}u_{F}^{2})/(\rho_{O}u_{O}^{2})=33.3}$, and the wall effusivity
is $\bm{e_{w}=1.9\times 10^{4}}$ (SI).

One-seventh power-law velocity profiles are imposed at both reactants inlets
through NSCBC boundary conditions [224] adapted for real-gas thermodynamics.
All walls are resolved and a no-slip velocity condition is imposed. Only the
temperature field in the lip of the central injector is coupled with the fluid
domain, such that the only non-adiabatic solid boundaries are the interfaces
between the splitter plate and the fluid domain. The strategy to compute
thermal quantities at this non-adiabatic boundary is detailed in the following
paragraphs.

At this point it is interesting to define some quantities relevant to
conjugate heat transfer. The fluid effusivity is defined as $e_{f}=(\rho
C_{p}\lambda)^{1/2}$, with $\lambda$ the thermal conductivity and $C_{p}$ the
specific heat capacity. We also define the fluid effusivity ratio by
$\kappa_{f}=e_{w}/e_{f}$. The effusivity allows to evaluate the ability of a
body to exchange heat with its surrounding. Thus, for $\kappa_{f}\gg 1$,
temperature at the solid/fluid interface is mostly determined by the wall
temperature, and this one can be considered as isothermal. On the contrary, if
$\kappa_{f}\ll 1$ the interface temperature is imposed by the fluid and the
wall may be considered as adiabatic. Intermediate values of $\kappa_{f}$
correspond to boundaries that are neither adiabatic nor isothermal. Thermal
properties (Tab.8.1) show the dominant thermal influence of the cryogenic CH4
stream on the lip. None of the interfaces can be considered as adiabatic. Most
importantly, the wall conduction characteristic time-scale is
$\tau_{w}=(\delta_{l}^{2}\rho C_{p}/\lambda)=0.6\ \rm ms$, which is comparable
to that of unsteady motions in the fluid. This very short time scale-scale,
that differs from most usual FWI applications, justifies the need for a
strongly coupled unsteady conjugate heat transfer problem, where both the
solid and the fluid temperatures vary on similar time scales.

#####  8.1.2 The AVTP solver

Resolving the conjugate heat-transfer problem requires to solve for the
unsteady temperature field within the lip (white domain in Fig. 8.2). In this
work, this is performed thanks to the solid thermal solver AVTP [253, 254],
developed at CERFACS. Its structure is based on AVBP, and it solves the heat
equation on unstructured hybrid meshes:

$\displaystyle\rho^{S}C_{p}^{S}\partialderivative{T^{S}}{t}=\partialderivative{xi}\left(\lambda^{S}\partialderivative{T^{S}}{xi}\right)$
(8.1)

where superscripts S denote quantities related to the solid. In Eq. (8.1) the
diffusion term is numerical integrated thanks to the $2\Delta$ operator [255].
The time-integration is here carried out with a 1st-order forward Euler
scheme:

$\displaystyle\partialderivative{T^{S}}{t}\approx\dfrac{T_{n+1}^{S}-T_{n}^{S}}{dt_{S}}$
(8.2)

where $T_{n+1}^{S}$ and $T_{n+1}^{S}$ are the solid temperature fields at the
instants $n+1$ and $n$, respectively, and $dt_{S}$ is the solid time-step.
This scheme is stable if the Fourier number $Fo$ is small enough:
$Fo=(\lambda^{S}dt_{S})/(\Delta_{x}^{2}C_{p}^{S}\rho^{S})<0.5$, with
$\Delta_{x}$ the cube-root of the cell volume. In the present simulation, the
solid time-step is fixed to a constant value, such that $Fo=0.1$. Classical
heat conduction boundary conditions are implemented in AVTP. Those include for
instance the isothermal condition ($T_{w}^{S}=T_{ref}$), the adiabatic
condition ($q_{w}^{S}=0$), and the heat-loss condition ($q_{w}^{S}=q_{ref}$),
where $T_{w}^{S}$ and $q_{w}^{S}$ are the lip temperature and heat-flux at the
fluid-solid interface.

In the present simulation, the mesh used by the thermal solver is uniform,
composed of triangular cells of size $\delta_{l}/120$. A quick evaluation of a
characteristic diffusion length in the lip shows that it is considerably
shorter than the axial length of the solid domain considered here: the lip
basis is therefore close to adiabatic conditions.

#####  8.1.3 Multi-physics code coupling

The most crucial components in the resolution of the conjugate heat transfer
problem is the coupling between the fluid solver and the solid thermal solver.
In this work, the Parallel Coupling Strategy (PCS) designed by Duchaine et al.
[256], and based on the OpenPalm library[257] is employed. In this approach,
depicted in Fig. 8.3, both solvers run in parallel on distinct set of
processors.

Figure 8.3: Flowchart of the Parallel Coupling Strategy (PCS). Both fluid and
solid solvers run in parallel on different set of processors. Boundary data
ensuring the physical coupling are exchanged at the fluid/solid interface at
meeting points, that are defined thanks to the exchange frequencies
$\bm{\alpha_{F}}$ and $\bm{\alpha_{S}}$.

Some meeting points are defined and serve to schedule the exchange of boundary
data values necessary to enforce a physical coupling. At each meeting point,
both solvers update their respective thermal boundary conditions that are then
used to continue the temporal integration. There exist different possibilities
to define this exchange, however Giles et al. [258] proved that the conjugate
heat transfer problem is stable only if an isothermal condition is imposed to
the fluid domain, and a heat-loss (or Neumann) condition is imposed to the
solid domain. Thus, at the coupling event $n$, AVBP-RG and AVTP respectively
exchange their wall heat flux $q_{w}^{F(n)}$ and their wall temperature
$T_{w}^{S(n)}$. After this coupling event, AVBP-RG computes $\alpha_{F}$
iterations of constant time-step $dt_{F}$, with an isothermal boundary at the
injector lip defined by $T_{w}^{S(n)}$. Meanwhile, AVTP computes $\alpha_{S}$
iterations at constant time-steps $dt_{S}$ with a heat-loss condition
determined by $q_{w}^{F(n)}$. At the next coupling event $n+1$, both boundary
values are updated from the fields that have been integrated, and so forth.
Note that since an unsteady and strongly coupled conjugate heat transfer
problem is considered here, both solvers need to be synchronized in time, that
is $\alpha_{F}dt_{F}=\alpha_{S}dt_{S}$. This is enforced by fixing $dt_{S}$ to
a constant value $dt_{S}=0.1(\Delta_{x}^{2}C_{p}^{S}\rho^{S})/\lambda^{S}$,
and by setting a large exchange frequency with $\alpha_{S}=50$ and
$\alpha_{F}=5$; $dt_{F}$ is then fixed to a constant value
$dt_{F}=dt_{S}\alpha_{S}dt_{S}/\alpha_{F}$. This approach ensures a fast
convergence of the solid solver and avoids any aliasing, but also results in a
significant computational cost.

####  8.2 Flame root stabilization mechansims

Two simulations are performed, the first one (superscript C) using conjugate
heat transfer at the injector lip (Fig. 8.2) thanks to the coupling strategy
described above, and the second one (superscript A) with adiabatic boundary
conditions at the lip. Characteristic quantities of the one-dimensional
counterflow diffusion flame shown in Fig. 5.10 are used to nondimensionalize
some of the results: the average volumetric heat release rate
$\dot{\Omega}_{0}=3.5\times 10^{11}\rm\ W/m^{3}$ is used to normalize the heat
release computed in the numerical simulations (noted $HR^{*}$), and the heat
release per flame surface area $\Phi_{0}=40.4\ \rm MW/m^{2}$ (obtained by
integration over the flame thickness) is used to normalize the wall heat flux
(noted $\Phi_{w}^{*}$). Note also that the thickness of this one-dimensional
flame in the physical space is $\rm 125\ \mu m$, which is much larger than the
smallest cell size used in the present simulations.

#####  8.2.1 Flame root structure

Both simulations were run and averaged over 12 ms, which roughly corresponds
to 5 convective times relative to a fluid particle in the slow dense oxygen
core, or equivalently to 20 characteristic diffusion times in the solid.
Averaged fields of heat release rate, fluid temperature, wall temperature, and
heat flux are shown in Fig. 8.4.

Figure 8.4: (1) Averaged temperature field and velocity streamlines in the
fluid. Averaged heat flux magnitude and vector field in the wall. (2) Averaged
heat release rate and velocity vector field in the fluid. Averaged temperature
and heat flux streamlines in the wall. Superscripts A (resp. C) denote the
simulation with adiabatic conditions at the lip (resp. the conjugate heat
transfer problem).

The flame root structure obtained in the adiabatic case presents some
distinguishable features. In Fig. 8.4-(1A) two corotating recirculation zones
are visible: a first one ($R_{1}^{A}$) is located in the lip’s wake, while the
second one ($R_{2}^{A}$) generated by upstream flow separation lies beneath
the lip, in the oxygen injection stream. A smaller vortex appears to be
trapped and stretched between $R_{1}^{A}$ and $R_{2}^{A}$, where a pocket of
hot gases ($HG^{A}$) is trapped underneath the lip. The density isoline
$\rho=0.9\rho_{O2}^{inj}$ allows the visualization of the position of the
dense oxygen core, and this one appears to retract under the influence of
$HG^{A}$, which facilitates cryogenic O2 pseudo-boiling. In Fig. 8.4-(2A), two
distinct reaction zones can be identified: an intense primary reaction zone
($P_{R}^{A}$) is located in the lip’s wake, at a distance of approximately
$1.5\delta_{l}$ from the lip’s end, in a region of high shear between the two
streams. This stabilization mechanism is similar to that described in [191]
for a LO2/GH2 flame. However, unlike the LO2/GH2 flame, the present case also
displays a weaker secondary reaction zone ($S_{R}^{A}$), situated in the zone
of high shear between $R_{1}^{A}$ and $R_{2}^{A}$. This secondary reaction
zone, anchored onto the lip within the oxygen stream, is fueled through the
pocket of hot gases $HG^{A}$, which facilitates cryogenic reactants
vaporization, and therefore their mixing. This effect is noticeable as
$S_{R}^{A}$ lies close to the stoichiometry line $Z=0.2$: gaseous CH4 invades
the dense O2 stream and is responsible for the anchoring of the secondary
flame front beneath the lip. Compared to the LO2/GH2 flame studied in [191],
this stabilization mechanism is probably due to a significantly higher
momentum flux ratio $J=(\rho_{F}u_{F}^{2})/(\rho_{O}u_{O}^{2})$ ($J=33.3$
here, while $J=0.24$ in [191]). The stabilization mechanism in the adiabatic
case seems extremely dependent on the formation and the expansion of the
pocket of hot gases $HG^{A}$. The temperatures reached by the lip extremity
for this adiabatic case range from approximately $300$ K to $3500$ K, and are
not realistic.

On the contrary (see Fig. 8.4-(2C)), the fully coupled problem yields a
relatively cold time-averaged lip temperature, which does not exceed $450$ K
at its extremity. It is also worth noting the large temperature gradient
across the lip thickness, varying from $450$ K at the lower corner of the lip
($L_{low}$) to $300$ K at the point $L_{up}$ near its upper corner. The cold
lip induces significant heat losses at the wall (up to $\Phi_{w}^{*}=4.5$, see
Fig. 8.4-(1C)) and prevents the formation of the region of hot gases $HG^{A}$
seen for the adiabatic case, thereby changing the entire stabilization
mechanism. The two main recirculation zones $R_{1}^{C}$ and $R_{2}^{C}$ still
exist, but are slightly more separated than in the adiabatic case, so that no
area of high shear is found at their common frontier. Mixing is also notably
affected by wall heat losses: because of lower temperatures at the lip,
vaporization of reactants is less efficient, and gaseous CH4 is not
transported within the oxygen stream, as in the adiabatic case. Only one
primary region of reaction ($P_{R}^{C}$) is observed in the coupled
simulation. It is located close to the lip extremity, and extends over a long
narrow area spreading from the lower corner $L_{low}$ to the point $L_{up}$
closer to the upper corner. This suggests periodic oscillations of the flame,
as discussed below. A wider downstream flame brush also tends to indicate
large amplitude fluttering of the flame front induced by oscillations of the
flame root.

Detailed flame structures (Fig. 8.4-(3A) and (3C)) in the mixture fraction
space show that with adiabatic boundary conditions, the flame mostly burns in
a pure diffusion mode, as its structure is very close to that of the one-
dimensional flame presented in Fig. 5.10.

Figure 8.5: Scatter plots of the flame structure in mixture fraction space,
compared to the flame displayed in Fig. 5.10. Superscripts A (resp. C) denote
the simulation with adiabatic conditions at the lip (resp. the conjugate heat
transfer problem).

Conversely, the coupled problem gives a flame that burns in diffusion as well
as partially-premixed modes: reactions do not occur close to the wall due to
low temperatures, and reactants can mix before reacting.

#####  8.2.2 Flame root dynamics

In the adiabatic case, the primary reaction zone $P_{R}^{A}$ does not exhibit
significant oscillations, and similarly to the LO2/GH2 flame of [191], it is
stabilized in the highly sheared mixing layer between the two reactant
streams. The weaker secondary reaction zone $S_{R}^{A}$, anchored beneath the
lip is more sensitive to perturbations in the cryogenic oxidizer stream, and
undergoes slight intermittentcy. The flame root obtained with FWI effects is
more unstable and displays strong periodic self-sustained fluctuations. The
primary reaction zone $P_{R}^{C}$ oscillates between two positions $L_{low}$
located near the lower lip corner, and $L_{up}$ closer to the upper lip
corner. Meanwhile, the lip temperature field changes in phase with the flame
root motion: at the lower corner $L_{low}$, temperature varies between $345$ K
and $505$ K, while it varies between $278$ K and $355$ K at the upper corner.
Snapshots of temperature, heat release, wall heat flux, and wall temperature
fields at two instants $t_{up}$ and $t_{low}$, corresponding to the two
distinct flame root locations, are provided in Fig. 8.6.

Figure 8.6: (1) Instantaneous temperature field and heat release contours in
the fluid. The dark line represents the $\bm{\rho=0.9\rho_{O2}^{inj}}$
isoline. Instantaneous heat flux magnitude and vector field in the wall. (2)
Instantaneous heat release rate and velocity vector field in the fluid. The
blue line represents the stoichiometry line $\bm{Z=0.2}$. Instantaneous
temperature field with heat flux streamlines in the wall. Subscripts u (resp.
l) denote an instant when the flame root is located at the point $\bm{L_{up}}$
(resp. $\bm{L_{low}}$).

When $P_{R}^{C}$ is close to the point $L_{up}$, the flame is subjected to a
strong shear due to the proximity of the high momentum CH4 stream. As a
consequence, the flame root is a relatively straight and unperturbed segment,
followed by large scale vortices downstream. At the same time, wall heat
losses contribute to rapid heat-up of the upper corner of the lip
($\Phi_{w}^{*}=8.5$), resulting in a relatively homogeneous temperature field
across the lip thickness ($T_{Low}=360$ K and $T_{up}=340$ K). Conversely,
when the main reaction zone $P_{R}^{C}$ is located near the lower corner
$L_{low}$, shear is weaker due to low velocities in the oxygen stream. The
flame root is stretched and rolled up by large scale vortical structures in
the mixing layer, and it undergoes a flapping motion. In the meantime, the
wall temperature gradient across the lip thickness reaches its maximal value
($T_{Low}=490$ K and $T_{up}=280$ K), due to the proximity of the reaction
zone to the lip’s lower corner.

In order to further investigate coupling mechanisms responsible for the self-
sustained oscillations affecting simultaneously the flame root and the wall
temperature, two quantities of interest are recorded over time. The first one
is the radial location $y_{fr}(t)$ of the flame root (noted fr), and the
second one is the radial location $y_{wh}(t)$ of the wall hotspot (noted wh),
defined as the point of maximum temperature in the lip (Fig. 8.6 right). Both
locations are normalized by the lip thickness $\delta_{l}$. Coordinates origin
is taken at the point $L_{up}$ indicated in Fig. 8.6. Temporal evolution and
Fourier spectra of the two signals are shown in Fig. 8.7.

Figure 8.7: Left: time-series of the flame root location $\bm{y_{fr}(t)}$ and
the wall hotspot location $\bm{y_{wh}(t)}$. Right: Fourier spectra associated
to these two signals.

A strong correlation between flame root and wall hotspot fluctuations appears.
Both are dominated by two modes, with a first fundamental mode at roughly
$450$Hz and its harmonic at approximately $900$Hz. Thus, self-sustained
oscillations affecting the flame root seem to result from a complex flame/heat
transfer coupling at the wall. Pressure fluctuations were also recorded at
several locations in the domain, and as resulting spectra did not exhibit any
of the previous frequencies, the influence of acoustics was ruled out as a
potential cause of the unsteady flame anchoring. The wall hotspot and flame
root coupled dynamics can be further studied by considering the difference
$\Delta y(t)=y_{wh}(t)-y_{fr}(t)$ (still normalized by the lip thickness
$\delta_{l}$). One-dimensional Dynamic Mode Decomposition [243] is used to
isolate the most energetic components of the signals $y_{fr}(t)$ and
$y_{wh}(t)$, corresponding to the peaks in Fig. 8.7. The resulting phase-
portrait in the plane ($\Delta y(t),\dot{\Delta}y(t)$) is displayed in Fig.
8.8.

Figure 8.8: (a) Phase-portrait of the coupled system in the plane ($\bm{\Delta
y(t),\dot{\Delta}y(t)}$). Red circles indicate unstable equilibrium points at
cusps in the trajectory. Green squares represent unstable attracting points.
(b) Schematic view of an analogous forced double-well oscillator.

The origin ($0,0$) in the phase-plane represents instants where wall hotspot
and flame root face each other and move at the same speed. A quick examination
of the phase-portrait allows to rule out the existence of any distinct closed
limit-cycle, as the trajectory is chaotic and alternates between small orbits
around the origin (i.e. fr and wh are close to each other and have similar
motion), and larger orbits (i.e. fr and wh are far from each other). Most
importantly, 4 cusps located near the origin are clearly visible on the
trajectory showing that $y_{fr}\approx y_{wh}$ is an unstable equilibrium
position. After reaching a cusp, the trajectory is sent back to a larger
orbit. These dynamics are very similar to that of a particle in a forced
double-well oscillator, the central bump representing the unstable equilibrium
situation where wall hotspot and flame root face each other. The particle
cannot reach exactly the top of the central bump, and slight overshoots in its
displacement cause it to fall back on the other side and oscillate on a larger
orbit. Analogously, flame root and wall hotspot are never exactly in phase,
and overshoots in the flame displacement cause this equilibrium situation to
diverge to large values of $\Delta y$. Thus, the coupled dynamics consists of
a combination of small-amplitude oscillations around the equilibrium state
$y_{fr}\approx y_{wh}$, and larger-amplitude oscillations with the flame root
position $y_{fr}$ and the wall hotspot position $y_{wh}$ decorrelated from
each other. Physically, this unsteady coupling can be explained by unsteady
heat convection phenomena occurring in the fluid. (1) When the flame is
stabilized near the point $L_{up}$ large vortices downstream of the straight
flame root bring hot gases close to the lower corner $L_{low}$ (see in Fig.
8.6). The lip hotspot then moves from the point $L_{up}$ to $L_{low}$, and the
flame root follows its motion. (2) Conversely, when the flame is stabilized
near $L_{low}$, vortices convect hot gases towards the upper corner, and the
reverse mechanism takes place.

####  8.3 Conclusion

Adiabatic boundary conditions are often used in numerical simulations of flame
dynamics, especially in LRE conditions where little experimental data is
available. As a result, the effects of Flame-Wall Interactions on both the
anchoring of the flame root, and on its forced and intrinsic dynamics is a
problem that was mostly left unexplored. The present study circumvented the
lack of experimentally measured wall temperatures by resolving an unsteady
conjugate heat-transfer problem, where the temperature field within the
injector lip was strongly coupled to a two-dimensional high-fidelity
simulation of the doubly transcritical LO2/LCH4 coaxial flame studied
throughout this thesis. The use of the detailed kinetic scheme for CH4
oxycombustion in LRE conditions, derived in Sec. 5.3, was also a key feature
of this simulation.

Significantly different stabilization mechanisms are found between simulations
using adiabatic walls and the fully coupled computations. More precisely, in
the former case, the flame is mainly stabilized in a sheared mixing layer
between reactants streams, similarly to LO2/GH2 flames. However, adiabatic
walls appear to facilitate cryogenic reactants vaporization and therefore
their mixing. This effect results in the formation of a secondary reaction
front within the O2 injection stream. Conversely, the coupled problem exhibits
no stable mixing layer and strong self-sustained oscillations of the flame
root and of the temperature field in the lip. Both were found to be bimodal
and to oscillate at similar frequencies, indicating a strong influence of heat
transfer in the lip on the flame root location. As a result, successive flame
root locations were observed (1) in the proximity of the CH4 stream, where the
flame is stabilized by high shear, (2) and in the proximity of the O2 stream,
where the flame root is stretched and rolled-up by vortical structures. The
chaotic coupled dynamics of this system can be compared to that of a particle
sliding in a forced double-well oscillator. These self-sustained oscillations
are expected to be important for the flame response to external forcing, but
also for the overall flame stabilization process. Thus, the three-dimensional
flame dynamics LES performed in Chapter 6 and 7 would need to be recomputed by
adding the resolution of a conjugate heat transfer problem. However, as this
task requires immense computational resources, it was deemed unfeasible in the
present thesis.

Finally, note that this study assumes perfectly smooth injector walls.
Although this point was not investigated, flow features taking place at such
small scale are expected to be highly sensitive to wall roughness, which could
potentially alter the reported flame root stabilization mechanisms.

### Conclusions and Perspectives

This PhD thesis is dedicated to the development of methods and models intended
to enable the numerical prediction of thermoacoustic instabilities in complex
industrial systems, such as Liquid Rocket Engines or gas turbines. A wide
variety of computational tools, ranging from cheap low-order models to costly
high-fidelity simulations, were utilized and perfected. This work was carried
out following two principal directions: the first one consisted in designing
more accurate and efficient low-order models for the resolution of unstable
thermoacoustic modes in complex combustors, while the second one was centered
on the use of high-fidelity Large Eddy Simulations and Direct Numerical
Simulations to explore the complex dynamics in an academic high-pressure
cryogenic burner, representative of a Liquid Rocket Engine.

###### Low-Order Modeling of thermoacoustic instabilities in complex systems

The first axis of this PhD thesis lies in the vein of previous studies
conducted at CERFACS, most notably by Bauerheim et al. [106, 259, 107, 52]. It
appeared that usual low-order modeling methods rely on drastic simplifications
of the combustors’ geometries, which significantly alter their applicability
to industrial systems. The present work therefore aimed at improving these
existing methods by designing a novel low-order strategy that would allow for
the prediction of unstable thermoacoustic modes in combustion systems that are
of particular interest for CERFACS industrial partners. This was achieved
through the implementation of the code STORM, based on the following
fundamental principles:

  * •

The number of Degrees of Freedom is reduced as mush as possible, such that a
low computational cost is ensured. This was achieved by making use of modal
expansions, which generally yield lower-order problems than other strategies
commonly encountered.

  * •

The LOM is highly modular and flexible, which eases the straightforward
modification of physical or geometrical parameters. In this matter, the LOM
code STORM relies on a divide and conquer strategy, where complex combustors
are split into acoustic networks of simpler elements, thanks to the state-
space formalism.

  * •

Non-rigid-wall frontiers, such as complex impedance boundaries that are
ubiquitous in industrial combustors, can be accurately modeled. This
particular point was addressed thanks to the frame modal expansion, a
generalization of the classical Galerkin expansion, that unlike the latter, is
able to precisely resolve the acoustic fields in the vicinity of non-trivial
boundaries.

  * •

Geometrically complex boundaries, which usually consist of wide and curved
panels, can be accounted for. This is enabled by the combination of the frame
expansion with the newly introduced surface modal expansion method. This
latter approach is analogous to a Galerkin expansion, in the case of a
curvilinear Helmholtz problem formulated in the vicinity of a topologically
complex surface.

The major contribution of this thesis to the low-order modeling of
thermoacoustic instabilities is undoubtedly the frame expansion method, which
constitutes the backbone of STORM. The modularity of flexibility of the
approach was further enhanced thanks to the design of an Application
Programming Interface (API) intended to generate and manage the complex
acoustic networks that STORM can resolve. The implementation of this API was
initiated during this PhD thesis, and will be continued in future ones. Even
though this work did not present any application of STORM to actual industrial
combustors, such results were nonetheless recently obtained on a Safran
Aircraft Engines annular combustor. A strong emphasis was however put on
empirically assessing the convergence and the precision of the proposed
numerical methods. The author believes that the strategies introduced during
this PhD have the potential to drastically improve state-of-the-art
thermoacoustic low-order models. Future development at CERFACS will include
the use of STORM for Uncertainty Quantification, adjoint optimization, robust
design, as well as temporal nonlinear simulations of thermoacoustic limit-
cycles and other chaotic phenomena.

###### High-fidelity simulations of a doubly-transcritical LO2/LCH4 coaxial
jet-flame

The second research direction followed during this PhD was guided by the
urgent need to characterize the doubly-transcritical combustion of methane and
oxygen, a regime that is mostly unexplored, but that is nonetheless likely to
occur in future reusable Liquid Rocket Engines. In this matter, several high-
fidelity simulations of the Mascotte academic test rig, operating in
conditions relevant to that of an industrial LRE, were performed. The
specificities of these simulations lied in the non-ideal thermodynamics that
are resolved, as well as in the use of a complex kinetic scheme that was
purposely derived for high-pressure methane oxicombustion. The allocation of
considerable computational resources through the PRACE, GENCI, and MeteoFrance
infrastructures allowed several facets of doubly-transcritical combustion to
be explored:

  * •

The first LES of a doubly-transcritical LO2/ LCH4 coaxial jet flame was used
to investigate its structure as well as its intrinsic dynamics. These latter
displayed several distinctive features never reported before. Most notably,
break-up modes of the oxidizer dense core, as well as a vortex-shedding mode
in the fuel stream were evidenced.

  * •

The first Flame Transfer Function of a doubly-transcritical LO2/LCH4 coaxial
jet flame was computed over a wide frequency range. This was achieved thanks
to a series of LES where small-amplitude harmonic modulations were imposed at
the fuel inlet to trigger linear flame fluctuations. Frequency-dependent
preferential response regions were identified, and were observed to be
strongly correlated to the generation of vortex-ring pairs in the methane
annular jet. In addition, an analysis of the mechanisms contributing to the
heat-release oscillations proved the prevalence of the species diffusivity
variations.

  * •

The nonlinear response of a LO2/LCH4 coaxial jet flame was assessed by
performing another set of similar LES, at a larger modulation amplitude. The
first Flame Describing Function of such flame was computed over a large
frequency range. Saturation mechanisms were identified, as well as nonlinear
production of higher harmonics.

  * •

Finally, the first DNS of a doubly-transcritical LO2/LCH4 flame utilizing both
complex chemistry and a conjugate heat transfer problem at the injector lip
was conducted. This high-fidelity simulation aimed at assessing the effects of
Flame-Wall Interactions on the flame root anchoring mechanism. Remarkably
different stabilization modes were observed between this non-adiabatic flame
and an adiabatic one. In particular, the strongly coupled unsteady conjugate
heat transfer problem yielded coherent self-sustained oscillations affecting
both the flame root location and the temperature field within the injector
lip. Such complex anchoring mechanism is expected to have first order effects
on the overall flame topology and on its response to acoustic pulsations,
which questions the use of adiabatic boundary conditions in simulations of LRE
combustion.

High-fidelity simulations’ results discussed in this work are expected to
provide a new and deeper insight into the dynamics of LRE cryogenic flames.
The proposed analyzes of these results have the potential to guide the future
development of theoretical models for the response of such flames to acoustic
perturbations. In addition, quantitative transfer functions computed during
this PhD thesis can be embedded into low-order acoustic networks for the
prediction of thermoacoustic instabilities in multi-injectors LREs.

### Publications

###### Ranked A papers

This PhD thesis led to the following publications:

  1. 1.

Laurent, C., Bauerheim, M., Poinsot, T., & Nicoud, F. (2019). A novel modal
expansion method for low-order modeling of thermoacoustic instabilities in
complex geometries. Combustion and Flame, 206, 334-348.

  2. 2.

Laurent, C., Badhe, A., Nicoud, F. (2020). Including complex boundary
conditions in low-order modeling of thermoacoustic instabilities. (Submitted
to Combustion and Flame)

  3. 3.

Laurent, C., Esclapez, L., Maestro, D., Staffelbach, G., Cuenot, B., Selle,
L., Schmitt, T., Duchaine, F., Poinsot, T. (2019). Flame wall interaction
effects on the flame root stabilization mechanisms of a doubly-transcritical
LO2/LCH4 cryogenic flame. Proceedings of the Combustion Institute, 37(4),
5147-5154.

  4. 4.

Laurent, C., Staffelbach, G., Nicoud, F., Poinsot, T. (2020). Heat-release
dynamics in a doubly-transcritical LO2/LCH4 cryogenic coaxial jet flame
subjected to fuel inflow acoustic modulation. (Proceedings of the Combustion
Institute, Accepted)

  5. 5.

Laurent, C., Poinsot, T. (2020). Nonlinear dynamics in the response of a
doubly-transcritical LO2/LCH4 coaxial jet flame to fuel inflow acoustic
perturbations. (Combustion and Flame, Submitted)

###### Conferences

  * •

Detailed Chemistry effects on the flame root stabilization mechanisms of a
doubly-transcritical LO2/LCH4 flame, 4eme Colloque INCA, Paris, October 2017.

  * •

Flame wall interaction effects on the flame root stabilization mechanisms of a
doubly-transcritical LO2/LCH4 cryogenic flame, Journée François Lacas,
Orléans, January 2018.

  * •

Modal expansions for accurate thermoacoustic low-order models, CRCT, Toulouse,
March 2018.

  * •

Flame wall interaction effects on the flame root stabilization mechanisms of a
doubly-transcritical LO2/LCH4 cryogenic flame, 37th International Symposium on
Combustion, Dublin, August 2018.

###  Appendix A A library of state-space realizations

####  A.1 Acoustic subdomain

The state-space representation for a subdomain $\Omega_{i}$ belonging to an
acoustic network is given in this section. It is assumed that the connection
boundary $S_{ci}$ consists of $M_{S}$ sub-surfaces $(S_{ci}^{(m)})_{1\leq
m\leq M_{S}}$ of surface area $\Delta S_{ci}^{(m)}$, and located at the
center-points $\vec{x}_{sm}$. It is also supposed that these sub-surfaces are
small enough, such that acoustic variables can be approximated by their values
at the points $\vec{x}_{sm}$. In addition, $\Omega_{i}$ contains $M_{H}$ heat-
sources $(H_{i}^{(l)})_{1\leq l\leq M_{H}}$. The block structure of the
dynamics equation, derived from Eq. (3.28), writes:

$\displaystyle\begin{aligned}
&\dfrac{d}{dt}\underbrace{\begin{pmatrix}\Gamma_{1}(t)\\\
\dot{\Gamma}_{1}(t)\\\ \vdots\\\ \Gamma_{N}(t)\\\ \dot{\Gamma}_{N}(t)\\\
\end{pmatrix}}_{\mathbf{X}^{\Omega_{i}}(t)}=\underbrace{\left(\begin{array}[]{ccccc}0&1&&&\\\
-\omega_{1}^{2}&-\delta&&&\\\ &&\ddots&&\\\ &&&0&1\\\
&&&-\omega_{N}^{2}&-\delta\\\
\end{array}\right)}_{\mathbf{A}^{\Omega_{i}}}\begin{pmatrix}\Gamma_{1}(t)\\\
\dot{\Gamma}_{1}(t)\\\ \vdots\\\ \Gamma_{N}(t)\\\ \dot{\Gamma}_{N}(t)\\\
\end{pmatrix}+\\\\[5.0pt]
&\underbrace{\left(\begin{array}[]{cccccc}\mathbf{B}_{1}^{S_{ci}^{(1)}}&\ldots&\mathbf{B}_{1}^{S_{ci}^{(M_{S})}}&\mathbf{B}_{1}^{H_{i}^{(1)}}&\ldots&\mathbf{B}_{1}^{H_{i}^{(M_{H})}}\\\
\vdots&\ddots&\vdots&\vdots&\ddots&\vdots\\\
\mathbf{B}_{N}^{S_{ci}^{(1)}}&\ldots&\mathbf{B}_{N}^{S_{ci}^{(M_{S})}}&\mathbf{B}_{N}^{H_{i}^{(1)}}&\ldots&\mathbf{B}_{N}^{H_{i}^{(M_{H})}}\\\
\end{array}\right)}_{\mathbf{B}^{\Omega_{i}}}\underbrace{\begin{pmatrix}\mathbf{U}^{S_{ci}^{(1)}}\\\
\vdots\\\ \mathbf{U}^{S_{ci}^{(M_{S})}}\\\ Q_{1}(t)\\\ \vdots\\\
Q_{M_{H}}(t)\\\ \end{pmatrix}}_{\mathbf{U}^{\Omega_{i}}(t)}\end{aligned}$
(A.1)

where the blocks $\mathbf{B}_{n}^{S_{ci}^{(m)}}$, $\mathbf{U}^{S_{ci}^{(m)}}$
(used to compute the surface source terms from the boundaries $S_{ci}^{(m)}$),
and $\mathbf{B}_{n}^{H_{i}^{(l)}}$ (used to compute the volume source terms
from the heat-sources $H_{i}^{(l)}$) are given by:

$\displaystyle\begin{aligned}
&\mathbf{B}_{n}^{S_{ci}^{(m)}}=\rho_{0}c_{0}^{2}\Delta
S_{ci}^{(m)}\left(\begin{array}[]{cc}0&0\\\
-\left[\mathbf{\Lambda}^{-1}\bm{\phi}(\vec{x}_{sm})\right]_{n}&\left[\mathbf{\Lambda}^{-1}\bm{\nabla_{s}\phi}(\vec{x}_{sm})\right]_{n}\\\
\end{array}\right)\\\ &\mathbf{B}_{n}^{H_{i}^{(l)}}=\begin{pmatrix}0\\\
(\gamma-1)\left[\mathbf{\Lambda}^{-1}\langle\bm{\phi},\mathcal{H}_{i}^{(l)}\rangle\right]_{n}\\\
\end{pmatrix}\ ,\
\mathbf{U}^{S_{ci}^{(m)}}=\begin{pmatrix}u_{s}^{S_{ci}^{(m)}}(\vec{x}_{sm},t)\\\
\varphi^{S_{ci}^{(m)}}(\vec{x}_{sm},t)\\\ \end{pmatrix}\end{aligned}$ (A.2)

In this equation, the state vector $\mathbf{X}^{\Omega_{i}}(t)$ is of size
$2N$, where $N$ is the number of eigenmodes used in the modal expansion. The
dynamics matrix $\mathbf{A}^{\Omega_{i}}$ is block-diagonal of size $2N\times
2N$. The first $2M_{S}$ columns of the input matrix $\mathbf{B}^{\Omega_{i}}$
and $2M_{S}$ elements of the input vector $\mathbf{U}^{\Omega_{i}}(t)$
correspond to surface source terms imposed at adjacent boundaries
$S_{ci}^{(m)}$. The $M_{H}$ last columns of $\mathbf{B}^{\Omega_{i}}$ and the
$M_{H}$ last elements of $\mathbf{U}^{\Omega_{i}}(t)$ correspond to the
volumetric heat release forcing.

In the state-space framework that is used here, any subsystem $\Omega_{i}$
outputs both the normal velocity
$u_{s}(\vec{x}_{sm})=\vec{u}(\vec{x}_{sm}).\vec{n}_{s}$ and the pressure
$p(\vec{x}_{sm})$ at each one of the $M_{S}$ surface elements. In addition, it
is also possible to incorporate in the output vector (not detailed here)
pressure and velocity at any point within $\Omega_{i}$, such that those can
then be passed as reference pressure/velocity to an active flame. Thus, the
following equation is used to compute the output vector for the subdomain
$\Omega_{i}$:

$\displaystyle\underbrace{\begin{pmatrix}\mathbf{Y}^{S_{ci}^{(1)}}\\\
\vdots\\\ \mathbf{Y}^{S_{ci}^{(M_{S})}}\\\
\end{pmatrix}}_{\mathbf{Y}^{\Omega_{i}}(t)}=\underbrace{\left(\begin{array}[]{ccc}\mathbf{C}_{1}^{S_{ci}^{(1)}}&\ldots&\mathbf{C}_{N}^{S_{ci}^{(1)}}\\\
\vdots&\ddots&\vdots\\\
\mathbf{C}_{1}^{S_{ci}^{(M_{S})}}&\ldots&\mathbf{C}_{N}^{S_{ci}^{(M_{S})}}\\\
\end{array}\right)}_{\mathbf{C}^{\Omega_{i}}}\underbrace{\begin{pmatrix}\Gamma_{1}(t)\\\
\dot{\Gamma}_{1}(t)\\\ \vdots\\\ \Gamma_{N}(t)\\\ \dot{\Gamma}_{N}(t)\\\
\end{pmatrix}}_{\mathbf{X}^{\Omega_{i}}(t)}$ (A.6)

where the blocks $\mathbf{C}_{n}^{S_{ci}^{(m)}}$ and
$\mathbf{Y}^{S_{ci}^{(m)}}$ are expressed as:

$\displaystyle\mathbf{Y}^{S_{ci}^{(m)}}=\begin{pmatrix}u_{s}(\vec{x}_{sm})\\\
p(\vec{x}_{sm})\\\ \end{pmatrix}\ ,\
\mathbf{C}_{n}^{S_{ci}^{(m)}}=\left(\begin{array}[]{cc}-\dfrac{1}{\rho_{0}}\nabla_{s}\phi_{n}(\vec{x}_{sm})&0\\\
0&\phi_{n}(\vec{x}_{sm})\\\ \end{array}\right)$ (A.9)

Note that the feedthrough matrix $\mathbf{D}$ is zero. At the end, the state-
space representation of the acoustics in the subdomain $\Omega_{i}$ is defined
by Eq. (A.1) to Eq. (A.9). It is a $2M_{S}$-input $2M_{S}$-output system (in
the absence of heat release source terms), whose dynamics are described by a
$2N\times 2N$ matrix.

####  A.2 Acoustic subdomain with a complex connection boundary

The state-space representation for the acoustics in a subdomain $\Omega_{i}$
needs to be adapted to account for complex impedance or conductivity on
geometrically elaborate boundaries. For conciseness, $\Omega_{i}$ is assumed
to be adjacent to only one complex boundary $S_{ci}^{(m)}$. Since other types
of state-space interconnections (simple point-wise boundaries, heat sources,
etc.) do not require any adaptation, they are not considered here. The
reformulated dynamics equation, formulated thanks to Eq. (4.2) and Eq. (4.32),
has the following block structure:

$\displaystyle\begin{aligned}
&\dfrac{d}{dt}\underbrace{\begin{pmatrix}\Gamma_{1}(t)\\\
\dot{\Gamma}_{1}(t)\\\ \vdots\\\ \Gamma_{N}(t)\\\ \dot{\Gamma}_{N}(t)\\\
\end{pmatrix}}_{\mathbf{X}^{\Omega_{i}}(t)}=\underbrace{\left(\begin{array}[]{ccccc}0&1&&&\\\
-\omega_{1}^{2}&-\alpha&&&\\\ &&\ddots&&\\\ &&&0&1\\\
&&&-\omega_{N}^{2}&-\alpha\\\
\end{array}\right)}_{\mathbf{A}^{\Omega_{i}}}\begin{pmatrix}\Gamma_{1}(t)\\\
\dot{\Gamma}_{1}(t)\\\ \vdots\\\ \Gamma_{N}(t)\\\ \dot{\Gamma}_{N}(t)\\\
\end{pmatrix}+\\\\[5.0pt]
&\underbrace{\left(\begin{array}[]{ccc}\mathbf{B}_{1,1}^{S_{ci}^{(m)}}&\ldots&\mathbf{B}_{1,K_{S}}^{S_{ci}^{(m)}}\\\
\vdots&\ddots&\vdots\\\
\mathbf{B}_{N,1}^{S_{ci}^{(m)}}&\ldots&\mathbf{B}_{N,K_{S}}^{S_{ci}^{(m)}}\\\
\end{array}\right)}_{\mathbf{B}^{\Omega_{i}}}\underbrace{\begin{pmatrix}\mathbf{U}_{1}^{S_{ci}^{(m)}}\\\
\vdots\\\ \mathbf{U}_{K_{S}}^{S_{ci}^{(m)}}\\\
\end{pmatrix}}_{\mathbf{U}^{\Omega_{i}}(t)}\\\ \end{aligned}$ (A.10)

where the blocks $\mathbf{B}_{n,k}^{S_{ci}^{(m)}}$ and
$\mathbf{U}_{k}^{S_{ci}^{(m)}}$, used to compute the surface forcing from the
complex boundary $S_{ci}^{(m)}$, are given by:

$\displaystyle\begin{aligned}
&\mathbf{B}_{n,k}^{S_{ci}^{(m)}}=\rho_{0}c_{0}^{2}\left(\begin{array}[]{cc}0&0\\\
-\mathpzc{s}_{i}\left[\mathbf{\Lambda}^{-1}\left(\bm{\phi}|\mathscr{K}_{k}^{(m)}\right)\right]_{n}&\left[\mathbf{\Lambda}^{-1}\left(\bm{\nabla_{s}\phi}|\mathscr{K}_{k}^{(m)}\right)\right]_{n}\\\
\end{array}\right)\\\
&\mathbf{U}_{k}^{S_{ci}^{(m)}}=\dfrac{1}{\lambda_{k}}\begin{pmatrix}\mathpzc{s}_{i}\left(u_{s}^{S_{ci}^{(m)}}|\mathscr{K}_{k}\right)\\\
\left(\varphi^{S_{ci}^{(m)}}|\mathscr{K}_{k}\right)\end{pmatrix}=\begin{pmatrix}\mu_{k}(t)\\\
\nu_{k}(t)\end{pmatrix}\end{aligned}$ (A.11)

where $\mathpzc{s}_{i}=\vec{u}_{s}^{\Omega_{i}}.\vec{e}_{\xi}=\pm 1$ defines
the orientation of the subdomain surface normal with respect to the complex
boundary $S_{ci}^{(m)}$. The state vector $\mathbf{X}^{\Omega_{i}}(t)$ is of
size $2N$, with $N$ the number of eigenmodes used in the subdomain frame
expansion. The input matrix $\mathbf{B}^{\Omega_{i}}$ ($2N\times 2K_{S}$)
serves to compute the surface source terms from the $2K_{S}$ elements of the
input vector $\mathbf{U}^{\Omega_{i}}$, which comprises the surface modal
amplitudes of both the normal velocity and the acoustic potential from the
complex boundary $S_{ci}^{(m)}$.

The subsystem $\Omega_{i}$ outputs both the projections of its normal velocity
$u_{s}(\vec{x}_{s})$ and pressure $p(\vec{x}_{s})$ onto each one of the
$K_{S}$ surface modes of $S_{ci}^{(m)}$. This is achieved thanks to the output
equation:

$\displaystyle\underbrace{\begin{pmatrix}\mathbf{Y}_{1}^{S_{ci}^{(m)}}\\\
\vdots\\\ \mathbf{Y}_{K_{S}}^{S_{ci}^{(m)}}\\\
\end{pmatrix}}_{\mathbf{Y}^{\Omega_{i}}(t)}=\underbrace{\left(\begin{array}[]{ccc}\mathbf{C}_{1,1}^{S_{ci}^{(m)}}&\ldots&\mathbf{C}_{1,N}^{S_{ci}^{(m)}}\\\
\vdots&\ddots&\vdots\\\
\mathbf{C}_{K_{S},1}^{S_{ci}^{(m)}}&\ldots&\mathbf{C}_{K_{S},N}^{S_{ci}^{(m)}}\\\
\end{array}\right)}_{\mathbf{C}^{\Omega_{i}}}\underbrace{\begin{pmatrix}\Gamma_{1}(t)\\\
\dot{\Gamma}_{1}(t)\\\ \vdots\\\ \Gamma_{N}(t)\\\ \dot{\Gamma}_{N}(t)\\\
\end{pmatrix}}_{\mathbf{X}^{\Omega_{i}}(t)}$ (A.15)

where the blocks $\mathbf{C}_{k,n}^{S_{ci}^{(m)}}$ and
$\mathbf{Y}_{k}^{S_{ci}^{(m)}}$ are expressed as:

$\displaystyle\mathbf{Y}_{k}^{S_{ci}^{(m)}}=\mathpzc{s}_{i}\begin{pmatrix}\left(u_{s}|\mathscr{K}_{k}\right)\\\
\left(p|\mathscr{K}_{k}\right)\end{pmatrix}\ ,\
\mathbf{C}_{k,n}^{S_{ci}^{(m)}}=\mathpzc{s}_{i}\left(\begin{array}[]{cc}-\dfrac{1}{\rho_{0}}\left(\nabla_{s}\phi_{n}|\mathscr{K}_{k}\right)&0\\\
0&\left(\phi_{n}|\mathscr{K}_{k}\right)\\\ \end{array}\right)$ (A.18)

Finally, the state-space realization of the subdomain $\Omega_{i}$ with one
adjacent complex boundary $S_{ci}^{(m)}$ is a 2$K_{S}$-inputs-2$K_{S}$-outputs
systems, whose intrinsic dynamics are governed by a $2N\times 2N$ matrix.

####  A.3 Cross-section change between two long ducts

In this section, the state-space representation for an acoustically compact
subdomain $\Omega_{sc}$ enclosing a cross-section change between two long
ducts $\Omega_{1}$ and $\Omega_{2}$ of respective cross-section area $S_{1}$
and $S_{2}$ is derived. Although the acoustic field can be considered purely
longitudinal in both ducts $\Omega_{1}$ and $\Omega_{2}$, it is locally
multidimensional in the neighborhood of the section change. The subsystem
$\Omega_{sc}$ is therefore defined such that it encloses the region where the
acoustic flow is multidimensional. In practice, this translates into
$L_{sc}\simeq\max(\sqrt{D_{1}},\sqrt{D_{2}})$. However, since the cross-
section change is assumed acoustically compact (since $L_{sc}\ll
L_{1},L_{2}$), the detailed topological structure of the acoustic flow in
$\Omega_{sc}$ has very little effect on the global eigenfrequencies and
eigenmodes of the whole geometry. For this reason, the acoustic pressure in
$\Omega_{sc}$ is not expanded onto an eigenmodes family, but instead a
specific treatment relying on a more global description is applied. The
volume-averaged pressure $\overline{p}(t)$ and velocity $\overline{u}(t)$ in
$\Omega_{sc}$ are introduced. Volume-averaging the linearized Euler equations
in $\Omega_{sc}$ yields (acoustic losses due to vorticity conversion being
neglected):

$\displaystyle\left\\{\begin{aligned}
&\dfrac{d\overline{u}}{dt}=\dfrac{1}{\rho_{0}L_{sc}}(p^{\Omega_{1}}(t)-p^{\Omega_{2}}(t))\\\
&\dfrac{d\overline{p}}{dt}=\dfrac{2\rho_{0}c_{0}^{2}}{L_{sc}(S_{1}+S_{2})}(S_{1}u_{s}^{\Omega_{1}}(t)+S_{2}u_{s}^{\Omega_{2}}(t))\end{aligned}\right.$
(A.19)

where $p^{\Omega_{1}}(t)$ and $u_{s}^{\Omega_{1}}(t)$ (resp.
$p^{\Omega_{2}}(t)$ and $u_{s}^{\Omega_{2}}(t)$) are the pressure and the
normal acoustic velocity imposed by the subdomain $\Omega_{1}$ (resp.
$\Omega_{2}$) onto $\Omega_{sc}$. The + sign in the second equation is due to
the fact that the outer normals at $S_{c1}$ and $S_{c2}$ are defined with
respect to $\Omega_{1}$ and $\Omega_{2}$ respectively, and are therefore
pointing in opposite directions. For low frequencies, these conservation
equations reduce to the classical quasi-static jump relations
$p^{\Omega_{1}}=p^{\Omega_{2}}$ and
$S_{1}u_{s}^{\Omega_{1}}=-S_{2}u_{s}^{\Omega_{2}}$. Pressures and normal
velocities imposed by $\Omega_{sc}$ onto the two ducts $\Omega_{1}$ and
$\Omega_{2}$ can then be calculated thanks to 0th-order approximations (i.e.
piece-wise constant functions):

$\displaystyle\left\\{\begin{aligned}
&u_{s}(x_{sc}=0,t)=\dfrac{1}{2}(1+\dfrac{S_{2}}{S_{1}})\overline{u}(t)\ ,\
p(x_{sc}=0,t)=\overline{p}(t)\\\
&u_{s}(x_{sc}=L_{sc},t)=-\dfrac{1}{2}(1+\dfrac{S_{1}}{S_{2}})\overline{u}(t)\
,\ p(x_{sc}=L_{sc},t)=\overline{p}(t)\end{aligned}\right.$ (A.20)

Utilizing these relations, the state-space realization of the subdomain
$\Omega_{sc}$ then writes:

$\displaystyle\begin{aligned}
\dfrac{d}{dt}&\underbrace{\begin{pmatrix}\overline{u}(t)\\\ \overline{p}(t)\\\
\overline{\varphi}(t)\end{pmatrix}}_{\mathbf{X}^{(sc)}(t)}=\underbrace{\begin{pmatrix}0&0&0\\\
0&0&0\\\
0&-1/\rho_{0}&0\end{pmatrix}}_{\mathbf{A}^{(sc)}}\begin{pmatrix}\overline{u}(t)\\\
\overline{p}(t)\\\ \overline{\varphi}(t)\end{pmatrix}\\\
&+\underbrace{\begin{pmatrix}0&\dfrac{1}{\rho_{0}L_{sc}}&0&-\dfrac{1}{\rho_{0}L_{sc}}\\\
\dfrac{2S_{1}c_{0}^{2}\rho_{0}}{L_{sc}(S_{1}+S_{2})}&0&\dfrac{2S_{2}c_{0}^{2}\rho_{0}}{L_{sc}(S_{1}+S_{2})}&0\\\
0&0&0&0\end{pmatrix}}_{\mathbf{B}^{(sc)}}\underbrace{\begin{pmatrix}u_{s}^{\Omega_{1}}(t)\\\
p^{\Omega_{1}}(t)\\\ u_{s}^{\Omega_{2}}(t)\\\
p^{\Omega_{2}}(t)\end{pmatrix}}_{\mathbf{U}^{(sc)}(t)}\end{aligned}$ (A.21)

The state-vector $\mathbf{X}^{(sc)}(t)$ contains the volume-averaged velocity
$\overline{u}(t)$, pressure $\overline{p}(t)$, and acoustic potential
$\overline{\varphi}(t)=-\int\overline{p}(t^{\prime})dt^{\prime}/\rho_{0}$. The
first two components of the input vector $\mathbf{U}^{(sc)}(t)$ are imposed by
the duct $\Omega_{1}$, while the last two are imposed by the duct
$\Omega_{2}$. We then note $\mathbf{B}_{1}^{(sc)}$ (resp.
$\mathbf{B}_{2}^{(sc)}$) the matrix consisting of the first two columns (resp.
last two columns) of $\mathbf{B}^{(sc)}$. In Eq. (A.21), the first two lines
essentially impose acoustic momentum and acoustic volume flux conservation.
For low frequencies, these conservation relations reduce to the classical
quasi-static jump relations $p^{\Omega_{1}}=p^{\Omega_{2}}$ and
$S_{1}u_{s}^{\Omega_{1}}=-S_{2}u_{s}^{\Omega_{2}}$. The third equation of the
dynamical system is a time-integrator that facilitates the computation of the
output vector, as the two ducts $\Omega_{1}$ and $\Omega_{2}$ require normal
acoustic velocity and acoustic potential as inputs. Outputs for the state-
space representation of $\Omega_{sc}$ are computed thanks to first-order
approximations:

$\displaystyle\underbrace{\begin{pmatrix}u_{s}(0,t)\\\ \varphi(0,t)\\\
u_{s}(L_{sc},t)\\\ \varphi(L_{sc},t)\end{pmatrix}}_{\mathbf{Y}^{(sc)}(t)}$
$\displaystyle=\underbrace{\begin{pmatrix}\dfrac{1}{2}(1+\dfrac{S_{2}}{S_{1}})&0&0\\\
0&0&1\\\ -\dfrac{1}{2}(1+\dfrac{S_{1}}{S_{2}})&0&0\\\
0&0&1\end{pmatrix}}_{\mathbf{C}^{(sc)}}\begin{pmatrix}\overline{u}(t)\\\
\overline{p}(t)\\\ \overline{\varphi}(t)\end{pmatrix}$ (A.22)

In Eq. (A.22), the feedthrough matrix $\mathbf{D}$ is zero. The first two
lines are outputs that are to be imposed to the first duct $\Omega_{1}$, while
the last two lines are outputs that are to be imposed to the second duct
$\Omega_{2}$. We then note $\mathbf{C}_{1}^{(sc)}$ (resp.
$\mathbf{C}_{2}^{(sc)}$) the matrix formed with the first two rows (resp. last
two rows) of $\mathbf{C}^{(sc)}$.

The specific treatment applied to the subdomain $\Omega_{sc}$ enclosing the
neighboring region of the section change allows to derive a generic state-
space representation for acoustically compact cross-section variations,
independently of the exact geometry of this region. Indeed the use of volume-
average acoustic variables and 0th-order approximations, permits to avoid
dealing with the topological complexity of the multidimensional acoustic flow
in $\Omega_{sc}$. It is also worth noting that even though the derivation of
the state-space representation for $\Omega_{sc}$ was here conducted in the
absence of acoustic sources, and without conversion from acoustic to entropy
or vorticity waves in the neighborhood of the section change, these
assumptions could be omitted to obtain more general representations.

####  A.4 Cross-section change between a long duct a large cavity

In the case of a small domain of thickness $L_{sc}$ enclosing a cross-section
change $\Omega_{sc}$ located between a long duct $\Omega_{1}$ giving on a
large cavity $\Omega_{2}$, the flux conservation equation
$S_{1}u_{s}^{\Omega_{1}}=S_{2}u_{s}^{\Omega_{2}}$ is degenerate (since
$S_{1}\ll S_{2}$ and $u_{s}^{\Omega_{2}}\ll u_{s}^{\Omega_{1}}$). The state-
space realization of Appendix A.3 therefore requires some adaptation:

$\displaystyle\dfrac{d}{dt}$
$\displaystyle\underbrace{\begin{pmatrix}\overline{u}(t)\end{pmatrix}}_{\mathbf{X}^{(sc)}(t)}=\underbrace{\begin{pmatrix}0\\\
\end{pmatrix}}_{\mathbf{A}^{(sc)}}\begin{pmatrix}\overline{u}(t)\end{pmatrix}+\underbrace{\begin{pmatrix}0&\dfrac{1}{\rho_{0}L_{sc}}&0&-\dfrac{1}{\rho_{0}L_{sc}}\end{pmatrix}}_{\mathbf{B}^{(sc)}}\underbrace{\begin{pmatrix}u_{s}^{\Omega_{1}}(t)\\\
p^{\Omega_{1}}(t)\\\ u_{s}^{\Omega_{2}}(t)\\\
p^{\Omega_{2}}(t)\end{pmatrix}}_{\mathbf{U}^{(sc)}(t)}$ (A.23)

In the low-frequency limit, this dynamics equation only imposes the pressure
continuity between both subdomains ($p^{\Omega_{1}}(t))=p^{\Omega_{2}}(t)$),
without affecting the velocity. It is therefore necessary to account for the
acoustic flux exchange in the output equation, such that the two subdomains
are fully coupled. This is achieved by considering a non-zero feedthrough
matrix:

$\displaystyle\underbrace{\begin{pmatrix}u_{s}(0,t)\\\ \varphi(0,t)\\\
u_{s}(L_{sc},t)\\\ \varphi(L_{sc},t)\end{pmatrix}}_{\mathbf{Y}^{(sc)}(t)}$
$\displaystyle=\underbrace{\begin{pmatrix}1\\\ 0\\\ 0\\\
0\end{pmatrix}}_{\mathbf{C}^{(sc)}}\begin{pmatrix}\overline{u}(t)\end{pmatrix}+\underbrace{\begin{pmatrix}0&0&0&0\\\
0&0&0&0\\\ -1&0&0&0\\\
0&0&0&0\end{pmatrix}}_{\mathbf{D}^{(sc)}}\begin{pmatrix}u_{s}^{\Omega_{1}}(t)\\\
p^{\Omega_{1}}(t)\\\ u_{s}^{\Omega_{2}}(t)\\\ p^{\Omega_{2}}(t)\end{pmatrix}$
(A.24)

####  A.5 Active flame approximated by a Multi-Pole expansion

For a given active flame, fluctuations of heat release rate $Q(t)$ are
governed by a FTF, that is approximated by a PBF expansion such as the one of
Eq. (3.45). In Chapter 3, that is applied to a constant time-delay FTF
$e^{-j\omega\tau}$, but this type of approximation is more generic, and its
coefficients can be adjusted to fit FTF data points obtained through
experiments or through pulsed LES. Recasting each Pole Base Function into the
time-domain leads to the state-space realization of the flame:

$\displaystyle\begin{aligned}
&\dfrac{d}{dt}\underbrace{\begin{pmatrix}Z_{1}(t)\\\ \dot{Z}_{1}(t)\\\
\vdots\\\ Z_{M_{PBF}}(t)\\\ \dot{Z}_{M_{PBF}}(t)\\\
\end{pmatrix}}_{\mathbf{X}^{(FTF)}(t)}=\underbrace{\left(\begin{array}[]{ccccc}0&1&&&\\\
-\omega_{01}^{2}&2c_{1}&&&\\\ &&\ddots&&\\\ &&&0&1\\\
&&&-\omega_{0M_{PBF}}^{2}&2c_{M_{PBF}}\\\
\end{array}\right)}_{\mathbf{A}^{(FTF)}}\begin{pmatrix}Z_{1}(t)\\\
\dot{Z}_{1}(t)\\\ \vdots\\\ Z_{M_{PBF}}(t)\\\ \dot{Z}_{M_{PBF}}(t)\\\
\end{pmatrix}+\underbrace{\begin{pmatrix}0\\\ -1/\overline{u}\\\ \vdots\\\
0\\\ -1/\overline{u}\\\
\end{pmatrix}}_{\mathbf{B}^{(FTF)}}\underbrace{\begin{pmatrix}u^{(ref)}(t)\\\
\end{pmatrix}}_{\mathbf{U}^{(FTF)}(t)}\\\ \end{aligned}$ (A.25)

In Eq. (A.25), the flame input vector $\mathbf{U}^{(FTF)}(t)$ has a single
entry, the reference fluctuating velocity $u^{(ref)}(t)$. The input matrix
$\mathbf{B}^{(FTF)}$ is simply used to normalize the reference velocity
fluctuations by the mean flow at the reference point. Note that other
formulations including a pressure value as reference are also possible. The
state-vector $\mathbf{X}^{(FTF)}(t)$ contains abstract variables that serve as
intermediates in the calculation of heat release rate. The state variable
$\dot{Z}_{k}(t)$ can be interpreted as the proportion of normalized heat
release fluctuating at frequencies contained in the band of width $2c_{k}$
centered around $\omega_{0k}$. This state-space realization is completed by
the following output equation:

$\displaystyle\begin{aligned} \underbrace{\begin{pmatrix}Q(t)\\\
\end{pmatrix}}_{\mathbf{Y}^{(FTF)}(t)}=\underbrace{\left(\begin{array}[]{ccccc}0&-2\overline{Q}a_{1}&\ldots&0&-2\overline{Q}a_{M_{PBF}}\\\
\end{array}\right)}_{\mathbf{C}^{(FTF)}}\begin{pmatrix}Z_{1}(t)\\\
\dot{Z}_{1}(t)\\\ \vdots\\\ Z_{M_{PBF}}(t)\\\ \dot{Z}_{M_{PBF}}(t)\\\
\end{pmatrix}\\\ \end{aligned}$ (A.26)

The flame output vector $\mathbf{Y}^{(FTF)}(t)$ only comprises the heat
release rate $Q(t)$, which is reconstructed from a linear combination of the
individual components $\dot{Z}_{k}(t)$ contained in the state-vector.

###### Remark

Once a state-space realization (based on a PBF expansion) has been obtained
for a FTF with a given time-delay $\tau$, deriving a new state-space
representation for another FTF with a time-delay $\tau^{\prime}$ is
straightforward and does not require to repeat the PBF expansion fitting
procedure. Instead, the following relations can be used:

$\displaystyle\exp\left(-j\omega\tau^{\prime}\right)=\exp\left(-j\left(\omega\dfrac{\tau^{\prime}}{\tau}\right)\tau\right)\approx\sum_{k=1}^{M_{PBF}}\dfrac{-2a_{k}j\left(\omega\dfrac{\tau^{\prime}}{\tau}\right)}{\left(\omega\dfrac{\tau^{\prime}}{\tau}\right)^{2}+2c_{k}j\left(\omega\dfrac{\tau^{\prime}}{\tau}\right)-\omega_{0k}^{2}}$
(A.27)

This directly yields a new PBF expansion for the FTF with a constant time-
delay $\tau^{\prime}$:

$\displaystyle
e^{-j\omega\tau^{\prime}}\approx\sum_{k=1}^{M_{PBF}}\dfrac{-2a_{k}^{\prime}j\omega}{\omega^{2}+2c_{k}^{\prime}j\omega-\omega_{0k}^{\prime
2}}$ (A.28)

where the PBF coefficients are all rescaled in a similar fashion:
$a_{k}^{\prime}=a_{k}\tau/\tau^{\prime}$,
$c_{k}^{\prime}=c_{k}\tau/\tau^{\prime}$, and
$\omega_{0k}^{\prime}=\omega_{0k}\tau/\tau^{\prime}$. Note, however, that the
FTF cutoff frequency and the width of the cutoff region are also rescaled by
this same factor $\tau/\tau^{\prime}$.

####  A.6 Non-reflective boundary at the end of a long duct

Similarly to the case of a cross-section change between two long ducts
(Appendix A.3), the state-space realization of an impedance
$Z(j\omega)=1/(r+j\omega/\omega_{c})$ at the end of a duct $\Omega_{1}$ is
derived by performing a volume averaging of the linearized Euler equations in
a small control volume of thickness $L_{b}$ enclosing the boundary. The
resulting state-space representation writes:

$\displaystyle\dfrac{d}{dt}\underbrace{\begin{pmatrix}\overline{p}(t)\\\
\overline{\varphi}(t)\\\ p^{\Omega_{1}}(t)\\\
\dot{p}^{\Omega_{1}}(t)\end{pmatrix}}_{\mathbf{X}^{(NR)}(t)}=\underbrace{\begin{pmatrix}0&0&0&-\dfrac{c_{0}}{L_{b}\omega_{c}}\\\
-1/\rho_{0}&0&0&0\\\ 0&0&0&1\\\
0&0&-\omega_{0}^{2}&-\omega_{0}\end{pmatrix}}_{\mathbf{A}^{(NR)}}\begin{pmatrix}\overline{p}(t)\\\
\overline{\varphi}(t)\\\ p^{\Omega_{1}}(t)\\\
\dot{p}^{\Omega_{1}}(t)\end{pmatrix}+\underbrace{\begin{pmatrix}\dfrac{\rho_{0}c_{0}^{2}}{L_{b}}&-\dfrac{rc_{0}}{L_{b}}\\\
0&0\\\ 0&0\\\
0&\omega_{0}^{2}\end{pmatrix}}_{\mathbf{B}^{(NR)}}\underbrace{\begin{pmatrix}u_{s}^{\Omega_{1}}(t)\\\
p^{\Omega_{1}}(t)\end{pmatrix}}_{\mathbf{U}^{(NR)}(t)}$ (A.29)

where $\omega_{0}$ is set to a very large value to build a time-derivator in
the last two lines. The non-reflective condition appears as a particular case
of this state-space equation: setting a large value for $\omega_{c}$ reduces
the first line of Eq. (A.29) to the quasi-static equality
$\rho_{0}c_{0}u_{s}^{\Omega_{1}}(t)=rp^{\Omega_{1}}(t)$, and choosing a value
of $r$ very close to unity yields the zero-reflection relation. The state-
space realization is completed by the following output equation:

$\displaystyle\underbrace{\begin{pmatrix}u_{s}(x_{0s},t)\\\
\varphi(x_{0s},t)\end{pmatrix}}_{\mathbf{Y}^{(NR)}(t)}$
$\displaystyle=\underbrace{\begin{pmatrix}0&0&0&0\\\
0&1&0&0\end{pmatrix}}_{\mathbf{C}^{(NR)}}\begin{pmatrix}\overline{p}(t)\\\
\overline{\varphi}(t)\end{pmatrix}$ (A.30)

####  A.7 Complex connection boundary

The surface modal amplitudes $\nu_{k}(t)$ and $\mu_{k}(t)$ (Eq. (4.30))
entirely characterize the acoustics dynamics on the complex boundary
$S_{ci}^{(m)}$, and are therefore used to build its state-space representation
from Eq. (4.31):

$\displaystyle\dfrac{d}{dt}\underbrace{\begin{pmatrix}\mathbf{X}^{S_{ci}^{(m)}}_{1}\\\
\vdots\\\ \mathbf{X}^{S_{ci}^{(m)}}_{K_{S}}\\\
\end{pmatrix}}_{\mathbf{X}^{S_{ci}^{(m)}}(t)}=\underbrace{\left(\begin{array}[]{ccc}\mathbf{A}_{1}^{S_{ci}^{(m)}}&&\\\
&\ddots&\\\ &&\mathbf{A}_{K_{S}}^{S_{ci}^{(m)}}\\\
\end{array}\right)}_{\mathbf{A}^{S_{ci}^{(m)}}}\begin{pmatrix}\mathbf{X}^{S_{ci}^{(m)}}_{1}\\\
\vdots\\\ \mathbf{X}^{S_{ci}^{(m)}}_{K_{S}}\\\
\end{pmatrix}+\underbrace{\left(\begin{array}[]{ccc}\mathbf{B}_{1}^{\Omega_{i,j}}&&\\\
&\ddots&\\\ &&\mathbf{B}_{K_{S}}^{\Omega_{i,j}}\\\
\end{array}\right)}_{\mathbf{B}^{S_{ci}^{(m)}}}\underbrace{\begin{pmatrix}\mathbf{U}_{1}^{\Omega_{i,j}}\\\
\vdots\\\ \mathbf{U}_{K_{S}}^{\Omega_{i,j}}\\\
\end{pmatrix}}_{\mathbf{U}^{S_{ci}^{(m)}}(t)}$ (A.37)

where the blocks $\mathbf{X}^{S_{ci}^{(m)}}_{k}$,
$\mathbf{A}_{k}^{S_{ci}^{(m)}}$, $\mathbf{B}_{k}^{\Omega_{i,j}}$, and
$\mathbf{U}_{k}^{\Omega_{i,j}}$ are defined as:

$\displaystyle\begin{aligned}
&\mathbf{X}^{S_{ci}^{(m)}}_{k}=\begin{pmatrix}\nu_{k}(t)\\\
\dot{\nu}_{k}(t)\\\ \mu_{k}(t)\\\
\mathbf{X}_{\mathbf{K_{R}},k}(t)\end{pmatrix}\ ,\
\mathbf{A}_{k}^{S_{ci}^{(m)}}=\left(\begin{array}[]{cccc}0&1&0&\mathbf{0}\\\
-\omega_{k}^{2}&0&0&\mathbf{0}\\\
0&0&0&\dfrac{1}{\lambda_{k}L_{\mathscr{D}}\overline{\rho_{0}}}\mathbf{C_{K_{R}}}\\\
0&0&0&\mathbf{A_{K_{R}}}\end{array}\right)\\\
&\mathbf{B}_{k}^{\Omega_{i,j}}=\left(\begin{array}[]{cccc}0&0&0&0\\\
-\dfrac{\overline{c_{0}}^{2}}{\lambda_{k}L_{\mathscr{D}}}&0&-\dfrac{\overline{c_{0}}^{2}}{\lambda_{k}L_{\mathscr{D}}}&0\\\
0&\dfrac{1}{\overline{\rho_{0}}\lambda_{k}L_{\mathscr{D}}}&0&\dfrac{1}{\overline{\rho_{0}}\lambda_{k}L_{\mathscr{D}}}\\\
\mathbf{B_{K_{R}}}&0&0&0\\\
\end{array}\right),\mathbf{U}_{k}^{\Omega_{i,j}}=\begin{pmatrix}\mathpzc{s}_{i}\left(u_{s}^{\Omega_{i}}|\mathscr{K}_{k}\right)\\\
\mathpzc{s}_{i}\left(\varphi^{\Omega_{i}}|\mathscr{K}_{k}\right)\\\
\mathpzc{s}_{j}\left(u_{s}^{\Omega_{j}}|\mathscr{K}_{k}\right)\\\
\mathpzc{s}_{j}\left(\varphi^{\Omega_{j}}|\mathscr{K}_{k}\right)\\\
\end{pmatrix}\end{aligned}$ (A.38)

In this equation, the subscripts ${}_{\mathbf{K_{R}}}$ refer to the SISO
state-space realization associated to the inverse Fourier transform of the
complex conductivity $K_{R}(j\omega)$. The blocks
$\mathbf{B}_{k}^{\Omega_{i,j}}$ and $\mathbf{U}_{k}^{\Omega_{i,j}}$ are
employed to compute the source terms from the two adjacent subdomains
$\Omega_{i}$ and $\Omega_{j}$. To be consistent with the subdomains input
vectors, $S_{ci}^{(m)}$ must output the modal amplitudes $\nu_{k}(t)$ and
$\mu_{k}(t)$, which is achieved thanks to:

$\displaystyle\underbrace{\begin{pmatrix}\mathbf{Y}_{1}^{\Omega_{i,j}}\\\
\vdots\\\ \mathbf{Y}_{K_{S}}^{\Omega_{i,j}}\\\
\end{pmatrix}}_{\mathbf{Y}^{S_{ci}^{(m)}}(t)}=\underbrace{\left(\begin{array}[]{ccc}\mathbf{C}_{1}^{\Omega_{i,j}}&&\\\
&\ddots&\\\ &&\mathbf{C}_{K_{S}}^{\Omega_{i,j}}\\\
\end{array}\right)}_{\mathbf{C}^{S_{ci}^{(m)}}}\underbrace{\begin{pmatrix}\mathbf{X}^{S_{ci}^{(m)}}_{1}\\\
\vdots\\\ \mathbf{X}^{S_{ci}^{(m)}}_{K_{S}}\\\
\end{pmatrix}}_{\mathbf{X}^{S_{ci}^{(m)}}(t)}$ (A.42)

with:

$\displaystyle\mathbf{Y}_{k}^{\Omega_{i,j}}=\begin{pmatrix}\mu_{k}(t)\\\
\nu_{k}(t)\\\ \mu_{k}(t)\\\ \nu_{k}(t)\\\ \end{pmatrix}\ ,\
\mathbf{C}_{k}^{\Omega_{i,j}}=\left(\begin{array}[]{cccc}0&0&1&\mathbf{0}\\\
1&0&0&\mathbf{0}\\\ 0&0&1&\mathbf{0}\\\ 1&0&0&\mathbf{0}\\\
\end{array}\right)$ (A.47)

####  A.8 Howe’s model Rayleigh conductivity

For a Rayleigh conductivity $K_{R}$ expressed through the approximation to
Howe’s model given in Eq. (4.37), the inverse Fourier transform of the term
$-j\omega\rho_{u}d^{2}/K_{R}(j\omega)$ in Eq. (4.25) can be translated into a
SISO state-space realization that is directly embedded into the complex
boundary state-space realization provided in Eq. (A.37)-(A.38). Its expression
writes:

$\displaystyle\derivative{t}\underbrace{\left(z(t)\right)}_{\mathbf{X_{K_{R}}}(t)}=\underbrace{\left(-K_{R}^{A}/K_{R}^{B}\right)}_{\mathbf{A_{K_{R}}}}\left(z(t)\right)+\underbrace{\left(-1/K_{R}^{B}\right)}_{\mathbf{B_{K_{R}}}}\left(\mathpzc{s}_{i}\left(u_{s}^{\Omega_{i}}|\mathscr{K}_{k}\right)\right)$
(A.48)

and the output $Y_{K_{R}}$ is computed through:

$\displaystyle\left(Y_{K_{R}}(t)\right)=\underbrace{\left(-\rho_{u}d^{2}\right)}_{\mathbf{C_{K_{R}}}}\left(z(t)\right)$
(A.49)

###  Appendix B Reference solution for a liner in a cylindrical domain

Let us consider an annular and a cylindrical subdomains separated by a ribbon-
like multi-perforated liner characterized by the Rayleigh conductivity
$K_{R}(j\omega)$, such as the system presented in Sec. 4.3. In both the inner
cylinder and the outer annulus, the pressure is solution of the classical
homogeneous Helmholtz equation. In these two domains, the separation of
variables $p(r,\theta)=R(r)\Theta(\theta)$ is performed, which yields a set of
Bessel differential equations governing $R(r)$. Solutions of these equations
write:

$\displaystyle R(r\leq R_{1})=AJ_{n}(kr)\ ,\ R(R_{1}<r\leq
R_{2})=BJ_{n}(kr)+CY_{n}(kr)$ (B.1)

where $A,B,C$, and $k=\omega/c_{0}$ are complex constants to be determined,
and $J_{n}$ and $Y_{n}$ are the Bessel functions of the 1st and 2nd kind,
respectively. The azimuthal component $\Theta(\theta)$ is a simple linear
combination of $\cos(n\theta)$ and $\sin(n\theta)$. The rigid-wall conditions
at $r=R_{2}$ combined with the flux continuity at $r=R_{1}$ as well as
pressure jump through the multi-perforated liner leads to a linear system
$\mathbf{M}\mathbf{X}=\mathbf{0}$ where $\mathbf{X}={}^{t}(A\ B\ C)$ and
$\mathbf{M}$ is given by:

$\displaystyle\left(\begin{array}[]{ccc}0&J_{n}^{\prime}(kR_{2})&Y_{n}^{\prime}(kR_{2})\\\
\dfrac{d^{2}}{K_{R}(j\omega)}kJ_{n}^{\prime}(kR_{1})+J_{n}(kR_{1})&-J_{n}(kR_{1})&-Y_{n}^{\prime}(kR_{1})\\\
J_{n}(kR_{1})&-J_{n}(kR_{1})+\dfrac{d^{2}}{K_{R}(j\omega)}kJ_{n}^{\prime}(kR_{1})&-Y_{n}(kR_{1})+\dfrac{d^{2}}{K_{R}(j\omega)}kY_{n}^{\prime}(kR_{1})\\\
\end{array}\right)$ (B.5)

This linear system admits non-trivial solutions if and only if
$det(\mathbf{M})=0$, which results in a dispersion relation $F(\omega)=0$ in
the complex plane. The resolution of this equation yields the
eigenfrequencies, their associated growth rates, and ultimately the mode
shapes thanks to Eq. (B.1).

### References

  * [1] M. Brebion, Joint numerical and experimental study of thermoacoustic instabilities.  PhD thesis, Institut National Polytechnique de Toulouse, 2017.
  * [2] D. Dunn-Rankin, Lean combustion: technology and control.  Academic Press, 2011.
  * [3] D. G. Crighton, A. P. Dowling, J. E. F. Williams, M. Heckl, and F. Leppington, Modern methods in analytical acoustics.  Lecture Notes, New-York: Springer Verlag, 1992.
  * [4] S. Candel, “Combustion instabilities coupled by pressure waves and their active control.,” in 24th Symp. (Int.) on Combustion, pp. 1277–1296, The Combustion Institute, Pittsburgh, 1992.
  * [5] T. Poinsot and D. Veynante, Theoretical and Numerical Combustion.  R.T. Edwards, 2nd edition, 2005.
  * [6] T. Lieuwen and V. Yang, Combustion Instabilities in Gas Turbine Engines. Operational Experience, Fundamental Mechanisms and Modeling.  Prog. in Astronautics and Aeronautics AIAA Vol 210, 2005.
  * [7] T. Poinsot, “Prediction and control of combustion instabilities in real engines,”  Proc. Combust. Inst. , pp. 1–28, 2017.
  * [8] J. L. Rayleigh, “The explanation of certain acoustical phenomena,” Nature, vol. 18, pp. 319–321, 1878.
  * [9] P. L. Rijke, “Notiz uber eine neue art, die in einer an beiden enden offenen rohre enthaltene luft in schwingungen zu ersetzen,” Annalen der Physik, vol. 107, p. 339, 1859.
  * [10] T. J. Poinsot, A. C. Trouve, D. P. Veynante, S. M. Candel, and E. J. Esposito, “Vortex-driven acoustically coupled combustion instabilities,” Journal of fluid mechanics, vol. 177, pp. 265–292, 1987.
  * [11] L. Crocco and S. Cheng, “Theory of combustion instability in liquid propellant rocket motors,(1956).”
  * [12] D. T. Harrje, “Liquid propellant rocket combustion instability,” NASA-SP-194, 1972.
  * [13] J. C. Oefelein and V. Yang, “Comprehensive Review of Liquid-Propellant Combustion Instabilities in F-1 Engines,” Journal of Propulsion and Power, vol. 9, no. 5, pp. 657–677, 1993.
  * [14] F. E. C. Culick and P. Kuentzmann, Unsteady Motions in Combustion Chambers for Propulsion Systems.  NATO Research and Technology Organization, 2006.
  * [15] T. Lieuwen, “Modeling premixed combustion-acoustic wave interactions: A review,”  J. Prop. Power , vol. 19, no. 5, pp. 765–781, 2003.
  * [16] T. Lieuwen and A. Banaszuk, “Background noise effects on combustor stability,”  J. Prop. Power , vol. 21, no. 1, pp. 25–31, 2005.
  * [17] S. Ducruix, T. Schuller, D. Durox, and S. Candel, “Combustion dynamics and instabilities: elementary coupling and driving mechanisms,”  J. Prop. Power , vol. 19, no. 5, pp. 722–734, 2003.
  * [18] S. Ducruix, T. Schuller, D. Durox, and S. Candel, “Combustion instability mechanisms in premixed combustors,”  Prog. Energy Combust. Sci. , vol. 210, p. 179, 2005.
  * [19] J. O’Connor, V. Acharya, and T. Lieuwen, “Transverse combustion instabilities: Acoustic, fluid mechanic, and flame processes,”  Prog. Energy Combust. Sci. , vol. 49, pp. 1–39, Aug. 2015.
  * [20] F. E. Marble and S. Candel, “Acoustic disturbances from gas nonuniformities convected through a nozzle,”  J. Sound Vib. , vol. 55, pp. 225–243, 1977.
  * [21] G. E. Abouseif, J. A. Keklak, and T. Y. Toong, “Ramjet rumble: The low-frequency instability mechanism in coaxial dump combustors,”  Combust. Sci. Technol. , vol. 36, no. 1-2, pp. 83–108, 1984.
  * [22] J. J. Keller, W. Egli, and J. Hellat, “Thermally induced low-frequency oscillations,” Zeitschrift für angewandte Mathematik und Physik ZAMP, vol. 36, pp. 250–274, 1985.
  * [23] C. S. Goh and A. S. Morgans, “The influence of entropy waves on the thermoacoustic stability of a model combustor,”  Combust. Sci. Technol. , vol. 185, no. 2, pp. 249–268, 2013.
  * [24] E. Motheau, F. Nicoud, and T. Poinsot, “Mixed acoustic - entropy combustion instabilities in gas turbines,”  J. Fluid Mech. , vol. 749, pp. 542–576, 6 2014.
  * [25] F. E. C. Culick, “Some recent results for nonlinear acoustic in combustion chambers,”  AIAA J. , vol. 32, no. 1, pp. 146–169, 1994.
  * [26] L. Crocco, “Aspects of combustion instability in liquid propellant rocket motors. part II.,”  J. American Rocket Society , vol. 22, pp. 7–16, 1952.
  * [27] N. Noiray and B. Schuermans, “On the dynamic nature of azimuthal thermoacoustic modes in annular gas turbine combustion chambers,” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 469, p. 20120535, Dec. 2012.
  * [28] K. Kedia, H. Altay, and A. Ghoniem, “Impact of flame-wall interaction on premixed flame dynamics and transfer function characteristics,”  Proc. Combust. Inst. , vol. 33, pp. 1113–1120, 2011.
  * [29] P. Palies, T. Schuller, D. Durox, and S. Candel, “Modeling of premixed swirling flames transfer functions,”  Proc. Combust. Inst. , vol. 33, no. 2, pp. 2967–2974, 2011.
  * [30] C. F. Silva, T. Emmert, S. Jaensch, and W. Polifke, “Numerical study on intrinsic thermoacoustic instability of a laminar premixed flame,” Combustion and Flame, vol. 162, no. 9, pp. 3370–3378, 2015.
  * [31] T. Kaiser, G. Öztarlik, L. Selle, and T. Poinsot, “Impact of symmetry breaking on the flame transfer function of a laminar premixed flame,” Proceedings of the Combustion Institute, vol. 37, no. 2, pp. 1953–1960, 2019\.
  * [32] H. Pitsch, “Large eddy simulation of turbulent combustion,”  Ann. Rev. Fluid Mech. , vol. 38, pp. 453–482, 2006.
  * [33] P. Wolf, G. Staffelbach, L. Gicquel, J.-D. Muller, and T. Poinsot, “Acoustic and large eddy simulation studies of azimuthal modes in annular combustion chambers,”  Combust. Flame, vol. 159, pp. 3398–3413, Nov. 2012.
  * [34] A. Urbano, Q. Douasbin, L. Selle, G. Staffelbach, B. Cuenot, T. Schmitt, S. Ducruix, and S. Candel, “Study of flame response to transverse acoustic modes from the LES of a 42-injector rocket engine,”  Proc. Combust. Inst. , vol. 36, no. 2, pp. 2633–2639, 2017.
  * [35] L. Y. M. Gicquel, G. Staffelbach, and T. Poinsot, “Large eddy simulations of gaseous flames in gas turbine combustion chambers,”  Prog. Energy Combust. Sci. , vol. 38, no. 6, pp. 782 – 817, 2012.
  * [36] L. Chong, T. Komarek, R. Kaess, S. Foller, and W. Polifke, “Identification of flame transfer functions from les of a premixed swirl burner,” in ASME Turbo expo (A. P. GT2010-22769, ed.), (Glasgow, UK), June 2010.
  * [37] M. Merk, R. Gaudron, C. Silva, M. Gatti, C. Mirat, T. Schuller, and W. Polifke, “Prediction of combustion noise of an enclosed flame by simultaneous identification of noise source and flame dynamics,” Proceedings of the Combustion Institute, vol. 37, no. 4, pp. 5263–5270, 2019.
  * [38] T. Schuller, D. Durox, and S. Candel, “A unified model for the prediction of laminar flame transfer functions: comparisons between conical and v-flames dynamics,”  Combust. Flame, vol. 134, pp. 21–34, 2003.
  * [39] M. Bauerheim, A. Ndiaye, P. Constantine, S. MOREAU, and F. Nicoud, “Symmetry breaking of azimuthal thermoacoustic modes: the UQ perspective,”  J. Fluid Mech. , vol. 789, pp. 534–566, Jan. 2016.
  * [40] L. Magri, “Adjoint methods as design tools in thermoacoustics,” Applied Mechanics Reviews, vol. 71, no. 2, 2019.
  * [41] L. Vingert, M. Habiballah, P. Vuillermoz, and S. Zurbach, “Mascotte, a test facility for cryogenic combustion research at high pressure,” 51st International Astronautical Congress, Rio de Janeiro, Brazil, 2000, 2000.
  * [42] A. Urbano, L. Selle, G. Staffelbach, B. Cuenot, T. Schmitt, S. Ducruix, and S. Candel, “Exploration of combustion instability triggering using Large Eddy Simulation of a multiple injector liquid rocket engine,”  Combust. Flame, vol. 169, pp. 1–12, May 2016.
  * [43] M. T. Migliorino and C. Scalo, “Real-fluid effects on standing-wave thermoacoustic instability,” Journal of Fluid Mechanics, vol. 883, 2020\.
  * [44] F. E. C. Culick, “Non linear behavior of acoustic waves in combustion chambers,”  Acta Astronautica , vol. 3, pp. 714–757, 1976.
  * [45] V. Yang and F. CULICK, “Third-order nonlinear acoustic waves and triggering of pressure oscillations in combustion chambers. i-longitudinal modes,” in 23rd Joint Propulsion Conference, p. 1873, 1987.
  * [46] F. Nicoud and K. Wieczorek, “About the zero mach number assumption in the calculation of thermoacoustic instabilitie,” Int. J. Spray and Combustion Dynamic, vol. 1, pp. 67–112, 2009.
  * [47] F. Nicoud, L. Benoit, C. Sensiau, and T. Poinsot, “Acoustic modes in combustors with complex impedances and multidimensional active flames,”  AIAA J. , vol. 45, pp. 426–441, 2007.
  * [48] E. Motheau, L. Selle, and F. Nicoud, “Accounting for convective effects in zero-Mach-number thermoacoustic models,”  J. Sound Vib. , vol. 333, no. 1, pp. 246–262, 2013.
  * [49] L. S. Chen, S. Bomberg, and W. Polifke, “Propagation and generation of acoustic and entropy waves across a moving flame front,” Combustion and Flame, vol. 166, pp. 170–180, 2016.
  * [50] M. Howe, “On the theory of unsteady high reynolds number flow through a circular aperture,” Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences, vol. 366, no. 1725, pp. 205–223, 1979.
  * [51] D. Yang and A. S. Morgans, “The acoustics of short circular holes opening to confined and unconfined spaces,” Journal of Sound and Vibration, vol. 393, pp. 41–61, 2017.
  * [52] M. Bauerheim, A. Ndiaye, P. Constantine, S. Moreau, and F. Nicoud, “Symmetry breaking of azimuthal thermoacoustic modes: The UQ perspective,” Journal of Fluid Mechanics, vol. 789, pp. 534–566, 2016.
  * [53] A. Avdonin, S. Jaensch, C. F. Silva, M. Češnovar, and W. Polifke, “Uncertainty quantification and sensitivity analysis of thermoacoustic stability with non-intrusive polynomial chaos expansion,” Combustion and Flame, vol. 189, pp. 300–310, 2018.
  * [54] J. G. Aguilar and M. P. Juniper, “Adjoint methods for elimination of thermoacoustic oscillations in a model annular combustor via small geometry modifications,” in ASME Turbo Expo 2018: Turbomachinery Technical Conference and Exposition, American Society of Mechanical Engineers, 2018.
  * [55] L. Crocco and S. I. Cheng, Theory of combustion instability in liquid propellant rocket motors, vol. Agardograph No 8.  Butterworths Science, 1956.
  * [56] A. Cuquel, D. Durox, and T. Schuller, “Impact of flame base dynamics on the nonlinear frequency response of conical flames,”  C. R. Acad. Sci. Mecanique, vol. 341, pp. 171–180, 2013.
  * [57] M. Fleifil, A. M. Annaswamy, Z. A. Ghoneim, and A. F. Ghoniem, “Response of a laminar premixed flame to flow oscillations: A kinematic model and thermoacoustic instability results,”  Combust. Flame, vol. 106, no. 4, pp. 487–510, 1996.
  * [58] N. Magina, D. H. Shin, V. Acharya, and T. Lieuwen, “Response of non-premixed flames to bulk flow perturbations,” Proceedings of the Combustion Institute, vol. 34, no. 1, pp. 963–971, 2013.
  * [59] N. A. Magina and T. C. Lieuwen, “Effect of axial diffusion on the response of diffusion flames to axial flow perturbations,” Combustion and Flame, vol. 167, pp. 395–408, 2016.
  * [60] S. Candel, D. Durox, T. Schuller, J.-F. Bourgouin, and J. P. Moeck, “Dynamics of Swirling Flames,”  Ann. Rev. Fluid Mech. , vol. 46, pp. 147–173, Jan. 2014.
  * [61] D. Mejia, M. Miguel-Brebion, A. Ghani, T. Kaiser, F. Duchaine, L. Selle, and T. Poinsot, “Influence of flame-holder temperature on the acoustic flame transfer functions of a laminar flame,”  Combust. Flame, vol. 188, pp. 5 – 12, 2018\.
  * [62] M. Gatti, R. Gaudron, C. Mirat, L. Zimmer, and T. Schuller, “Impact of swirl and bluff-body on the transfer function of premixed flames,” Proceedings of the Combustion Institute, vol. 37, no. 4, pp. 5197–5204, 2019\.
  * [63] S. Hermeth, G. Staffelbach, L. Y. Gicquel, V. Anisimov, C. Cirigliano, and T. Poinsot, “Bistable swirled flames and influence on flame transfer functions,”  Combust. Flame, vol. 161, no. 1, pp. 184 – 196, 2014.
  * [64] N. Noiray, D. Durox, T. Schuller, and S. Candel, “A unified framework for nonlinear combustion instability analysis based on the flame describing function,”  J. Fluid Mech. , vol. 615, pp. 139–167, 2008.
  * [65] Preetham, H. Santosh, and T. Lieuwen, “Dynamics of laminar premixed flames forced by harmonic velocity disturbances,” Journal of Propulsion and Power, vol. 24, no. 6, pp. 1390–1402, 2008.
  * [66] G. Ghirardo, B. Ćosić, M. P. Juniper, and J. P. Moeck, “State-space realization of a describing function,” Nonlinear Dynamics, vol. 82, no. 1-2, pp. 9–28, 2015.
  * [67] A. Orchini and M. P. Juniper, “Flame Double Input Describing Function analysis,” Combustion and Flame, 2016.
  * [68] M. Haeringer, M. Merk, and W. Polifke, “Inclusion of higher harmonics in the flame describing function for predicting limit cycles of self-excited combustion instabilities,” Proceedings of the Combustion Institute, vol. 37, no. 4, pp. 5255–5262, 2019.
  * [69] K. Kashinath, I. C. Waugh, and M. P. Juniper, “Nonlinear self-excited thermoacoustic oscillations of a ducted premixed flame : bifurcations and routes to chaos,” Journal of Fluid Mechanics, pp. 1–26, 2014.
  * [70] A. Orchini and M. P. Juniper, “Heat release response to forced flow oscillations of a low-order modelled laboratory scale dump combustor,” in ASME turbo expo 2015: turbine technical conference and exposition, American Society of Mechanical Engineers Digital Collection, 2015.
  * [71] a. Orchini, S. J. Illingworth, and M. P. Juniper, “Frequency domain and time domain analysis of thermoacoustic oscillations with wave-based acoustics,” Journal of Fluid Mechanics, vol. 775, pp. 387–414, 2015.
  * [72] Z. T. Ren, S. H. Li, and M. Zhu, “Dynamic responses of anchored ducted flames to harmonic velocity forcing,” International Journal of Spray and Combustion Dynamics, vol. 10, no. 1, pp. 72–85, 2018.
  * [73] B. Semlitsch, A. Orchini, A. P. Dowling, and M. P. Juniper, “G-equation modelling of thermoacoustic oscillations of partially premixed flames,” International Journal of Spray and Combustion Dynamics, vol. 9, no. 4, pp. 260–276, 2017.
  * [74] K. Balasubramanian and R. I. Sujith, “Thermoacoustic instability in a Rijke tube: Non-normality and nonlinearity,”  Phys. Fluids , vol. 20, no. 4, p. 044103, 2008\.
  * [75] S. J. Illingworth, I. C. Waugh, and M. P. Juniper, “Finding thermoacoustic limit cycles for a ducted Burke-Schumann flame,” Proceedings of the Combustion Institute, vol. 34, no. 1, pp. 911–920, 2013.
  * [76] L. Magri and M. P. Juniper, “Global modes, receptivity, and sensitivity analysis of diffusion flames coupled with duct acoustics,” Journal of Fluid Mechanics, 2014.
  * [77] A. Avdonin, M. Meindl, and W. Polifke, “Thermoacoustic analysis of a laminar premixed flame using a linearized reactive flow solver,” Proceedings of the Combustion Institute, vol. 37, no. 4, pp. 5307–5314, 2019.
  * [78] H. Yu, T. Jaravel, M. Ihme, M. P. Juniper, and L. Magri, “Data assimilation and optimal calibration in nonlinear models of flame dynamics,” Journal of Engineering for Gas Turbines and Power, vol. 141, no. 12, 2019.
  * [79] E. Gullaud and F. Nicoud, “Effect of perforated plates on the acoustics of annular combustors,” AIAA journal, vol. 50, no. 12, pp. 2629–2642, 2012\.
  * [80] T. Sayadi, V. Le Chenadec, P. J. Schmid, F. Richecoeur, and M. Massot, “Thermoacoustic instability - A dynamical system and time domain analysis,” Journal of Fluid Mechanics, vol. 753, pp. 448–471, 2014.
  * [81] M. R. Bothien, J. P. Moeck, A. Lacarelle, and C. O. Paschereit, “Time domain modelling and stability analysis of complex thermoacoustic systems,” Proceedings of the Institution of Mechanical Engineers, Part A: Journal of Power and Energy, vol. 221, no. 5, pp. 657–668, 2007.
  * [82] B. Schuermans, V. Bellucci, and C. O. Paschereit, “Thermoacoustic modeling and control of multi burner combustion systems,” in ASME Turbo Expo 2003, collocated with the 2003 International Joint Power Generation Conference, pp. 509–519, American Society of Mechanical Engineers, 2003.
  * [83] B. Schuermans, Modeling and control of thermoacoustic instabilities.  PhD thesis, EPFL, 2003.
  * [84] T. Emmert, M. Meindl, S. Jaensch, and W. Polifke, “Linear state space interconnect modeling of acoustic systems,” Acta Acustica united with Acustica, vol. 102, no. 5, pp. 824–833, 2016.
  * [85] B. Friedland, Control system design: an introduction to state-space methods.  Courier Corporation, 2012.
  * [86] R. Redheffer, “On a certain linear fractional transformation,” Studies in Applied Mathematics, vol. 39, no. 1-4, pp. 269–286, 1960.
  * [87] F. Nicoud, L. Benoit, C. Sensiau, and T. Poinsot, “Acoustic modes in combustors with complex impedances and multidimensional active flames,” AIAA journal, vol. 45, no. 2, pp. 426–441, 2007.
  * [88] S. Camporeale, B. Fortunato, and G. Campa, “A finite element method for three-dimensional analysis of thermo-acoustic combustion instability,” Journal of Engineering for Gas Turbines and Power, vol. 133, no. 1, p. 011506, 2011.
  * [89] G. A. Mensah, G. Campa, and J. P. Moeck, “Efficient computation of thermoacoustic modes in industrial annular combustion chambers based on bloch-wave theory,” Journal of Engineering for Gas Turbines and Power, vol. 138, no. 8, p. 081502, 2016.
  * [90] W. Krebs, G. Walz, P. Flohr, and S. Hoffmann, “Modal analysis of annular combustors: Effect of burner impedance,” in ASME Turbo Expo 2001: Power for Land, Sea, and Air, American Society of Mechanical Engineers Digital Collection, 2001.
  * [91] A. Andreini, C. Bianchini, B. Facchini, F. Simonetti, and A. Peschiulli, “Assessment of numerical tools for the evaluation of the acoustic impedance of multi-perforated plates,” in ASME 2011 Turbo Expo: Turbine Technical Conference and Exposition, pp. 1065–1077, American Society of Mechanical Engineers Digital Collection, 2011.
  * [92] A. Giusti, A. Andreini, B. Facchini, I. Vitale, and F. Turrini, “Thermoacoustic analysis of a full annular aero-engine lean combustor with multi-perforated liners,” in 19th AIAA/CEAS Aeroacoustics Conference, p. 2099, 2013.