2.1. Common general equations

The 1-D modelling approach considered in this work treats the primary and secondary streams as two parallel quasi-1-D domains along the axial coordinate $x$ , in an axisymmetric geometry with radial coordinate as shown schematically in figure 1. The domains span the full ejector: single-stream inlets merge at $x=0$ and continue as a coupled, double quasi-1-D mixing pipe. Cross-sections $A_i(x)$ , with $i \in [p,s]$ , vary smoothly along $x$ and reveal a local discontinuity for the secondary at $x=0$ , where inlet inclination and finite primary-nozzle thickness are treated as in § 3.2.

Figure 1.

The ejector is modelled with two 1-D domains that are coupled in the mixing pipe. The geometry is assumed to be axisymmetric. The boundary conditions consist of the total pressures $p_{t,i}$ and total temperatures $T_{t,i}$ at the inlets, and the static back pressure $p_b$ . The dividing streamline, indicated by the dotted line, is computed with (2.9) from the compound choking theory or imposed from the RANS simulations. Momentum exchange is accounted for through friction forces (see (2.17)) or through imposed total pressure gradients from RANS. The schematics shown here are intended for illustration purposes only and do not represent the geometry used in the RANS simulations (see figure 6).

Figure 1. Long description

A diagram illustrating different models of ejector performance under choked conditions. Panel A: Model 1 shows compound choking with calibrated closure. The diagram features a cross-sectional view of an ejector with labeled parameters such as total pressure and temperature at the primary and secondary inlets, static back pressure, and the dividing streamline. The relationship between these parameters is depicted through various lines and annotations, including the force Fi and the area derivative dAi/dx. Panel B: Model 2 illustrates compound choking with closure from computational fluid dynamics (CFD). Similar to Model 1, it shows the ejector’s cross-sectional view with labeled parameters and annotations. The total pressure gradient dp_ti/dx is derived from RANS simulations. Panel C: Model 3 depicts Fabri choking with closure from CFD. This model includes the Mach number Mas equaling 1 and the density ratio rd(x) derived from RANS simulations. The total pressure gradient dp_ti/dx is also derived from RANS simulations. Panel D: Model 4 shows Fabri choking with closure from compound choking theory. It includes the Mach number Mas equaling 1 and the force Fi. The area derivative dAi/dx is given by equation 2.9.

The framework describes the axial evolution of cross-section-averaged flow quantities in each stream. Radial flow variations are averaged out, with their effect represented through variable cross-sectional areas and exchange terms at the stream interfaces and the wall. The governing equations are the conservation laws of mass, momentum and energy in variable-area ducts, with ideal-gas behaviour and no inter-stream mass transfer – hence, the interface in figure 1 is a dividing streamline. Written in terms of static pressure $p$ , the total pressure $p_t$ and the total temperature $T_t$ , these read

(2.1) \begin{align} \frac {1}{p_i} \frac {\text{d}p_i}{\text{d}x} = &\left [\frac {\gamma {Ma}_i^2}{1-{Ma}_i^2}\right ]\frac {1}{A_i}\frac {\text{d} A_i}{\text{d}x} + \left [\frac {1 + \left (\gamma - 1\right ){Ma}_i^2}{1-{Ma}_i^2}\right ] \frac {F_i}{A_i p_i} - \left [\frac {\left (\gamma -1\right ){Ma}_i^2}{1-{Ma}_i^2}\right ]\frac {Q_i}{p_i A_i u_i}, \\[-28pt] \nonumber \end{align}

(2.2) \begin{align} \frac {1}{p_{t,i}}\frac {\text{d}p_{t,i}}{\text{d}x} = &\frac {F_i}{A_i p_i} - \left [\frac {\left (\gamma -1\right ){Ma}_i^2}{2 + \left (\gamma -1\right ){Ma}_i^2}\right ]\frac {Q_i}{p_i A_i u_i}, \\[-28pt] \nonumber \end{align}

(2.3) \begin{align} \frac {1}{T_{t,i}}\frac {\text{d}T_{t,i}}{\text{d}x} = &\left [1 + \frac {\gamma -1}{2}{Ma}_i^2\right ]^{-1}\frac {Q_i}{p_i A_i u_i}, \\[0pt] \nonumber \end{align}

with $\gamma$ the specific-heat ratio, $F_i$ and $Q_i$ the momentum and energy exchange received by stream $i$ per unit length ( $\text{d}x$ ), $u_i$ the axial velocity, and $ \textit{Ma}_i$ the Mach numbers defined as

(2.4) \begin{align} \textit{Ma}_i = u_i / a_i, \quad \mbox{where } \ a_i = \sqrt {\gamma R T_i}, \end{align}

with $T_i$ the static temperature for each stream and $R$ the specific gas constant. The relations (2.1)–(2.3) apply both in the separate inlets and in the mixing pipe with the practical difference that the definitions of the momentum and energy exchanges $F_i$ and $Q_i$ depend on the domain. For example, wall friction acts where a stream is in contact with a wall, whereas inter-stream shear appears only in the mixing pipe. Details are given in § 3. The reader is referred to Van den Berghe et al. (Reference Van den Berghe, Mendez and Bartosiewicz2024) and to Shapiro (Reference Shapiro1953) for more details on the derivation of these equations.

