# A quasi-one-dimensional critical-throat acoustic boundary condition for thermally choked dual-mode ramjet nozzles

Frédéric Olivon, Aurelien Genot, Jean-Étienne Durand, Avraham Hirschberg, Estelle Piot

► To cite this version:

Frédéric Olivon, Aurelien Genot, Jean-Étienne Durand, Avraham Hirschberg, Estelle Piot. A quasi-one-dimensional critical-throat acoustic boundary condition for thermally choked dual-mode ramjet nozzles. 2025. hal-05197093

HAL Id: hal-05197093

https://hal.science/hal-05197093v1

Preprint submitted on 2 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution - ShareAlike 4.0 International License
A quasi-one-dimensional critical-throat acoustic boundary condition
for thermally choked dual-mode ramjet nozzles

Frédéric Olivon¹,²), Aurelien Genot²), Jean-Étienne Durand¹), Avraham Hirschberg³)
and Estelle Piot²)

1) DMPE, ONERA, Université Paris-Saclay, Palaiseau, 91120, France.

2) DMPE, ONERA, Université de Toulouse, Toulouse, 31000, France.

3) Dept. of Applied Physics, Technische Universiteit Eindhoven, Eindhoven,
5600 MB, The Netherlands.

Summary

Thermally choked nozzles, where choking is induced
by heat addition rather than a geometrical throat,
are a promising solution for dual-mode ramjet transi-
tions to hypersonic speed. Despite their relevance,
thermal throat boundary conditions for quasi-one-
dimensional acoustic modelling are not derived in the
state-of-the-art. This study introduces a generalized
critical-throat acoustic boundary condition applica-
ble to both geometrically and thermally choked flow
configurations. A dedicated one-dimensional linear
acoustic solver is formulated to incorporate this con-
dition and is validated against two-dimensional Euler
simulations. Particular attention is paid to the impact
of entropy and acoustic waves at the critical throat.
The results show that the new boundary condition
improves the prediction of acoustic reflection, entropy
noise production, and transmission coefficients, espe-
cially under thermally choked conditions where the
commonly used quasi-steady assumption fails. For
both the thermal-throat and geometric-throat con-
figurations, the deviation in the acoustic transmis-
sion coefficient between the linear acoustic model us-
ing the proposed boundary condition and the simu-
lations remains below 2%, while the deviation in the
entropy-noise transmission coefficient remains below
5%, demonstrating the robustness of the proposed
boundary condition.

1 Introduction

Ramjets and scramjets are key propulsion systems for
high-speed flight due to their higher propulsive effi-
ciency compared to turbojets and turbofans at equiv-
alent velocities. Ramjets operate with subsonic com-
bustion and are effective up to moderate flight Mach
numbers (Mach 5), but their performance deteriorates
at higher speeds due to substantial pressure losses and
the effects of chemical dissociation and recombina-
tion [1]. In contrast, scramjets use supersonic com-

<figure><img src="image_2.png" /><figcaption>Figure 1: (Colour online) Sketch of an idealized dual-mode ramjet incorporating a thermal throat.</figcaption></figure>

bustion, achieving greater efficiency at Mach num-
bers above 5, making them suitable for hypersonic
air breathing applications [2]. To extend operability
across a wider range of Mach numbers in flight, dual-
mode ramjets have been developed. These engines
can seamlessly shift between the subsonic and super-
sonic combustion regimes, offering adaptability in a
wide range of flight conditions [3, 4]. However, this
versatility introduces significant design complexity.

In conventional ramjets, a geometric throat is used
to choke the flow, hereafter referred to as a geomet-
rically choked-flow nozzle. Although effective in low-
flight Mach number combustion regimes, this configu-
ration becomes unsuitable under scramjet conditions.
An alternative is the use of a thermal throat (fig. 1),
produced by the volumetric heat release from com-
bustion [5, 6]. Unlike the geometric configuration,
this design uses a purely divergent nozzle, enabling
a continuous transition between combustion regimes
without altering the geometry of the engine.

Quasi-one-dimensional linear acoustic models are powerful tools for analysing and designing the acoustic behaviour of combustion chambers due to their numerical efficiency. However, they involve a one-dimensional acoustic boundary condition at the critical-throat location. For example, Lin et al. [7] proposed a model that describes the acoustic feedback loop in a scramjet combustor featuring a thermal throat, but treated the flame and critical throat as a single acoustic boundary, thus neglecting their distinct physical roles.

In quasi-one-dimensional acoustic models, it is com-
monly assumed that Mach number fluctuations vanish
at the critical throat [8–13]. While this assumption
holds under quasi-steady, isentropic flow conditions
for geometrically choked nozzles, Olivon et al. [14, 15]
demonstrated its failure under unsteady inlet condi-
tions for both geometrical and thermal throats. In
particular, this quasi-steady critical Mach number as-
sumption induces a singularity at the critical-throat
location [15]. Building on the work of Marble and
Candel [8], a generalized boundary condition has been
proposed for isentropic flows [15]. The key idea be-
hind this theory is to impose a regular solution for the
acoustic field at the throat. This was originally pro-
posed by Tsien [16], following a suggestion by Crocco,
and later used by Marble and Candel [8] for a particu-
lar nozzle geometry. Considering a series expansion of
the solution in powers of the distance from the crit-
ical throat, Tsien [16] discarded the solution whose
leading-order term is inversely proportional to that
distance. A generalization of this approach for ther-
mally choked configurations is proposed in the present
study.

The effect of heat addition on flow dynamics has
been extensively studied in subsonic flows, particu-
larly in the context of acoustic wave propagation.
These studies show that steady heating can induce
acoustic absorption [17, 18] and generate entropy
waves as acoustic waves propagated through the
heated region [19–21]. Beyond the subsonic regime,
investigations have focused on geometrically choked
flows, providing important insights into aero-acoustic
responses and wave dynamics [22–25]. However, none
of these studies have addressed heat addition as the
primary choking mechanism.

All the aforementioned works assume a calorically
perfect gas, in which the specific heat capacities are
considered temperature-independent and the flow is
assumed to be isentropic. However, in thermally-
choked flows, characterized by significant temperature
variations, theses assumptions are inherently ques-
tionable [26–28]. Olivon et al. [14] numerically in-
vestigated the impact of temperature-dependent spe-
cific heat capacities and demonstrated that both the
acoustic-reflection and entropy-noise-production coef-
ficients can vary by as much as 20%. These findings
indicate that the heat capacity temperature depen-
dency must be taken into account to model thermally-
choked flows using linear acoustic models.

In the present work, the inviscid assumption is used
which is common in the literature on geometrical-
throat nozzles [8, 10, 13], because boundary layers
are thin. Although these layers modify the mean cen-
tral flow, their influence on thermal throat dynamics
has been shown to be limited [29]. The volumetric
heat source term P<sub>v</sub> is assumed to be steady, and
fluctuations in heat release are not considered. This
implies no feedback loop between the heat release and
flow field, as adopted in previous studies [17, 22, 30].
Yeddula et al. [31] highlighted this model limitation;

nevertheless, it allows for a first step before introduc-
ing more complex phenomena. Finally, the analysis is
limited to acoustic and entropy waves; indirect noise
arising from gas mixture non-uniformity, as investi-
gated by Gentil et al. [13, 32], is not taken into ac-
count.

The objective of this study is to develop and val-
idate a boundary condition that describes the be-
haviour of the critical throat subjected to both acous-
tic and entropy wave forcing in a thermally choked
flow configuration. The quasi-one-dimensional lin-
ear acoustic model employed is validated against two-
dimensional Euler numerical simulations.

The paper is organized as follows. The quasi-one-
dimensional linear acoustic model is first presented in
section 2, along with the new critical throat boundary
condition model in section 3. The methodology and
numerical parameters used for the validation of the
proposed boundary condition are then introduced in
section 4. Results and discussion follow in section 5.
Conclusions are summarized in section 6.

## 2 Quasi-one-dimensional linear acoustic model

The following assumptions define the limits of va-
lidity of the present model. First, the flow is mod-
elled as quasi-one-dimensional, and only plane waves
are considered. Second, the Euler equations are em-
ployed, neglecting viscous and thermal diffusion ef-
fects. Third, a non-uniform, steady volumetric heat
source term P<sub>v</sub> is introduced to model the energy re-
lease from the flame. The influence of the acoustic
perturbations on this heat release are neglected.

2.1 Flow model

For an ideal gas, assuming no losses due to viscosity,
thermal diffusion, or external force fields, the stan-
dard quasi-one-dimensional equations for mass, mo-
mentum, and energy conservation are expressed as
follows [8]:

