© The Author(s), Published by EDP Sciences, 2025

1 Introduction

Ramjets and scramjets are key propulsion systems for high-speed flight due to their higher propulsive efficiency compared to turbojets and turbofans at equivalent velocities. Ramjets operate with subsonic combustion 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 recombination [1]. In contrast, scramjets use supersonic combustion, achieving greater efficiency at Mach numbers 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 supersonic 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 geometrically choked-flow nozzle. Although effective in low-flight Mach number combustion regimes, this configuration becomes unsuitable under scramjet conditions. An alternative is the use of a thermal throat (Fig. 1), produced by the volumetric heat release from combustion [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.

Figure 1. Sketch of an idealized dual-mode ramjet incorporating a thermal throat.

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 commonly 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 conditions for both geometrical and thermal throats. In particular, this quasi-steady critical Mach number assumption 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 behind this theory is to impose a regular solution for the acoustic field at the throat. This was originally proposed by Tsien [16], following a suggestion by Crocco, and later used by Marble and Candel [8] for a particular nozzle geometry. Considering a series expansion of the solution in powers of the distance from the critical throat, Tsien [16] discarded the solution whose leading-order term is inversely proportional to that distance. A generalization of this approach for thermally choked configurations is proposed in the present study.

The effect of heat addition on flow dynamics has been extensively studied in subsonic flows, particularly 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 questionable [26–28]. Olivon et al. [14] numerically investigated the impact of temperature-dependent specific heat capacities and demonstrated that both the acoustic-reflection and entropy-noise-production coefficients can vary by as much as 20%. These findings indicate that the heat capacity temperature dependency 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 central flow, their influence on thermal throat dynamics has been shown to be limited [29]. The volumetric heat source term P _v 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 introducing more complex phenomena. Finally, the analysis is limited to acoustic and entropy waves; indirect noise arising from gas mixture non-uniformity, as investigated by Gentil et al. [13, 32], is not taken into account.

The objective of this study is to develop and validate a boundary condition that describes the behaviour of the critical throat subjected to both acoustic and entropy wave forcing in a thermally choked flow configuration. The quasi-one-dimensional linear 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 discussions follow in Section 5. Conclusions are summarized in Section 6.

2 Quasi-one-dimensional linear acoustic model

The following assumptions define the limits of validity of the present model. First, the flow is modelled as quasi-one-dimensional, and only plane waves are considered. Second, the Euler equations are employed, neglecting viscous and thermal diffusion effects. Third, a non-uniform, steady volumetric heat source term P _v is introduced to model the energy release 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 standard quasi-one-dimensional equations for mass, momentum, and energy conservation are expressed as follows [8]:

$$\begin{array}{c}\hfill \{\begin{array}{cc}& \frac{\mathrm{D}\rho }{\mathrm{D}t}+\rho \frac{\partial u}{\partial x}+\rho u\alpha =0,\hfill \\ \hfill [1ex]& \frac{\mathrm{D}u}{\mathrm{D}t}+\frac{1}{\rho }\frac{\partial p}{\partial x}=0,\hfill \\ \hfill [1ex]& \frac{\mathrm{D}p}{\mathrm{D}t}-c^2\frac{\mathrm{D}\rho }{\mathrm{D}t}=(\gamma -1)P_v,\hfill \end{array}\end{array}$$(1)

where D/Dt = ∂/∂t + u ∂/∂x is the material derivative, ρ is the density, u is the longitudinal velocity, α = (dA/dx)/A is the relative variation of the nozzle cross-sectional area A, p is the static pressure, c is the speed of sound, γ = c _p /c _v is the ratio of specific heats at constant pressure c _p and constant volume c _v , and P _v 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:

$$\begin{array}{c}\hfill p=\rho rT,\end{array}$$(2)

where T is the static temperature, r = ℛ/𝒲 is the specific gas constant, with ℛ = 8.3145 J mol⁻¹K⁻¹ the universal gas constant, and 𝒲 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:

$$\begin{array}{c}\hfill c_p=\sum _{n=0}^7\frac{{\alpha }_n}{\mathcal{W}}{T}^n,\end{array}$$(3)

where α _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:

$$\begin{array}{c}\hfill \frac{\mathrm{d}s}{c_p}=\frac{\mathrm{d}T}{T}-(\gamma -1)\frac{\mathrm{d}p}{\gamma p},\end{array}$$(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)=\overline{y}(x)+y\prime (x,t)$, where $\overline{y}$ denotes the steady component and y′ the small perturbation of y = [ρ, u, p, T, c, γ]. 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 \prime /\overline{\rho }$, $U=u\prime /\overline{u}$, $P=p\prime /(\overline{\gamma }\overline{p})$ yields the following linearized system of differential equations:

$$\begin{array}{c}\hfill \{\begin{array}{cc}& \frac{\partial D}{\partial t}+\overline{u}\frac{\partial }{\partial x}(U+D)=0,\hfill \\ \hfill [1ex]& \frac{\partial U}{\partial t}+\overline{u}\frac{\partial U}{\partial x}+\frac{{\overline{c}}^2}{\overline{u}}\frac{\partial P}{\partial x}\hfill \\ \hfill & +2\frac{\mathrm{d}\overline{u}}{\mathrm{d}x}U+(\frac{\kappa {\overline{c}}^2}{\overline{u}}\frac{\mathrm{d}\overline{T}}{\mathrm{d}x}-\overline{\gamma }\frac{\mathrm{d}\overline{u}}{\mathrm{d}x})P+\frac{\mathrm{d}\overline{u}}{\mathrm{d}x}D=0,\hfill \\ \hfill [1ex]& \frac{\partial P}{\partial t}+\overline{u}\frac{\partial }{\partial x}(U+P)+(\overline{\gamma }-1)\frac{{\overline{P}}_v}{\overline{\gamma }\overline{p}}U\hfill \\ \hfill & +(\frac{{\overline{P}}_v}{\overline{p}}((\overline{\gamma }-1)-\kappa \overline{\gamma }\overline{T})\hfill \\ \hfill & +\kappa \overline{u}(\frac{\mathrm{d}\overline{T}}{\mathrm{d}x}-\frac{\overline{\gamma }\overline{T}}{\overline{\rho }}\frac{\mathrm{d}\overline{\rho }}{\mathrm{d}x}))P\hfill \\ \hfill & +\kappa \overline{T}(\frac{{\overline{P}}_v}{\overline{p}}+\frac{\overline{u}}{\overline{\rho }}\frac{\mathrm{d}\overline{\rho }}{\mathrm{d}x})D=0,\hfill \end{array}\end{array}$$(5)

where $\kappa =(\mathrm{d}\overline{\gamma }/\mathrm{d}T)/\overline{\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 equation (5) is rewritten using the convention $y\prime (x,t)=\widehat{y}(x)exp(+j\omega t)$, where $\widehat{y}$ is the complex amplitude of the Fourier transform of the quantity y, and ω = 2π f is the angular frequency. Introducing the unknown vector $\widehat{X}={(\widehat{U},\widehat{P},\widehat{D})}^T$ leads to the following matrix equation:

$$\begin{array}{c}\hfill \frac{\mathrm{d}\widehat{X}}{\mathrm{d}x}+{\overline{K}}_1^{-1}{\overline{K}}_2\widehat{X}=0,\end{array}$$(6)

with

$$\begin{array}{c}\hfill {\overline{K}}_1^{-1}=\frac{1}{\overline{u}(1-{\overline{\mathrm{M}}}^2)}\left[\begin{array}{ccc}0& -{\overline{\mathrm{M}}}^2& 1\\ 0& {\overline{\mathrm{M}}}^2& -{\overline{\mathrm{M}}}^2\\ 1-{\overline{\mathrm{M}}}^2& {\overline{\mathrm{M}}}^2& -1\end{array}\right]\end{array}$$

where $\overline{\text{M}}=\overline{u}/\overline{c}$ is the steady-flow Mach number, and

$$\begin{array}{cc}& {\overline{K}}_2=\hfill \\ \hfill & \left[\begin{array}{ccc}0& 0& j\omega \\ 2\frac{\mathrm{d}\overline{u}}{\mathrm{d}x}+j\omega & \frac{\kappa {\overline{c}}^2}{\overline{u}}\frac{\mathrm{d}\overline{T}}{\mathrm{d}x}-\overline{\gamma }\frac{\mathrm{d}\overline{u}}{\mathrm{d}x}& \frac{\mathrm{d}\overline{u}}{\mathrm{d}x}\\ \frac{(\overline{\gamma }-1){\overline{P}}_v}{\overline{\gamma }\overline{p}}& {\overline{K}}_2^{32}& \kappa \overline{T}(\frac{{\overline{P}}_v}{\overline{p}}+\frac{\overline{u}}{\overline{\rho }}\frac{\mathrm{d}\overline{\rho }}{\mathrm{d}x})\end{array}\right]\hfill \end{array}$$(7)

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

$$\begin{array}{cc}\hfill {\overline{K}}_2^{32}& =\frac{{\overline{P}}_v}{\overline{p}}((\overline{\gamma }-1)-\kappa \overline{\gamma }\overline{T})\hfill \\ \hfill & +\kappa \overline{u}(\frac{\mathrm{d}\overline{T}}{\mathrm{d}x}-\frac{\overline{\gamma }\overline{T}}{\overline{\rho }}\frac{\mathrm{d}\overline{\rho }}{\mathrm{d}x})+j\omega .\hfill \end{array}$$(8)

It is noteworthy that ${\overline{K}}_1$ is invertible if and only if $\overline{u}\ne 0$ and $\overline{\mathrm{M}}\ne 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 $\overline{C}=-{\overline{K}}_1^{-1}{\overline{K}}_2$ and assuming a discretization of the nozzle over j ∈ ⟦0, …, l − 1, l, l + 1, …, n⟧ elements (the index l denotes the critical-throat position), each of size Δx, and using an explicit fourth-order Runge–Kutta (RK4) numerical scheme [34] for integration, the following relation is obtained:

$$\begin{array}{c}\hfill {\widehat{X}}_{j+1}={\widehat{X}}_j+\frac{\mathrm{\Delta }x}{6}(k_1+2k_2+2k_3+k_4),\end{array}$$(9)

with the intermediate slopes defined as, ∀j ≠ l:

$$\begin{array}{cc}\hfill k_1& ={\overline{C}}_j{\widehat{X}}_j,\hfill \\ \hfill k_2& =\frac{{\overline{C}}_j+{\overline{C}}_{j+1}}{2}({\widehat{X}}_j+\frac{\mathrm{\Delta }x}{2}{k}_1),\hfill \\ \hfill k_3& =\frac{{\overline{C}}_j+{\overline{C}}_{j+1}}{2}({\widehat{X}}_j+\frac{\mathrm{\Delta }x}{2}{k}_2),\hfill \\ \hfill k_4& ={\overline{C}}_{j+1}({\widehat{X}}_j+\mathrm{\Delta }x{k}_3).\hfill \end{array}$$

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 Sect. 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 Sect. 3) and are used to solve the system equation (6). However, although the system is solved for the fluctuations $\widehat{U}(x)$, $\widehat{P}(x)$, and $\widehat{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- (${\widehat{P}}^{+}$) and upstream-propagating (${\widehat{P}}^{-}$) wave components using the following relations:

$$\begin{array}{c}\hfill {\widehat{P}}^{\pm }=\frac{1}{2}(\widehat{P}\pm \overline{\mathrm{M}}\widehat{U}).\end{array}$$(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 Figure 2.

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.

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

$$\begin{array}{c}\hfill \widehat{\sigma }=\frac{\widehat{s}}{{\overline{c}}_p}=\widehat{P}-\widehat{D}.\end{array}$$(11)

Since entropy disturbances are convected with the mean flow $\overline{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”):

$$\begin{array}{c}\hfill R_a=\frac{{\widehat{P}}_{\mathrm{in}}^{-}}{{\widehat{P}}_{\mathrm{in}}^{+}}·\end{array}$$(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:

$$\begin{array}{c}\hfill T_a=\frac{{\widehat{P}}_{\mathrm{out}}^{+}}{{\widehat{P}}_{\mathrm{in}}^{+}}·\end{array}$$(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 (${\widehat{P}}_{\mathrm{in}}^{+}=0$):

$$\begin{array}{c}\hfill R_s=\frac{{\widehat{P}}_{\mathrm{in}}^{-}}{{\widehat{\sigma }}_{\mathrm{in}}},\end{array}$$(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:

$$\begin{array}{c}\hfill T_s=\frac{{\widehat{P}}_{\mathrm{out}}^{+}}{{\widehat{\sigma }}_{\mathrm{in}}}·\end{array}$$(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 equation (6), three boundary conditions must be applied at the system boundaries. Following the methodology of Duran and Moreau [10], the unknown vector $\widehat{X}(x)$ is computed in two stages (see Fig. 2).

In the subsonic region, two boundary conditions must be applied at the inlet (x _in), and one at the upstream side of the critical-throat (denoted x _cu).

To ensure that the system equation (6) remains invertible ($\overline{\text{M}}\ne 1$), the locations where boundary conditions are applied near the critical throat are shifted slightly away from the sonic point by an infinitesimal distance ϵ, as illustrated in Figure 2, leading to x _cu = x _* − ϵ and x _cd = x _* + ϵ. This distance ϵ is set equal to the spatial discretization Δ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 ${\widehat{X}}_{\mathrm{in}}$. For instance, ${\widehat{P}}_{\mathrm{in}}$ and ${\widehat{D}}_{\mathrm{in}}$ are prescribed as input values, while the third component $\widehat{U}$ is initially guessed, forming the inlet unknown vector ${\widehat{X}}_{\mathrm{in},g}={({\widehat{U}}_{\mathrm{guess}},{\widehat{P}}_{\mathrm{in}},{\widehat{D}}_{\mathrm{in}})}^T$. The linear acoustic model equation (9) is then used to compute the spatial profiles of $\widehat{P}(x)$, $\widehat{D}(x)$, and $\widehat{U}(x)$.

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

$$\begin{array}{c}\hfill {\widehat{X}}_{\mathrm{in},g+1}={\widehat{X}}_{\mathrm{in},g}-\frac{h\left({\widehat{X}}_{\mathrm{in},g}\right)}{h^{\prime }\left({\widehat{X}}_{\mathrm{in},g}\right)},\end{array}$$(16)

where the derivative $h\prime ({\widehat{X}}_{\mathrm{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 × 10⁻⁸. 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 ${\widehat{X}}_{\mathrm{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 ${\overline{x}}_{\mathrm{out}}$. The resulting solution satisfies the continuity condition $\widehat{X}(x_{\mathrm{cd}})=\widehat{X}(x_{\mathrm{cu}})$.

3.1 Inlet boundary conditions

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

$$\begin{array}{c}\hfill \{\begin{array}{cc}\hfill {\widehat{P}}_{\mathrm{in}}& =\eta ,\hfill \\ \hfill {\widehat{D}}_{\mathrm{in}}& ={\widehat{P}}_{\mathrm{in}}.\hfill \end{array}\end{array}$$(17)

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

$$\begin{array}{c}\hfill \{\begin{array}{cc}\hfill {\widehat{D}}_{\mathrm{in}}& ={\widehat{P}}_{\mathrm{in}}-\eta ,\hfill \\ \hfill {\widehat{U}}_{\mathrm{in}}& =-{\widehat{P}}_{\mathrm{in}}/{\overline{\mathrm{M}}}_{\mathrm{in}}.\hfill \end{array}\end{array}$$(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 ()_* hereafter denotes quantities evaluated at the steady critical-throat location. For isentropic gas flows with constant heat capacity ratio γ, this quasi-steady critical-throat boundary condition (denoted QSS) can be expressed as:

$$\begin{array}{c}\hfill 2{\widehat{U}}_{\ast }=\gamma {\widehat{P}}_{\ast }-{\widehat{D}}_{\ast }.\end{array}$$(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:

$$\begin{array}{c}\hfill 2{\widehat{U}}_{\ast }=(\gamma -1){\widehat{P}}_{\ast }+{\widehat{\sigma }}_{\ast }.\end{array}$$(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, $\overline{u}(x)$, is locally approximated at the throat as $\overline{u}(\tilde{x})={(\mathrm{d}\overline{u}/\mathrm{d}x)}_{\ast }\tilde{x}$, where ${(\mathrm{d}\overline{u}/\mathrm{d}x)}_{\ast }$ denotes the velocity gradient at the critical-throat. The origin of the coordinate $\tilde{x}$ is chosen such that the extrapolated velocity profile $\overline{u}(\tilde{x})$ vanishes at $\tilde{x}=0$. The throat position is ${\tilde{x}}_{\ast }={\overline{c}}_{\ast }/{(\mathrm{d}\overline{u}/\mathrm{d}x)}_{\ast }$. Based on this definition, the following time-space transformation is introduced as performed by Marble and Candel [8]:

$$\begin{array}{cc}\hfill \tau & =\frac{{\overline{c}}_{\ast }}{{\tilde{x}}_{\ast }}t;\xi ={\left(\frac{\tilde{x}}{{\tilde{x}}_{\ast }}\right)}^2.\hfill \end{array}$$

This boundary condition is based on the linearized quasi-one-dimensional Euler equations (5), rewritten using the above time-space transformation. Introducing the dimensionless frequency $\mathrm{\Omega }=2\pi f/{(\mathrm{d}\overline{u}/\mathrm{d}x)}_{\ast }$, the system equation (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 $\widehat{U}$, $\widehat{P}$, and $\widehat{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 − ξ). Such a condition would lead to a singularity at the critical throat. Under these assumptions, and because at the throat ξ = 1, ${\overline{c}}_{\ast }=\overline{c}$ and terms in ξ − 1 vanish, the generalized boundary condition at the critical-throat position (denoted GEN) becomes:

$$\begin{array}{c}\hfill {\overline{g}}_U{\widehat{U}}_{\ast }={\overline{g}}_P{\widehat{P}}_{\ast }+{\overline{g}}_D{\widehat{D}}_{\ast },\end{array}$$(21)

where the coefficients are defined as:

$$\begin{array}{c}\hfill \{\begin{array}{cc}\hfill {\overline{g}}_U& =2+j\mathrm{\Omega }-\frac{({\overline{\gamma }}_{\ast }-1)}{{(\mathrm{d}\overline{u}/\mathrm{d}x)}_{\ast }}\frac{{\overline{P}}_{v,\ast }}{{\overline{\gamma }}_{\ast }{\overline{p}}_{\ast }},\hfill \\ \hfill [1ex]{\overline{g}}_P& ={\overline{\gamma }}_{\ast }+j\mathrm{\Omega }+\frac{{\overline{P}}_{v,\ast }}{{\overline{p}}_{\ast }{(\mathrm{d}\overline{u}/\mathrm{d}x)}_{\ast }}\hfill \\ \hfill & [({\overline{\gamma }}_{\ast }-1)-{\kappa }_{\ast }{\overline{\gamma }}_{\ast }{\overline{T}}_{\ast }]-\frac{2{\kappa }_{\ast }{\overline{\gamma }}_{\ast }{\overline{T}}_{\ast }}{{\overline{\rho }}_{\ast }}{\frac{\mathrm{d}\overline{\rho }}{\mathrm{d}\xi }|}_{\ast },\hfill \\ \hfill {\overline{g}}_D& ={\kappa }_{\ast }{\overline{T}}_{\ast }(\frac{{\overline{P}}_{v,\ast }}{{\overline{p}}_{\ast }{(\mathrm{d}\overline{u}/\mathrm{d}x)}_{\ast }}+\frac{2}{{\overline{\rho }}_{\ast }}{\frac{\mathrm{d}\overline{\rho }}{\mathrm{d}\xi }|}_{\ast })-1.\hfill \end{array}\end{array}$$(22)

This critical-throat boundary condition is an extension 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 condition is also of interest. At very high frequencies (Ω ≫ 1), the relation ${\widehat{U}}_{\ast }={\widehat{P}}_{\ast }$ emerges, indicating that in a heated medium, the generation of entropy due to propagating acoustic waves becomes negligible. This is because ${\widehat{D}}_{\ast }={\widehat{P}}_{\ast }-{\widehat{\sigma }}_{\ast }$ becomes negligible compared to ${\widehat{P}}_{\ast }$. 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 x _cu 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:

$$\begin{array}{cc}\hfill \text{QSS:}& 2{\widehat{U}}_{\ast }-\gamma {\widehat{P}}_{\ast }+{\widehat{D}}_{\ast }=0,\hfill \end{array}$$(23)

$$\begin{array}{cc}\hfill \text{GEN:}& {\overline{g}}_U{\widehat{U}}_{\ast }-{\overline{g}}_P{\widehat{P}}_{\ast }-{\overline{g}}_D{\widehat{D}}_{\ast }=0.\hfill \end{array}$$(24)

4 Validation methodology

To validate the GEN boundary condition proposed in equation (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 simulations 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 configuration and a thermally-choked-flow configuration. The results comparisons are performed on the acoustic coefficients (R _a , R _s , T _a and T _s ) previously defined in Section 2.3.

In the following, only the main steps of the methodology are outlined.

4.1 Thermally choked flow nozzle

For the thermally-choked flow configuration, a simplified dual-mode ramjet geometry, based on the works of [6, 14, 26, 28], is used, as shown in Figure 3. The nozzle geometry consists of three sections: a straight-duct upstream section (0 ≤ x/L <  x_u/L = 0.42), a smooth transition section (x_u/L = 0.42 ≤ x/L <  x_d/L = 0.58) formed by a wall in circular arc with angle α = 3°, and a divergent duct with a constant divergence angle α = 3° (x/L >  x_d/L = 0.58). The total length and inlet height of the nozzle are L = 2.4 m and y_in = 25 × 10⁻³ m, respectively. The circular arc radius of the transition region is R_t/y_in = 305. A slip-adiabatic boundary condition is imposed on the upper nozzle wall, while a symmetry boundary condition is applied along the lower edge.

Figure 3. Simplified two-dimensional planar sketch of a dual-mode ramjet geometry (not to scale).

The incoming airflow at the inlet is subsonic, as detailed in Table 1. The inlet conditions are representative 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 _v (x) is divided into two regions. In the straight duct (Fig. 3), the flow is not heated, P _v is set to zero. In the divergent section (x >  x _u), the volumetric power distribution P _v follows a log-normal profile along x, as proposed by Wang et al. [4] and Olivon et al. [14].

To determine P _v (x), a quasi-one-dimensional steady-state model is used, based on the Mach number 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 _v (x). It determines the gradient of the heat release per unit mass, dq/dx:

$$\begin{array}{c}\hfill P_v(x)=\frac{\dot{m}}{\mathrm{A}(x)}\frac{\mathrm{d}q}{\mathrm{d}x},\end{array}$$(25)

where $\dot{m}=\overline{\rho }\overline{u}\text{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 = δ _*/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 γ _in:

$$\begin{array}{c}\hfill {\mathrm{A}}_{\mathrm{is}}(x)=\frac{{\mathrm{A}}_{\mathrm{in}}{\mathrm{M}}_{\mathrm{in}}}{{\mathrm{M}}_{\mathrm{tht}}}{\left(\frac{2+({\gamma }_{\mathrm{in}}-1){\mathrm{M}}_{\mathrm{tht}}^2}{2+({\gamma }_{\mathrm{in}}-1){\mathrm{M}}_{\mathrm{in}}^2}\right)}^{\frac{{\gamma }_{\mathrm{in}}+1}{2({\gamma }_{\mathrm{in}}-1)}}.\end{array}$$(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 γ = γ _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 × 10⁻⁴. 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 Δx/L = 4.1 × 10⁻⁴. This resolution ensures adequate spatial discretization of acoustic waves, providing at least 100 points per wavelength for frequencies up to $\mathrm{\Omega }=2\pi f/{(\mathrm{d}\overline{u}/\mathrm{d}x)}_{\ast }=30$. Since the flow is primarily longitudinal, the same discretization (Δx ≃ Δ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 Δt = 1 × 10⁻⁶ s, ensuring a Courant–Friedrichs–Lewy (CFL) number of $({\overline{u}}_{\mathrm{in}}+{\overline{c}}_{\mathrm{in}})\mathrm{\Delta }t/\mathrm{\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 Sect. 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 Figure 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 equation (21).

Figure 4. 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 ($\overline{\text{M}}$, (a)); the static pressure ($\overline{p}$, (b)) and the static temperature ($\overline{T}$, (c)) are displayed. Quantities are measured on the symmetry axis of the computational domain.

For both thermally and geometrically choked flow configurations, the steady Mach number profiles are similar, with a maximum relative deviation of approximately 5 × 10⁻⁴ 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 _* ≈ 1740 K at the throat. This significant temperature increase induces a variation in the heat capacity ratio, from γ _in = 1.379 to γ _* = 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)).

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 equation (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, ${\widehat{\text{M}}}_{\ast }/{\widehat{P}}_{\ast }$, resulting from the isentropic acoustic forcing defined in equation (17), plotted as a function of the dimensionless frequency $\mathrm{\Omega }=2\pi f/{(\mathrm{d}\overline{u}/\mathrm{d}x)}_{\ast }$. 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 equations (21) and (22), the analytical expression for the normalised critical-throat Mach number fluctuations is given by:

Figure 5. Modulus (a) and phase angle (b) of the normalized critical Mach number fluctuations due to isentropic acoustic forcing at the inlet ((Eq. 17)), ${\widehat{\text{M}}}_{\ast }/{\widehat{P}}_{\ast }$, plotted against the dimensionless frequency $\mathrm{\Omega }=2\pi f/{(\mathrm{d}\overline{u}/\mathrm{d}x)}_{\ast }$. 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.

$$\begin{array}{cc}\hfill \frac{{\widehat{\mathrm{M}}}_{\ast }}{{\widehat{P}}_{\ast }}& =[\frac{{\overline{g}}_P}{{\overline{g}}_U}-{\overline{\gamma }}_{\ast }\frac{(1+{\kappa }_{\ast }{\overline{T}}_{\ast })}{2}]\hfill \\ \hfill & +[\frac{{\overline{g}}_D}{{\overline{g}}_U}-\frac{(1+{\kappa }_{\ast }{\overline{T}}_{\ast })}{2}]\frac{{\widehat{D}}_{\ast }}{{\widehat{P}}_{\ast }}·\hfill \end{array}$$(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-{({\widehat{\text{M}}}_{\ast }/{\widehat{P}}_{\ast })}^{\mathrm{L}}/{({\widehat{\text{M}}}_{\ast }/{\widehat{P}}_{\ast })}^{\mathrm{CFD}}|$. The forcing amplitude is set to $|{\widehat{P}}_{\mathrm{in}}|=\eta ={10}^{-2}$, ensuring that non-linear effects remain negligible. To verify the linearity of the response, three frequencies Ω = [3, 8.5, 22] were tested with two isentropic forcing amplitudes (|P′_in|=1 × 10⁻² and 1 × 10⁻³). The relative deviations between the responses were on the order of 1 × 10⁻³, confirming linear behaviour.

For the geometric-throat configuration, the QSS boundary condition ${\widehat{\text{M}}}_{\ast }=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 Ω <  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 $\underset{\mathrm{\Omega }\to 0}{lim}|{\widehat{\text{M}}}_{\ast }/{\widehat{P}}_{\ast }|\simeq 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 ${\widehat{\text{M}}}_{\ast }/{\widehat{P}}_{\ast }$ peaks. This maximum, occurring at Ω = 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 ${\widehat{D}}_{\ast }/{\widehat{P}}_{\ast }=1-{\widehat{\sigma }}_{\ast }/{\widehat{P}}_{\ast }=1$ across all frequencies, resulting in the absence of a local maximum, as shown in Figure 5. In contrast, for a heated flow, the ratio ${\widehat{\text{M}}}_{\ast }/{\widehat{P}}_{\ast }$ becomes strongly dependent on the entropy generation term, ${\widehat{D}}_{\ast }/{\widehat{P}}_{\ast }=1-{\widehat{\sigma }}_{\ast }/{\widehat{P}}_{\ast }$, since the density fluctuations ${\widehat{D}}_{\ast }$ are directly linked to entropy fluctuations via equation (11). At low frequencies (0 <  Ω <  3.1), the entropy waves generated by the propagation of acoustic waves induce additional pressure fluctuations, thereby amplifying the Mach number fluctuations at the critical throat. At moderate frequencies (3.1 <  Ω <  10), the generation of entropy waves induced by acoustic forcing decreases, leading to reduced Mach number fluctuations at the throat. At high frequencies, Ω = 𝒪(10), the ratio ${\widehat{D}}_{\ast }/{\widehat{P}}_{\ast }$ 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 equation (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, ${\widehat{\text{M}}}_{\ast }/{\widehat{\sigma }}_{\ast }$, are plotted as a function of the dimensionless frequency $\mathrm{\Omega }=2\pi f/{(\mathrm{d}\overline{u}/\mathrm{d}x)}_{\ast }$.

Figure 6. Modulus (a) and phase angle (b) of the normalized critical Mach number fluctuations, ${\widehat{\text{M}}}_{\ast }/{\widehat{\sigma }}_{\ast }$, resulting from entropy forcing at the inlet ((Eq. 18)), plotted against the dimensionless frequency $\mathrm{\Omega }=2\pi f/{(\mathrm{d}\overline{u}/\mathrm{d}x)}_{\ast }$. 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.

The relative deviation between the quasi-one-dimensional linear acoustic model (L) and the CFD simulations is quantified as $|1-{({\widehat{\text{M}}}_{\ast }/{\widehat{\sigma }}_{\ast })}^{\mathrm{L}}/{({\widehat{\text{M}}}_{\ast }/{\widehat{\sigma }}_{\ast })}^{\mathrm{CFD}}|$. The amplitude of the incoming entropy wave, ${\widehat{\sigma }}_{\mathrm{in}}$, is set to the order of 10⁻², 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 Ω >  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 ${\widehat{\text{M}}}_{\ast }=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={\widehat{P}}_{\mathrm{in}}^{-}/{\widehat{P}}_{\mathrm{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 $\mathrm{\Omega }=2\pi f/{(\mathrm{d}\overline{u}/\mathrm{d}x)}_{\ast }$.

Figure 7. Modulus (a) and phase angle (b) of the acoustic reflection coefficient, $R_a={\widehat{P}}_{\mathrm{in}}^{-}/{\widehat{P}}_{\mathrm{in}}^{+}$, plotted as a function of the dimensionless frequency $\mathrm{\Omega }=2\pi f/{(\mathrm{d}\overline{u}/\mathrm{d}x)}_{\ast }$. 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.

For both choking configurations, a fair agreement is observed between the quasi-one-dimensional linear acoustic model using the GEN boundary condition (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 (Ω <  10) when using the GEN boundary condition. At higher frequencies (Ω >  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 (Ω <  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={\widehat{P}}_{\mathrm{in}}^{-}/{\widehat{\sigma }}_{\mathrm{in}}$, is plotted in Figure 8 as a function of the dimensionless frequency $\mathrm{\Omega }=2\pi f/{(\mathrm{d}\overline{u}/\mathrm{d}x)}_{\ast }$. This coefficient represents the ratio between the upstream-propagating acoustic wave ${\widehat{P}}_{\mathrm{in}}^{-}$ and the incoming entropy wave ${\widehat{\sigma }}_{\mathrm{in}}$.

Figure 8. Modulus (a) and phase angle (b) of the entropy-noise production coefficient, $R_s={\widehat{P}}_{\mathrm{in}}^{-}/{\widehat{\sigma }}_{\mathrm{in}}$, plotted as a function of the dimensionless frequency $\mathrm{\Omega }=2\pi f/{(\mathrm{d}\overline{u}/\mathrm{d}x)}_{\ast }$. 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.

Fair agreement is observed between the quasi-one-dimensional linear acoustic model using the GEN boundary condition (…. and - - -) and CFD solver (⃞. and ⃘) 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 Ω = 5. In all cases, the deviation remains below 25% for the thermally-choked-flow configuration (⃘) 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 Ω <  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 condition (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={\widehat{P}}_{\mathrm{out}}^{+}/{\widehat{P}}_{\mathrm{in}}^{+}$, as a function of the dimensionless frequency $\mathrm{\Omega }=2\pi f/{(\mathrm{d}\overline{u}/\mathrm{d}x)}_{\ast }$ for both configurations, as detailed in Figure 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).

Figure 9. Modulus (a) and phase angle (b) of the acoustic-transmission coefficient $T_a={\widehat{P}}_{\mathrm{out}}^{+}/{\widehat{P}}_{\mathrm{in}}^{+}$, plotted against the dimensionless frequency $\mathrm{\Omega }=2\pi f/{(\mathrm{d}\overline{u}/\mathrm{d}x)}_{\ast }$. 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.

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: ⃘, geometric throat: ⃞) across the entire frequency range considered (Ω ≤ 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 (Ω <  5), a significant portion of the acoustic wave energy is converted into entropy waves. At higher frequencies (Ω >  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 Ω >  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={\widehat{P}}_{\mathrm{out}}^{+}/{\widehat{\sigma }}_{\mathrm{in}}$, plotted as a function of the dimensionless frequency $\mathrm{\Omega }=2\pi f/{(\mathrm{d}\overline{u}/\mathrm{d}x)}_{\ast }$, 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.

Figure 10. Modulus (a) and phase angle (b) of the entropy-noise-transmission coefficient $T_s={\widehat{P}}_{\mathrm{out}}^{+}/{\widehat{\sigma }}_{\mathrm{in}}$, plotted against the dimensionless frequency $\mathrm{\Omega }=2\pi f/{(\mathrm{d}\overline{u}/\mathrm{d}x)}_{\ast }$. 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.

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 (${\overline{\text{M}}}_{\ast }=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.

Conflicts of interest

The authors declare no conflicts of interest in regards to this article.

Data availability statement

The research data associated with this article are included within the article.

References

Appendix 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)=\overline{y}(x)+y\prime (x,t)$, where $\overline{y}$ denotes the steady component and y′ the small perturbation of y = [ρ, u, p, T, c, γ]:

$$\begin{array}{c}\hfill \{\begin{array}{cc}& \frac{\partial U}{\partial t}+\overline{u}\frac{\partial U}{\partial x}+\frac{{\overline{c}}^2}{\overline{u}}\frac{\partial P}{\partial x}\hfill \\ \hfill & +2\frac{\mathrm{d}\overline{u}}{\mathrm{d}x}U+(\frac{\kappa {\overline{c}}^2}{\overline{u}}\frac{\mathrm{d}\overline{T}}{\mathrm{d}x}-\overline{\gamma }\frac{\mathrm{d}\overline{u}}{\mathrm{d}x})P\hfill \\ \hfill & +\frac{\mathrm{d}\overline{u}}{\mathrm{d}x}D=0,\hfill \\ \hfill [1ex]& \frac{\partial P}{\partial t}+\overline{u}\frac{\partial }{\partial x}(U+P)+(\overline{\gamma }-1)\frac{{\overline{P}}_v}{\overline{\gamma }\overline{p}}U\hfill \\ \hfill & +(\frac{{\overline{P}}_v}{\overline{p}}((\overline{\gamma }-1)-\kappa \overline{\gamma }\overline{T})\hfill \\ \hfill & +\kappa \overline{u}(\frac{\mathrm{d}\overline{T}}{\mathrm{d}x}-\frac{\overline{\gamma }\overline{T}}{\overline{\rho }}\frac{\mathrm{d}\overline{\rho }}{\mathrm{d}x}))P\hfill \\ \hfill & +\kappa \overline{T}(\frac{{\overline{P}}_v}{\overline{p}}+\frac{\overline{u}}{\overline{\rho }}\frac{\mathrm{d}\overline{\rho }}{\mathrm{d}x})D=0,\hfill \end{array}\end{array}$$(A.1)

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

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

$$\begin{array}{cc}\hfill \tau & =\frac{{\overline{c}}_{\ast }}{{\tilde{x}}_{\ast }}t;\xi ={\left(\frac{\tilde{x}}{{\tilde{x}}_{\ast }}\right)}^2.\hfill \end{array}$$

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

$$\begin{array}{c}\hfill \{\begin{array}{cc}& \frac{\partial U}{\partial \tau }+2\xi \frac{\partial U}{\partial \xi }+2\frac{{\overline{c}}^2}{{\overline{c}}_{\ast }^2}\frac{\partial P}{\partial \xi }\hfill \\ \hfill & +2U+(\frac{2\kappa {\overline{c}}^2}{{\overline{c}}_{\ast }^2}\frac{\mathrm{d}\overline{T}}{\mathrm{d}\xi }-\overline{\gamma })P+D=0,\hfill \\ \hfill [1ex]& \frac{\partial P}{\partial \tau }+2\xi \frac{\partial }{\partial \xi }(U+P)+\frac{(\overline{\gamma }-1)}{{(\mathrm{d}\overline{u}/\mathrm{d}x)}_{\ast }}\frac{{\overline{P}}_v}{\overline{\gamma }\overline{p}}U\hfill \\ \hfill & +(\frac{{\overline{P}}_v}{\overline{p}{(\mathrm{d}\overline{u}/\mathrm{d}x)}_{\ast }}((\overline{\gamma }-1)-\kappa \overline{\gamma }\overline{T})\hfill \\ \hfill & +2\kappa \xi (\frac{\mathrm{d}\overline{T}}{\mathrm{d}\xi }-\frac{\overline{\gamma }\overline{T}}{\overline{\rho }}\frac{\mathrm{d}\overline{\rho }}{\mathrm{d}\xi }))P\hfill \\ \hfill & +\kappa \overline{T}(\frac{{\overline{P}}_v}{\overline{p}{(\mathrm{d}\overline{u}/\mathrm{d}x)}_{\ast }}+\frac{2\xi }{\overline{\rho }}\frac{\mathrm{d}\overline{\rho }}{\mathrm{d}\xi })D=0.\hfill \end{array}\end{array}$$(A.2)

Introducing the dimensionless frequency $\mathrm{\Omega }=2\pi f/{(\mathrm{d}\overline{u}/\mathrm{d}x)}_{\ast }$, the system equation (A.2) is further expressed in the frequency domain. Applying the Fourier transform with the convention $y(\xi ,\tau )=\widehat{y}(\xi )exp(+j\mathrm{\Omega }\tau )$, the above system of equations writes:

$$\begin{array}{c}\hfill \{\begin{array}{cc}& 2\xi \frac{\mathrm{d}\widehat{U}}{\mathrm{d}\xi }+2\frac{{\overline{c}}^2}{{\overline{c}}_{\ast }^2}\frac{\mathrm{d}\widehat{P}}{\mathrm{d}\xi }+(2+j\mathrm{\Omega })\widehat{U}\hfill \\ \hfill & +(\frac{2\kappa {\overline{c}}^2}{{\overline{c}}_{\ast }^2}\frac{\mathrm{d}\overline{T}}{\mathrm{d}\xi }-\overline{\gamma })\widehat{P}+\widehat{D}=0,\hfill \\ \hfill [1ex]& 2\xi \frac{\mathrm{d}\widehat{U}}{\mathrm{d}\xi }+2\xi \frac{\mathrm{d}\widehat{P}}{\mathrm{d}\xi }+\frac{(\overline{\gamma }-1)}{{(\mathrm{d}\overline{u}/\mathrm{d}x)}_{\ast }}\frac{{\overline{P}}_v}{\overline{\gamma }\overline{p}}\widehat{U}\hfill \\ \hfill & +(\frac{{\overline{P}}_v}{\overline{p}{(\mathrm{d}\overline{u}/\mathrm{d}x)}_{\ast }}((\overline{\gamma }-1)-\kappa \overline{\gamma }\overline{T})\hfill \\ \hfill & +2\kappa \xi (\frac{\mathrm{d}\overline{T}}{\mathrm{d}\xi }-\frac{\overline{\gamma }\overline{T}}{\overline{\rho }}\frac{\mathrm{d}\overline{\rho }}{\mathrm{d}\xi })+j\mathrm{\Omega })\widehat{P}\hfill \\ \hfill & +\kappa \overline{T}(\frac{{\overline{P}}_v}{\overline{p}{(\mathrm{d}\overline{u}/\mathrm{d}x)}_{\ast }}+\frac{2\xi }{\overline{\rho }}\frac{\mathrm{d}\overline{\rho }}{\mathrm{d}\xi })\widehat{D}=0.\hfill \end{array}\end{array}$$(A.3)

Subtracting the momentum equation from the energy equation yields:

$$\begin{array}{c}\hfill \begin{array}{cc}& 2(\xi -\frac{{\overline{c}}^2}{{\overline{c}}_{\ast }^2})\frac{\mathrm{d}\widehat{P}}{\mathrm{d}\xi }-(2+j\mathrm{\Omega }-\frac{(\overline{\gamma }-1)}{{(\mathrm{d}\overline{u}/\mathrm{d}x)}_{\ast }}\frac{{\overline{P}}_v}{\overline{\gamma }\overline{p}})\widehat{U}\hfill \\ \hfill & +\left(\begin{array}{cc}& \overline{\gamma }+j\mathrm{\Omega }+\frac{{\overline{P}}_v}{\overline{p}{(\mathrm{d}\overline{u}/\mathrm{d}x)}_{\ast }}((\overline{\gamma }-1)-\kappa \overline{\gamma }\overline{T})\hfill \\ \hfill & -\frac{2\kappa \xi \overline{\gamma }\overline{T}}{\overline{\rho }}\frac{\mathrm{d}\overline{\rho }}{\mathrm{d}\xi }+2\kappa (\xi -\frac{{\overline{c}}^2}{{\overline{c}}_{\ast }^2})\frac{\mathrm{d}\overline{T}}{\mathrm{d}\xi }\hfill \end{array}\right)\widehat{P}\hfill \\ \hfill & +(\kappa \overline{T}(\frac{{\overline{P}}_v}{\overline{p}{(\mathrm{d}\overline{u}/\mathrm{d}x)}_{\ast }}+\frac{2\xi }{\overline{\rho }}\frac{\mathrm{d}\overline{\rho }}{\mathrm{d}\xi })-1)\widehat{D}=0.\hfill \end{array}\end{array}$$(A.4)

This equation is now evaluated at the steady critical-throat position. The assumption of Tsien [16], as proposed by Crocco, is used such that $\widehat{U}$, $\widehat{P}$, and $\widehat{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 − ξ), as such a condition would lead to a singularity at the critical-throat.

Under these assumptions and because at the throat ξ = 1 and ${\overline{c}}_{\ast }=\overline{c}$, hence the term proportional to ${(\mathrm{d}\widehat{P}/\mathrm{d}\xi )}_{\ast }$ vanishes and the boundary condition at the critical-throat position becomes equation (21).

Appendix B

CEDRE solver

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

$$\begin{array}{c}\hfill \{\begin{array}{cc}& \frac{\partial \rho }{\partial t}+\mathrm{\nabla }·\left(\rho \mathbf{u}\right)=0,\hfill \\ \hfill & \frac{\partial }{\partial t}\left(\rho \mathbf{u}\right)+\mathrm{\nabla }·(\rho \mathbf{u}\otimes \mathbf{u})+\mathrm{\nabla }p=0,\hfill \\ \hfill & \frac{\partial }{\partial t}\left(\rho e_0\right)+\mathrm{\nabla }·\left((\rho e_0+p)\mathbf{u}\right)=P_v,\hfill \end{array}\end{array}$$(B.1)