Energy-exchange terms are retained in (2.1)–(2.3) for completeness, but are neglected in the present study ( $Q_i=0)$ . The walls are treated as adiabatic and the inlet total temperatures are identical (consistent with the RANS dataset; see § 6). Analysis of the reference data shows that inter-stream energy exchange is limited and localised, and that the corresponding contributions in (2.1) and (2.2) are at least two orders of magnitude smaller than those associated with area variation and momentum exchange $F_i$ . This is consistent with the large-eddy simulation of Debroeyer et al. (Reference Debroeyer, Toulorge, Rasquin, Winckelmans and Bartosiewicz2025), which reports only a short region of total enthalpy exchange in the converging portion of the mixing pipe. The transfer of kinetic energy from the primary to the secondary stream is offset by heat transfer in the opposite direction, because of the low temperature of the strongly expanded primary stream. Under this assumption, (2.3) implies uniform total temperature, while (2.2) shows that total-pressure variation is governed solely by the force terms $F_i$ . Therefore, neglecting $Q_i$ has a negligible impact.

Energy-exchange closures may be required when inlet total temperatures differ significantly, but this is left for future work.

The radial position $r_d(x)$ of the dividing streamline determines the local cross-sections of both domains and is defined implicitly from the primary mass flow rate as

(2.5) \begin{align} \dot m_p(x)=\int _0^{r_d(x)} \rho (x,r)\,u(x,r)\,2\pi r\,\text{d}r, \end{align}

so that, at the exit of the primary nozzle $(x=0)$ , the streamline originates at the nozzle lip (figure 1).

The cross-sectional areas $A_i$ and their axial gradients are prescribed by geometry in the inlet sections, but remain unknown in the mixing region where the position of the dividing streamline must be determined. These unknown areas satisfy the geometric constraint $A_p(x) + A_s(x) = A(x)$ , with $A(x)$ the total cross-section of the mixing duct, leaving a single degree of freedom. A unique solution can then be obtained by closing the system in one of two ways: by enforcing equal static pressure in both streams – yielding the compound-choking-based models (e.g. Pearson et al. Reference Pearson, Holiday and Smith1958; Bernstein et al. Reference Bernstein, Heiser and Hevenor1967) – or by prescribing the dividing streamline from an external source, which defines the Fabri-choking-based models.

In this work, both approaches were implemented and evaluated, highlighting their respective advantages and limitations, as well as their potential common ground. Furthermore, two alternatives are considered for the exchange terms in (2.1)–(2.3): they were either defined through correlations or extracted from the axisymmetric RANS simulations presented in § 6. The test cases in § 7 are organised according to these two aspects and are introduced later.

2.2. Selected models

The four model configurations investigated in this work combine two treatments of momentum exchange with two choking mechanisms, as illustrated in figure 1. The compound-based models include a fully predictive case using empirical friction correlations (Model 1, § 2.2.1) and a non-predictive one where total pressure gradients are imposed from the RANS simulations (Model 2, § 2.2.2). Comparing these two isolates the role of the closure terms and reveals the intrinsic limitations of the compound formulation.

The Fabri-based models follow the same logic: the dividing streamline and forces are either imposed from RANS data (Model 3, § 2.2.3) or taken from the compound model while retaining its closure correlations (Model 4, § 2.2.4). Comparing these highlights how the choking mechanism itself influences the predicted flow evolution and helps identify conditions where the Fabri and compound frameworks may converge.

2.2.1. Model 1: compound choking with calibrated closures

Parallel streams choke according to the compound flow theory when a sonic condition involving all streams is met. This is commonly modelled through the constraint of uniform static pressure at any streamwise position $x$ , as proposed by Pearson et al. (Reference Pearson, Holiday and Smith1958) and Bernstein et al. (Reference Bernstein, Heiser and Hevenor1967) among others. Summing (2.1) for the primary and secondary streams yields the pressure and setting $Q_i=0$ ,

(2.6) \begin{align} \frac {1}{p} \frac {\text{d}p}{\text{d}x} = \dfrac {1}{\beta }\left( \dfrac {\text{d}A}{\text{d}x} + \sum \limits _{i\in [p,s]}\left [\dfrac {1 + \left (\gamma - 1\right ){Ma}_i^2}{\gamma {Ma}_i^2}\right ] \dfrac {F_i}{p}\right )\!, \end{align}

where $\beta$ is the compound choking indicator:

(2.7) \begin{align} \beta = \sum _{i\in [p,s]} A_i\dfrac {1-{Ma}_i^2}{\gamma {Ma}_i^2}. \end{align}

This plays a similar role as $1 - {Ma}^2$ in (2.1): if $\beta$ equals zero, the system of equations is singular, which corresponds to the sonic point of the parallel streams or the so-called ‘compound’ flow. Note that the definition (2.7) of $\beta$ is identical to that in the original isentropic theory by Pearson et al. (Reference Pearson, Holiday and Smith1958) and Bernstein et al. (Reference Bernstein, Heiser and Hevenor1967). The non-isentropic effects are handled through the exchange terms $F_i$ . Alternatively, these have been included in the definition of $\beta$ by Debroeyer et al. (Reference Debroeyer, Winckelmans and Bartosiewicz2026). The reader is referred to Van den Berghe et al. (Reference Van den Berghe, Mendez and Bartosiewicz2024) for a detailed analysis of the compound flow theory in a 1-D formulation. Hedges & Hill (Reference Hedges and Hill1974) defined an equivalent Mach number $ \textit{Ma}_{\textit{eq}}$ related to $\beta$ as follows:

(2.8) \begin{align} {Ma}_{\textit{eq}} = \left (\gamma \frac {\beta }{A}+1\right )^{-\tfrac {1}{2}}. \end{align}