$$
\left\{
\begin{array}{l}
\displaystyle \frac{\mathrm{D}\rho}{\mathrm{D}t} + \rho \frac{\partial u}{\partial x} + \rho u \alpha = 0, \\[1em]
\displaystyle \frac{\mathrm{Du}}{\mathrm{D}t} + \frac{1}{\rho} \frac{\partial p}{\partial x} = 0, \\[1em]
\displaystyle \frac{\mathrm{D}p}{\mathrm{D}t} - c^2 \frac{\mathrm{D}\rho}{\mathrm{D}t} = (\gamma - 1) P_v,
\end{array}
\right.
\qquad (1)
$$

where D/Dt = ∂/∂t + u∂/∂x is the material deriva-
tive, ρ is the density, u is the longitudinal velocity,
α = (dA/dx)/A is the relative variation of the noz-
zle cross-sectional area A, p is the static pressure, c
is the speed of sound, γ = c<sub>p</sub>/c<sub>v</sub> is the ratio of spe-
cific heats at constant pressure c<sub>p</sub> and constant vol-
ume c<sub>v</sub>, and P<sub>v</sub> is the volumetric heat-rate addition
source term. For clarity and conciseness, the system of equations (1) will hereafter be referred to as the quasi-one-dimensional Euler equations.

To close the quasi-one-dimensional Euler equations (1), the ideal gas law is used:

$$p = \rho r T, \qquad (2)$$

where T is the static temperature, $r = \mathcal{R}/\mathcal{W}$ is the specific gas constant, with $\mathcal{R} = 8.3145 \text{ J mol}^{-1} \text{ K}^{-1}$ the universal gas constant, and $\mathcal{W}$ the molar mass of the gas.

The gas considered here is the air, described as an ideal gas with a specific isobaric heat capacity $c_p$, function of the temperature. It is computed using a seventh-order polynomial function of the static temperature T:

$$c_p = \sum_{n=0}^{7} \frac{\alpha_n}{\mathcal{W}} T^n, \qquad (3)$$

where $\alpha_n$ are the heat capacity coefficients provided by Lemmon et al. [33].

For a uniform gas composition, the Gibbs energy equation provides a relation between the variations of entropy $s$, temperature $T$ and pressure $p$ such as:

$$\frac{ds}{c_p} = \frac{dT}{T} - (\gamma - 1) \frac{dp}{\gamma p}, \qquad (4)$$

for a reversible process.

## 2.2 Linear acoustic propagation model

To study the propagation of acoustic and entropy waves through an heated nozzle with non-uniform steady flow field, quasi-one-dimensional Euler equations (1) are linearized using the perturbation convention $y(x,t) = \bar{y}(x) + y'(x,t)$, where $\bar{y}$ denotes the steady component and $y'$ the small perturbation of $y = [\rho, u, p, T, c, \gamma]$. The volumetric heat-rate source term $P_v$ being prescribed as steady, there are no fluctuations ($P'_v = 0$). Neglecting second-order perturbations and introducing the normalized quantities $D = \rho'/\bar{\rho}$, $U = u'/\bar{u}$, $P = p'/(\gamma\bar{p})$ yields the following linearized system of differential equations:

$$\left\{
\begin{aligned}
& \frac{\partial D}{\partial t} + \bar{u} \frac{\partial}{\partial x} (U + D) = 0, \\
& \frac{\partial U}{\partial t} + \bar{u} \frac{\partial U}{\partial x} + \frac{\bar{c}^2}{\bar{u}} \frac{\partial P}{\partial x} \\
& \qquad + 2 \frac{d\bar{u}}{dx} U + \left( \frac{\kappa \bar{c}^2}{\bar{u}} \frac{d\bar{T}}{dx} - \bar{\gamma} \frac{d\bar{u}}{dx} \right) P + \frac{d\bar{u}}{dx} D = 0, \\
& \frac{\partial P}{\partial t} + \bar{u} \frac{\partial}{\partial x} (U + P) + (\bar{\gamma} - 1) \frac{\bar{P}_v}{\bar{\gamma}\bar{p}} U \\
& \qquad + \left( \frac{\bar{P}_v}{\bar{p}} \left( (\bar{\gamma} - 1) - \kappa\bar{\gamma}\bar{T} \right) + \kappa\bar{u} \left( \frac{d\bar{T}}{dx} - \frac{\bar{\gamma}\bar{T}}{\bar{p}} \frac{d\bar{\rho}}{dx} \right) \right) P \\
& \qquad + \kappa\bar{T} \left( \frac{\bar{P}_v}{\bar{p}} + \frac{\bar{u}}{\bar{p}} \frac{d\bar{\rho}}{dx} \right) D = 0,
\end{aligned}
\right.
\qquad (5)$$

where $\kappa = (d\bar{\gamma}/dT)/\bar{\gamma}$ is the relative isentropic ratio variation. This system is an extension to heated and non-isentropic flow of the Linearized Euler Equations (LEE) system described in [8, 10, 13].

Assuming a mono-frequency forcing, the linearized quasi-one-dimensional Euler system eq. (5) is rewritten using the convention $y'(x,t) = \hat{y}(x) \exp(+j\omega t)$, where $\hat{y}$ is the complex amplitude of the Fourier transform of the quantity $y$, and $\omega = 2\pi f$ is the angular frequency. Introducing the unknown vector $\hat{X} = (\hat{U}, \hat{P}, \hat{D})^T$ leads to the following matrix equation:

$$\frac{d\hat{X}}{dx} + \bar{K}_1^{-1} \bar{K}_2 \hat{X} = 0, \qquad (6)$$

with

$$\bar{K}_1^{-1} = \frac{1}{\bar{u}(1 - \bar{M}^2)} \begin{bmatrix} 0 & -\bar{M}^2 & 1 \\ 0 & \bar{M}^2 & -\bar{M}^2 \\ 1 - \bar{M}^2 & \bar{M}^2 & -1 \end{bmatrix}$$

where $\bar{M} = \bar{u}/\bar{c}$ is the steady-flow Mach number, and
$\bar{K}_2 =$

$$\bar{K}_2 = \begin{bmatrix} 0 & 0 & j\omega \\ 2 \frac{d\bar{u}}{dx} + j\omega & \frac{\kappa\bar{c}^2}{\bar{u}} \frac{d\bar{T}}{dx} - \bar{\gamma} \frac{d\bar{u}}{dx} & \frac{d\bar{u}}{dx} \\ \frac{(\bar{\gamma} - 1)\bar{P}_v}{\bar{\gamma}\bar{p}} & \bar{K}_2^{32} & \kappa\bar{T} \left( \frac{\bar{P}_v}{\bar{p}} + \frac{\bar{u}}{\bar{p}} \frac{d\bar{\rho}}{dx} \right) \end{bmatrix} \qquad (7)$$

where the matrix coefficient $\bar{K}_2^{32}$ is given by:

$$\bar{K}_2^{32} = \frac{\bar{P}_v}{\bar{p}} ((\bar{\gamma} - 1) - \kappa\bar{\gamma}\bar{T}) + \kappa\bar{u} \left( \frac{d\bar{T}}{dx} - \frac{\bar{\gamma}\bar{T}}{\bar{p}} \frac{d\bar{\rho}}{dx} \right) + j\omega. \qquad (8)$$

It is noteworthy that $\bar{K}_1$ is invertible if and only if $\bar{u} \neq 0$ and $\bar{M} \neq 1$ as pointed out by Duran and Moreau [10]. This implies that the present acoustic model is not valid at the sonic (critical) throat.

Introducing the matrix $\bar{C} = -\bar{K}_1^{-1}\bar{K}_2$ and assuming a discretization of the nozzle over $j \in [0, ..., l - 1, l, l + 1, ..., n]$ elements (the index $l$ denotes the critical-throat position), each of size $\Delta x$, and using an explicit fourth-order Runge-Kutta (RK4) numerical scheme [34] for integration, the following relation is obtained:

$$\hat{X}_{j+1} = \hat{X}_j + \frac{\Delta x}{6} (k_1 + 2k_2 + 2k_3 + k_4), \qquad (9)$$

with the intermediate slopes defined as, $\forall j \neq l$:

$$\begin{aligned}
k_1 &= \bar{C}_j \hat{X}_j, \\
k_2 &= \frac{\bar{C}_j + \bar{C}_{j+1}}{2} \left( \hat{X}_j + \frac{\Delta x}{2} k_1 \right), \\
k_3 &= \frac{\bar{C}_j + \bar{C}_{j+1}}{2} \left( \hat{X}_j + \frac{\Delta x}{2} k_2 \right), \\
k_4 &= \bar{C}_{j+1} \left( \hat{X}_j + \Delta x k_3 \right).
\end{aligned}$$
To ensure convergence of the solution, a Newton-Raphson root-finding method is employed to impose a continuity condition at the critical-throat location (see section 3).

## 2.3 Waves decomposition and acoustic coefficients

The acoustic and entropy waves propagating through the flow are introduced via the inlet boundary conditions (see section 3) and are used to solve the system eq. (6). However, although the system is solved for the fluctuations $\hat{U}(x)$, $\hat{P}(x)$, and $\hat{D}(x)$ along the spatial domain $x$, these variables do not directly correspond to pure acoustic or entropy plane waves. Therefore, a wave decomposition is required.

In a region of uniform flow, the acoustic pressure field can be decomposed into downstream- ($\hat{P}^+$) and upstream-propagating ($\hat{P}^-$) wave components using the following relations:

$$ \hat{P}^{\pm} = \frac{1}{2} (\hat{P} \pm \bar{M} \hat{U}), \quad (10) $$

Such a uniform flow region is assumed both upstream and downstream of the nozzle (i.e., at the inlet and the outlet), as detailed in fig. 2.

The linearized entropy fluctuation is defined as $\sigma = s'/\bar{c}_p$, which, in the frequency domain and according to the Gibbs energy equation eq. (4), becomes:

$$ \hat{\sigma} = \frac{\hat{s}}{\bar{c}_p} = \hat{P} - \hat{D}. \quad (11) $$

Since entropy disturbances are convected with the mean flow $\bar{u}$, they propagate in the downstream direction.

Based on the decomposition into downstream- and upstream-propagating plane waves, the reflection and transmission coefficients can be defined. For isentropic acoustic waves, the inlet acoustic reflection coefficient $R_a$ is defined as the ratio of the upstream- to downstream-propagating acoustic wave amplitudes at the inlet (subscript "in"):

$$ R_a = \frac{\hat{P}_{\text{in}}^{-}}{\hat{P}_{\text{in}}^{+}}. \quad (12) $$

The acoustic transmission coefficient $T_a$ is defined as the ratio of the downstream-propagating isentropic acoustic wave amplitudes at the outlet (subscript "out") and inlet:

$$ T_a = \frac{\hat{P}_{\text{out}}^{+}}{\hat{P}_{\text{in}}^{+}}. \quad (13) $$

Similarly, the entropy-noise production coefficient $R_s$ is defined as the ratio of the upstream-propagating

isentropic acoustic wave to the incoming entropy wave for a vanishing incoming acoustic wave ($\hat{P}_{\text{in}}^{+} = 0$):

$$ R_s = \frac{\hat{P}_{\text{in}}^{-}}{\hat{\sigma}_{\text{in}}}, \quad (14) $$

and the entropy-noise transmission coefficient $T_s$ is the ratio of the downstream-propagating isentropic wave at the outlet to the entropy wave at the inlet:

$$ T_s = \frac{\hat{P}_{\text{out}}^{+}}{\hat{\sigma}_{\text{in}}}. \quad (15) $$

It is noteworthy that this definition of the transmission coefficient assumes that the flow is uniform and non-heated at the outlet of the geometry. Consequently, the coefficient depends on the specific location where the outlet quantities are measured. This coefficient is primarily used to validate the supersonic region of the quasi-one-dimensional linear acoustic model.

## 3 Boundary conditions

To solve eq. (6), three boundary conditions must be applied at the system boundaries. Following the methodology of Duran and Moreau [10], the unknown vector $\hat{X}(x)$ is computed in two stages (see fig. 2).

In the subsonic region, two boundary conditions must be applied at the inlet ($x_{\text{in}}$), and one at the upstream side of the critical-throat (denoted $x_{\text{cu}}$).

To ensure that the system eq. (6) remains invertible ($\bar{M} \neq 1$), the locations where boundary conditions are applied near the critical throat are shifted slightly away from the sonic point by an infinitesimal distance $\epsilon$, as illustrated in fig. 2, leading to $x_{\text{cu}} = x_{*} - \epsilon$ and $x_{\text{cd}} = x_{*} + \epsilon$. This distance $\epsilon$ is set equal to the spatial discretization $\Delta x$.

As only two boundary conditions are imposed at the inlet, a root-finding method must be employed to determine the third component of the unknown vector $\hat{X}_{\text{in}}$. For instance, $\hat{P}_{\text{in}}$ and $\hat{D}_{\text{in}}$ are prescribed as input values, while the third component $\hat{U}$ is initially guessed, forming the inlet unknown vector $\hat{X}_{\text{in},g} = (\hat{U}_{\text{guess}}, \hat{P}_{\text{in}}, \hat{D}_{\text{in}})^T$. The linear acoustic model eq. (9) is then used to compute the spatial profiles of $\hat{P}(x)$, $\hat{D}(x)$, and $\hat{U}(x)$.

At the critical-throat position $x_{\text{cu}}$, the corresponding boundary condition (which will be detailed in section 3.2) is evaluated as a root function, denoted $h(\hat{X}_{\text{in},g}) = 0$. Based on the outcome, the guessed inlet condition is updated using the Newton-Raphson method:

$$ \hat{X}_{\text{in},g+1} = \hat{X}_{\text{in},g} - \frac{h(\hat{X}_{\text{in},g})}{h'(\hat{X}_{\text{in},g})}, \quad (16) $$
<figure><img src="image_3.png" /><figcaption>Figure 2: Sketch of the boundary condition treatment for a thermally choked flow configuration. A similar decomposition can be applied to a nozzle with a geometric throat.</figcaption></figure>

where the derivative $h'(\hat{X}_{in,g})$ is approximated using a first-order explicit finite-difference scheme.

To ensure convergence of the root-finding method, the relative error on the critical-throat boundary condition is set to $1 \times 10^{-8}$. The initial guess of the unknown vector is based on an isentropic downstream-propagating acoustic wave for acoustic forcing and downstream-propagating entropy wave for entropy forcing. It has been verified that this initial guess does not affect the final value of $\hat{X}_{in}$.

Subsequently, the supersonic region is solved by imposing three boundary conditions just downstream of the critical-throat position, denoted $x_{cd}$ and no conditions at the outlet $\bar{x}_{out}$. The resulting solution satisfies the continuity condition $\hat{X}(x_{cd}) = \hat{X}(x_{cu})$.

## 3.1 Inlet boundary conditions

For an isentropic acoustic forcing of angular frequency $\omega = 2\pi f$, the inlet acoustic pressure fluctuation $\hat{P}_{in}$ is imposed with an amplitude $\eta$. As the forcing is assumed to be perfectly isentropic, i.e., $\hat{\sigma}_{in} = \hat{P}_{in} - \hat{D}_{in} = 0$, the following inlet boundary conditions apply:

$$ \left\{ \begin{array}{l} \hat{P}_{in} = \eta, \\ \hat{D}_{in} = \hat{P}_{in}. \end{array} \right. \qquad (17) $$

For entropy forcing of angular frequency $\omega$, an entropy amplitude $\eta$ is applied at the inlet, and the downstream-propagating acoustic wave is set to zero, i.e., $\hat{P}_{in}^{+} = \hat{P}_{in} + \bar{M}_{in} \hat{U}_{in} = 0$, resulting in the following two inlet boundary conditions:

$$ \left\{ \begin{array}{l} \hat{D}_{in} = \hat{P}_{in} - \eta, \\ \hat{U}_{in} = -\hat{P}_{in}/\bar{M}_{in}. \end{array} \right. \qquad (18) $$

## 3.2 Critical-throat acoustic-boundary condition

In the literature [8–13], a commonly used approach is the quasi-steady critical Mach number assumption, which states that Mach number fluctuations vanish at the steady critical-throat position, i.e., $M'_{*} = 0$.

The subscript $\ast$ hereafter denotes quantities evaluated at the steady critical-throat location. For isentropic gas flows with constant heat capacity ratio $\gamma$, this quasi-steady critical-throat boundary condition (denoted QSS) can be expressed as:

$$ 2\hat{U}_{\ast} = \gamma \hat{P}_{\ast} - \hat{D}_{\ast}. \qquad (19) $$

Using the linearized Gibbs equation (eq. (11)), the entropy-based form of this condition, as given by Marble and Candel [8] and Duran and Moreau [10], becomes:

$$ 2\hat{U}_{\ast} = (\gamma - 1)\hat{P}_{\ast} + \hat{\sigma}_{\ast}. \qquad (20) $$

However, it has been shown in previous works [15, 26] that the critical-throat position is not steady under acoustic and entropy mono-harmonic plane wave forcing, and that the quasi-steady assumption does not hold, even for isentropic nozzles (except at very low frequencies). In what follows, a new boundary condition is proposed for thermally-choked-flows.

As performed by Tsien [16] and Marble and Candel [8], the steady flow velocity profile, $\bar{u}(x)$, is locally approximated at the throat as $\bar{u}(\tilde{x}) = (d\bar{u}/dx)_{\ast} \tilde{x}$, where $(d\bar{u}/dx)_{\ast}$ denotes the velocity gradient at the critical-throat. The origin of the coordinate $\tilde{x}$ is chosen such that the extrapolated velocity profile $\bar{u}(\tilde{x})$ vanishes at $\tilde{x} = 0$. The throat position is $\tilde{x}_{\ast} = \bar{c}_{\ast}/(d\bar{u}/dx)_{\ast}$. Based on this definition, the following time-space transformation is introduced as performed by Marble and Candel [8]:

$$ \tau = \frac{\bar{c}_{\ast}}{\tilde{x}_{\ast}} t \quad ; \quad \xi = \left( \frac{\tilde{x}}{\tilde{x}_{\ast}} \right)^2. $$

This boundary condition is based on the linearized quasi-one-dimensional Euler equations eq. (5), rewritten using the above time-space transformation. Introducing the dimensionless frequency $\Omega = 2\pi f/(d\bar{u}/dx)_{\ast}$, the system eq. (5) is further expressed in the frequency domain. This dimensionless frequency characterizes how the critical throat interacts with the propagating waves through the local velocity gradient. The condition is derived by subtracting the energy equation from the momentum equation (see appendix A for details). The assumption of Tsien [16], as proposed by Crocco, is adopted such that $\hat{U}$, $\hat{P}$, and $\hat{D}$, as well as their derivatives, are continuous at the critical-throat location. It is further assumed that none of the leading-order terms in the expansion of these quantities is proportional to $1/(1 - \xi)$. Such a condition would lead to a singularity at the critical throat. Under these assumptions, and because at the throat $\xi = 1$, $\bar{c}_{\ast} = \bar{c}$ and terms in $\xi - 1$ vanish, the generalized boundary condition at the critical-throat position (denoted GEN) becomes:

$$ \bar{g}_{U} \hat{U}_{\ast} = \bar{g}_{P} \hat{P}_{\ast} + \bar{g}_{D} \hat{D}_{\ast}, \qquad (21) $$
where the coefficients are defined as:

$$
\left\{
\begin{aligned}
\bar{g}_U &= 2 + j\Omega - \frac{(\bar{\gamma}_* - 1)}{(d\bar{u}/dx)_*} \frac{\bar{P}_{v,*}}{\bar{\gamma}_*\bar{p}_*}, \\
\bar{g}_P &= \bar{\gamma}_* + j\Omega + \frac{\bar{P}_{v,*}}{\bar{p}_*(d\bar{u}/dx)_*} \left[ (\bar{\gamma}_* - 1) - \kappa_* \bar{\gamma}_* \bar{T}_* \right] \\
&\quad - \frac{2\kappa_* \bar{\gamma}_* \bar{T}_*}{\bar{\rho}_*} \left. \frac{d\bar{\rho}}{d\xi} \right|_* , \\
\bar{g}_D &= \kappa_* \bar{T}_* \left( \frac{\bar{P}_{v,*}}{\bar{p}_*(d\bar{u}/dx)_*} + \frac{2}{\bar{\rho}_*} \left. \frac{d\bar{\rho}}{d\xi} \right|_* \right) - 1.
\end{aligned}
\right.
\quad (22)
$$

This critical-throat boundary condition is an ex-
tension of the boundary condition for isentropic flow
proposed by Olivon et al. [15], to heated and non-
calorically perfect gas flows. Further details on the
derivation of this model are provided in appendix A.

The asymptotic behaviour of this boundary con-
dition is also of interest. At very high frequencies
(Ω ≫ 1), the relation ̂U* = P* emerges, indicating
that in a heated medium, the generation of entropy
due to propagating acoustic waves becomes negligible.
This is because ̂D* = P* - ̂σ* becomes negligible com-
pared to P*. Furthermore, this relation suggests that
acoustic reflection at the critical throat vanishes, and
the boundary condition can be approximated at high
frequencies by an isentropic, non-reflective boundary
condition.

As explained previously, at the critical-upstream
position xcu the boundary condition is evaluated as
a root function. Depending on whether the quasi-
steady (QSS) or the generalized formulation (GEN)
is chosen, it writes :

$$
\text{QSS: } 2\hat{U}_{*} - \gamma\hat{P}_{*} + \hat{D}_{*} = 0, \quad (23)
$$

$$
\text{GEN: } \bar{g}_U \hat{U}_* - \bar{g}_P \hat{P}_* - \bar{g}_D \hat{D}_* = 0. \quad (24)
$$

4 Validation methodology

To validate the GEN boundary condition proposed
in eq. (21), results from the quasi-one-dimensional
acoustic model using this boundary condition (eq. (6))
are compared with two-dimensional Computational
Fluid Dynamics (CFD) simulations. These simula-
tions are conducted at low amplitudes, with η =
1 × 10⁻², ensuring that a linear unsteady-perturbation
propagation is expected. Two choking configurations
are investigated: an isentropic geometric-throat con-
figuration and a thermally-choked-flow configuration.
The results comparisons are performed on the acous-
tic coefficients (Rₐ, Rₛ, Tₐ and Tₛ) previously defined
in section 2.3.

In the following, only the main steps of the method-
ology are outlined.

<figure><img src="image_4.png" /><figcaption>Figure 3: Simplified two-dimensional planar sketch of a dual-mode ramjet geometry (not to scale).</figcaption></figure>

<table><caption>Table 1: Nozzle-inlet Mach number M<sub>in</sub>, total tem-
perature T<sub>0,in</sub>, and static pressure p<sub>in</sub>, corresponding
to a typical flight Mach number M = 3.</caption>
  <thead>
    <tr>
      <th>Parameter</th>
      <th>Value</th>
      <th>Unit</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>Mach number: M<sub>in</sub></td>
      <td>0.6</td>
      <td>-</td>
    </tr>
    <tr>
      <td>Total temperature: T<sub>0,in</sub></td>
      <td>600</td>
      <td>K</td>
    </tr>
    <tr>
      <td>Static pressure: p<sub>in</sub></td>
      <td>2</td>
      <td>bar</td>
    </tr>
  </tbody>
</table>

4.1 Thermally choked flow nozzle

For the thermally-choked flow configuration, a simpli-
fied dual-mode ramjet geometry, based on the works
of [6, 14, 26, 28], is used, as shown in fig. 3. The
nozzle geometry consists of three sections: a straight-
duct upstream section (0 ≤ x/L < x<sub>u</sub>/L = 0.42),
a smooth transition section (x<sub>u</sub>/L = 0.42 ≤ x/L <
x<sub>d</sub>/L = 0.58) formed by a wall in circular arc with
angle α = 3°, and a divergent duct with a constant di-
vergence angle α = 3° (x/L > x<sub>d</sub>/L = 0.58). The to-
tal length and inlet height of the nozzle are L = 2.4 m
and y<sub>in</sub> = 25 × 10<sup>-3</sup> m, respectively. The circular arc
radius of the transition region is R<sub>t</sub>/y<sub>in</sub> = 305. A slip-
adiabatic boundary condition is imposed on the upper
nozzle wall, while a symmetry boundary condition is
applied along the lower edge.

The incoming airflow at the inlet is subsonic, as de-
tailed in table 1. The inlet conditions are representa-
tive of typical dual-mode ramjet engine operation [6].

To establish a steady thermal throat deep within
the divergent nozzle via non-uniform heat release, the
approach from prior works [14, 26, 28] is adopted.
The spatial distribution of P<sub>v</sub>(x) is divided into two
regions. In the straight duct (fig. 3), the flow is not
heated, P<sub>v</sub> is set to zero. In the divergent section
(x > x<sub>u</sub>), the volumetric power distribution P<sub>v</sub> follows
a log-normal profile along x, as proposed by Wang
et al. [4] and Olivon et al. [14].

To determine P<sub>v</sub>(x), a quasi-one-dimensional
steady-state model is used, based on the Mach num-
ber and static temperature formulations by Shapiro
[35] for an ideal gas, and the steady-state thermal
throat relation of Heiser and Pratt [36]. This model
does not directly define P<sub>v</sub>(x). It determines the gra-
dient of the heat release per unit mass, dq/dx:
$$P_v(x) = \frac{\dot{m}}{A(x)} \frac{dq}{dx}, \qquad (25)$$

where $\dot{m} = \bar{\rho} \bar{u} A$ is the quasi-one-dimensional steady mass flow rate, calculated using isentropic relations and the inlet conditions provided in table 1. The model is tuned to achieve a steady-state thermal throat located at $x_*/L = 0.8$, corresponding to a flame length of $(x_* - x_u)/L = \delta_*/L = 0.38$. This location aligns with the mean spatial position of the thermal throat in Case C4 of the JAPHAR study [6], where the feasibility of a thermally-choked nozzle was experimentally demonstrated for a hydrogen inflow at Mach 5.

## 4.2 Geometrically choked flow nozzle

To validate the critical-throat acoustic boundary condition, an isentropic flow through a nozzle with a geometric throat is also considered. This nozzle is constructed using the same spatial Mach number profile, $M(x)$, as that obtained from the thermally choked flow configuration described in section 4.1. The geometry of the isentropic flow nozzle, denoted $A_{is}$, is determined using the following expression, assuming a constant specific heat ratio $\gamma_{in}$:

$$A_{is}(x) = \frac{A_{in}M_{in}}{M_{tht}} \left( \frac{2 + (\gamma_{in} - 1)M_{tht}^2}{2 + (\gamma_{in} - 1)M_{in}^2} \right) \frac{\gamma_{in} + 1}{2(\gamma_{in} - 1)}. \qquad (26)$$

The steady Mach number profile $M_{tht}(x)$ is extracted along the symmetry axis from the CFD simulation of the thermally choked flow and is used to compute the equivalent isentropic area distribution $A_{is}(x)$. Since the cross-section $A_{is}(x)$ is calculated using a constant Poisson ratio $\gamma = \gamma_{in}$, and further isentropic-flow calculations account for temperature variation (see eq. (3)), a small deviation arises between $M_{is}(x)$ and $M_{tht}(x)$. This deviation remains below $5 \times 10^{-4}$. This has no impact on the conclusions drawn from comparing CFD solver results with the linear quasi-one-dimensional flow model.

## 4.3 Computational setup

The CFD simulations are performed using the ONERA unstructured code CEDRE [37], which solves the Euler equations [14] (see appendix B for the full equations).

The computational setup (see fig. 3) consists of a two-dimensional planar domain, as described in section 4.1, meshed with approximately 60 000 quadrilateral cells of longitudinal size $\Delta x/L = 4.1 \times 10^{-4}$. This resolution ensures adequate spatial discretization of acoustic waves, providing at least 100 points per wavelength for frequencies up to $\Omega = 2\pi f/(d\bar{u}/dx)_* = 30$. Since the flow is primarily longitudinal, the

same discretization ($\Delta x \approx \Delta y$) is applied in both the longitudinal and transverse directions, resulting in approximately 2400 points along the longitudinal axis and 25 points along the transverse axis. Spatial discretization employs a second-order multi-slope MUSCL scheme [38] and an HLLC flux scheme [39].

Temporal integration is carried out using an implicit second-order Runge-Kutta scheme [40], with a time step of $\Delta t = 1 \times 10^{-6}$ s, ensuring a Courant-Friedrichs-Lewy (CFL) number of $(\bar{u}_{in} + \bar{c}_{in})\Delta t/\Delta x = 0.78 < 1$. More details on mesh and time independence can be found in appendix B.

Unsteady solutions corresponding to acoustic forcing are obtained using a modulated post-processing time step based on the inlet forcing frequency. To ensure satisfactory spectral analysis, 75 periods are computed, and 100 points per period are stored. Spectral analysis (Fast Fourier Transform) is performed on the last 35 periods to eliminate transient effects.

For the geometric-throat configuration (see section 4.2), the same numerical parameters are used to ensure mesh and time-step independence.

## 4.4 Steady base-flow fields

Once the steady-state is reached, the flow fields (Mach number M, static pressure $p$, and static temperature $T$) are extracted from the CFD solver and displayed in fig. 4. As expected, the thermal throat is located at the reduced position $x_*/L = 0.8$. Due to the non-uniform volumetric heat source term $P_v$, the steady Mach number increases continuously. This increase remains smooth around the critical-throat position, validating the assumption of a locally linear velocity profile at the throat, as required by the critical-throat acoustic boundary condition framework proposed in eq. (21).

For both thermally and geometrically choked flow configurations, the steady Mach number profiles are similar, with a maximum relative deviation of approximately $5 \times 10^{-4}$ between the two cases. As the velocity increases, the static pressure decreases accordingly. In the thermally choked configuration, energy addition causes the static temperature to rise from $T_{in} = 560$ K at the inlet to $T_* \approx 1740$ K at the throat. This significant temperature increase induces a variation in the heat capacity ratio, from $\gamma_{in} = 1.379$ to $\gamma_* = 1.304$, highlighting the importance of accounting for this variation in the quasi-one-dimensional linear acoustic model (eq. (6)).

This steady flow-field is used to as an input to quasi-one-dimensional linear acoustic model (eq. (6)).
<figure><img src="image_5.png" /><figcaption>Figure 4: (Colour online) Steady-state flow fields as a function of the longitudinal length $x/L$ for the geometric (· · ·) and thermal throat (---) choking cases obtained from the CFD solver. The steady Mach number ($\bar{M}$, (a)); the static pressure ($\bar{p}$, (b)) and the static temperature ($\bar{T}$, (c)) are displayed. Quantities are measured on the symmetry axis of the computational domain.</figcaption></figure>

# 5 Discussion and results

## 5.1 Validation of the critical-throat acoustic-boundary condition

The GEN boundary condition (eqs. (21) and (22)) is validated using thermally and geometrically choked flow configurations. To this end, the Euler equations are solved using both the quasi-one-dimensional linear acoustic model eq. (6) and the CFD solver described in section 4.3. To highlight the impact of the new critical-throat boundary condition, results from the quasi-one-dimensional linear acoustic model using the QSS boundary condition (eq. (19)) are also presented.

Figure 5 shows the normalized critical-throat Mach number fluctuations, $\hat{M}_*/\hat{P}_*$, resulting from the isentropic acoustic forcing defined in eq. (17), plotted as a function of the dimensionless frequency $\Omega = 2\pi f/(d\bar{u}/dx)_*$. Both thermally choked and geometrically choked configurations are considered. Results are shown for the CFD solver: thermal throat (○) and geometric throat (□); and for the quasi-one-dimensional linear acoustic model using either the GEN or the QSS boundary conditions, respectively: thermal throat (---, · · ·) and geometric throat (· · ·, ---).

Based on the GEN boundary condition detailed in eqs. (21) and (22), the analytical expression for the normalised critical-throat Mach number fluctuations

is given by:

$$ \frac{\hat{M}_*}{\hat{P}_*} = \left[ \frac{\bar{g}_P}{\bar{g}_U} - \bar{\gamma}_* \frac{(1 + \kappa_* \bar{T}_*)}{2} \right] + \left[ \frac{\bar{g}_D}{\bar{g}_U} - \frac{(1 + \kappa_* \bar{T}_*)}{2} \right] \frac{\hat{D}_*}{\hat{P}_*} \quad (27) $$

The relative deviation between the quasi-one-dimensional linear acoustic model (L) results (with GEN or QSS boundary condition) and the CFD results is quantified as $|1 - (\hat{M}_*/\hat{P}_*)^L/(\hat{M}_*/\hat{P}_*)^{CFD}|$. The forcing amplitude is set to $|\hat{P}_{in}| = \eta = 10^{-2}$, ensuring that non-linear effects remain negligible. To verify the linearity of the response, three frequencies $\Omega = [3, 8.5, 22]$ were tested with two isentropic forcing amplitudes ($|P'_{in}| = 1 \times 10^{-2}$ and $1 \times 10^{-3}$). The relative deviations between the responses were on the order of $1 \times 10^{-3}$, confirming linear behaviour.

For the geometric-throat configuration, the QSS boundary condition $\hat{M}_* = 0$ (---) is only valid in the very low-frequency limit and is generally inaccurate. A relative error of 100% is obtained when Mach number fluctuations are imposed to vanish at the throat. In contrast, the linear model using the GEN boundary condition (· · ·) closely matches the CFD results (□), with a maximum deviation of approximately 2% at very low frequencies, which decreases monotonically as the frequency increases.

In the thermally choked configuration, the GEN boundary condition also performs well, with deviations of about 10% for $\Omega < 7$, dropping below 0.1% at higher frequencies, demonstrating the accuracy of the proposed linear acoustic model (---). At very low frequencies, a finite value of the normalized Mach number fluctuations $\lim_{\Omega \to 0} |\hat{M}_*/\hat{P}_*| \approx 0.2$ is observed in the thermal throat configuration due to the presence of the volumetric power addition $P_v$, contradicting the commonly used QSS assumption (---). Note that applying the QSS condition, $M'_* = 0$, leads to a local singularity in the acoustic field at the critical-throat, which contradicts the physical assumption of a continuous and smooth solution.

Under the GEN boundary condition, a local maximum in the deviation is observed proximate to the frequency at which the normalized critical Mach number fluctuation $\hat{M}_*/\hat{P}_*$ peaks. This maximum, occurring at $\Omega = 3.1$, is primarily attributed to the generation of entropy waves induced by the propagation of acoustic waves through a heated medium. For an unheated flow, i.e., in the geometric throat case, the density and pressure fluctuations satisfy $\hat{D}_*/\hat{P}_* = 1 - \hat{\sigma}_*/\hat{P}_* = 1$ across all frequencies, resulting in the absence of a local maximum, as shown in fig. 5. In contrast, for a heated flow, the ratio $\hat{M}_*/\hat{P}_*$ becomes strongly dependent on the entropy generation term, $\hat{D}_*/\hat{P}_* = 1 - \hat{\sigma}_*/\hat{P}_*$, since the density fluctuations $\hat{D}_*$ are directly linked to entropy fluctuations via eq. (11). At low frequencies ($0 < \Omega < 3.1$), the entropy waves generated by the propagation of acoustic waves induce
<figure><img src="image_6.png" /><figcaption>Figure 5: (Colour online) Modulus (a) and phase angle (b) of the normalized critical Mach number fluctuations due to isentropic acoustic forcing at the inlet (eq. (17)), $\hat{M}_*/\hat{P}_*$, plotted against the dimensionless frequency $\Omega = 2\pi f/(d\bar{u}/dx)_*$. Symbols: (○) and (□) denote CFD results for the thermal and geometric throats, respectively; (---, --) and (···, ···) represent the corresponding results from the linear acoustic model (L) using the GEN (eqs. (21) and (22)) or QSS (eq. (19)) boundary conditions, respectively. Panel (c) displays the relative deviation between the CFD and linear models results.</figcaption></figure>

additional pressure fluctuations, thereby amplifying the Mach number fluctuations at the critical throat. At moderate frequencies ($3.1 < \Omega < 10$), the generation of entropy waves induced by acoustic forcing decreases, leading to reduced Mach number fluctuations at the throat. At high frequencies, $\Omega = O(10)$, the ratio $\hat{D}_*/\hat{P}_*$ asymptotically approaches unity, indicating that acoustic wave propagation no longer generates significant entropy fluctuations.

Figure 6 shows the response at the critical throat under entropy-wave forcing defined by eq. (18) for both geometric-throat (···, ---) and thermal-throat (----, --) choking configurations. These are computed using the quasi-one-dimensional linear acoustic model with the GEN (eqs. (21) and (22)) and QSS (eq. (19)) boundary conditions, respectively. The corresponding CFD results are also shown for the thermal-throat (○) and geometric-throat (□) cases. The normalized critical Mach number fluctuations, $\hat{M}_*/\hat{\sigma}_*$, are plotted as a function of the dimensionless frequency $\Omega = 2\pi f/(d\bar{u}/dx)_*$.

The relative deviation between the quasi-one-dimensional linear acoustic model (L) and the CFD simulations is quantified as $|1 - (\hat{M}_*/\hat{\sigma}_*)^L/(\hat{M}_*/\hat{\sigma}_*)^{CFD}|$. The amplitude of the incoming entropy wave, $\hat{\sigma}_{in}$, is set to the order of $10^{-2}$, ensuring that non-linear effects remain negligible. For the geometric-throat case (□), the deviation between CFD and the linear model with the GEN boundary condition remains within a few percent at low frequencies and decreases rapidly as frequency increases. In the thermally choked configuration (○), a larger deviation, up to 61%, is observed at very low frequencies, but it drops below 1% for $\Omega > 5$. In both choking cases, the QSS boundary condition yields a constant relative deviation of 100%, due to the erroneous prescription of a vanishing critical Mach number fluctuations $\hat{M}_* = 0$.

For both types of inlet forcing (acoustic and entropy), the GEN boundary condition therefore demonstrates a fair accuracy in capturing the Mach number fluctuations at the critical throat for both geometrically and thermally choked configurations.

## 5.2 Acoustic-reflection coefficients

Figure 7 presents the acoustic-reflection coefficient, $R_a = \hat{P}_{in}^-/\hat{P}_{in}^+$, defined as the ratio of the upstream- to downstream-propagating acoustic waves at the inlet, for both geometric-throat and thermal-throat configurations. Results from both the quasi-one-dimensional linear acoustic model using both GEN (eqs. (21) and (22)) or QSS (eq. (19)) boundary condition (lines) and CFD solver (markers) are shown as a function of the dimensionless frequency $\Omega = 2\pi f/(d\bar{u}/dx)_*$.

For both choking configurations, a fair agreement is observed between the quasi-one-dimensional linear acoustic model using the GEN boundary condition
<figure><img src="image_7.png" /><figcaption>Figure 6: (Colour online) Modulus (a) and phase angle (b) of the normalized critical Mach number fluctuations, $\hat{M}_*/\hat{\sigma}_*$, resulting from entropy forcing at the inlet (eq. (18)), plotted against the dimensionless frequency $\Omega = 2\pi f/(d\bar{u}/dx)_*$. Symbols: (○) and (□) denote CFD results for the thermal and geometric throats, respectively; (---, ---) and (···, ···) indicate the corresponding results from the linear acoustic model (L) using the GEN (eqs. (21) and (22)) or QSS (eq. (19)) boundary conditions, respectively. Panel (c) presents the relative deviation between the CFD and linear models results.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7: (Colour online) Modulus (a) and phase angle (b) of the acoustic reflection coefficient, $R_a = \hat{P}_{in}^-/\hat{P}_{in}^+$, plotted as a function of the dimensionless frequency $\Omega = 2\pi f/(d\bar{u}/dx)_*$. Symbols: (○) and (□) denote CFD results for the thermal and geometric throats, respectively; (---, ---) and (···, ···) indicate the corresponding results from the linear acoustic model (L) using the GEN (eqs. (21) and (22)) or QSS (eq. (19)) boundary conditions, respectively. Panel (c) presents the relative deviation between the CFD and linear models results.</figcaption></figure>

(thermal throat: ---, geometric throat: ···) and the CFD solver (thermal throat: ○, geometric throat: □), as indicated by the relative deviation $|1 - R_a^L/R_a^{CFD}|$. The linear acoustic model accurately predicts the reflection coefficient for both the geometric-throat (···) and thermal-throat (---) configurations, with deviations remaining within 1% at low frequencies ($\Omega < 10$) when using the GEN boundary condition. At higher frequencies ($\Omega > 15$), however, discrepancies in both the modulus and phase angle become apparent, with deviations reaching up to 36%. This increase in deviation has been verified to be independent of the mesh resolution, time step, and forcing amplitude.

The deviation between CFD and the quasi-one-dimensional linear acoustic model using the QSS boundary condition is also shown for the geometric (---) and thermal throat (---) configurations. At low frequencies ($\Omega < 10$), the QSS boundary condition leads to deviations of up to 50% for the thermal throat case, and only a few percent for the geometric case. This highlights the importance of the GEN boundary condition in accurately capturing the behaviour of thermally choked flows, where the QSS assumption fails. For the geometric case, the improvement brought by the new boundary condition is less significant.

The entropy-noise-production coefficient, $R_s = 729$
<figure><img src="image_9.png" /></figure>

Figure 8: (Colour online) Modulus (a) and phase angle (b) of the entropy-noise production coefficient, $R_s = \hat{P}_{in}^-/\hat{\sigma}_{in}$, plotted as a function of the dimensionless frequency $\Omega = 2\pi f/(d\bar{u}/dx)_*$. Symbols: (o) and (□) denote CFD results for the thermal and geometric throats, respectively; (---, --) and (···, ···) indicate the corresponding results from the linear acoustic model (L) using the GEN (eqs. (21) and (22)) or QSS (eq. (19)) boundary conditions, respectively. Panel (c) presents the relative deviation between the CFD and linear models results.

$\hat{P}_{in}^-/\hat{\sigma}_{in}$, is plotted in fig. 8 as a function of the dimensionless frequency $\Omega = 2\pi f/(d\bar{u}/dx)_*$. This coefficient represents the ratio between the upstream-propagating acoustic wave $\hat{P}_{in}^-$ and the incoming entropy wave $\hat{\sigma}_{in}$.

Fair agreement is observed between the quasi-one-dimensional linear acoustic model using the GEN boundary condition (··· and ---) and CFD solver (□ and o) for both the geometric-throat and thermal-throat configurations, respectively. The relative deviation, defined as $|1 - R_s^L/R_s^{CFD}|$, remains below 2% for frequencies below $\Omega = 5$. In all cases, the deviation remains below 25% for the thermally-choked-flow configuration (o) and below 8% for the geometric-throat nozzle (□).

The deviation between the CFD results and the quasi-one-dimensional linear acoustic model using the QSS boundary condition is higher than that obtained with the GEN boundary condition for both the thermal (---) and geometric (····) throat configurations. For frequencies below $\Omega < 10$, the QSS boundary condition results in a deviation of a few percent for the geometric case and up to 40% for the thermal throat configuration.

Here again, the linear acoustic model based on the newly proposed critical-throat boundary condi-

tion (eqs. (21) and (22)) shows good agreement with the two-dimensional CFD results. For the thermal throat configuration, the GEN boundary condition significantly reduces the deviation between the CFD and the linear acoustic model. In contrast, the benefits of this new boundary condition are less pronounced for the geometric case.

## 5.3 Acoustic-transmission coefficients

The transmission coefficient is now analyzed for both geometric-throat and thermal-throat configurations. Figure 9 shows the transmission coefficient, $T_a = \hat{P}_{out}^+/\hat{P}_{in}^+$, as a function of the dimensionless frequency $\Omega = 2\pi f/(d\bar{u}/dx)_*$ for both configurations, as detailed in fig. 2. Results are provided for the quasi-one-dimensional linear acoustic model using either the GEN- or the QSS-boundary conditions (lines), and for the CFD solver (markers).

The quasi-one-dimensional linear model using the GEN boundary condition (thermal throat: ---, geometric throat: ···) accurately predicts the acoustic transmission coefficient $T_a$, closely matching the CFD results (thermal throat: o, geometric throat: □) across the entire frequency range considered ($\Omega \le 25$). The relative deviation, $|1 - T_a^L/T_a^{CFD}|$, quantifies the difference between the linear acoustic model and CFD solver results. For both the geometric-throat (···) and thermal-throat (---) configurations, the deviation remains below 2% at all frequencies.

Contrary to the behaviour observed in the geometric throat configuration, the acoustic transmission coefficient $T_a$ in the thermally choked flow remains below unity across all frequencies. This indicates that, at very low frequencies ($\Omega < 5$), a significant portion of the acoustic wave energy is converted into entropy waves. At higher frequencies ($\Omega > 10$), as the thermal throat approaches a purely isentropic acoustic radiation behaviour, the transmission coefficient converges toward an asymptotic value.

The deviation between the results of the CFD and the linear acoustic model using the QSS boundary condition is also detailed for both choking configurations. For the geometric case (---), the deviation is comparable to that obtained with the GEN boundary condition, except at very low and very high frequencies, suggesting that the new boundary condition has limited benefits on the prediction of the acoustic transmission coefficient $T_a$. In contrast, for the thermally choked case (···), the QSS boundary condition induces a significant deviation, up to 140% at very low frequencies and around 10% for frequencies above $\Omega > 7.5$. This highlights that the new boundary condition is essential for accurately predicting the acoustic behaviour in such thermally choked configurations.

Figure 10 presents the entropy-noise transmission coefficient, $T_s = \hat{P}_{out}^+/\hat{\sigma}_{in}$, plotted as a function of
<figure><img src="image_10.png" /><figcaption>Figure 9: (Colour online) Modulus (a) and phase angle (b) of the acoustic-transmission coefficient $T_a = \hat{P}_{\text{out}}^+/\hat{P}_{\text{in}}^+$, plotted against the dimensionless frequency $\Omega = 2\pi f/(d\bar{u}/dx)_*$. Symbols: (○) and (□) denote CFD results for the thermal and geometric throats, respectively; (---, --) and (···, ···) indicate the corresponding results from the linear acoustic model (L) using the GEN (eqs. (21) and (22)) or QSS (eq. (19)) boundary conditions, respectively. Panel (c) presents the relative deviation between the CFD and linear models results.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 10: (Colour online) Modulus (a) and phase angle (b) of the entropy-noise-transmission coefficient $T_s = \hat{P}_{\text{out}}^+/\hat{\sigma}_{\text{in}}$, plotted against the dimensionless frequency $\Omega = 2\pi f/(d\bar{u}/dx)_*$. Symbols: (○) and (□) denote CFD results for the thermal and geometric throats, respectively; (---, --) and (···, ···) indicate the corresponding results from the linear acoustic model (L) using the GEN (eqs. (21) and (22)) or QSS (eq. (19)) boundary conditions, respectively. Panel (c) presents the relative deviation between the CFD and linear model results.</figcaption></figure>

the dimensionless frequency $\Omega = 2\pi f/(d\bar{u}/dx)_*$, for both geometric-throat and thermal-throat configurations. Results from the quasi-one-dimensional linear acoustic model (lines) and the CFD solver (markers) are shown.

Fair agreement is observed between the quasi-one-dimensional linear acoustic model using the GEN boundary condition (L) and the CFD solver, with deviations of only a few percent at low frequencies for both choking configurations. The deviation increases to approximately 5% at very high frequencies. In contrast, the QSS boundary condition results in more significant discrepancies: up to 5% for the geometric-throat configuration (····), and around 22% for the thermal-throat case (---).

Once again, while the GEN boundary condition has a limited impact on the entropy-noise transmission coefficient $T_s$ in the geometric-throat configuration, it proves essential for accurate predictions in the thermally choked case.
## 6 Conclusions

Quasi-one-dimensional linear acoustic solvers are powerful tools for investigating and predicting acoustic oscillations and their sources, particularly direct and indirect combustion noise within the combustion chambers of ramjet engines. In ramjet mode, the transition from subsonic to supersonic flow is achieved through the formation of a thermal throat, driven by heat release from combustion.

In the literature, the boundary condition at the critical throat ($\bar{M}_* = 1$) is often based on the quasi-steady assumption that Mach number fluctuations vanish, i.e., $M'_* = 0$. However, the validity of the condition is questionable, especially for thermally choked flow configurations.

In the present study, an acoustic critical-throat boundary condition applicable to both geometrically- and thermally-choked throats is proposed for non-isentropic (heated) and ideal-gas flows with a uniform composition. To validate this boundary condition, a quasi-one-dimensional linear acoustic model is developed under the same framework. Its predictions are compared against results from two-dimensional Euler flow simulations. To assess the model across a range of choking scenarios, two configurations are considered: one featuring a geometric-throat, and another exhibiting a thermal-throat induced by a steady volumetric non-uniform heat source term.

The results of the quasi-one-dimensional linear acoustic model show good agreement with those of the two-dimensional CFD simulations when the new acoustic boundary condition is applied at the critical throat. In contrast, significant deviations are observed when using the commonly adopted quasi-steady critical-throat boundary condition, $M'_* = 0$, particularly for the transmission coefficients in the case of a thermal throat. It should also be noted that enforcing the QSS boundary condition at the critical throat introduces a local singularity in the acoustic field, which lacks physical justification. The proposed boundary condition eliminates this singularity and accurately captures the acoustic behaviour of thermally choked flow configurations, highlighting its importance for modelling such systems.

Further validation is required for additional thermally choked flow configurations, including variations in geometry and steady-state volumetric heat release profiles. To better capture the acoustic behaviour of such engines, the assumption of steady volumetric heat addition should be relaxed by explicitly modelling the thermoacoustic coupling. Finally, the single-species flow assumption will need to be lifted in order to investigate the effects of indirect compositional noise.

## A Quasi-one-dimensional critical-throat boundary condition

To obtain the quasi-one-dimensional critical-throat boundary condition model for a thermally choked flow nozzle, the quasi-one-dimensional linearized Euler equations for momentum and energy are considered using the perturbation convention $y(x, t) = \bar{y}(x) + y'(x, t)$, where $\bar{y}$ denotes the steady component and $y'$ the small perturbation of $y = [\rho, u, p, T, c, \gamma]$:

$$ 
\left\{
\begin{aligned}
\frac{\partial U}{\partial t} + \bar{u} \frac{\partial U}{\partial x} + \frac{\bar{c}^2}{\bar{u}} \frac{\partial P}{\partial x} \\
+ 2 \frac{d\bar{u}}{dx} U + \left( \frac{\kappa \bar{c}^2}{\bar{u}} \frac{dT}{dx} - \bar{\gamma} \frac{d\bar{u}}{dx} \right) P + \frac{d\bar{u}}{dx} D &= 0, \\
\frac{\partial P}{\partial t} + \bar{u} \frac{\partial}{\partial x} (U + P) + (\bar{\gamma} - 1) \frac{\bar{P}_v}{\bar{\gamma} \bar{p}} U \\
+ \left( \frac{\bar{P}_v}{\bar{p}} \left( (\bar{\gamma} - 1) - \kappa \bar{\gamma} \bar{T} \right) + \kappa \bar{u} \left( \frac{dT}{dx} - \frac{\bar{\gamma} \bar{T}}{\bar{\rho}} \frac{d\bar{\rho}}{dx} \right) \right) P \\
+ \kappa \bar{T} \left( \frac{\bar{P}_v}{\bar{p}} + \frac{\bar{u}}{\bar{\rho}} \frac{d\bar{\rho}}{dx} \right) D &= 0,
\end{aligned}
\right.
\quad (28) $$

where $P = p'/(\bar{\gamma}\bar{p})$, $U = u'/(\bar{u})$, and $D = \rho'/(\bar{\rho})$ are the normalized pressure, velocity, and density fluctuations, respectively.

The steady flow velocity profile, $\bar{u}(x)$, is locally approximated at the throat as $\bar{u}(\tilde{x}) = (d\bar{u}/dx)_* \tilde{x}$, where $(d\bar{u}/dx)_*$ denotes the velocity gradient at the critical-throat. The origin of the coordinate $\tilde{x}$ is chosen such that the extrapolated velocity profile $\bar{u}(\tilde{x})$ vanishes at $\tilde{x} = 0$. The throat position is $\tilde{x}_* = \bar{c}_*/(d\bar{u}/dx)_*$. Based on this definition, the following time-space transformation is introduced as performed by Marble and Candel [8]:

$$ \tau = \frac{\bar{c}_*}{\tilde{x}_*} t \quad ; \quad \xi = \left( \frac{\tilde{x}}{\tilde{x}_*} \right)^2. $$

The momentum and total energy equations of the linearized Euler system (eq. (28)) can be rewritten using the time-space transformation as follows:

$$ 
\left\{
\begin{aligned}
\frac{\partial U}{\partial \tau} + 2\xi \frac{\partial U}{\partial \xi} + 2 \frac{\bar{c}^2}{\bar{c}_*^2} \frac{\partial P}{\partial \xi} \\
+ 2U + \left( \frac{2\kappa \bar{c}^2}{\bar{c}_*^2} \frac{dT}{d\xi} - \bar{\gamma} \right) P + D &= 0, \\
\frac{\partial P}{\partial \tau} + 2\xi \frac{\partial}{\partial \xi} (U + P) + \frac{(\bar{\gamma} - 1) \bar{P}_v}{(d\bar{u}/dx)_* \bar{\gamma} \bar{p}} U \\
+ \left( \frac{\bar{P}_v}{\bar{p}(d\bar{u}/dx)_*} \left( (\bar{\gamma} - 1) - \kappa \bar{\gamma} \bar{T} \right) + 2\kappa \xi \left( \frac{dT}{d\xi} - \frac{\bar{\gamma} \bar{T}}{\bar{\rho}} \frac{d\bar{\rho}}{d\xi} \right) \right) P \\
+ \kappa \bar{T} \left( \frac{\bar{P}_v}{\bar{p}(d\bar{u}/dx)_*} + \frac{2\xi}{\bar{\rho}} \frac{d\bar{\rho}}{d\xi} \right) D &= 0.
\end{aligned}
\right.
\quad (29) $$
Introducing the dimensionless frequency $\Omega = 2\pi f/(d\bar{u}/dx)_*$, the system eq. (29) is further expressed in the frequency domain. Applying the Fourier transform with the convention $y(\xi, \tau) = \hat{y}(\xi) \exp(+j\Omega\tau)$, the above system of equations writes:

$$ egin{cases} 
2\xi \frac{d\hat{U}}{d\xi} + 2\frac{\bar{c}_*^2}{\bar{c}_*^2} \frac{d\hat{P}}{d\xi} \\ 
+ (2 + j\Omega) \hat{U} + \left( \frac{2\kappa\bar{c}_*^2}{\bar{c}_*^2} \frac{d\bar{T}}{d\xi} - \bar{\gamma} \right) \hat{P} + \hat{D} = 0, \\ 
2\xi \frac{d\hat{U}}{d\xi} + 2\xi \frac{d\hat{P}}{d\xi} \\ 
+ \frac{(\bar{\gamma} - 1)}{(d\bar{u}/dx)_*} \frac{\bar{P}_v}{\bar{\gamma}\bar{\rho}\bar{P}} \hat{U} \\ 
+ \left( \frac{\bar{P}_v}{\bar{\rho}(d\bar{u}/dx)_*} \right) \left( (\bar{\gamma} - 1) - \kappa\bar{\gamma}\bar{T} \right) + 2\kappa\xi \left( \frac{d\bar{T}}{d\xi} - \frac{\bar{\gamma}\bar{T}}{\bar{\rho}} \frac{d\bar{\rho}}{d\xi} \right) + j\Omega \right) \hat{P} \\ 
+ \kappa\bar{T} \left( \frac{\bar{P}_v}{\bar{\rho}(d\bar{u}/dx)_*} + \frac{2\xi}{\bar{\rho}} \frac{d\bar{\rho}}{d\xi} \right) \hat{D} = 0. 
\end{cases} 
\quad (30) $$

Subtracting the momentum equation from the energy equation yields:

$$ egin{aligned} & 2(\xi - \frac{\bar{c}_*^2}{\bar{c}_*^2}) \frac{d\hat{P}}{d\xi} - \left( 2 + j\Omega - \frac{(\bar{\gamma} - 1)}{(d\bar{u}/dx)_*} \frac{\bar{P}_v}{\bar{\gamma}\bar{\rho}\bar{P}} \right) \hat{U} \\ & + \left( \bar{\gamma} + j\Omega + \frac{\bar{P}_v}{\bar{\rho}(d\bar{u}/dx)_*} \right) \left( (\bar{\gamma} - 1) - \kappa\bar{\gamma}\bar{T} \right) \hat{P} \\ & - \frac{2\kappa\xi\bar{\gamma}\bar{T}}{\bar{\rho}} \frac{d\bar{\rho}}{d\xi} + 2\kappa(\xi - \frac{\bar{c}_*^2}{\bar{c}_*^2}) \frac{d\bar{T}}{d\xi} \\ & + \left( \kappa\bar{T} \left( \frac{\bar{P}_v}{\bar{\rho}(d\bar{u}/dx)_*} + \frac{2\xi}{\bar{\rho}} \frac{d\bar{\rho}}{d\xi} \right) - 1 \right) \hat{D} = 0. \end{aligned} 
\quad (31) $$

This equation is now evaluated at the steady critical-throat position. The assumption of Tsien [16], as proposed by Crocco, is used such that $\hat{U}$, $\hat{P}$, and $\hat{D}$, as well as their derivatives, are continuous at the critical-throat location. It is further assumed that none of the derivatives are proportional to $1/(1 - \xi)$, as such a condition would lead to a singularity at the critical-throat.

Under these assumptions and because at the throat $\xi = 1$ and $\bar{c}_* = \bar{c}$, hence the term proportional to $(d\hat{P}/d\xi)_*$ vanishes and the boundary condition at the critical-throat position becomes eq. (21).

## B CEDRE solver

The unstructured Computational Fluid Dynamics (CFD) ONERA code CEDRE [37] solves the classical two-dimensional Euler equations:

$$ egin{cases} 
\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) = 0, \\ 
\frac{\partial}{\partial t} (\rho \mathbf{u}) + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) + \nabla p = 0, \\ 
\frac{\partial}{\partial t} (\rho e_0) + \nabla \cdot ((\rho e_0 + p) \mathbf{u}) = P_v, 
\end{cases} 
\quad (32) $$

where $\rho$ denotes the mass density, $\mathbf{u} = (u, v)$ the velocity vector, $p$ the static pressure, $e_0 = e + |\mathbf{u}|^2/2$ the total specific energy ($e$ being the specific internal energy), $P_v$ the volumetric heat source term, $t$ the time, $\nabla$ the Nabla operator, and $\otimes$ denotes the outer product.

The potential influence of computational resolution, both spatial and temporal, has been assessed by varying the computational mesh size and time step.

Four distinct computational grids and four time-step values have been tested. The mesh size $\Delta x$ was varied from 0.5 mm (235 000 grid points) to 4 mm (3000 grid points), while the time step $\Delta t$ ranged from $1 \times 10^{-7}$ s to $5 \times 10^{-6}$ s.

The results show that the pressure as well as the velocity fluctuations (so the acoustic waves) become independent of spatial resolution for mesh sizes smaller than $\Delta x = 1$ mm at a tested frequency $\Omega = 30$. Consequently, the selected mesh size and time step, $(\Delta x, \Delta t) = (1 \text{ mm}, 1 \times 10^{-6} \text{ s})$, corresponding to a CFL number of $(\bar{u}_{in} + \bar{c}_{in})\Delta t/\Delta x = 0.78 < 1$ $(\bar{u}_{in}$ and $\bar{c}_{in}$ being respectively the inlet steady-state longitudinal velocity and speed of sound), offer the best compromise between numerical accuracy and computational cost, yielding an absolute relative error of only $1 \times 10^{-3}$ compared to the finest refinement.

To validate the linear acoustic response of the flow, simulations were performed at the same frequency ($\Omega = 2\pi f/(d\bar{u}/dx)_* = 30$) for acoustic forcing amplitudes of 1% and 0.1%. A relative difference of less than $1 \times 10^{-3}$ was observed in the acoustic-reflection coefficient, confirming the validity of the linear acoustics assumption.

## References

[1] I. Smith, The effect of dissociation in rockets, ramjets, and wind tunnels, Proceedings of the Institution of Mechanical Engineers, Conference Proceedings 178 (1963) 27-32. doi:10.1243/PIME_CONF_1963_178_193_02.

[2] R. Fry, A century of ramjet propulsion technology evolution, J. Propul. Power 20 (2004) 27-58. doi:10.2514/1.9178.

[3] F. Falempin, Ramjet and dual mode operation, Technical Report RTP-EN-AVT-150, NATO, 2008.

[4] H. Wang, Q. Yang, X. Xu, Effect of thermal choking on ejection process in a rocket-based combined cycle engine, App. Therm. Eng. 116 (2017) 197-204. doi:10.1016/j.applthermaleng.2017.01.059.

[5] W. Knuth, P. Gloyer, J. Goodman, R. Litchfor, Thermally-choked combustor technology, Technical 19960009786, NASA, 1993. URL: https://ntrs.nasa.gov/citations/19960009786.

[6] O. Dessornes, D. Scherrer, Tests of the JAPHAR dual-mode ramjet engine, Aerosp. Sci. Technol. 9 (2005) 211-221. doi:10.1016/j.ast.2005.01.007.
[7] K.-C. Lin, K. Jackson, R. Behdadnia, T. Jackson, F. Ma, V. Yang, Acoustic characterization of an ethylene-fueled scramjet combustor with a cavity flameholder, J. Propul. Power 26 (2010) 1161–1170. doi:10.2514/1.43338.

[8] F. Marble, S. Candel, Acoustic disturbance from gas non-uniformities convected through a nozzle, J. Sound Vib. 55 (1977) 225–243. doi:10.1016/0022-460X(77)90596-X.

[9] W. Moase, M. Brear, C. Manzie, The forced response of choked nozzles and supersonic diffusers, J. Fluid Mech. 585 (2007) 281–304. doi:10.1017/S0022112007006647.

[10] I. Duran, S. Moreau, Solution of the quasi-one-dimensional linearized Euler equations using flow invariants and the Magnus expansion, J. Fluid Mech. 723 (2013) 190–231. doi:10.1017/jfm.2013.118.

[11] L. Magri, On indirect noise in multicomponent nozzle flows, J. Fluid Mech. 828 (2017) R2. doi:10.1017/jfm.2017.591.

[12] M. Huet, A. Emmanuelli, T. Le Garrec, Entropy noise modelling in 2D choked nozzle flows, J. Sound Vib. 488 (2020) 115637. doi:10.1016/j.jsv.2020.115637.

[13] Y. Gentil, G. Daviller, S. Moreau, Combustion noise modelling of thermally perfect and multi-species gas flow in nozzles, J. Fluid Mech. 1014 (2025) A39. doi:10.1017/jfm.2025.10330.

[14] F. Olivon, A. Genot, J.-E. Durand, A. Hirschberg, E. Piot, On the linear aero-acoustic response of a thermally-choked-flow nozzle to acoustic and entropy plane waves, Submitted to J. Sound. Vib. (2025). doi:10.2139/ssrn.5297823.

[15] F. Olivon, A. Genot, L. Hirschberg, S. Moreau, A. Hirschberg, On critical-throat boundary-condition in quasi-one-dimensional linearized-Euler equations models, Submitted to J. Fluid Mech. (2025).

[16] H. S. Tsien, The transfer functions of rocket nozzles, J. Am. Rocket Soc. 22 (1952) 139–143. doi:10.2514/8.4448.

[17] N. Karimi, M. Brear, W. Moase, Acoustic and disturbance energy analysis of a flow with heat communication, J. Fluid Mech. 597 (2008) 67–89. doi:10.1017/S0022112007009573.

[18] S. R. Yeddula, R. Gaudron, A. S. Morgans, Acoustic absorption and generation in ducts of smoothly varying area sustaining a mean flow and a mean temperature gradient, J. Sound Vib. 515 (2021) 116437. doi:10.1016/j.jsv.2021.116437.

[19] N. Karimi, M. J. Brear, W. H. Moase, On the interaction of sound with steady heat communicating flows, J. Sound Vib. 329 (2010) 4705–4718. doi:10.1016/j.jsv.2010.05.009.
1098

1099

1100

1101

1102

1103

[31] S. Yeddula, J. Guzman Inigo, A. Morgans, Effect of steady and fluctuating heat transfer on acous- tic and entropy transfer functions of a nozzle, in: 28th AIAA/CEAS Aeroacoustics (2022) Conference, American Institute of Aeronautics and Astronautics, 2022, pp. 1–12. doi:10.2514/6.2022-2978.

1104

1105

1106

1107

1108

[32] Y. Gentil, G. Daviller, S. Moreau, N. Treleaven, T. Poinsot, Theoretical analysis and numerical val- idation of the mechanisms controlling composition noise, J. Sound Vib. 585 (2024) 118463. doi:10.1016/ j.jsv.2024.118463.

1109

1110

1111

1112

1113
1114
1115
1116

[33] E. W. Lemmon, R. T. Jacobsen, S. G. Penoncello, D. G. Friend, Thermodynamic properties of air and mixtures of nitrogen, argon, and oxygen from 60 to 2000 K at pressures to 2000 MPa, J. Phys. Chem. Ref. Dat. 29 (2000) 331–385. doi:10.1063/1.1285884.

[34] E. Süli, D. F. Mayers, An Introduction to Numerical Analysis, 1 ed., Cambridge University Press, 2003. doi:10.1017/CB09780511801181.

1117

1118

[35] A. H. Shapiro, The dynamics and thermodynamics of compressible fluid flow. 1, 23. print ed., Wiley, 1976.

1119

1120

[36] W. Heiser, D. Pratt, Hypersonic Airbreathing Propulsion, AIAA, Washington, D.C, 1994.

1121

1122

1123

1124

1125

[37] A. Refloch, B. Courbet, A. Murrone, P. Villedieu, C. Laurent, P. Gilbank, J. Troyes, L. Tessé, G. Chain- eray, J. Dargaud, E. Quémerais, F. Vuillot, CE- DRE software, AerospaceLab J. (2011). URL: https: //hal.science/hal-01182463.

1126

1127

1128

1129

[38] C. Le Touze, A. Murrone, H. Guillard, Multislope MUSCL method for general unstructured meshes, J. Comp. Phys. 284 (2015) 389–418. doi:10.1016/j. jcp.2014.12.032.

1130

1131

1132

[39] E. Toro, M. Spruce, W. Speares, Restoration of the contact surface in the HLL-Riemann solver, Shock Waves J. 4 (1994) 25–34. doi:10.1007/BF01414629.

1133

1134

1135

[40] J. Butcher, Implicit Runge-Kutta processes, Math. Comp. J. 18 (1964) 50–64. doi:10.1090/ S0025-5718-1964-0159424-9.