The evolution of the dividing streamline is computed explicitly by inverting (2.1) to obtain the equation for the individual cross-sections $A_i$ :

(2.9) \begin{align} \frac {\text{d}A_i}{\text{d}x} = \left [A_i\frac {1-{Ma}_i^2}{\gamma {Ma}_i^2}\right ]\frac {1}{p} \frac {\text{d}p}{\text{d}x} - \left [\frac {1 + \left (\gamma - 1\right ){Ma}_i^2}{\gamma {Ma}_i^2}\right ] \frac {F_i}{p} , \end{align}

in which the pressure gradient is computed with (2.6). The compound choking model therefore provides a direct means to compute the evolution of the dividing streamline under the assumption of equal static pressures between the streams. Moreover, the formulation naturally accommodates the exchange terms, which appear in the governing equations in a generic form. In particular, the pressure constraint remains satisfied whether the friction forces are prescribed from predictive correlations or derived from the total pressure gradients imposed by the RANS simulations.

The constraint of uniform static pressure eliminates the degree of freedom in (2.1)–(2.3) introduced by the dividing streamline. This leaves only the exchange terms undefined. The momentum exchange is embodied in this configuration by friction forces between the streams, and between the secondary stream and the wall as shown by Van den Berghe et al. (Reference Van den Berghe, Mendez and Bartosiewicz2024). The wall friction is modelled with a classic friction coefficient $f_w$ :

(2.10) \begin{align} F_{w,i} = \frac {1}{2} f_w \gamma p_i {Ma}_i^2 l_{w,i}, \end{align}

where $l_{w,i}$ denotes the perimeter of the wall(s) at a given position $x$ . Note that the inner and the outer walls are considered for the secondary inlet. Wall friction is the only active force in the inlets ( $x \lt 0$ ), so $F_i = -F_{w,i}$ in (2.1) and (2.2). It only acts on the secondary stream in the mixing pipe, so the local static pressure and Mach number of the secondary stream and the perimeter of the entire cross-section $A$ are used in (2.10). The friction coefficient is computed using the correlation of Van Driest (Reference Van Driest1951):

(2.11) \begin{align} \dfrac {0.242}{\sqrt {f_w}} \sqrt {1 - \lambda ^2} \dfrac {\arcsin (\lambda )}{\lambda } = 0.41 + \log _{10} \left (f_w Re_x\right ) + \log _{10} \left ((1 - \lambda ^2)\left (1 - \dfrac {\theta \lambda ^2}{1 + \theta }\right )\right )\!, \end{align}

where

(2.12) \begin{align} 1-\lambda ^2 = \left (1 + \frac {\gamma - 1}{2} {Ma}_s^2\right )^{-1}\!, \quad Re_x = \dfrac {\rho _s u_s x}{\mu _s}\, \quad \text{and} \quad \theta = \dfrac {S}{T_s}, \end{align}

where $\rho$ denotes the density, $S = 110.4$ K and the law of Sutherland (Reference Sutherland1893) is used for the dynamic viscosity $\mu$ , with a reference viscosity $\mu _{\textit{ref}}$ of $1.716 \ 10^{-5}$ Pa s at $T_{\textit{ref}} = 273.15$ K. The friction between the streams is defined following Papamoschou (Reference Papamoschou1993):

(2.13) \begin{align} F_{\textit{ps}} = \frac {1}{2} f_{\textit{ps}} \dfrac {\rho _p+\rho _s}{2} \left (u_p - u_s\right ) \left |u_p - u_s\right | l_{\textit{ps}}, \end{align}

where $f_{\textit{ps}}$ denotes a friction coefficient and $l_{\textit{ps}}$ denotes the perimeter of the primary stream:

(2.14) \begin{align} l_{\textit{ps}} = 2 \sqrt {\pi A_p}. \end{align}

The inter-stream friction coefficient $f_{\textit{ps}}$ is computed as by Papamoschou (Reference Papamoschou1993):

(2.15) \begin{align} f_{\textit{ps}} = 0.013 \dfrac {(1 + \zeta ) (1 + \eta )}{1 + \zeta \eta } \big (0.25 + 0.75 \exp \big(\!-\!3 {Ma}_c^2 \big)\big ), \end{align}

where $\zeta =u_s/u_p$ , $\eta =\sqrt {\rho _s/\rho _p}$ and the local convective Mach number $ \textit{Ma}_c$ is defined as follows:

(2.16) \begin{align} {Ma}_c = \dfrac {u_p - u_s}{a_p + a_s}. \end{align}

The inter-stream friction is defined to be positive if $u_p \gt u_s$ , in which case it counteracts the primary stream. The forces $F_p$ and $F_s$ are thus defined as follows in the mixing pipe:

(2.17) \begin{align} F_p = -F_{\textit{ps}} \quad \text{and} \quad F_s = F_{\textit{ps}} - F_{w,s}. \end{align}

The effect of these forces on the sonic point is discussed in detail by Van den Berghe et al. (Reference Van den Berghe, Mendez and Bartosiewicz2024).

3.1. Treatment of the inlets

For all the models discussed in the previous section, the friction coefficients in the inlets are taken using the empirical correlations in (2.10)–(2.12).

The boundary conditions specify the total pressure $p_{t,i}$ and total temperature $T_{t,i}$ at each inlet. The governing equations (2.1)–(2.3) form an explicit set of ordinary differential equations (ODEs) that are integrated axially from the inlets, where the total quantities are known. A shooting method is used to adjust the inlet static pressure until the target sonic condition is met (see § 4).

The same formulation applies to both streams. However, since the primary jet is choked in its converging–diverging nozzle, it can be solved independently. The secondary flow, however, depends on the back pressure $p_b$ at the outlet of the diffuser and must therefore be solved in conjunction with the mixing section, as discussed later.

Figure 2.

The jump conditions from the inlets to the mixing pipe. The primary cross-sections and, hence, the flow states are identical on both sides. The secondary stream requires a force balance, where the momentum fluxes and pressure forces are projected on the x-axis with their respective angles. The result is a double flow state at the inlet of the mixing pipe $x=0$ , which generally has different static pressures. Schematics are intended for illustration purposes only and do not represent the geometry used in the RANS simulations (see figure 6).

3.2. Jump relations at the inlet of the mixing pipe

The full flow states are known at the outlets of the separate domains ( $x=0$ ) after integrating the respective ODEs in the inlets. The dividing streamline originates at the bottom of the primary nozzle lip, so the primary cross-sections and variables match between the inlet, denoted with the subscript $L$ , and in the mixing pipe, with subscript $R$ (see figure 2). This is not the case for the secondary stream due to its inclined inlet and a potential finite thickness $A_w$ of the nozzle. The secondary flow state on the right is found through conservation of mass, momentum and energy. The control volume is infinitesimally thin ( $\Delta x \rightarrow 0$ in figure 2). Since $A_{s,R}+A_{p,R}=A(x=0)$ , three unknowns remain (e.g. density, velocity and temperature). Two forces intervene in the axial momentum balance: the force on the bottom wall $\bar {F}_{w,b}$ and the force on the nozzle lip $\bar {F}_{w,ps}$ . It is assumed that friction acting on these surfaces is negligible compared with the pressure force, which has been confirmed during the post-processing of the RANS simulations in this work. These forces are therefore orthogonal to the surfaces on which they act. The inclinations of the bottom and top wall are generally different and are denoted here with $\alpha _b\gt 0$ and $\alpha _t\gt 0$ . The boundary of the secondary inlet is defined as the line originating from the top wall at the exit of the primary nozzle ( $x=0$ ), with a central angle $\alpha _c = 0.5(\alpha _b + \alpha _t)$ . Projecting the momentum fluxes and the forces on the x-axis with their respective angles leads to the following balance:

(3.1) \begin{align} \left (\dot {m}_s u_{s,L} + p_{s,L} A_{s,L}\right ) \cos {\alpha _c} + F_{w,b} \sin {\alpha _b} + F_{w,ps} = \dot {m}_s u_{s,R} + p_{s,R} A_{s,R}, \end{align}

where $u$ denotes the average velocity as described in the 1-D equations. The magnitude of the forces $F_{w,b}$ and $F_{w,ps}$ need to be estimated from the averaged quantities of the left and right states. The expansion of the secondary stream from the inclined section on the left to $x=0$ is generally small, hence the force on the bottom wall is approximated as $F_{w,b} = p_{s,L} A_{w,b}$ , where the area can be computed from the geometry. Similarly, the force on the nozzle lip is defined as $F_{w,ps} = p_{s,R} A_{w,ps}$ , where the (unknown) static pressure on the right side is applied since it is in direct contact with this wall.

Conservation of energy implies that the total temperature is conserved, since there is no source term and the walls are assumed to be adiabatic. The remaining unknowns in (3.1) are the velocity $u_{s,R}$ and the static pressure $p_{s,R}$ . These are linked through the known mass flow rate $\dot {m}_s$ , the known total temperature $T_{t,s}$ and the ideal gas law:

(3.2) \begin{align} \dot {m}_s = \frac {p_{s,R}}{R T_{s,R}} A_{s,R} u_{s,R} \, \quad \text{and} \quad T_{s,R} = T_{t,s,R} - \frac {u_{s,R}^2}{2c_P}. \end{align}

The two equations in (3.2) allow elimination of $p_{s,R}$ from the momentum balance (3.1), leading to a nonlinear equation in $u_{s,R}$ , which can be solved iteratively. The velocity $u_{s,L}$ from the secondary inlet is a good initial guess. The secondary flow state can then be computed from $\dot {m}_s$ , $T_{t,s,R}$ and $u_{s,R}$ .

The resulting right primary and secondary states generally have a different static pressure because the primary stream is over- or under-expanded. A uniform static pressure $p(x)=p_{{\kern-1pt}p}(x)=p_s(x)$ is required for the compound model (see § 2.2.1), so first, a pressure equalisation mechanism is defined later. This is not required when using the Fabri-based model, where the static pressures are left free in each stream.

Figure 3.

An under- or over-expanded primary stream is brought to a uniform static pressure over a finite distance using Prandtl–Meyer expansion theory or oblique shocks, respectively. The procedure results in the angle $\theta$ of the dividing streamline, which suffices to complement the conservation equations (2.1)–(2.3).

3.3. Pressure equalisation mechanism for the compound theory

At the inlet of the mixing pipe, the primary stream is typically either under- or over-expanded relative to the surrounding secondary stream, as evidenced by the shock-train structures observed in numerical studies (Bartosiewicz et al. Reference Bartosiewicz, Aidoun, Desevaux and Mercadier2005; Hemidi et al. Reference Hemidi, Henry, Leclaire, Seynhaeve and Bartosiewicz2009). This conflicts with the assumption of equal static pressures in the compound flow theory and thus with model 1 (§ 2.2.1) and model 2 (§ 2.2.2). To address this limitation, the present work introduces a pressure equalisation mechanism that brings the two streams to a uniform static pressure over a finite distance downstream of the primary nozzle outlet, as illustrated schematically in figure 3(a).

The proposed mechanism is inspired by the first cell of the shock train and reproduces the local pressure-adjustment process up to the point where the average static pressures of the two streams first become equal. The analysis is local in nature and provides the orientation of the dividing streamline as a function of the local static-pressure ratio, as shown in figure 3. Separate formulations are adopted for under-expanded and over-expanded primary streams, though the underlying principle remains the same.

For an under-expanded primary stream, further expansion is modelled through a diverging cross-section computed using Prandtl–Meyer expansion theory, following the approach of del Valle et al. (Reference del Valle, Jabardo, Ruiz and Alonso2012). This formulation implicitly assumes that the primary jet is supersonic at the nozzle outlet. Since practical ejectors typically employ convergent–divergent primary nozzles, this condition is generally satisfied in practical operating conditions and does not significantly restrict the applicability of the model. The process is assumed isentropic, ensuring conservation of total pressure from the current state $(p_{{\kern-1pt}p}, \textit{Ma}_{{\kern-1pt}p})$ to a virtual state $(p_s, \widetilde {\textit{Ma}}_p)$ characterised by the same static pressure as the secondary stream:

(3.3) \begin{align} p_{t,p} = p_s \left (1 + \frac {\gamma - 1}{2}\,\widetilde {{Ma}}_p^{\,2} \right )^{({\gamma }/({\gamma - 1}))}, \end{align}

from which the expanded Mach number $\widetilde {{Ma}}_p$ is determined. The deviation angle $\theta$ is then defined as the difference between the Prandtl–Meyer angles $\nu$ and $\tilde {\nu }$ corresponding to the local and virtual primary Mach numbers $ \textit{Ma}_{{\kern-1pt}p}$ and $\widetilde {{Ma}}_p$ , respectively. For completeness, the Prandtl–Meyer angle is given by (Shapiro Reference Shapiro1953)

(3.4) \begin{align} \nu ({Ma}) \;=\; \sqrt {\frac {\gamma + 1}{\gamma - 1}}\, \arctan \!\left (\sqrt {\frac {\gamma - 1}{\gamma + 1}\,\big ({Ma}^{2}-1\big )}\right ) \;-\; \arctan \!\left (\sqrt {{Ma}^{2}-1}\right )\!. \end{align}

If the primary and secondary static pressures were equal, the dividing streamline would be tangent to the wall at the outlet of the primary nozzle (without sudden change of angle). The computed deviation angle $\tilde {\nu }-\nu$ should thus be added to the angle $\theta _p$ of the primary nozzle at its outlet to obtain the angle of the dividing streamline:

(3.5) \begin{align} \theta = \theta _p + \left (\tilde {\nu }-\nu \right )\!. \end{align}

An oblique shock is computed for the pressure ratio $p_s/p_{{\kern-1pt}p}$ and a primary Mach number $ \textit{Ma}_{{\kern-1pt}p}$ in case of an over-expanded primary stream. The shock angle $\alpha$ can be computed from the following equation (Shapiro Reference Shapiro1953):

(3.6) \begin{align} \frac {p_s}{p_{{\kern-1pt}p}} = 1 + \frac {2\gamma }{\gamma +1} \left ( {Ma}_p^2 \sin ^2(\alpha ) - 1 \right )\!. \end{align}

The deviation angle $\theta$ is directly related to the shock angle $\alpha$ and the Mach number $ \textit{Ma}_{{\kern-1pt}p}$ :

(3.7) \begin{align} \tan (\theta ) = -\frac {2 \cot (\alpha ) \left ( {Ma}_p^2 \sin ^2(\alpha ) - 1 \right )}{{Ma}_p^2 \left ( \gamma + \cos (2\alpha ) \right ) + 2}. \end{align}

Note the minus sign for the inward deviation of the primary flow. In the axisymmetric geometry, the axial derivative of the primary cross-section can be decomposed in terms of the deviation angle $\theta$ :

(3.8) \begin{align} \frac {\text{d}A_p}{\text{d}x} = 2 \pi r_p \frac {\text{d}r_p}{\text{d}x} = 2 \sqrt {\pi A_p} \tan (\theta ). \end{align}

The axial derivative of the secondary cross-section follows from the geometry $A(x)$ of the mixing pipe. These local gradients of the cross-sections allow to integrate the governing equations (2.1)–(2.3) (which remain valid in the mixing pipe) to the next grid point, using the closure relations (2.17). The overall process is not isentropic due to the friction forces; the isentropic pressure equalisation mechanism is used only to determine the local angle of the dividing streamline. The above-mentioned procedure is repeated until the static pressure difference drops below an arbitrary threshold (10 Pa in this work).

This approach reproduces only the first cell of the shock train and cannot capture downstream behaviour, as reflections of Mach lines on the centreline and dividing streamline are inherently two-dimensional and thus excluded from the present one-dimensional model. Nevertheless, the proposed pressure equalisation mechanism offers two main advantages: (i) its local formulation makes it directly compatible with standard ODE integrators; and (ii) the pressure equalisation occurs over a short distance, enabling the compound flow theory to be applied as far upstream as possible to accurately locate the sonic point.

3.4. Normal compound shocks

An oblique shock train typically forms in the diffuser under on-design conditions, decelerating the flow to subsonic speed and re-establishing the back pressure imposed by the boundary conditions (Bartosiewicz et al. Reference Bartosiewicz, Aidoun, Desevaux and Mercadier2005; Hemidi et al. Reference Hemidi, Henry, Leclaire, Seynhaeve and Bartosiewicz2009). In simplified models, this system is often represented as a normal shock in 0-D formulations (Huang et al. Reference Huang, Chang, Wang and Petrenko1999; Chen et al. Reference Chen, Liu, Chong, Yan, Little and Bartosiewicz2013; Metsue et al. Reference Metsue, Debroeyer, Poncet and Bartosiewicz2021) and it also emerges naturally in the unsteady single-stream model of Van den Berghe et al. (Reference Van den Berghe, Dias, Bartosiewicz and Mendez2023).

In the compound-flow framework of model 1 and model 2, however, a normal-shock representation cannot be directly applied because the compound supersonic region may include a subsonic component, rendering standard shock relations invalid. Furthermore, shocks of unequal strength would develop in the two streams due to their differing velocities, producing a downstream pressure imbalance and violating the assumption of uniform static pressure. The inherently two-dimensional nature of the oblique shock train and its interaction with boundary layers further preclude an accurate one-dimensional description under the constraint of pressure uniformity.

Since the 1-D formulation cannot resolve the flow structure downstream of the shock, the mixing between the two streams is also undefined. To apply normal-shock relations in this context, it is therefore necessary to define an equivalent single-stream state that represents the combined effect of the primary and secondary streams at the shock location. This equivalent, or ‘fully mixed’, state is a notional construct used solely for the computation of the normal shock and is defined to have the same total mass flow rate, momentum and energy as the actual double-stream configuration. The normal-shock relations are then applied between this equivalent upstream state, denoted by the subscript $L$ , and the downstream state, denoted by $R$ . The properties of the equivalent state upstream of the shock are obtained from the conservation of mass and energy as

(3.9) \begin{align} \dot {m}_L = \dot {m}_p + \dot {m}_s \quad \text{and} \quad T_{t,L} = \frac {\dot {m}_p T_{t,p} + \dot {m}_s T_{t,s}}{\dot {m}_p + \dot {m}_s}. \end{align}

The mass flow rate and the momentum flux can be defined in terms of the static pressure $p$ , the total temperature $T_t$ , the cross-section $A$ and the Mach number $ \textit{Ma}$ :

(3.10) \begin{align} \dot {m} = p A {Ma} \sqrt {\frac {\gamma }{R T_t}}\sqrt {1+\frac {\gamma -1}{2}{Ma}^2}\, \quad \text{and} \quad \dot {m} u + p A = pA \big (1 + \gamma {Ma}^2\big ). \end{align}

The static pressure and the cross-section can be eliminated by taking the ratio of the momentum flux and the mass flow rate, which becomes a function of the Mach number and the total temperature:

(3.11) \begin{align} \frac {\dot {m} u + p A}{\dot {m}} = \frac {1 + \gamma {Ma}^2}{{Ma} \sqrt {1+(({\gamma -1})/{2}){Ma}^2}} \sqrt {\frac {R T_t}{\gamma }}. \end{align}

Hence, conservation of momentum between the double state and the equivalent state $L$ implies

(3.12) \begin{align} \frac {\sum _{i \in [p,s]}\left (\dot {m}_i u_i + p_i A_i\right )}{\sum _{i \in [p,s]}\dot {m}_i} = \left [\frac {1 + \gamma {Ma}_L^2}{{Ma}_L \sqrt {1+(({\gamma -1})/{2}){Ma}_L^2}}\right ] \sqrt {\frac {R T_{t,L}}{\gamma }}, \end{align}

where the total temperature $T_{t,L}$ is known from (3.9). The term between the brackets is nonlinear in $ \textit{Ma}_L$ , and admits a subsonic and a supersonic solution for a given ratio of momentum to mass flow rate. The minimum corresponds to a sonic flow, indicating a minimal momentum for a given mass flow rate (or a maximal mass flow rate for a given amount of momentum). The supersonic solution is retained, followed by the normal shock relations to compute the right state $R$ (see for example Shapiro Reference Shapiro1953):

(3.13) \begin{align} M_R^2 = \frac {( \gamma - 1 ) M_L^2 + 2}{2\gamma M_L^2 - ( \gamma - 1 )}, \quad p_R = p_L \left (1 + \frac {2 \gamma }{\gamma + 1} \left ( M_L^2 - 1 \right )\right )\, \quad \text{and} \quad T_{t,R} = T_{t,L}. \end{align}

The flow downstream is computed as an isentropic single stream, which has an analytical solution considering that the total pressure and total temperature are conserved and that the mass flow rate is known. The Mach number $ \textit{Ma}_m$ of the single stream respects the following (nonlinear) equation in any point downstream of the shock:

(3.14) \begin{align} \dot {m}_L = p_{t,R} A(x) \sqrt {\frac {\gamma }{R T_{t,R}}} {Ma}_m(x) \left [1 + \frac {\gamma - 1}{2} {Ma}_m^2(x)\right ]^{-( {1}/{2})(( {\gamma +1})/({\gamma -1}))}. \end{align}

3.5. Filtering, differentiation and interpolation

Imposing cross-section or total-pressure gradients from the RANS simulations requires computing and interpolating these quantities so they can be evaluated by the solver at any axial position $x$ . The primary and secondary cross-sections and total pressures are obtained directly from the geometry, the dividing streamline, and the flow field (see § 6). However, their derivatives in (2.1), obtained by numerical differentiation, are sensitive to noise due to interpolation and cross-stream averaging.

The distributions of the total pressures and the primary cross-section are therefore first filtered before applying finite differences with Numpy’s gradient function. The cross-section of the secondary stream and its derivatives follow from the known geometry and the primary stream. The filter used in this work is based on the filter of Chambolle (Reference Chambolle2004) available in the denoise_tv_chambolle function from Scikit-image (Van der Walt et al. Reference Van der Walt, Schönberger, Nunez-Iglesias, Boulogne, Warner, Yager, Gouillart and Yu2014). The total variation denoising formulation avoids spurious overshoot near sharp transitions, but it is subject to boundary effects, which are problematic in the current work because the dividing streamline originates from the lip of the primary nozzle. Therefore, the signal used in the model is a combination of the original and the filtered signal:

(3.15) \begin{align} \tilde {u}[j] = w[j] u[j] + \hat {w}[j] \hat {u}[j], \end{align}

where $\tilde {u}[j]$ denotes the $j{\text{th}}$ entry of the final discrete signal, composed of the original and the filtered signals $u$ and $\hat {u}$ . These are blended via a partition of unity that balances the unfiltered and the filtered signal such that the sum of the weights equals one at each point, that is,

(3.16) \begin{align} w[j] = \begin{cases} 0.5 + 0.5 \cos \left (\pi \frac {i}{60}\right ) &\text{ if } j \leqslant 60,\\[12pt] 0.5 - 0.5 \cos \left (\pi \frac {i - (N-61)}{60}\right ) &\text{ if } j \gt N-60,\\[4pt] 0 &\text{ otherwise} \end{cases} \end{align}

and

(3.17) \begin{align} \hat {w}[j] = \begin{cases} 0.5 - 0.5 \cos \left (\pi \frac {i}{60}\right ) &\text{ if } j \leqslant 60,\\[12pt] 0.5 + 0.5 \cos \left (\pi \frac {i - (N-61)}{60}\right ) &\text{ if } j \gt N-60,\\[4pt] 1 &\text{ otherwise.} \end{cases} \end{align}

The original signal is thus kept near the boundaries, transitioning over 60 entries to the filtered signal which is kept elsewhere. The width of 60 entries has been chosen arbitrarily as a compromise between filter intrusiveness near the boundary and noise tolerance.

Finally, since the full solver requires the evaluation of the governing equations (2.1)–(2.3) at any axial position $x$ , an interpolation within the filtered grid was applied. This is a cubic interpolation using Scipy’s ‘interp1d’ function (Virtanen et al. Reference Virtanen2020).

4.1. Primary inlet

Given the inlet totals $(p_{t,p},T_{t,p})$ , the primary mass flow is obtained by a bracketed shooting on the inlet static pressure $p_{{\kern-1pt}p}(x_0)$ to enforce $M_p=1$ (and $N_p=0$ in (4.3)). By definition, the upper bound is $p_{{\kern-1pt}p}(x_0)\leqslant p_{t,p}$ , while the lower bound is the isentropic sonic pressure

(4.1) \begin{align} p_{{sonic}} = p_{t,p}\left (\frac {\gamma +1}{2}\right )^{-\gamma /(\gamma -1)}, \end{align}

since the flow is subsonic in the converging section. The search therefore brackets $p_{{\kern-1pt}p}(x_0)\in [p_{{sonic}},\,p_{t,p}]$ and proceeds by bisection: at each iteration, a trial $p_{{\kern-1pt}p}^{(k)}$ (the bracket midpoint) is used to integrate (2.1)–(2.3) from $x_0$ to either the point where $M_p=1$ or the nozzle exit. Reaching the exit while still subsonic indicates $p_{{\kern-1pt}p}^{(k)}$ is too high and hence the upper bound is set to $p_{{\kern-1pt}p}^{(k)}$ . Otherwise, the bracket is tightened around the $M_p=1$ solution. This monotone update exploits the fact that decreasing $p_{{\kern-1pt}p}(x_0)$ increases the mass flow towards choking (see the flow chart in figure 4).

As the sonic condition is approached, the denominator in (2.1) tends to zero. To obtain a finite ratio and allow the flow to continuously expand to supersonic conditions, also the numerator must tend to zero (see Shapiro Reference Shapiro1953). Isolating the common denominator and defining the numerator $N_p$ :

(4.2) \begin{align} \frac {1}{p_{{\kern-1pt}p}} \frac {\text{d}p_{{\kern-1pt}p}}{\text{d}x} = \frac {N_p}{1-{Ma}_p^2}, \end{align}

where

(4.3) \begin{align} N_p = \gamma {Ma}_p^2\frac {1}{A_p}\frac {\text{d} A_p}{\text{d}x} + \left (1 + \left (\gamma - 1\right ){Ma}_p^2\right ) \frac {F_p}{A_p p_{{\kern-1pt}p}}. \end{align}

A finite pressure gradient is obtained if

(4.4) \begin{align} \lim _{(p,{Ma}_p)\to (p^\ast ,1)} \frac {1}{p_{{\kern-1pt}p}} \frac {\text{d}p_{{\kern-1pt}p}}{\text{d}x}= \frac {\text{d}N_p/\text{d}x}{-{\rm d}{Ma}_p^2/\text{d}x} \end{align}

following de l’Hôpital’s rule. The derivative of the Mach number can be analytically evaluated from the conservation equation (see (A3)), which introduces the derivative of the static pressure in the denominator. This leads to the quadratic equation (A16) for the static pressure gradient in the sonic point, which is derived in Appendix A (see also Shapiro Reference Shapiro1953; Van den Berghe et al. Reference Van den Berghe, Mendez and Bartosiewicz2024). The static pressure gradient can thus be calculated in all admissible conditions.

Figure 4.

Flow chart for computing the choked flow in the primary nozzle. The static pressure at the inlet is found iteratively to satisfy the sonic condition $ \textit{Ma}_{{\kern-1pt}p} \rightarrow 1, N_p \rightarrow 0$ in (4.2). The approximations in Appendix A allow to integrate the governing equations through the sonic point without numerical issues related to division by zero (see also Restrepo & Simões-Moreira Reference Restrepo and Simões-Moreira2022; Van den Berghe et al. Reference Van den Berghe, Mendez and Bartosiewicz2024).

If the procedure in figure 4 runs into the sonic point, the numerator $N_p$ is computed with (4.3). If its absolute value drops below an arbitrary threshold ( $10^{-6}$ in this work), the correct sonic condition is found and the calculation proceeds. Otherwise, the bracket is adjusted. A negative value of the numerator indicates that the flow would expand if it is subsonic ( $\text{d}p/\text{d}x \lt 0$ in (4.2)). This occurs for example in the convergent section for an isentropic flow (see (4.3) with $F=0$ ). The fact that the flow already reached $ \textit{Ma}_{{\kern-1pt}p}=1$ indicates that the flow rate is too high: the Mach number should be lower, allowing the flow to expand further without running into the singularity of (4.2). In the isentropic example, the sonic point should be reached further downstream at the throat. The guessed static pressure is thus too low if $N_p(x) \lt 0$ and $M_p(x)=1$ . Therefore, the used average pressure $p(x_0)$ becomes the new lower bound $p_{\textit{min}}$ . The inverse is true if $N_p\gt 0$ . The friction force $F_p\lt 0$ interferes in the numerator $N_p$ , but it does not alter the mechanism discussed earlier. Note that $\text{d}A/\text{d}x$ must be positive to compensate the negative term containing the force $F_p$ in (4.3): the wall friction force pushes the sonic section downstream in a divergent section of the nozzle. The algorithm iterates until the tolerance on the numerator is reached or until the bracket reaches the following relative tolerance:

(4.5) \begin{align} \frac {p_{\textit{max}}-p_{\textit{min}}}{0.5 (p_{\textit{min}} + p_{\textit{max}})} \leqslant 10^{-6}. \end{align}

The above-mentioned procedure provides the solution up to the sonic point. At this point, the static pressure gradient is undetermined as (2.1) approaches the limit $0/0$ , so a second-order approximation of the pressure gradient near the sonic point is computed using Taylor expansions (see (A18) in Appendix A) if $M_p \leqslant 1.01$ . Once the primary stream expands beyond $ \textit{Ma}_{{\kern-1pt}p} = 1.01$ , the pressure is computed from (4.2). The threshold of $1.01$ was found to provide a good compromise between accuracy and smoothness of the transition from the approximation (A18) to the governing equation (4.2). In practice, the approximation (A18) is blended with (4.2) in a single function that automatically switches between both depending on the primary Mach number $ \textit{Ma}_{{\kern-1pt}p}$ , as indicated in figure 4. The interested reader is referred to Van den Berghe et al. (Reference Van den Berghe, Mendez and Bartosiewicz2024) for a more detailed description of this approach.

4.2. On-design operation

The procedure for the rest of the ejector is similar to that of the primary nozzle, as shown in figure 5: the static pressure is guessed at the secondary inlet, the governing equations are integrated to the compound- or Fabri-sonic point ( $\beta =0$ in (2.6) or $ \textit{Ma}_s=1$ in (2.1)) and the solution is accepted if the numerator of (2.6) or (2.1) equals zero.

An approximation of the static pressure gradient is then computed near the sonic point similarly as in the primary nozzle. The required derivatives for the compound case can be found in the work of Van den Berghe et al. (Reference Van den Berghe, Mendez and Bartosiewicz2024) and those for the Fabri case in Appendix A. The relevant equations are identical to those of the single stream in the primary nozzle, with the difference that the force term also includes the inter-stream friction. The blended equations are integrated to the outlet, switching back to the direct equations once $ \textit{Ma}_{\textit{eq}} \geqslant 1.01$ or $ \textit{Ma}_s \geqslant 1.01$ .

The flow might become sonic in the secondary inlet, in which case, the guessed static pressure is increased to reduce the mass flow rate and reach a sonic condition in the mixing pipe instead. For the geometry analysed in this work, the sonic point is always located in the mixing pipe, which is the most typical scenario. The analysis of the sign of the numerator in the computed sonic point is identical to that above-mentioned for the primary nozzle due to the similarity of (2.1) and (2.6).

The resulting solution is fully supersonic in the diffuser, so the normal shock from § 3.4 is required to match the back pressure. Its axial position is found iteratively to respect the outlet condition.

Figure 5.

The solution procedure iterates on the static pressure at the secondary inlet to find the correct compound- or Fabri-sonic point. A sonic flow in the secondary inlet is not accepted as a valid solution since it is atypical for ejectors. The procedure differs in the compound and Fabri case through the sonic condition and the corresponding approximations of the pressure gradient.