1. Introduction
Parallel compressible streams are at the heart of many applications, the most classic example being jet engines with a high-speed core and a slower surrounding jet. Examples of internal flows with multiple streams include ejectors or (sc)ramjets. Such heterogeneous flows are referred to as compound flows. This study investigates the choking of internal compound flows.
Recently, Lamberts et al. (Reference Lamberts, Chatelain, Bourgeois and Bartosiewicz2018) observed that compound flows could choke upstream of the geometrical throat, in contrast to the classic downstream displacement of the sonic section due to friction. By numerically increasing the wall roughness, Kracik & Dvorak (Reference Kracik and Dvorak2023) showed that friction has this effect on compound flows as well. Lamberts et al. (Reference Lamberts, Chatelain, Bourgeois and Bartosiewicz2018) attributed the upstream shift to non-isentropic effects, but the inherent cause remained unknown. While both Lamberts et al. (Reference Lamberts, Chatelain, Bourgeois and Bartosiewicz2018) and Kracik & Dvorak (Reference Kracik and Dvorak2023) addressed this problem through axisymmetric Reynolds-averaged Navier–Stokes (RANS) simulations, to the best of the authors’ knowledge, no theoretical analysis of the sonic section displacement has been provided.
The first theoretical treatment of compound flows was proposed by Pearson, Holiday & Smith (Reference Pearson, Holiday and Smith1958). These authors introduced the idea of splitting the domain into parallel one-dimensional (1-D) streams bounded by dividing streamlines and walls. The position of these streamlines can be identified by solving mass, momentum and energy conservation in each stream, together with a force balance at the dividing streamlines. The authors assumed adiabatic and isentropic conditions for the compound flow. Hence, only a pressure force is exerted between the streams, and the force balance requires equal static pressure for the entire cross-section. This formulation leads to a closed system of equations linking the evolution of mass, momentum and energy in each stream to the evolution of the respective cross-sections and streamwise pressure distribution. Pearson et al. (Reference Pearson, Holiday and Smith1958) state that the assumption of equal pressure is reasonable if the curvature of the streamlines remains limited.
Through some manipulations of the conservation equations, Bernstein, Heiser & Hevenor (Reference Bernstein, Heiser and Hevenor1967) showed that a single equation can be obtained for the uniform static pressure. Similarly to a single 1-D stream (see e.g. Shapiro Reference Shapiro1953), this equation is singular with a division of zero by zero if the flow is choked. The major result of the study by Pearson et al. (Reference Pearson, Holiday and Smith1958) lies in the fact that a compound flow can be choked with a combination of subsonic and supersonic streams. They showed that information cannot travel upstream in the subsonic streams due to the assumption of uniform static pressure. The derivation by Bernstein et al. (Reference Bernstein, Heiser and Hevenor1967) confirms the findings of Pearson et al. (Reference Pearson, Holiday and Smith1958) identifying the compound choking from the singularity of the pressure equation. The same compound choking condition was found by Hoge & Segars (Reference Hoge and Segars1965) by minimising the impulse function defined in Chapter 4 of Shapiro (Reference Shapiro1953). An alternative approach by Agnone (Reference Agnone1971) identifies the compound choking from the singularity of the governing equations, written in a matrix form, while Fage (Reference Fage1976) defines it through the so-called ‘elasticity of the fluid boundary’ between the streams. Bernstein et al. (Reference Bernstein, Heiser and Hevenor1967) gave the name ‘compound flow theory’, which we also adopt. These works did not consider exchanges between the streams; hence, the flow remained isentropic.
Later work investigated the effect of forces on the compound flow. Hoge & Segars (Reference Hoge and Segars1965) explored the effects of friction between the streams and the walls, assuming isentropic and adiabatic conditions but without deriving an explicit model to account for their impact on the flow. These authors postulated that the exchange of momentum between the streams would have a minimal impact on the compound choking mechanism. A more explicit approach is presented by Clark (Reference Clark1995), Papamoschou (Reference Papamoschou1996) and Grazzini, Mazzelli & Milazzo (Reference Grazzini, Mazzelli and Milazzo2016), who used the correlation proposed by Papamoschou (Reference Papamoschou1993) to compute the shear stress between the streams. The authors used the resulting compound flow equations with momentum exchange in a 1-D ejector model but did not investigate the choking mechanism. Earlier work by Otis (Reference Otis1976) treats the mixing and choking of a compound flow for constant-area ducts. An important conclusion of that work is that additional effects such as forces or heat exchange do not affect the choking condition as derived in the original compound theory. To the best of the author's knowledge, a general treatment of compound choking with momentum exchange in variable area ducts has not yet been proposed in the literature.
This work proposes such an extension and investigates the effect of friction between the streams and the walls on the choking mechanism. It is analytically shown that the sonic section moves upstream or downstream depending on the relative magnitude of the interstream and wall friction. The analysis is similar to ‘generalised 1-D flow’ (Shapiro Reference Shapiro1953, Chapter 8), which describes quasi-1-D compressible flow of a single stream with any combination of area change, friction, heat exchange and/or mass injection. Similarly, the isentropic compound flow theory, which only includes area change, is generalised in this work by including momentum exchanges. Our contribution is both theoretical and numerical: (i) we propose a simple theoretical model that captures the position of the sonic section, and (ii) we propose a quasi-third-order approximation of the pressure gradient in the vicinity of the sonic section to circumvent the numerical issues near the singular sonic point.
The governing equations are derived in § 2 and converted into an explicit system of ordinary differential equations (ODEs). Section 3 goes into more detail on the singularity of the equation for the static pressure at the sonic section. The numerical solution of the model is discussed in § 4. Numerical test cases are presented in § 5, to compare the results with the 1-D predictions in § 6. Section 7 closes the paper with conclusions and perspectives.
2. Model definition
The governing equations are derived in their general form in § 2.1. Section 2.2 introduces the treatment of the friction forces in the equation governing the streamwise pressure gradient. The corresponding choking mechanism is analysed in § 2.3.
2.1. Governing equations
Consider an internal flow consisting of two parallel streams (cf. figure 1) of the same ideal gas. The theory can be extended to any number of streams, but this generalisation is left for future work. The streams are modelled as uniform 1-D streams, i.e. in terms of cross-stream averaged quantities. The main assumption for compound-compressible flows is that the static pressure is uniform in all cross-sections, hence
\begin{equation} p_p(x) = p_s(x) = p(x), \end{equation}
where 
$p$ denotes the static pressure, 
$x$ is the axial coordinate and the subscripts 
$p$ and 
$s$ refer to the primary and secondary stream. The distribution of the cross-sections 
$A_i(x)$, with 
$i \in [p,s]$, is initially unknown and should be calculated. These are subject to a geometrical constraint imposed by the cross-section of the channel 
$A(x)$:
\begin{equation} \sum_{i\in[p,s]} A_i(x) = A(x). \end{equation}The governing equations consist of the conservation of mass, momentum and energy in each stream. For a steady, quasi-1-D flow, these read as follows:
$$\begin{gather} \frac{{\rm d}}{{\rm d}\kern0.7pt x} (\rho_i A_i u_i) = 0, \end{gather}$$
$$\begin{gather}\frac{{\rm d}}{{\rm d}\kern0.7pt x} (\rho_i A_i u_i^2) ={-}A_i \frac{{\rm d}p_i}{{\rm d}\kern0.7pt x} + F_i, \end{gather}$$
$$\begin{gather}\frac{{\rm d}}{{\rm d}\kern0.7pt x} (\rho_i A_i u_i h_{t,i}) = 0, \end{gather}$$
where 
$\rho$ denotes the density, 
$u$ denotes the axial velocity, 
$F_i$ is the net force per unit length exerted on stream 
$i$ and 
$h_t$ denotes the total enthalpy. The force 
$F_i$ accounts for the friction between the streams and at the wall, and is defined in § 2.2. The equations imply that the streams do not exchange mass and that the separating line between the streams is the dividing streamline. Heat exchange is not considered.
Figure 1. The considered axisymmetric nozzle has a contour 
$A(x)$ and contains two streams, which have the same total temperature 
$T_t$, but can have a different total pressure 
$p_t$. They exchange momentum through a friction force 
$F_{ps}$ at the dividing streamline and a wall friction force 
$F_w$. The static pressure 
$p$ is assumed to be uniform in each cross-section. The solution consists of the streamwise distribution of the flow variables 
$p$, 
$p_{t,p}$, 
$p_{t,s}$, 
$T_{t,p}$, 
$T_{t,s}$ and the individual cross-sections 
$A_p$ and 
$A_s$.
This system of equations can be converted to an explicit system of ODEs for numerical integration. The derivation can be found in classic textbooks; for example in Chapter 8 of Shapiro (Reference Shapiro1953). The equations in terms of the static pressure 
$p$, the total pressure 
$p_t$ and the total temperature 
$T_t$ read as follows:
$$\begin{gather} \frac{1}{p_i} \frac{{\rm d}p_i}{{\rm d}\kern0.7pt x} = \left[\frac{1 + \left(\gamma - 1\right)\textit{Ma}_i^2}{1-\textit{Ma}_i^2}\right] \frac{F_i}{A_i p_i} + \left[\frac{\gamma \textit{Ma}_i^2}{1-\textit{Ma}_i^2}\right]\frac{1}{A_i}\frac{{\rm d} A_i}{{\rm d}\kern0.7pt x}, \end{gather}$$
$$\begin{gather}\frac{1}{p_{t,i}}\frac{{\rm d}p_{t,i}}{{\rm d}\kern0.7pt x} = \frac{F_i}{A_i p_i}, \end{gather}$$
$$\begin{gather}\frac{1}{T_{t,i}}\frac{{\rm d}T_{t,i}}{{\rm d}\kern0.7pt x} = 0, \end{gather}$$
having introduced the ratio of specific heats 
$\gamma$ and the Mach numbers 
$\textit {Ma}_i$,
\begin{equation} \textit{Ma}_i = u_i / a_i \quad \mbox{and} \quad a_i = \sqrt{\gamma R T_i} , \end{equation}
where 
$R$ denotes the specific gas constant. Equations (2.6)–(2.8) do not form a closed system of equations because the cross-sections 
$A_i(x)$ are unknown (as opposed to a single stream where 
$A_i = A$). Nonetheless, their sum should equal the known channel's cross-section 
$A(x)$ (2.2). The derivative 
${\rm d}A/{{\rm d}\kern0.7pt x}$ can be computed from the channel's profile, so the unknown terms 
${\rm d}A_i/{{\rm d}\kern0.7pt x}$ can be eliminated from the system by summing over the streams 
$i$. The gradient of the cross-section 
$A_i$ is first isolated from (2.6),
\begin{equation} \frac{{\rm d}A_i}{{\rm d}\kern0.7pt x} = \left[A_i\frac{1-\textit{Ma}_i^2}{\gamma \textit{Ma}_i^2}\right]\frac{1}{p} \frac{{\rm d}p}{{\rm d}\kern0.7pt x} - \left[\frac{1 + \left(\gamma - 1\right)\textit{Ma}_i^2}{\gamma \textit{Ma}_i^2}\right] \frac{F_i}{p}, \end{equation}
where the index for the pressure is omitted because of (2.1). Note that (2.10) can be integrated to compute the cross-section of a single stream if the pressure gradient is known. Summing over the streams 
$i$ gives
\begin{equation} \sum_{i\in[p,s]} \frac{{\rm d}A_i}{{\rm d}\kern0.7pt x} = \frac{{\rm d}A}{{\rm d}\kern0.7pt x} = \sum_{i\in[p,s]} \left[A_i\frac{1-\textit{Ma}_i^2}{\gamma \textit{Ma}_i^2}\right]\frac{1}{p} \frac{{\rm d}p}{{\rm d}\kern0.7pt x} - \sum_{i\in[p,s]}\left[\frac{1 + \left(\gamma - 1\right)\textit{Ma}_i^2}{\gamma \textit{Ma}_i^2}\right] \frac{F_i}{p}. \end{equation}The pressure gradient is the only remaining unknown, for which an explicit equation is obtained from (2.11),
\begin{equation} \frac{1}{p} \frac{{\rm d}p}{{\rm d}\kern0.7pt x} = \dfrac{1}{\beta}\left( \dfrac{{\rm d}A}{{\rm d}\kern0.7pt x} + \sum_{i\in[p,s]}\left[\dfrac{1 + \left(\gamma - 1\right)\textit{Ma}_i^2}{\gamma \textit{Ma}_i^2}\right] \dfrac{F_i}{p}\right), \end{equation}
where 
$\beta$ denotes the compound choking indicator,
\begin{equation} \beta = \sum_{i\in[p,s]} A_i\dfrac{1-\textit{Ma}_i^2}{\gamma \textit{Ma}_i^2}. \end{equation} Bernstein et al. (Reference Bernstein, Heiser and Hevenor1967) showed that 
$\beta > 0$ corresponds to a compound-subsonic flow, 
$\beta =0$ to a compound-sonic flow and 
$\beta < 0$ to a compound-supersonic flow. As a result, the pressure gradient and the area gradient of a compound-subsonic flow share the same sign in the absence of forces, leading to an expansion in a convergent channel and a compression in a divergent channel. Recently, Metsue et al. (Reference Metsue, Debroeyer, Poncet and Bartosiewicz2021) extended the compound choking condition to real gases and showed that the condition 
$\beta =0$ corresponds to the maximal combined mass flow rate.
A convenient variable to identify the compound flow regime is the so-called equivalent Mach number 
$\textit {Ma}_{eq}$, introduced by Hedges & Hill (Reference Hedges and Hill1974),
\begin{equation} \textit{Ma}_{eq} = \left(\gamma \frac{\beta}{A}+1\right)^{-{1}/{2}}. \end{equation} This equals one for a compound choked flow. Section 2.3 analyses choking in more detail. It is important to note that the forces do not affect 
$\beta$ and hence the compound choking criterion, as pointed out by Otis (Reference Otis1976). Moreover, as shown by Bernstein et al. (Reference Bernstein, Heiser and Hevenor1967), compound waves can no longer travel upstream of a location 
$x$ if 
$\beta (x) \leqslant 0$, which constitutes the fundamental principle of choking as in the case of a single stream.
Equations (2.7), (2.8), (2.10) and (2.12) form a closed system for the (uniform) static pressure 
$p$ and the total pressure 
$p_{t,i}$, the total temperature 
$T_{t,i}$ and the cross-section 
$A_i$ in each stream. These variables fully define the state of the flow and other variables can easily be calculated from them using classic gas dynamic relations. The explicit nature of the equations facilitates numerical integration with a shooting method. More details on the numerical solution are provided in § 4.
2.2. Model closure
The compound flow is assumed to be axisymmetric, with the primary stream surrounded by the secondary stream. Planar flows are left as technical extensions for future studies. Hence, wall friction acts only on the secondary stream. A friction coefficient 
$f_w$ is defined such that the friction force per unit length 
$F_w$ equals
\begin{equation} F_w = \tfrac{1}{2} f_w \rho_s u_s^2 l_w = \tfrac{1}{2} f_w \gamma p \textit{Ma}_s^2 l_w, \end{equation}
where 
$l_w$ denotes the wetted perimeter of the cross-section. In a circular duct, 
$l_w$ is related to the cross-section 
$A$ as follows:
\begin{equation} l_w = 2 {\rm \pi}r = 2 \sqrt{{\rm \pi} A}, \end{equation}
where 
$r$ denotes the radius of the total cross-section 
$A$. Similarly, the friction between the streams is defined as (Papamoschou Reference Papamoschou1993)
\begin{equation} F_{ps} = \frac{1}{2} f_{ps} \dfrac{\rho_p+\rho_s}{2} \left(u_p - u_s\right) \left|u_p - u_s\right| l_{ps}, \end{equation}
where 
$f_{ps}$ denotes a friction coefficient and 
$l_{ps}$ denotes the perimeter of the primary stream,
\begin{equation} l_{ps} = 2 {\rm \pi}r_p = 2 \sqrt{{\rm \pi} A_p}. \end{equation} The friction force between the streams opposes the stream with the highest velocity and is defined to be positive by (2.17) if 
$u_p > u_s$. Therefore, 
$F_{ps}$ opposes the primary stream, while the secondary stream is entrained by the interstream friction and decelerated by the wall friction,
\begin{equation} F_p ={-} F_{ps} \quad \mbox{and} \quad F_s = F_{ps} - F_w. \end{equation} Two approaches were used in this work for both coefficients 
$f_w$ and 
$f_{ps}$ and tested to the RANS simulations in § 6.1.
The first approach was to leverage known correlations for similar flow configurations. The formula of Van Driest (Reference Van Driest1951) for the local wall friction coefficient 
$f_w$ under a turbulent boundary layer in a compressible flow gives
\begin{equation} \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{equation}where
\begin{equation} 1-\lambda^2 = \left(1 + \frac{\gamma - 1}{2} \textit{Ma}_s^2\right)^{{-}1}, \quad Re_x = \dfrac{\rho_s u_s x}{\mu_s} \quad \mbox{and}\quad \theta = \dfrac{S}{T_s}, \end{equation}
with the dynamic viscosity 
$\mu$ computed with the law of Sutherland (Reference Sutherland1893), a reference viscosity 
$\mu _{ref}$ of 
$1.716 \times 10^{-5}$ Pa s at 
$T_{ref} = 273.2$ K, and 
$S = 110.4$ K. For the interstream friction coefficient 
$f_{ps}$, the formula of Papamoschou (Reference Papamoschou1993) gives
\begin{equation} f_{ps} = 0.013 \dfrac{(1 + u_s/u_p) (1 + \sqrt{\rho_s / \rho_p})}{1 + (u_s/u_p) \sqrt{\rho_s / \rho_p}} \left(0.25 + 0.75 \exp({-}3 \textit{Ma}_c^2)\right), \end{equation}
with all quantities evaluated locally and the local convective Mach number 
$\textit {Ma}_c$ is defined as follows:
\begin{equation} \textit{Ma}_c = \dfrac{u_p - u_s}{a_p + a_s}. \end{equation} The second approach was to infer these coefficients from the RANS results using an inverse method. In this setting, both coefficients 
$f_w$ and 
$f_{ps}$ are assumed to be constant and are determined by an optimiser which seeks to minimise the 
$l_2$ norm of the discrepancy between the proposed model and the results of RANS simulations. This approach is expected to be more accurate due to the calibration, provided that the assumption of constant coefficient is adequate. On the other hand, these coefficients are not expected to generalise well outside the reference data, while the correlations (2.20) and (2.22) make the model self-sufficient and closed.
2.3. Choking mechanism
Substituting the forces 
$F_p$ and 
$F_s$ as defined in the section above in (2.12) leads to a rational equation for the static pressure gradient,
\begin{equation} \frac{1}{p} \frac{{\rm d}p}{{\rm d}\kern0.7pt x} = \dfrac{N}{\beta}, \end{equation}
with 
$\beta$ defined in (2.13) and the numerator,
\begin{equation} N = \dfrac{{\rm d}A}{{\rm d}\kern0.7pt x} + N_{ps} + N_w, \end{equation}where
$$\begin{gather} N_{ps} = \dfrac{f_{ps}l_{ps}}{2}\frac{\textit{Ma}_p^2 - \textit{Ma}_s^2}{\gamma \textit{Ma}_p^2 \textit{Ma}_s^2} \dfrac{\rho_p+\rho_s}{2} \left(u_p - u_s\right) \left|u_p - u_s\right| \dfrac{1}{p}, \end{gather}$$
$$\begin{gather}N_w ={-} \dfrac{f_{w} l_{w}}{2} \left(1 + \left(\gamma - 1\right)\textit{Ma}_s^2\right). \end{gather}$$ The term 
$N_w$ stems from wall friction ( 
$f_w$) and is negative by definition. The term 
$N_{ps}$ is due to friction between the streams and it is positive unless the terms 
$\textit {Ma}_p^2-\textit {Ma}_s^2$ and 
$u_p-u_s$ have a different sign. In the natural case 
$u_p>u_s$, this is possible only if 
$T_p/T_s \geqslant u^2_p/u^2_s$. This scenario is not encountered in Lamberts et al. (Reference Lamberts, Chatelain, Bourgeois and Bartosiewicz2018) and Kracik & Dvorak (Reference Kracik and Dvorak2023) and thus not discussed further. Therefore, the interstream friction has a positive contribution to the numerator 
$N$ in (2.24) in the considered conditions and the wall friction has a negative contribution. This leads to some fundamental implications.
A compound-subsonic flow (
$\beta > 0$) expands in a convergent section (
${\rm d}A/{{\rm d}\kern0.7pt x} < 0$) and under the action of wall friction (as in a classic Fanno flow). Conversely, it recovers pressure in a divergent section and under the action of interstream friction. This force opposes the stream with the highest velocity, seeking to equalise velocities between streams. Consequently, momentum exchange between the streams promotes flow uniformity and results in a net increase in static pressure. The opposite is true for a compound-supersonic flow.
Equation (2.24) is singular in the sonic section (
$\beta = 0$), similarly to a single stream where the denominator equals 
$1 - \textit {Ma}^2$ (cf. Shapiro Reference Shapiro1953; Restrepo & Simões-Moreira Reference Restrepo and Simões-Moreira2022). Nevertheless, classic nozzle flows and compound flows expand continuously through the sonic section, meaning that their pressure gradient is finite at 
$\beta =0$. This is possible according to (2.24) only if both the numerator 
$N$ and the denominator 
$\beta$ tend to zero, leading to an indeterminate form with a finite value. This observation is corroborated by (2.11), which simplifies to 
$N=0$ if 
$\beta =0$. In isentropic conditions ( 
$f_{ps} = f_w = 0$), a compound flow becomes sonic at the throat (
${\rm d}A/{{\rm d}\kern0.7pt x} = 0$), as shown by Bernstein et al. (Reference Bernstein, Heiser and Hevenor1967). In non-isentropic conditions, the relative magnitudes of 
$N_{ps}$ and 
$N_w$ in (2.25) determine the position of the sonic section.
Wall friction has a negative contribution to the numerator. Therefore, in the absence of interstream friction ( 
$f_{ps}=0$), one must have 
${\rm d}A/{{\rm d}\kern0.7pt x}>0$ to respect 
$N=0$. This implies that the sonic section is located in a divergent section, hence downstream of the throat. The wall friction has the same effect on a nozzle flow with a single stream (Shapiro Reference Shapiro1953; Beans Reference Beans1970; Restrepo & Simões-Moreira Reference Restrepo and Simões-Moreira2022). Likewise, in the absence of wall friction ( 
$f_w=0$), one must have 
${\rm d}A/{{\rm d}\kern0.7pt x} < 0$ to respect 
$N=0$. This implies that the interstream friction pushes the sonic section upstream of the throat. This analysis offers a plausible explanation for the observed displacement of the sonic section by Lamberts et al. (Reference Lamberts, Chatelain, Bourgeois and Bartosiewicz2018) and Kracik & Dvorak (Reference Kracik and Dvorak2023). The interplay between both forces is investigated further in § 6.
The indeterminate pressure gradient can be computed by applying l'Hôpital's rule on (2.24) (cf. Bernstein et al. Reference Bernstein, Heiser and Hevenor1967; Restrepo & Simões-Moreira Reference Restrepo and Simões-Moreira2022),
\begin{equation} \left(\frac{1}{p}\frac{{\rm d}p}{{\rm d}\kern0.7pt x}\right)^* = \frac{N^*}{\beta^*} = \frac{0}{0} = \frac{({\rm d}N/{{\rm d}\kern0.7pt x})^*}{({\rm d}\beta / {{\rm d}\kern0.7pt x})^*}, \end{equation}
where the superscript 
$*$ refers to the sonic section. The derivatives 
${\rm d}N/{{\rm d}\kern0.7pt x}$ and 
${\rm d}\beta /{{\rm d}\kern0.7pt x}$ give rise to the derivatives of Mach numbers 
$\textit {Ma}_i$, cross-sections 
$A_i$, forces 
$F_i$ and the static pressure 
$p$ (see (2.25) and (2.13)), and these can be all related to the local pressure gradient. The relevant equations are presented in Appendix A.
For example, the gradient of 
$\beta$ is given by
\begin{equation} \dfrac{{\rm d}\beta}{{\rm d}\kern0.7pt x} = \sum_{i\in[p,s]} \left(\dfrac{{\rm d}A_i}{{\rm d}\kern0.7pt x}\dfrac{1-\textit{Ma}_i^2}{\gamma \textit{Ma}_i^2} - \dfrac{A_i}{\gamma \textit{Ma}_i^4}\dfrac{{\rm d}\textit{Ma}_i^2}{{\rm d}\kern0.7pt x}\right). \end{equation} Writing the Mach numbers 
$\textit {Ma}_i$ as a function of the total and the static pressure,
\begin{equation} \textit{Ma}_i^2 = \frac{2}{\gamma-1}\left(\left(\frac{p_{t,i}}{p}\right)^{(\gamma-1)/{\gamma}}-1\right), \end{equation}
the derivatives 
${\rm d}\textit {Ma}_i^2/{{\rm d}\kern0.7pt x}$ in (2.29) can be substituted by the gradients of the static and total pressures. The resulting expression for the gradient of 
$\beta$ thus contains first derivatives of the static pressure, the cross-sections 
$A_i$ and the total pressures 
$p_{t,i}$. The latter two can be substituted using (2.7) and (2.10), which contain the static pressure gradient. The resulting expression for the gradient of 
$\beta$ is linear with respect to the static pressure gradient, with nonlinear coefficients depending on the Mach numbers, forces and static pressure,
\begin{equation} \dfrac{{\rm d}\beta}{{\rm d}\kern0.7pt x} = c_{\beta0}(\textit{Ma}_p,\textit{Ma}_s)\left(\dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right) + c_{\beta1}(\textit{Ma}_p,\textit{Ma}_s,F_p,F_s,p). \end{equation} The analytical expressions of 
$c_{\beta 0}$ and 
$c_{\beta 1}$ are given by (A26) and (A27). The same approach leads to a similar linear equation for the derivative of the numerator (see (A21)). These equations are derived for constant friction coefficients, i.e. the correlations (2.20) and (2.22) are not explicitly differentiated. This approximation is deemed acceptable since the distribution of 
$f_w$ and 
$f_{ps}$ are relatively flat, as later shown in § 6. The result is an implicit equation for the static pressure gradient in the sonic section,
\begin{equation} \left(\frac{1}{p}\frac{{\rm d}p}{{\rm d}\kern0.7pt x}\right)^* = \dfrac{c_{N0}^* \left(\dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right)^* + c_{N1}^*}{c_{\beta0}^*\left(\dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right)^* + c_{\beta1}^*}. \end{equation} Rearranging the equation above leads to the quadratic equation (A28) for the static pressure gradient. The two roots correspond to the subsonic and the supersonic solution downstream of the sonic point (cf. Shapiro Reference Shapiro1953). For an isentropic compound flow, the derivative of the numerator corresponds to the second derivative of the cross-section 
${\rm d}^2A/{{\rm d}\kern0.7pt x}^2$ ( 
$f_{ps}=f_w=0$), leading to the tractable expression reported by Bernstein et al. (Reference Bernstein, Heiser and Hevenor1967). The additional terms introduced by the friction forces preclude such a concise expression, which is why the coefficients and their analytical expressions are reported in Appendix A. Note, however, that these equations are exact and analytical.
3. Approximations to circumvent the singularity at the sonic section
The pressure gradient can be computed in any section using (2.32) if the flow is compound-sonic or by using (2.24) otherwise. However, (2.24) is numerically ill-conditioned near the sonic section. Therefore, approximations of the pressure gradient near the sonic section help stabilising the numerical integration. In this section, such approximations are developed analytically.
L’Hôpital's rule used in § 2.3 arises from a Taylor series expansion of the numerator and the denominator around the sonic section 
$x^*$, which is then evaluated at the sonic section 
$x^*$ itself. More generally, these expansions lead to the following equation in 
$x$, which is still exact,
\begin{equation} \frac{1}{p}\frac{{\rm d}p}{{\rm d}\kern0.7pt x} = \frac{N}{\beta} = \frac{N^* + (x-x^*) \dfrac{{\rm d}N}{{\rm d}\kern0.7pt x}^* + \dfrac{1}{2} (x-x^*)^2 \dfrac{{\rm d}^2N}{{{\rm d}\kern0.7pt x}^2}^* + \dfrac{1}{6} (x-x^*)^3 \dfrac{{\rm d}^3N}{{{\rm d}\kern0.7pt x}^3}^* + \dots}{\beta^* + (x-x^*) \dfrac{{\rm d}\beta}{{\rm d}\kern0.7pt x}^* + \dfrac{1}{2} (x-x^*)^2 \dfrac{{\rm d}^2\beta}{{{\rm d}\kern0.7pt x}^2}^* + \dfrac{1}{6} (x-x^*)^3 \dfrac{{\rm d}^3\beta}{{{\rm d}\kern0.7pt x}^3}^* + \dots}. \end{equation} Since 
$N^* = \beta ^* = 0$, the equation above simplifies to the following expression:
\begin{equation} \frac{1}{p}\frac{{\rm d}p}{{\rm d}\kern0.7pt x} = \frac{N}{\beta} = \frac{\dfrac{{\rm d}N}{{\rm d}\kern0.7pt x}^* + \dfrac{1}{2} (x-x^*) \dfrac{{\rm d}^2N}{{{\rm d}\kern0.7pt x}^2}^* + \dfrac{1}{6} (x-x^*)^2 \dfrac{{\rm d}^3N}{{{\rm d}\kern0.7pt x}^3}^* + \dots}{\dfrac{{\rm d}\beta}{{\rm d}\kern0.7pt x}^* + \dfrac{1}{2} (x-x^*) \dfrac{{\rm d}^2\beta}{{{\rm d}\kern0.7pt x}^2}^* + \dfrac{1}{6} (x-x^*)^2 \dfrac{{\rm d}^3\beta}{{{\rm d}\kern0.7pt x}^3}^* + \dots}. \end{equation}Truncating the Taylor series leads to approximations of arbitrary order in space. The first-order approximation only retains the first derivatives, which are constant since they are evaluated at the sonic section. Hence, the pressure gradient in the vicinity of the sonic section is approximated as a constant, equal to its value in the sonic section. Note that the resulting expression is identical to (2.28) and thus leads to the quadratic (2.32) in § 2.3.
Higher accuracy is obtained by retaining higher-order derivatives. The analytical expressions for the second derivatives are given by (B32) and (B42). The derivation could continue up to the third derivative, but to avoid overly cumbersome expressions, this work proceeds with a third-order approximation by numerically differentiating the analytical second derivatives.
Numerically, one has to decide when to switch from the direct equation (2.24) to the approximation (3.2). This can be achieved by imposing a tolerance on 
$\beta$ (or 
$\textit {Ma}_{eq}$), e.g. using the approximation when 
$-\epsilon < \beta < \epsilon$, with 
$\epsilon$ an arbitrary small positive number (cf. Restrepo & Simões-Moreira Reference Restrepo and Simões-Moreira2022). An alternative condition is based on an estimate of the truncation error in the Taylor series, which can be estimated as the magnitude of the first term after the truncation. The resulting truncation errors on the numerator and the denominator, denoted by 
$\delta N$ and 
$\delta \beta$, can be propagated to the pressure gradient as follows:
\begin{equation} \delta \left(\frac{1}{p}\frac{{\rm d}p}{{\rm d}\kern0.7pt x}\right) = \sqrt{\left(\frac{1}{\beta}\right)^2\left(\delta N\right)^2 + \left(-\frac{N}{\beta^2}\right)^2\left(\delta \beta\right)^2}. \end{equation}
Equation (3.3) is analytically available for the first-order approximation using the second derivatives for 
$\delta N$ and 
$\delta \beta$ (cf. Appendix B). The second-order approximation requires third derivatives to estimate its error, which can be computed with finite differences (as in the quasi-third-order approximation). The switch from the direct equation (2.24) to the approximation (3.2) is based on the relative error on the static pressure gradient,
\begin{equation} \frac{1}{p}\frac{{\rm d}p}{{\rm d}\kern0.7pt x} = \begin{cases} \textrm{(3.2)} & \textrm{if } \dfrac{\delta \left(\dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right)}{\left(\dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right)} \leqslant 5\,\% , \\[10pt] \textrm{(2.24)} & \textrm{otherwise.} \end{cases} \end{equation}Examples of the distribution of the pressure gradient using the approximations of increasing order are shown in figure 2. The higher-order approximations are more accurate in a wider region around the sonic section. The estimated truncation error for the quasi-third-order approximation is not shown. The compound equations have been integrated using the blended equation (3.4) for the pressure. More details on the numerical integration are provided in § 4. The blended function (3.4) matches the direct equation (2.24) at the inlet and switches to the approximation (3.2) as the truncation error, computed with (3.3), drops below 5 %. The flow then passes through the sonic section using the approximation (3.2), effectively avoiding the indeterminate form (2.28) of the direct equation. Downstream, the direct equation is integrated when the threshold on the truncation error is exceeded again.
Figure 2. The pressure gradient computed directly with (2.24), with the approximation (3.2) (truncated after the first, second and third derivatives) and with the blended function (3.4), the shaded region indicates the truncation error estimated by (3.3), which is unavailable for the quasi-third-order approximation: (a) first order; (b) first order (zoom); (c) second order; (d) second order (zoom); (e) quasi-third-order; ( f) quasi-third-order (zoom).
In this study, the quasi-third-order approximation is used with a backward second-order finite difference formula, with a threshold of 5 % on the second-order approximation, since the fourth derivatives needed for the truncation error are not available.
4. Numerical solution
The governing equations (2.7), (2.8), (2.10), (2.24), with the approximation (3.2) near the sonic point, form a system of ODEs in function of the axial coordinate 
$x$, which can thus be integrated in space. The fourth-order Runge–Kutta method is used in this study, which is readily available in Python with the ‘solve_ivp’ function in SciPy (https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html). In this study, the mesh consists of 100 equally distributed points, which has been found sufficient to accurately respect the conservation equations.
The boundary conditions consist of the total pressures 
$p_{t,p}, p_{t,s}$, total temperatures 
$T_{t,p}, T_{t,s}$ and the cross-sections 
$A_{p}, A_{s}$ at the inlet (
$x=0$) and of the static pressure 
$p_b$ at the outlet (
$x=L$), as indicated in figure 1. The distribution of the channel's cross-section 
$A(x)$ is prescribed. Integrating the system of ODEs requires an ‘initial’ condition at the inlet or at any point from which the integration proceeds forward in 
$x$. In a single stream computation, the sonic condition 
$N^*=0$ is used to identify the sonic section and split the integration for the Mach number from the inlet to the sonic section and from the sonic section to the outlet (as proposed by Beans (Reference Beans1970) and Restrepo & Simões-Moreira (Reference Restrepo and Simões-Moreira2022)).
This is not a viable option for compound flows since the sonic state is not uniquely defined (as 
$\textit {Ma}=1$), but can consist of various combinations of cross-sections 
$A_i$ and Mach numbers 
$\textit {Ma}_i$. The system state is available at the inlet except for the static pressure 
$p(0)$. This was found with a shooting method until the boundary condition at the outlet was respected. Firstly, bounds for the inlet pressure are established, considering that 
$p(0)$ cannot exceed the total pressure of any stream by definition, that is the upper bound of 
$p(0)$ is
\begin{equation} p_{max} = \min\left(p_{t,p}(0), p_{t,s}(0)\right). \end{equation} The compound flow is assumed to be subsonic at the inlet (as in a classic nozzle flow), so the second limiting condition is a sonic flow at the inlet (
$\beta (0) = 0)$. Using (2.13) and (2.30), the sonic condition at the inlet leads to the lower bound 
$p_{min}$ for the static pressure at the inlet,
\begin{equation} \beta(0) = 0 = \sum_{i\in[p,s]}\dfrac{A_{i}}{\gamma}\left(\left[\frac{2}{\gamma-1}\left(\left(\frac{p_{t,i}(0)}{p_{min}}\right)^{(\gamma-1)/{\gamma}}-1\right)\right]^{{-}1} - 1\right), \end{equation}
where 
$p_{min}$ can be found iteratively. The bounds 
$[p_{min}, p_{max}]$ allow us to find the correct inlet pressure 
$p(0)$ with a bisection method, which iteratively restricts the bounds of the interval as depicted in figure 3. In each iteration, the governing equations are integrated using the mean of the current interval, after which either its bottom or top half is used in the next iteration. This update depends on whether the computed flow is choked.
Figure 3. Flowchart of the numerical procedure to solve the compound nozzle flow. The static pressure 
$p(0)$ at the inlet is found with a bisection method. A subsonic solution is accepted if the computed back pressure 
$p(L)$ matches the boundary condition 
$p_b$. A choked flow requires both 
$\beta$ and the numerator 
$N$ to equal zero (see (2.28)). The approximation (3.2) is computed after convergence of the choked procedure, which is then repeated, replacing (2.24) by (3.4).
In the subsonic case, the inlet pressure 
$p(0)$ and the computed outlet pressure 
$p(L)$ are positively correlated: a lower inlet pressure leads to a lower outlet pressure. If the computed back pressure is too low with respect to the boundary condition (
$p(L) < p_b$), the inlet pressure should increase and the average pressure 
$p(0)$ becomes the new lower bound of the bracket in the next iteration. The opposite is true if 
$p(L) < p_b$. If the pressures 
$p(L)$ and 
$p_b$ match, the solution has been found and the flow is compound-subsonic.
If the computed flow becomes sonic (
$\beta (x) = 0$), the integration is stopped locally at the sonic point and the numerator 
$N$ is computed with (2.25). If the condition 
$N=0$ is not respected (see (2.28)), the inlet pressure 
$p(0)$ must be updated. A negative numerator indicates that a subsonic flow tends to expand. If this occurs in the sonic section, the sonic point is reached too far upstream, and hence the inlet pressure should be increased. Considering, for example, an isentropic compound flow ( 
$f_{ps}=f_w=0$), where the numerator 
$N$ is reduced to the gradient 
${\rm d}A/{{\rm d}\kern0.7pt x}$, it is known that the flow becomes sonic at the throat (
${\rm d}A/{{\rm d}\kern0.7pt x} = 0$). The case 
$N= {\rm d}A/{{\rm d}\kern0.7pt x} < 0$ at the sonic section implies that the sonic point is reached in a convergent section; hence, the inlet pressure should be increased to avoid having a non-physical mass flow rate exceeding the choked one. The same argument holds for a non-isentropic flow.
If 
$N=0$, the flow is sonic in the correct section, and the approximation (3.4) can be computed with the desired order of accuracy. This slightly modifies the governing equations, so the choked solution also changes. Therefore, the procedure to find the inlet pressure 
$p(0)$ is repeated, replacing (2.24) with (3.4).
The resulting choked solution is available from the inlet to the sonic section (
$\beta (x) = 0$). The governing equations are then integrated up to the outlet, using the positive root 
$({\rm d}p/{{\rm d}\kern0.7pt x})^*$ of (2.32) to obtain the subsonic solution. The complete solution is thus subsonic throughout, except at the sonic point. This limiting case corresponds to the highest back pressure which chokes the flow, which is commonly referred to as the critical back pressure 
$p_b^*$ in ejector modelling. The final step checks whether the choked solution respects the boundary condition 
$p_b$. The flow remains choked if the imposed back pressure 
$p_b$ is lower than the critical value 
$p_b^*$. Otherwise, the static pressure at the inlet should increase, leading to a subsonic solution.
The supersonic solution can be obtained with the negative root 
$({\rm d}p/{{\rm d}\kern0.7pt x})^*$ of (2.32). If the back pressure 
$p_b$ lies between the pressures corresponding to the two continuous choked solutions, a shock appears in the diffuser as in the case of a nozzle flow with a single stream. The shock position can be found iteratively by imposing the computed back pressure 
$p(L)$ to match 
$p_b$ (see e.g. Restrepo & Simões-Moreira Reference Restrepo and Simões-Moreira2022). However, to the authors’ knowledge, a compound shock model has not been presented in the literature.
5. Validation test cases with axisymmetric simulations
The predictions of the 1-D theoretical model are compared with postprocessed axisymmetric RANS simulations of a compound nozzle in § 6. The numerical set-up of these two-dimensional (2-D) simulations is detailed in § 5.1, and the test cases are described in § 5.2.
5.1. Numerical set-up for the RANS validation
The computational domain is shown in figure 4. It consists of an axisymmetric convergent–divergent nozzle with two distinct inlets. The simulation was carried out with a wedge-shaped domain of five degrees and a single cell in the tangential direction. The contour of the upper wall was defined as a cosine between 
${\rm \pi} /2$ and 
$2{\rm \pi}$ with appropriate scaling to reach the dimensions shown in figure 4:
\begin{equation} r(x) = \dfrac{R_o + R_{th}}{2} + \dfrac{R_o - R_{th}}{2} \cos\left(\dfrac{3{\rm \pi}}{2} \dfrac{x}{L} + \dfrac{\rm \pi}{2}\right). \end{equation}
Figure 4. The computational domain consists of an axisymmetric convergent–divergent nozzle with two distinct inlets. The primary stream (index 
$p$) is surrounded by the secondary stream (index 
$s$). The reference geometry is defined by 
$L = 187.5$ mm, 
$R_i = 9.5$ mm, 
$R_{th} = 9.0$ mm, 
$R_o = 10$ mm and 
$R_p = 0.5 R_i$. The model validation considered various throat 
$R_{th}$ and the primary inlet 
$R_p$ radii, as described in § 5.2, while keeping the same resolution at the wall. The coarse mesh is displayed for visibility.
Hence, the radius at the inlet (
$x=0$) follows from the outlet radius 
$R_o$, the throat radius 
$R_{th}$ and the length 
$L$. The throat is located at one-third of the domain. The analysis of the sonic section displacement was focused on a test case with nozzle lengths of 187.5 mm, outlet radius 
$R_o=10$ mm, inlet radius 
$R_i = 9.5$ mm, throat radios 
$R_{th} = 9$ mm and 
$R_p = 0.5 R_i$. For the model validation, several other geometries were considered as reported in more details in § 5.2.
The mesh was generated using Gmsh (https://gmsh.info/) coupled to Python to allow full control of the number of elements and of the refinement near the walls. The simulations were performed with an adapted version of the compressible transient solver rhoCentralFoam in OpenFOAM v9 (https://openfoam.org/version/9/) with the 
$k-\omega$ shear stress transport turbulence model of Menter (Reference Menter1993). The adapted solver does not solve for the tangential component of the velocity, which remains equal to zero (no swirl). The law of Sutherland (Reference Sutherland1893) is used to compute the viscosity. The simulations are wall-resolved as shown in § 5.2.
The boundary conditions at the inlets consist of the total pressure and total temperature. The primary inlet is chosen to be supersonic, so static pressure is also imposed. A wave transmissive boundary condition is imposed at the exit as the flow is supersonic. The wall is adiabatic and treated with a slip or no-slip condition, depending on the test case. The ‘wedge’ boundary condition is imposed on the sides of the wedge to inform OpenFOAM of the axisymmetry.
The resulting 2-D flow fields are postprocessed by identifying the dividing streamline between the streams and by averaging over the resulting cross-sections. The radial position 
$r_{d}$ of the dividing streamline is found at any point 
$x$ by integrating the axial mass flux,
\begin{equation} \dot{m}_p(x) = \int_0^{r_{d}(x)} \rho(x,r) u(x,r) 2 {\rm \pi}r \, {\rm d}r \quad \mbox{and} \quad \dot{m}_s(x) = \int_{r_{d}(x)}^{R} \rho(x,r) u(x,r) 2 {\rm \pi}r \, {\rm d}r, \end{equation}
where 
$u$ denotes the axial component of the velocity vector. Knowing the radius 
$r_d(0)$ of the dividing streamline at the primary inlet, the primary mass flow rate is calculated. The evolution of the dividing streamline at any position 
$x$ follows from (5.2a,b). The second step consists of averaging the relevant flow variables across the ejector's section, i.e.
$$\begin{gather} \hat{\rho}_{i} (x) = \frac{1}{A_{i}}\int_{A_{i}}\rho \, {\rm d}A_{i}, \end{gather}$$
$$\begin{gather}\hat{u}_{i} (x) = \frac{1}{\hat{\rho}_{i} A_{i}}\int_{A_{i}}\rho u \, {\rm d}A_{i}, \end{gather}$$
$$\begin{gather}\hat{e}_{t,i} (x) = \frac{1}{\hat{\rho}_{i} A_{i}}\int_{A_{i}}\rho e_t \, {\rm d}A_{i}, \end{gather}$$
for 
$i \in [p,s]$ and where 
$e_t = c_v T + u^2/2$ denotes the total specific internal energy and 
$c_v$ is the specific heat at constant volume. The density-averaged velocity and internal energy ensure that the 1-D mass flow rate and internal energy match their integrated 2-D counterpart. The averaged static temperature, static pressure and Mach number follow from (5.3)–(5.5) and the specific gas constant 
$R$:
\begin{equation} \hat{T}_{i} = \frac{1}{c_v} \left(\hat{e}_{t,i}-\frac{1}{2}\hat{u}_{i}^2\right), \quad \hat{p}_i = \hat{\rho}_i R \hat{T}_i \quad \mbox{and} \quad \widehat{\textit{Ma}}_{i} = \frac{\hat{u}_{i}}{\sqrt{\gamma R \hat{T}_{i}}}.\end{equation}5.2. Selected test cases
Three simulations were carried out where the wall and interstream friction forces were numerically (de)-activated. The boundary conditions are listed in table 1. At both inlets, the turbulent mixing length was set equal to 7 % of the hydraulic diameter, as recommended by Versteeg (Reference Versteeg2007). The turbulence intensity was set to 1 %, as this value proved to best match the friction coefficients predicted by (2.22).
Table 1. Boundary conditions of the simulations presented in § 6. The primary inlet is supersonic with 
$\textit {Ma}_p=1.22$. They cover (i) an isentropic case, (ii) a case with friction between the streams but without wall friction and (ii) a case with both forces, respectively.
The first simulation is inviscid and thus serves as a reference for the original compound flow theory without momentum exchange. The second simulation is viscous with a slip condition at the wall: only friction forces between the streams are thus active in this case. Finally, wall friction is included in the third case by imposing a no-slip condition along the upper boundary. A mesh convergence study was carried out on the viscous simulation with the no-slip condition. The number of cells was doubled in each direction, leading to meshes with 
$400 \times 25$, 
$800 \times 52$ and 
$1600 \times 106$ cells. The resulting distributions of 
$y^+$ near the wall and the Mach number at the centreline are shown for the three meshes in figure 5. The medium and fine meshes show good agreement and are both wall-resolved. The wall friction Reynolds number 
$Re_{\tau } = \rho u_{\tau } D_h / \mu$ equals 
$1.9\times 10^4$ at the inlet and 
$5.3 \times 10^3$ at the outlet in case 3. The results on the fine mesh are used in the next section to profit from its higher resolution.
Figure 5. The coarse, medium and fine mesh contain 
$400 \times 25$, 
$800 \times 52$ and 
$1600 \times 106$ elements, respectively. The medium and fine meshes lead to a wall-resolved simulation and to good agreement for the Mach number at the centreline. The fine mesh is retained for the other simulations and the postprocessing. Here (a) 
$y^+$ of the first cell near the wall; (b) Mach number at the centreline.
A parametric study was carried out to validate the model in a wider range of conditions. The reference case is the viscous simulation with the no-slip condition at the wall (case 3 in table 1), as it is the most representative for practical applications. Three parameters were varied independently: (i) the total pressure at the secondary inlet between 
$p_{t,s} = 0.5$ and 3 bar; (ii) the radius of the throat between 
$R_{th}/R_o = 0.6$ and 0.9; and (iii) the radius of the primary inlet between 
$R_p/R_i = 0.3$ and 0.7 (see figure 4). The correlations (2.20) and (2.22) were used in all cases, so the model was not tuned further.
The resulting distributions of the static pressure and the Mach numbers depend on the geometry and the inlet conditions, so the static pressure imposed at the supersonic primary inlet changes too. However, the static pressure at the inlet corresponding to a choked flow without a shock train is initially unknown. The 1-D model proved to be a useful tool to compute this pressure, as it provides a cross-stream averaged description of a choked flow with a uniform static pressure profile. The 1-D model was thus run first, to provide the boundary conditions for the RANS simulations afterwards (see table 2). The primary stream is supersonic in most cases, except if 
$p_{t,s} \geqslant 2$ bar. In that case, only the total pressure and total temperature were imposed. The other boundary conditions remained identical to those in table 1.
Table 2. Boundary conditions for the parametric study. The primary static pressure 
$p_p$ was computed a priori with the 1-D model and served a posteriori as a boundary condition for the RANS simulations. The primary stream is subsonic at the inlet for 
$p_{t,s} \geqslant 2$ bar, so the static pressure is not imposed.
6. Results
The two approaches to close the model from § 2.2 are compared in § 6.1, followed by a validation of the model in various operating conditions and geometries in § 6.2. The effect of the friction forces on compound choking is illustrated in § 6.3. Compound choking is shown to be a plausibly more general flow blockage mechanism than individual streams reaching a unitary Mach number in § 6.4. This is demonstrated through a choked compound flow with a fully subsonic secondary stream, both with the 1-D model and with an axisymmetric RANS simulation.
6.1. Model closure
Before proceeding with the analysis of the investigated test cases, figure 6 showcases the impact of the friction coefficients at the wall and at the dividing streamline for the test case 3. Figures 6(a) and 6(b) show the streamwise evolution of the friction coefficients (i) obtained with the correlations (2.20) and (2.22), (ii) identified as constants by the inverse method as described in § 2.2 and (iii) extracted from the viscous RANS simulations as reference. The correlation (2.20) for the wall friction is adequate and follows the correct trend. The calibrated constant matches the reference on average. The friction between the streams is mostly underpredicted by correlation (2.22), as also observed by Grazzini et al. (Reference Grazzini, Mazzelli and Milazzo2016). This was found to be sensitive to the turbulent intensity at the inlets, which directly impacts the turbulent viscosity and hence the friction between the streams. Analysing this effect in further details is beyond the scope of the current work.
Figure 6. Comparison of the two approaches for determining the friction coefficients (see § 2.2). The model predicts the static pressure accurately with the correlations (2.20) and (2.22), and with the calibrated constants ( 
$f_w = 0.00377, f_{ps}=0.0355$). (a) Wall friction coefficient 
$f_w$; (b) interstream friction coefficient 
$f_{ps}$; (c) prediction of the static pressure with the correlations; (d) prediction of the static pressure with calibrated constants.
More importantly, the distribution of the static pressure is predicted accurately by the 1-D model in both cases, as shown in figures 6(c) and 6(d). The averaged static pressure is computed from the RANS simulations with (5.2a,b)–(5.6a–c) and is displayed with diamonds as the reference. For plotting purposes, the spacing in the markers is much larger than the mesh resolution.
It is worth noticing that the static pressure in the two streams oscillates near the inlet due to a weak shock train. This effect is not captured in the 1-D model because of the assumption of uniform static pressure. The agreement between the model and the RANS simulations shows that the correlations (2.20) and (2.22) are adequate for the purposes of this work. Moreover, the agreement between predictions with both constant and varying coefficients indicates a weak sensitivity to these parameters, justifying the exclusion of their derivatives in the approximation in § 3. All the results presented for the 1-D model in the remaining of this section were obtained with the correlations. This approach has the advantage of not requiring calibration since the model is fully closed when using (2.20) and (2.22).
6.2. Model validation
Figure 7 illustrates the distribution of the Mach numbers and the dividing streamlines for the different secondary inlet pressures (see table 2). The primary Mach number decreases significantly with increasing total pressure 
$p_{t,s}$, due to an increased static pressure at the inlet. Similarly, the secondary Mach number increases with increasing 
$p_{t,s}$, but is less sensitive. When the total pressures are equal (uniform flow), the Mach numbers are equal in each stream at the inlet. Downstream, the secondary Mach number decreases relative to the primary due to wall friction, which is not fully counteracted by the weak shear force from the quasiuniform flow.
Figure 7. Comparison of the model predictions and the postprocessed RANS simulations with varying total pressure 
$p_{t,s}$ at the secondary inlet and constant total pressure 
$p_{t,p} = 3$ bar at the primary inlet. The reference case is highlighted in bold. The geometry is identical to the reference case in § 5.1. The primary Mach number decreases with increasing 
$p_{t,s}$ due to a higher static pressure at the inlet and eventually becomes subsonic. The agreement degrades with a larger difference in inlet pressures. Here (a) Mach numbers; (b) dividing streamlines.
The model agreement with RANS simulations worsens as the inlet pressure difference grows. Notably, in these cases, the primary Mach number and dividing streamline oscillate near the inlet, due to the formation of a mild shock train. The resulting oblique shocks introduce pressure losses that are not captured by the 1-D model, which assumes a uniform static pressure across the streams. Additionally, the maximal relative error of the model amounts to 10 %, which is comparable to the scatter of the data on which the correlation (2.22) is based (see Papamoschou Reference Papamoschou1993). Despite this, the overall agreement is satisfactory, especially considering that the model was not specifically calibrated for these conditions.
The throat radius 
$R_{th}$ was varied between 60 % and 90 % of the outlet radius 
$R_o$, while maintaining the same total pressures as in case of table 1. A narrower throat induces a stronger expansion, resulting in lower Mach numbers upstream and higher Mach numbers downstream of the sonic section (see figure 8). The agreement between the model and the RANS simulations is excellent, with discrepancies below 10 %. The overpredicted Mach numbers at the outlet are attributed to underpredicted friction coefficients by (2.20) and (2.22) with respect to the RANS simulations, as also observed in figures 6(a) and 6(b). It is worth noticing that the radius of the wall at the inlet changes with the throat due to the cosine-shaped profile between 
${\rm \pi} /2$ and 
$2 {\rm \pi}$, as defined by (5.1). The primary inlet changes accordingly and extends halfway to the total inlet in each configuration.
Figure 8. Comparison of the model predictions and the postprocessed RANS simulations with varying radius of the throat. The reference case is highlighted in bold. The total pressures at the inlet are identical to the reference case 3 in table 1. The Mach numbers decrease at the inlet and increase at the outlet with a narrowing throat. The curvature becomes important for the smallest throat, inducing 2-D effects and reducing the accuracy of a 1-D formulation. Here (a) Mach numbers; (b) dividing streamlines.
Finally, the radius of the primary inlet is varied between 30 % and 70 % of the total radius at the inlet, while keeping the same profile of the wall and the same total pressures as the reference case. Figure 9 shows that both Mach numbers decrease with an increasing radius of the primary inlet. This is due to a higher average total pressure at the inlet (a larger portion with 
$p_{t,p} = 3$ bar instead of 
$p_{t,s} = 1.5$ bar), which in turn leads to a higher static pressure of the choked flow. Despite the lower Mach numbers in both streams, the equivalent Mach number is comparable. This surprising effect is possible through the nonlinear relation between the Mach numbers and the cross-sections in the terms of 
$\beta$ in (2.13). Clearly, the Mach numbers alone do not suffice to conclude on the state of the flow and should be complemented with information on the cross-sections. The 1-D model captures these trends accurately, although the Mach numbers are overpredicted near the outlet as observed earlier in figure 8.
Figure 9. Comparison of the model predictions and the postprocessed RANS simulations with varying radius of the primary inlet. The reference case is highlighted in bold. The total pressures at the inlet are identical to the reference case 3 in table 1. A larger primary cross-section at the inlet increases the average total pressure at the inlet and hence increases the static pressure required for choking. Consequently, both Mach numbers decrease with an increasing radius of the primary inlet. Here (a) Mach numbers; (b) dividing streamlines.
6.3. The effect of friction on compound choking
Figure 10 shows the Mach number fields of the three simulations in table 1, along with the dividing streamline in white and the sonic line in black. The shear layer in the inviscid case is only due to numerical diffusion, so the two streams remain clearly distinct across the dividing streamline. The sonic line ends on the wall downstream of the geometrical throat (indicated by the arrows). The primary stream is thus supersonic at the throat and the secondary stream is still subsonic. The secondary stream is accelerated through shear in the viscous simulations, causing the sonic line to move radially outwards. The sonic line no longer reaches the wall when a no-slip condition is imposed. The development of the boundary layer tends to slow down the flow, leading to a lower maximal Mach number.
Figure 10. Mach number contours of the three simulations (cf. table 1). The black line indicates the sonic line and tends towards the wall as the compound flow becomes supersonic. The dividing streamline is indicated in white. The arrows indicate the geometrical throat. The numerical domain depicted in figure 4 has been mirrored for visibility.
The averaged Mach numbers extracted from the RANS simulations are displayed with diamonds in figure 11 as the reference. The predictions of the proposed model are shown in full lines. It is worth noticing that the static pressure at the inlet for the 1-D model prediction is not imposed to be equal to the one from the RANS simulations, but rather results from the iterative procedure laid out in § 3. The only difference in the three cases is the friction coefficients, which are used to obtain the 1-D predictions. In the inviscid case, both naturally equal 0. The friction coefficients are computed with (2.20) and (2.22) for the viscous simulations.
Figure 11. Comparison of the Mach number computed with the 1-D compound flow theory and from the postprocessed RANS results. The sonic section (dotted line) coincides with the throat (red triangle) in the inviscid case 1. Friction between the streams pushes the sonic section upstream (case 2), while wall friction pushes it downstream in case 3 (see (2.25)). The used friction coefficients are given by (2.20) and (2.22). Here (a) case 1, sonic at 
$x/L = 0.333$; (b) case 2, sonic at 
$x/L = 0.323$; (c) case 3, sonic at 
$x/L = 0.359$.
The key finding in these results is the displacement of the sonic section, which is highlighted by the vertical dashed line and commented below. The sonic section is located in the geometrical throat (indicated by the red triangle) in the inviscid case, as predicted by the original compound flow theory (see e.g. Bernstein et al. Reference Bernstein, Heiser and Hevenor1967). The viscous simulation, without wall friction but with shear forces between the streams, shows an upstream shift towards the convergent section, as expected in the analysis of § 2.3. Conversely, wall friction is dominant in the final simulation, overcoming the effects of the shear force and ultimately pushing the sonic section in the divergent portion of the nozzle.
Note that the primary Mach number is slightly overestimated in all cases. This is attributed to the weak shock train at the inlet (see figure 6c,d and 10), which introduces a pressure loss in the primary stream.
The 1-D predictions allow for computing each term in the equation for the static pressure gradient (cf. (2.24)) separately. The evolution of the components in the numerator is shown in figure 12 for the slip and no-slip cases. The inviscid case is trivial since the gradient of the cross-section is the only remaining term. The terms are normalised with the maximal gradient of the cross-section, as this is a fixed quantity for the three simulations. The evolution of the two friction terms supports the theoretical considerations in § 2.3, the first (interstream friction) always being positive and the second (wall friction) being negative. As a result, the zero-crossing of the numerator moves up or downstream. The gradient of the cross-section balances the effect of the forces in either direction by moving the sonic section in the convergent or the divergent section. Note that the interstream friction tends to decrease from left to right as the difference in Mach number decreases, while friction tends to increase as the flow accelerates through the nozzle.
Figure 12. Components of the numerator in (2.24), normalised by the maximal area gradient. The positive interstream friction term pushes the zero-crossing of the sum upstream, while the negative wall friction term pushes it downstream. Here (a) slip, sonic at 
$x/L = 0.323$; (b) no slip, sonic at 
$x/L = 0.359$.
6.4. Flow blockage mechanism
As observed in figure 11, the secondary Mach number is below unity in the compound-sonic section if the primary stream is supersonic. It reaches unity farther downstream, where the compound flow continues to expand in the diverging section of the nozzle. It can thus be stated that the compound flow is sonic in the first section and that the secondary stream is sonic, separately, in a section farther downstream. Both sonic sections are plausible explanations for the blockage of the secondary mass flow rate, e.g. in ejectors, these choking mechanisms are referred to as ‘compound choking’ and ‘Fabri choking’ (see Lamberts et al. Reference Lamberts, Chatelain, Bourgeois and Bartosiewicz2018). There is currently no consensus on which of these mechanisms limits the secondary mass flow rate.
The excellent match between the postprocessed RANS simulations and the 1-D predictions of the compound flow theory suggests that the compound choking condition is the more general mechanism. Furthermore, the secondary Mach number reaches unity necessarily downstream of the compound-sonic section since a combination of a sonic and a supersonic stream is necessarily compound-supersonic (cf. (2.13)). This suggests that cases identified as ‘Fabri choking’ also feature compound choking farther upstream (see Hoge & Segars Reference Hoge and Segars1965). Moreover, Lamberts et al. (Reference Lamberts, Chatelain, Bourgeois and Bartosiewicz2018) and Kracik & Dvorak (Reference Kracik and Dvorak2023) showed that the static pressure required to reach compound choking is higher than the one required to reach 
$\textit {Ma}_s = 1$. This implies that the compound condition is less restrictive (the flow needs to expand less to reach the compound condition), but still sufficient to choke the flow. As a definitive proof of the generality of the compound choking condition to explain the mass flow limitation (compared with Fabri choking), a compound choked flow is demonstrated with a fully subsonic secondary stream.
To lower the secondary Mach number below unity, the divergent section of the nozzle presented in § 5 is replaced by a constant cross-section having the same diameter as the throat. The first simulation was carried out on a nozzle of the same length, while the second is cut shorter at 
$x/L=0.4$ (so downstream of the throat at 
$x/L=0.333$). The settings for the mesh, the boundary conditions and the numerical solver remain identical to the previous simulations. The absence of a divergent section would severely impact the viscous case with the no-slip condition at the wall since the compound sonic section is located in the divergent section. The resulting flow field forms a Fanno flow, which becomes nearly sonic at the outlet of the constant-area pipe. Therefore, the viscous case with a slip condition at the wall is retained, which becomes compound-sonic upstream of the throat regardless of the geometry downstream.
Figure 13 shows the evolution of the postprocessed Mach number of the convergent–constant and cut (still convergent–constant) configurations. The compound choking condition is reached at 
$x/L = 0.323$ in both cases, and the secondary mass flow rates are identical to 
$0.0622\ {\rm kg}\ {\rm s}^{-1}$ (since the geometry and the boundary conditions at the inlet are also identical). The flow expands less in the constant area section compared with the convergent–divergent case, but the secondary Mach number still reaches unity at 
$x/L = 0.56$. However, the secondary cross-section does not reach a minimum in this section, suggesting that the secondary stream does not choke separately. This is confirmed by the simulation with the shorter nozzle, which coincides with the results of the complete nozzle. Despite the low pressure imposed as a boundary condition at the outlet, the secondary flow does not accelerate beyond 
$\textit {Ma}_s = 0.97$ (both in the RANS simulation and in the 1-D model). Hence, no larger secondary mass flow rate can be obtained, and the compound flow is choked with a fully subsonic secondary stream. These results suggest that compound choking is the more general choking condition in ejectors and, more generally, for parallel compressible streams with different stagnation pressures.
Figure 13. Evolution of the Mach number for the convergent–constant and cut nozzle with friction between the streams and a slip condition at the wall (see figure 11b). The results coincide in the convergent section as the flow is choked in all cases. The secondary Mach number does not reach unity in the shorter nozzle. Here (a) case 2, convergent–constant configuration; (b) case 2, cut configuration.
7. Conclusion
The isentropic compound flow theory has been extended with friction forces between the streams and at the wall. The resulting 1-D model captures the upstream and downstream displacement of the sonic section observed by Lamberts et al. (Reference Lamberts, Chatelain, Bourgeois and Bartosiewicz2018) and Kracik & Dvorak (Reference Kracik and Dvorak2023). Moreover, the 1-D predictions have been validated numerically to 2-D axisymmetric RANS simulations with appropriate boundary conditions to include or eliminate the effect of friction and a parametric study has demonstrated that the model generalises well across different operating conditions and different geometries. An extension of this work to the modelling of ejectors is foreseen in the future.
To the authors’ knowledge, our work is the first to propose an explanation for choking of a compound flow upstream of the geometrical throat. The Taylor expansion of the equation for the static pressure gradient allows for circumventing its singularity, avoiding numerical instabilities and increasing the accuracy of the numerical predictions. Analytical expressions have been presented for all required derivatives for a second-order approximation, with a quasi-third-order approximation obtained with finite differences. These expressions and the governing equations have first been derived in a general way in terms of a force 
$F_i$, facilitating future developments such as a compound flow with more than two streams or a planar geometry instead of the axisymmetrical one considered in this work. Another open path includes heat exchange and different total temperatures at the inlet, which may affect the shifting mechanism due to friction between the streams, as discussed in § 2.3. In particular, a high primary inlet temperature may invert the effect of interstream friction. Finally, we believe that the developed model not only provides more insight into the choking of compound flows but also opens perspectives in the reduced-order modelling of ejectors.
Funding
J.V.d.B. is supported by F.R.S.-FNRS FRIA grant no. 47455.
Declaration of interests
The authors report no conflict of interest.
Appendix A. First-order approximation of the static pressure gradient
The static pressure gradient in the sonic section has been calculated in § 2.3 using l'Hôpital's rule (see (2.28)). This operation requires the first derivatives of the numerator 
$N$ and the compound choking indicator 
$\beta$. This section provides analytical expressions for these derivatives and those of the forces (
${\rm d}F_i/{{\rm d}\kern0.7pt x}$), ultimately leading to a quadratic equation for the static pressure gradient in the sonic section.
A.1. Friction forces
As discussed in § 2.3, the derivatives of the numerator and 
$\beta$ give rise to derivatives of Mach numbers 
$\textit {Ma}_i$, cross-section 
$A_i$, forces 
$F_i$ and the static pressure 
$p$. These can be substituted by the governing equations (2.7), (2.8), (2.10), resulting in expressions that are linear with respect to the static pressure gradient. This is also true for the friction forces. The first derivative of the wall friction force is given by the following equation:
\begin{equation} \dfrac{{\rm d}F_w}{{\rm d}\kern0.7pt x} = c_{Fw0} \left(\dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right) + c_{Fw1}, \end{equation}where
$$\begin{gather} c_{Fw0} = f_w \sqrt{\rm \pi} \sqrt{A} p \left( \textit{Ma}_s^2-2\right), \end{gather}$$
$$\begin{gather}c_{Fw1} = \dfrac{\sqrt{\rm \pi}\left(f_w \gamma A_s \textit{Ma}_s^2 p \dfrac{{\rm d}A}{{\rm d}\kern0.7pt x} +\left( 2 f_w \gamma-2 f_w\right) A F_s \textit{Ma}_s^2+4 f_w A F_s\right)}{2 \sqrt{A} A_s}. \end{gather}$$
The interstream friction force explicitly depends on densities, velocities and the primary cross-section (see (2.17)). The vector 
$\boldsymbol {q} = [\rho _p, \rho _s, u_p, u_s, A_p, A_s, \textit {Ma}_p^2, \textit {Ma}_s^2]$ is introduced to ease the notation. The first derivative of the friction between the streams is given by the following equation:
\begin{equation} \dfrac{{\rm d} F_{ps}}{{\rm d}\kern0.7pt x} = \dfrac{\partial F_{ps}}{\partial \boldsymbol{q}} \boldsymbol{\cdot} \dfrac{{\rm d}\boldsymbol{q}}{{\rm d}\kern0.7pt x}, \end{equation}
where the dot operator indicates an inner product. The equations below are derived for the case where 
$u_p > u_s$ to avoid the absolute value operator in (2.17). The derivatives in the case 
$u_p < u_s$ have the same magnitude, but only differ in the sign. It thus suffices to add a minus sign in the equation above. The partial derivatives of 
$F_{ps}$ with respect to 
$\boldsymbol {q}$ are given by
\begin{equation} \dfrac{\partial F_{ps}}{\partial \boldsymbol{q}} = \begin{pmatrix} \dfrac{f_{ps} \sqrt{\rm \pi} \sqrt{A_p} {{\left( u_p-u_s\right) }^{2}}}{2}\\[8pt] \dfrac{f_{ps} \sqrt{\rm \pi} \sqrt{A_p} {{\left( u_p-u_s\right) }^{2}}}{2}\\[8pt] f_{ps} \sqrt{\rm \pi} \sqrt{A_p} \left( \rho_s+\rho_p\right) \left( u_p-u_s\right) \\[8pt] -f_{ps} \sqrt{\rm \pi} \sqrt{A_p} \left( \rho_s+\rho_p\right) \left( u_p-u_s\right) \\[8pt] \dfrac{f_{ps} \sqrt{\rm \pi} \left( \rho_s+\rho_p\right) {{\left( u_p-u_s\right) }^{2}}}{4 \sqrt{A_p}}\\[8pt] 0\\ 0\\ 0 \end{pmatrix}. \end{equation}
Note that the last three quantities of 
$\boldsymbol {q}$ do not appear directly in the definition of 
$F_{ps}$. These are included in 
$\boldsymbol {q}$ because they appear in the second derivatives (see Appendix B). The derivative 
${\rm d}\boldsymbol {q}/{{\rm d}\kern0.7pt x}$ follows from the governing equations (see § 2.1 and Shapiro (Reference Shapiro1953)),
\begin{equation} \dfrac{{\rm d}\boldsymbol{q}}{{\rm d}\kern0.7pt x} = \boldsymbol{c}_{q0} \left(\dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right) + \boldsymbol{c}_{q1p} F_p + \boldsymbol{c}_{q1s} F_s, \end{equation}
where 
$\boldsymbol {c}_{q0}$, 
$\boldsymbol {c}_{q1p}$ and 
$\boldsymbol {c}_{q1s}$ are vectors of the same length as 
$\boldsymbol {q}$, and themselves functions of 
$\boldsymbol {q}$,
\begin{gather} \boldsymbol{c}_{q0} = \begin{pmatrix} \dfrac{\rho_p}{\gamma}\\[13pt] \dfrac{\rho_s}{\gamma}\\[13pt] -\dfrac{u_p}{\gamma \textit{Ma}_p^2} \\[13pt] -\dfrac{u_s}{\gamma \textit{Ma}_s^2} \\[13pt] \dfrac{A_p}{\gamma} \left( \dfrac{1}{\textit{Ma}_p^2}-1\right)\\[14pt] \dfrac{A_s}{\gamma} \left( \dfrac{1}{\textit{Ma}_s^2}-1\right)\\[14pt] -\dfrac{2}{\gamma} \left( \dfrac{\left( \gamma-1\right)}{2} \textit{Ma}_p^2 +1\right)\\[14pt] -\dfrac{2}{\gamma} \left( \dfrac{\left( \gamma-1\right)}{2} \textit{Ma}_s^2 +1\right) \end{pmatrix}, \end{gather}
\begin{gather} \boldsymbol{c}_{q1p} = \!\begin{pmatrix} \dfrac{\left( \gamma-1\right) \rho_p}{\gamma A_p p}\\[7pt] 0\\[3pt] \dfrac{u_p}{\gamma A_p \textit{Ma}_p^2 p}\\[7pt] 0\\[3pt] \dfrac{-\left( \dfrac{1}{\gamma \textit{Ma}_p^2}\right) -\dfrac{\gamma-1}{\gamma}}{p}\\[7pt] 0\\[3pt] \dfrac{2 \left( \dfrac{\left( \gamma-1\right) \textit{Ma}_p^2}{2}+1\right) }{\gamma A_p p}\\[6pt] 0 \end{pmatrix} \, \mbox{and} \, \boldsymbol{c}_{q1s} = \!\begin{pmatrix} 0\\ \dfrac{\left( \gamma-1\right) \rho_s}{\gamma A_s p}\\[14pt] 0\\ \dfrac{u_s}{\gamma A_s \textit{Ma}_s^2 p}\\[14pt] 0\\ \dfrac{-\left( \dfrac{1}{\gamma \textit{Ma}_s^2}\right) -\dfrac{\gamma-1}{\gamma}}{p}\\[14pt] 0\\ \dfrac{2 \left( \dfrac{\left( \gamma-1\right) \textit{Ma}_s^2}{2}+1\right) }{\gamma A_s p} \end{pmatrix}. \end{gather}Combining (A4)–(A8a,b) leads to the following equation:
\begin{equation} \dfrac{{\rm d}F_{ps}}{{\rm d}\kern0.7pt x} = c_{Fps0} \left(\dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right) + c_{Fps1}, \end{equation}where
\begin{equation} c_{Fps0} = \dfrac{\partial F_{ps}}{\partial \boldsymbol{q}} \boldsymbol{\cdot} \boldsymbol{c}_{q0}, \end{equation}and
\begin{equation} c_{Fps1} = \dfrac{\partial F_{ps}}{\partial \boldsymbol{q}} \boldsymbol{\cdot} \left(\boldsymbol{c}_{q1p} F_p + \boldsymbol{c}_{q1s} F_s\right). \end{equation}
The derivatives of the forces 
$F_i$ follow from (A1), (A9) and their definition (2.19a,b),
\begin{equation} \dfrac{{\rm d}F_i}{{\rm d}\kern0.7pt x} = c_{Fi0} \left(\dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right) + c_{Fi1}, \end{equation}
where 
$i\in [p,s]$ and
\begin{equation} \left.\begin{gathered} c_{Fp0} ={-}c_{Fps0}, \quad {\rm and} \quad c_{Fs0} = c_{Fps0} - c_{Fw0},\\ c_{Fp1} ={-}c_{Fps1}, \quad {\rm and} \quad c_{Fs1} = c_{Fps1} - c_{Fw1}. \end{gathered} \right\}\end{equation}A.2. The numerator of the equation for the static pressure
The numerator is differentiated in its general form (2.12) for any definition of the forces 
$F_i$ to facilitate potential modifications to the forces in future work,
\begin{equation} N = \dfrac{{\rm d}A}{{\rm d}\kern0.7pt x} + \sum_{i\in[p,s]}\left[\dfrac{1 + \left(\gamma - 1\right)\textit{Ma}_i^2}{\gamma \textit{Ma}_i^2}\right] \dfrac{F_i}{p}. \end{equation}Introducing the concise notation for the second term,
\begin{equation} Z_i = \left[\dfrac{1 + \left(\gamma - 1\right)\textit{Ma}_i^2}{\gamma \textit{Ma}_i^2}\right] \dfrac{F_i}{p},\end{equation}
the gradient of the numerator gives rise to first derivatives of 
$Z_i$:
\begin{equation} \dfrac{{\rm d}N}{{\rm d}\kern0.7pt x} = \dfrac{{\rm d}^2A}{{{\rm d}\kern0.7pt x}^2} + \sum_{i\in[p,s]} \dfrac{{\rm d}Z_i}{{\rm d}\kern0.7pt x}. \end{equation}These are given by
\begin{equation} \dfrac{{\rm d}Z_i}{{\rm d}\kern0.7pt x} = c_{Zi0} \left( \dfrac{{\rm d}F_i}{{\rm d}\kern0.7pt x}\right) + c_{Zi1} \left(\dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right) + c_{Zi2}, \end{equation}where
$$\begin{gather} c_{Zi0} = \dfrac{1}{p}\left(\dfrac{1 + \left(\gamma - 1\right)\textit{Ma}_i^2}{\gamma \textit{Ma}_i^2}\right), \end{gather}$$
$$\begin{gather}c_{Zi1} = \dfrac{\left(\left(\gamma-\gamma^2\right)\textit{Ma}_i^4 - \textit{Ma}_i^2 +2\right)F_i}{\gamma^2 p \textit{Ma}_i^4}, \end{gather}$$
$$\begin{gather}c_{Zi2} ={-}\dfrac{2 F_i^2 \left( \dfrac{\left( \gamma-1\right) \textit{Ma}_i^2}{2}+1\right) }{\gamma^2 p^2 A_i \textit{Ma}_i^4}. \end{gather}$$The combination of (A16), (A17) and (A12) yields the final expression for the first derivative of the numerator,
\begin{equation} \dfrac{{\rm d}N}{{\rm d}\kern0.7pt x} = c_{N0} \left(\dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right) + c_{N1}, \end{equation}where
$$\begin{gather} c_{N0} = \sum_{i\in[p,s]} \left(c_{Zi0} c_{Fi0} + c_{Zi1}\right), \end{gather}$$
$$\begin{gather}c_{N1} = \dfrac{{\rm d}^2A}{{{\rm d}\kern0.7pt x}^2} + \sum_{i\in[p,s]} \left(c_{Zi0} c_{Fi1} + c_{Zi2}\right), \end{gather}$$A.3. The compound choking indicator
The compound choking indicator can be written as a sum of components 
$\beta _i$ (cf. (2.13)):
\begin{equation} \beta = \sum_{i\in[p,s]} \beta_i = \sum_{i\in[p,s]} A_i\dfrac{1-\textit{Ma}_i^2}{\gamma \textit{Ma}_i^2}. \end{equation}
The first derivative of 
$\beta$ is given by the following equation:
\begin{equation} \dfrac{{\rm d}\beta}{{\rm d}\kern0.7pt x} = c_{\beta 0} \left(\dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right) + c_{\beta 1}, \end{equation}where
$$\begin{gather} c_{\beta 0} = \sum_{i\in[p,s]} \dfrac{A_i\left({\textit{Ma}_i^4}+\left( \gamma-3\right) \textit{Ma}_i^2+3 \right) }{\gamma^2 {\textit{Ma}_i^4}}, \end{gather}$$
$$\begin{gather}c_{\beta 1} = \sum_{i\in[p,s]} \dfrac{\left( \gamma-1\right) {\textit{Ma}_i^4}+\left( 3-2 \gamma\right) \textit{Ma}_i^2-3 }{\gamma^2 {\textit{Ma}_i^4}}\left(\dfrac{F_i}{p}\right). \end{gather}$$A.4. Static pressure gradient
Using l'Hôpital's rule on (2.12) leads to (2.32), which is equivalent to the following quadratic equation:
\begin{equation} c_{\beta 0}^*\left(\dfrac{1}{p^*} \dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}^*\right)^2 + \left(c_{\beta1}^* - c_{N0}^*\right)\left(\dfrac{1}{p^*} \dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}^*\right) -c_{N1}^* = 0. \end{equation}The stars indicate that the variables are evaluated at the sonic section. The two roots correspond to a subsonic or supersonic solution downstream of the sonic section. The first-order approximation in the vicinity of the sonic section consists of the same constant gradient in the throat given by the equation above (see figure 2).
Appendix B. Second-order approximation of the static pressure gradient
The second-order approximation requires the first and the second derivatives of the numerator 
$N$ and 
$\beta$ (see (3.2)). These give rise to second derivatives of the quantities encountered in Appendix A. In this section, the first derivatives are not expanded since they can be calculated directly using the expressions in Appendix A.
B.1. Mach number
Here
\begin{equation} \dfrac{{\rm d}^2\textit{Ma}_i^2}{{{\rm d}\kern0.7pt x}^2} = c_{\textit{Ma}_i^20} \left( \dfrac{{\rm d}}{{\rm d}\kern0.7pt x} \left(\dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right) \right) + c_{\textit{Ma}_i^21}, \end{equation}where
\begin{gather} c_{\textit{Ma}_i^20} ={-} \dfrac{2}{\gamma} \left( \dfrac{\left( \gamma-1\right) \textit{Ma}_i^2}{2}+1\right), \end{gather}
\begin{gather}\begin{aligned} c_{\textit{Ma}_i^21} &={-}\dfrac{\left( \gamma-1\right)}{\gamma} \left(\dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right) \left( \dfrac{{\rm d} \textit{Ma}_i^2}{{\rm d}\kern0.7pt x}\right) +\dfrac{\left( \gamma-1\right)}{\gamma} \dfrac{F_i}{p A_i} \left( \dfrac{{\rm d} \textit{Ma}_i^2}{{\rm d}\kern0.7pt x}\right)\\ &\quad +\dfrac{2}{\gamma p A_i} \left( \dfrac{\left( \gamma-1\right) \textit{Ma}_i^2}{2}+1\right) \left( \dfrac{{\rm d}F_i}{{\rm d}\kern0.7pt x}\right) - \dfrac{2 F_i}{\gamma p A_i} \left( \dfrac{\left( \gamma-1\right) \textit{Ma}_i^2}{2}+1\right) \left( \dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right), \end{aligned}\end{gather}
where the derivatives 
$\textrm {d}\textit {Ma}_i^2/{\textrm {d}{\kern0.9pt}x}$ are the last two components of the vector 
$\textrm {d}\boldsymbol {q}/{\textrm {d}{\kern0.9pt}x}$ in (A6).
B.2. Friction forces
The second derivative of the wall friction force is given by the following equation:
\begin{equation} \dfrac{{\rm d}^2 F_w}{{{\rm d}\kern0.7pt x}^2} = c_{Fw^20} \left( \dfrac{{\rm d}}{{\rm d}\kern0.7pt x} \left(\dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right)\right) + c_{Fw^21}, \end{equation}where
\begin{gather} c_{Fw^20} = f_w \sqrt{\rm \pi} \sqrt{A} p \left(\textit{Ma}_s^2 - 2\right), \end{gather}
\begin{gather}\begin{aligned}[b] c_{Fw^21} &= f_w \gamma \sqrt{\rm \pi} \sqrt{A} p c_{\textit{Ma}_s^21} +\dfrac{f_w \gamma \sqrt{\rm \pi} p }{\sqrt{A}} \left( \dfrac{{\rm d}A}{{\rm d}\kern0.7pt x}\right) \left( \dfrac{{\rm d}\textit{Ma}_s^2}{{\rm d}\kern0.7pt x}\right) \\ &\quad +2 f_w \gamma \sqrt{\rm \pi} \sqrt{A} p \left(\dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right) \left( \dfrac{{\rm d}\textit{Ma}_s^2}{{\rm d}\kern0.7pt x}\right) +\dfrac{f_w \gamma \sqrt{\rm \pi} \textit{Ma}_s^2 p}{2 \sqrt{A}} \left( \dfrac{{\rm d}^2 A}{{{\rm d}\kern0.7pt x}^2}\right)\\ &\quad -\dfrac{f_w \gamma \sqrt{\rm \pi} \textit{Ma}_s^2 p}{4 {{A}^{{3}/{2}}}} \left( \dfrac{{\rm d}A}{{\rm d}\kern0.7pt x}\right)^2 +\dfrac{f_w \gamma \sqrt{\rm \pi} \textit{Ma}_s^2 p}{\sqrt{A}} \left(\dfrac{1}{p}\dfrac{dp}{{\rm d}\kern0.7pt x}\right) \left( \dfrac{{\rm d}A}{{\rm d}\kern0.7pt x}\right)\\ &\quad +f_w \gamma \sqrt{\rm \pi} \sqrt{A} \textit{Ma}_s^2 {{\left(\dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right)}^{2}} p. \end{aligned}\end{gather} The second derivative of the interstream friction force follows from (A9) (with the same comment regarding the sign in case 
$u_p < u_s$),
\begin{equation} \dfrac{{\rm d}^2 F_{ps}}{{{\rm d}\kern0.7pt x}^2} = c_{Fps0} \left(\dfrac{{\rm d}}{{\rm d}\kern0.7pt x}\left(\dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right)\right) + \dfrac{{\rm d}c_{Fps0}}{{\rm d}\kern0.7pt x} \left(\dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right) + \dfrac{{\rm d}c_{Fps1}}{{\rm d}\kern0.7pt x}. \end{equation}
The derivatives of the coefficients 
$c_{Fps0}$ and 
$c_{Fps1}$ follow from (A10) and (A11),
\begin{align} \dfrac{{\rm d}c_{Fps0}}{{\rm d}\kern0.7pt x} &= \left(\dfrac{\partial^2 F_{ps}}{\partial \boldsymbol{q}^2} \dfrac{{\rm d}\boldsymbol{q}}{{\rm d}\kern0.7pt x}\right) \boldsymbol{\cdot} \boldsymbol{c}_{q0} + \dfrac{\partial F_{ps}}{\partial \boldsymbol{q}} \boldsymbol{\cdot} \dfrac{{\rm d}\boldsymbol{c}_{q0}}{{\rm d}\kern0.7pt x},\end{align}
\begin{align} \dfrac{{\rm d}c_{Fps1}}{{\rm d}\kern0.7pt x} &= \left(\dfrac{\partial^2 F_{ps}}{\partial \boldsymbol{q}^2} \dfrac{{\rm d} \boldsymbol{q}}{{\rm d}\kern0.7pt x}\right) \boldsymbol{\cdot} \left(\boldsymbol{c}_{q1p} F_p + \boldsymbol{c}_{q1s} F_s\right) \nonumber\\ &\quad + \dfrac{\partial F_{ps}}{\partial \boldsymbol{q}} \boldsymbol{\cdot} \left(\dfrac{{\rm d}\boldsymbol{c}_{q1p}}{{\rm d}\kern0.7pt x} F_p + \boldsymbol{c}_{q1p}\dfrac{{\rm d}F_p}{{\rm d}\kern0.7pt x} + \dfrac{{\rm d}\boldsymbol{c}_{q1s}}{{\rm d}\kern0.7pt x} F_s + \boldsymbol{c}_{q1s} \dfrac{{\rm d}F_s}{{\rm d}\kern0.7pt x} \right). \end{align}
The first derivatives 
$\partial F_{ps}/\partial \boldsymbol {q}$, 
$\textrm {d}\boldsymbol {q}/{\textrm {d}{\kern0.9pt}x}$, 
$\textrm {d}F_p/{\textrm {d}{\kern0.9pt}x}$ and 
$\textrm {d}F_s/{\textrm {d}{\kern0.9pt}x}$ are given by (A5), (A6) and (A12). The second derivatives of the friction force 
$F_{ps}$ with respect to 
$\boldsymbol {q}$ are given by the following matrix:
\begin{equation} \dfrac{\partial^2 F_{ps}}{\partial \boldsymbol{q}^2} = f_{ps} \sqrt{\rm \pi} \sqrt{A_p} \left(u_p-u_s\right) \boldsymbol{M}, \end{equation}where
\begin{align} \boldsymbol{M} = \!\begin{pmatrix} 0 & 0 & 1 & -1 & \dfrac{u_p-u_s}{4 A_p} & 0 & 0 & 0\\[14pt] 0 & 0 & 1 & -1 & \dfrac{ u_p-u_s}{4 A_p} & 0 & 0 & 0 \\[14pt] 1 & 1 & \dfrac{\rho_p+\rho_s}{u_p-u_s} & -\dfrac{\rho_p+\rho_s}{u_p-u_s} & \dfrac{\rho_p+\rho_s}{2 A_p} & 0 & 0 & 0 \\[14pt] -1 & -1 & -\dfrac{\rho_p+\rho_s}{u_p-u_s} & \dfrac{\rho_p+\rho_s}{u_p-u_s} & - \dfrac{\rho_p+\rho_s}{2 A_p} & 0 & 0 & 0 \\[14pt] \dfrac{u_p-u_s}{4 A_p} & \dfrac{u_p-u_s}{4 A_p} & \dfrac{\rho_p+\rho_s}{2 A_p} & -\dfrac{\rho_p+\rho_s}{2 A_p} & - \dfrac{ \left( \rho_p+\rho_s\right) \left( u_p-u_s\right)}{8 A_p^2} & 0 & 0 & 0 \\[14pt] 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \end{pmatrix}\!. \end{align}
The derivative of 
$\boldsymbol {c}_{q0}$ follows from (A7),
\begin{equation} \dfrac{{\rm d}\boldsymbol{c}_{q0}}{{\rm d}\kern0.7pt x} = \dfrac{\partial\boldsymbol{c}_{q0}}{\partial \boldsymbol{q}} \dfrac{{\rm d}\boldsymbol{q}}{{\rm d}\kern0.7pt x}, \end{equation}where
\begin{align} \dfrac{\partial \boldsymbol{C}_{q0}}{\partial \boldsymbol{q}} = \!\begin{pmatrix} \dfrac{1}{\gamma} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\[14pt] 0 & \dfrac{1}{\gamma} & 0 & 0 & 0 & 0 & 0 & 0 \\[14pt] 0 & 0 & - \dfrac{1}{\gamma \textit{Ma}_p^2} & 0 & 0 & 0 & \dfrac{u_p}{\gamma \textit{Ma}_p^4} & 0 \\[14pt] 0 & 0 & 0 & - \dfrac{1}{\gamma \textit{Ma}_s^2} & 0 & 0 & 0 & \dfrac{u_s}{\gamma \textit{Ma}_s^4} \\[14pt] 0 & 0 & 0 & 0 & \dfrac{\dfrac{1}{\textit{Ma}_p^2}-1}{\gamma} & 0 & - \dfrac{A_p}{\gamma \textit{Ma}_p^4} & 0 \\[14pt] 0 & 0 & 0 & 0 & 0 & \dfrac{\dfrac{1}{\textit{Ma}_s^2}-1}{\gamma} & 0 & - \dfrac{A_s}{\gamma \textit{Ma}_s^4} \\[14pt] 0 & 0 & 0 & 0 & 0 & 0 & - \dfrac{\gamma-1}{\gamma} & 0 \\[14pt] 0 & 0 & 0 & 0 & 0 & 0 & 0 & - \dfrac{\gamma-1}{\gamma} \end{pmatrix}\!. \end{align} The coefficients 
$\boldsymbol {c}_{q1p}$ and 
$\boldsymbol {c}_{q1s}$ explicitly depend on 
$\boldsymbol {q}$ and the static pressure 
$p$ (see (A8a,b)), so their derivatives are given by the following equations:
$$\begin{gather} \dfrac{{\rm d}\boldsymbol{c}_{q1p}}{{\rm d}\kern0.7pt x} = \dfrac{\partial\boldsymbol{c}_{q1p}}{\partial \boldsymbol{q}} \dfrac{{\rm d}\boldsymbol{q}}{{\rm d}\kern0.7pt x} + \dfrac{\partial \boldsymbol{c}_{q1p}}{\partial p} \dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}, \end{gather}$$
$$\begin{gather}\dfrac{{\rm d}\boldsymbol{c}_{q1s}}{{\rm d}\kern0.7pt x} = \dfrac{\partial\boldsymbol{c}_{q1s}}{\partial \boldsymbol{q}} \dfrac{{\rm d}\boldsymbol{q}}{{\rm d}\kern0.7pt x} + \dfrac{\partial \boldsymbol{c}_{q1s}}{\partial p} \dfrac{{\rm d}p}{{\rm d}\kern0.7pt x} , \end{gather}$$where
\begin{align} \dfrac{\partial \boldsymbol{C}_{q1p}}{\partial \boldsymbol{q}} = \begin{pmatrix} \dfrac{\gamma-1}{\gamma A_p p} & 0 & 0 & 0 & -\dfrac{\left( \gamma-1\right) \rho_p}{\gamma A_p^{2} p} & 0 & 0 & 0 \\[14pt] 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \dfrac{1}{\gamma A_p \textit{Ma}_p^2 p} & 0 & -\dfrac{u_p}{\gamma A_p^{2} \textit{Ma}_p^2 p} & 0 & -\dfrac{u_p}{\gamma A_p \textit{Ma}_p^4 p} & 0 \\[14pt] 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & \dfrac{1}{\gamma \textit{Ma}_p^4 p} & 0 \\[14pt] 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & -\dfrac{2 \left( \dfrac{\left( \gamma-1\right) \textit{Ma}_p^2}{2}+1\right) }{\gamma A_p^{2} p} & 0 & \dfrac{\gamma-1}{\gamma A_p p} & 0 \\[14pt] 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix}, \end{align}
\begin{align} \dfrac{\partial \boldsymbol{C}_{q1s}}{\partial \boldsymbol{q}} = \begin{pmatrix} 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & \dfrac{\gamma-1}{\gamma A_s p} & 0 & 0 & 0 & -\dfrac{\left( \gamma-1\right) \rho_s}{\gamma A_s^2 p} & 0 & 0 \\[14pt] 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \dfrac{1}{\gamma A_s \textit{Ma}_s^2 p} & 0 & -\dfrac{u_s}{\gamma A_s^2 \textit{Ma}_s^2 p} & 0 & -\dfrac{u_s}{\gamma A_s \textit{Ma}_s^4 p} \\[14pt] 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & \dfrac{1}{\gamma \textit{Ma}_s^4 p} \\[14pt] 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & -\dfrac{2 \left( \dfrac{\left( \gamma-1\right) \textit{Ma}_s^2}{2}+1\right) }{\gamma A_s^2 p} & 0 & \dfrac{\gamma-1}{\gamma A_s p} \end{pmatrix}, \end{align}
\begin{equation} \dfrac{\partial \boldsymbol{C}_{q1p}}{\partial p} = \begin{pmatrix} -\dfrac{\left( \gamma-1\right) \rho_p}{\gamma A_p {p^2}}\\[14pt] 0\\[6pt] -\dfrac{u_p}{\gamma A_p \textit{Ma}_p^2 {p^2}}\\[14pt] 0\\[6pt] \dfrac{\left( \dfrac{1}{\gamma \textit{Ma}_p^2}\right) +\dfrac{\gamma-1}{\gamma}}{{p^2}}\\[14pt] 0\\[6pt] -\dfrac{2 \left( \dfrac{\left( \gamma-1\right) \textit{Ma}_p^2}{2}+1\right) }{\gamma A_p {p^2}}\\[14pt] 0 \end{pmatrix} \quad \mbox{and} \quad \dfrac{\partial \boldsymbol{C}_{q1s}}{\partial p} = \begin{pmatrix} 0\\[6pt] -\dfrac{\left( \gamma-1\right) \rho_s}{\gamma A_s {p^2}}\\[14pt] 0\\[6pt] -\dfrac{u_s}{\gamma A_s \textit{Ma}_s^2 {p^2}}\\[14pt] 0\\[6pt] \dfrac{\left( \dfrac{1}{\gamma \textit{Ma}_s^2}\right) + \dfrac{\gamma-1}{\gamma}}{{p^2}}\\[14pt] 0\\[6pt] -\dfrac{2 \left( \dfrac{\left( \gamma-1\right) \textit{Ma}_s^2}{2}+1\right) }{\gamma A_s {p^2}} \end{pmatrix}. \end{equation}The second derivative of the interstream friction is thus given by
\begin{equation} \dfrac{{\rm d}^2 F_{ps}}{{{\rm d}\kern0.7pt x}^2} = c_{Fps^20} \left(\dfrac{{\rm d}}{{\rm d}\kern0.7pt x}\left(\dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right)\right) + c_{Fps^21}, \end{equation}where
\begin{equation} c_{Fps^20} = c_{Fps0} \quad \mbox{and} \quad c_{Fps^21} = \dfrac{{\rm d}c_{Fps0}}{{\rm d}\kern0.7pt x} \left(\dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right) + \dfrac{{\rm d}c_{Fps1}}{{\rm d}\kern0.7pt x} \end{equation}
are found with (A10), (B8), (B9) and the local pressure gradient. The second derivatives of the forces 
$F_i$ follow from (B4), (B19) and their definition (2.19a,b),
\begin{equation} \dfrac{{\rm d}^2F_i}{{{\rm d}\kern0.7pt x}^2} = c_{Fi^20} \left(\dfrac{{\rm d}}{{\rm d}\kern0.7pt x}\left(\dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right)\right) + c_{Fi^21}, \end{equation}
where 
$i\in [p,s]$ and
\begin{equation} \left.\begin{gathered} c_{Fp^20} ={-}c_{Fps^20} \quad \mbox{and} \quad c_{Fs^20} = c_{Fps^20} - c_{Fw^20},\\ c_{Fp^21} ={-}c_{Fps^21} \quad \mbox{and} \quad c_{Fs^21} = c_{Fps^21} - c_{Fw^21}. \end{gathered}\right\}\end{equation}B.3. The numerator of the equation for the static pressure
The second derivative of the numerator gives rise to the second derivatives of 
$Z_i$:
\begin{equation} \dfrac{{\rm d}^2N}{{{\rm d}\kern0.7pt x}^2} = \dfrac{{\rm d}^3A}{{{\rm d}\kern0.7pt x}^3} + \sum_{i\in[p,s]} \dfrac{{\rm d}^2Z_i}{{{\rm d}\kern0.7pt x}^2}. \end{equation}
The second derivative of the force term 
$Z_i$ (cf. (A15)) is given by the following equation:
\begin{align} \dfrac{{\rm d}^2Z_i}{{{\rm d}\kern0.7pt x}^2} &= c_{Zi^20} \left( \dfrac{{\rm d}}{{\rm d}\kern0.7pt x} \left(\dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right)\right) +c_{Zi^21} \left(\dfrac{{{\rm d}^2F_i}}{{{\rm d}\kern0.7pt x}^2}\right) +c_{Zi^22} \left( \dfrac{{\rm d}F_i}{{\rm d}\kern0.7pt x}\right)\nonumber\\ &\quad + c_{Zi^23} \left(\dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right)^{2} + c_{Zi^24} \left(\dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right) + C_{Z^25,i}, \end{align}where
$$\begin{gather} c_{Zi^20} ={-}\dfrac{\left( \gamma^2-\gamma\right) F_i \textit{Ma}_i^4+ F_i \textit{Ma}_i^2-2 F_i}{\gamma^2 p \textit{Ma}_i^4}, \end{gather}$$
$$\begin{gather}c_{Zi^21} = \dfrac{\left( \gamma-1\right) \textit{Ma}_i^2+1}{\gamma p \textit{Ma}_i^2}, \end{gather}$$
$$\begin{gather}c_{Zi^22} ={-}\dfrac{\left( 2 \gamma^2-2 \gamma\right) \textit{Ma}_i^4 + 2 \textit{Ma}_i^2 - 4}{\gamma^2 p \textit{Ma}_i^4}\left(\dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right) - \dfrac{\left(3 \gamma-3 \right) F_i \textit{Ma}_i^2+6 F_i}{\gamma^2 p^2 A_i \textit{Ma}_i^4}, \end{gather}$$
$$\begin{gather}c_{Zi^23} ={-}\dfrac{\left( \gamma^2-\gamma^3\right) F_i \textit{Ma}_i^6- F_i \textit{Ma}_i^4 +\left( 6 - 2\gamma\right) F_i \textit{Ma}_i^2-8 F_i}{\gamma^3 p \textit{Ma}_i^6}, \end{gather}$$
$$\begin{gather}c_{Zi^24} ={-}\dfrac{ \left(1-\gamma^2\right)F_i^2 \textit{Ma}_i^4+\left( 5 \gamma -9\right) F_i^2 \textit{Ma}_i^2+14 F_i^2}{\gamma^3 p^2 A_i \textit{Ma}_i^6}, \end{gather}$$
$$\begin{gather}c_{Zi^25} = \dfrac{\left( 3 \gamma-3\right) {{F_i}^{3}} \textit{Ma}_i^2+6 {{F_i}^{3}}}{\gamma^3 p^3 A_i^2 \textit{Ma}_i^6}. \end{gather}$$Grouping the terms excluding the second derivatives gives
\begin{equation} c_{Zi^22-5} = c_{Zi^22} \left( \dfrac{{\rm d}F_i}{{\rm d}\kern0.7pt x}\right) + c_{Zi^23} \left(\dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right)^{2} + c_{Zi^24} \left(\dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right) + C_{Zi^25}. \end{equation}The combination of (B23), (B24) and (B21) yields the final expression for the second derivative of the numerator,
\begin{equation} \dfrac{{\rm d}^2N}{{{\rm d}\kern0.7pt x}^2} = c_{N^20} \left(\dfrac{{\rm d}}{{\rm d}\kern0.7pt x}\left(\dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right)\right) + c_{N^21}, \end{equation}where
$$\begin{gather} c_{N^20} = \sum_{i\in[p,s]} \left(c_{Zi^20} + c_{Zi^21} \, c_{Fi^20}\right), \end{gather}$$
$$\begin{gather}c_{N^21} = \dfrac{{\rm d}^3A}{{{\rm d}\kern0.7pt x}^3} + \sum_{i\in[p,s]} \left(c_{Zi^21} c_{Fi^21} + c_{Zi^22-5}\right), \end{gather}$$B.4. The compound choking indicator
The second derivative of 
$\beta _i$ is given by the following equation:
\begin{align} \dfrac{{\rm d}^2\beta_i}{{{\rm d}\kern0.7pt x}^2} = c_{\beta i^2 0} \left( \dfrac{{\rm d}}{{\rm d}\kern0.7pt x} \left(\dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right)\right) +c_{\beta i^2 1} \left( \dfrac{{\rm d}F_i}{{\rm d}\kern0.7pt x}\right) + c_{\beta i^2 2}\left(\dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right)^{2} + c_{\beta i^2 3} \left(\dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right) + c_{\beta i^2 4} , \end{align}where
$$\begin{gather} c_{\beta i^2 0} = \dfrac{A_i {\textit{Ma}_i^4}+\left( \gamma-3 \right) A_i \textit{Ma}_i^2+3 A_i}{\gamma^2 \textit{Ma}_i^4}, \end{gather}$$
$$\begin{gather}c_{\beta i^2 1} = \dfrac{\left( \gamma-1\right) {\textit{Ma}_i^4}+\left( 3 -2 \gamma\right) {\textit{Ma}_i^2}-3 }{\gamma^2 p \textit{Ma}_i^4}, \end{gather}$$
$$\begin{gather}c_{\beta i^2 2} = \dfrac{- A_i \textit{Ma}_i^6 +\left( \gamma^2-5 \gamma+7\right) A_i {\textit{Ma}_i^4}+\left( 9 \gamma-18\right) A_i \textit{Ma}_i^2+15 A_i}{\gamma^3 \textit{Ma}_i^6}, \end{gather}$$
$$\begin{gather}c_{\beta i^2 3} = \dfrac{\left( 1-\gamma^2\right) F_i \textit{Ma}_i^6 -\left( 2 \gamma^2 -10 \gamma+10\right) F_i \textit{Ma}_i^4+\left( 30-19 \gamma\right) F_i \textit{Ma}_i^2-27 F_i}{\gamma^3 p \textit{Ma}_i^6}, \end{gather}$$
$$\begin{gather}c_{\beta i^2 4} = \dfrac{\left( 2 {\gamma^2}-5 \gamma+3\right) F_i^2 \textit{Ma}_i^4+\left( 10 \gamma-12\right) F_i^2 \textit{Ma}_i^2+12 F_i^2}{\gamma^3 p^2 A_i \textit{Ma}_i^6}. \end{gather}$$Grouping the terms excluding the second derivative of the pressure gives
\begin{equation} c_{\beta i^2 1-4} = c_{\beta i^2 1} \left( \dfrac{{\rm d}F_i}{{\rm d}\kern0.7pt x}\right) + c_{\beta i^2 2}\left(\dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right)^{2} + c_{\beta i^2 3} \left(\dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right) + c_{\beta i^2 4}, \end{equation}
the second derivative of 
$\beta$ becomes
\begin{equation} \dfrac{{\rm d}^2\beta}{{{\rm d}\kern0.7pt x}^2} = c_{\beta^20} \left(\dfrac{{\rm d}}{{\rm d}\kern0.7pt x}\left(\dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right)\right) + c_{\beta^21}, \end{equation}where
$$\begin{gather} c_{\beta^20} = \sum_{i\in[p,s]} c_{\beta i^20}, \end{gather}$$
$$\begin{gather}c_{\beta^21} = \sum_{i\in[p,s]} c_{\beta i^2 1-4}, \end{gather}$$B.5. Static pressure
The second derivative of the static pressure in the sonic section is given by
\begin{equation} \dfrac{{\rm d}}{{\rm d}\kern0.7pt x} \left(\dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right) = \dfrac{{\rm d}}{{\rm d}\kern0.7pt x}\left(\dfrac{N}{\beta}\right) = \dfrac{\beta ({\rm d}N/{{\rm d}\kern0.7pt x}) - N ({\rm d}\beta/{{\rm d}\kern0.7pt x})}{\beta^2}. \end{equation}
The equation above can be evaluated using (A21) and (A25) if 
$\beta \neq 0$. Otherwise, l'Hôpital's rule yields the following equation:
\begin{equation} \left(\dfrac{{\rm d}}{{\rm d}\kern0.7pt x} \left(\dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right)\right)^* = \dfrac{0}{0} = \dfrac{({\rm d}^2N/{{\rm d}\kern0.7pt x}^2)^*}{2 ({\rm d}\beta/{{\rm d}\kern0.7pt x})^*} - \dfrac{N^*}{\beta^*}\dfrac{({\rm d}^2\beta/{{\rm d}\kern0.7pt x}^2)^*}{2({\rm d}\beta/{{\rm d}\kern0.7pt x})^*}, \end{equation}
where 
$N^*/\beta ^*$ equals the static pressure gradient in the sonic section (see (2.28) and (A28)). The second derivatives 
$\textrm {d}^2N/{\textrm {d}{\kern0.9pt}x}^2$ and 
$\textrm {d}^2\beta /{\textrm {d}{\kern0.9pt}x}^2$ are given by (B32) and (B42), and are linear with respect to the second derivative of the static pressure,
\begin{align} \left(\dfrac{{\rm d}}{{\rm d}\kern0.7pt x} \left(\dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right)\right)^* &= \dfrac{1}{2 ({\rm d}\beta/{{\rm d}\kern0.7pt x})^*} \left(c_{N^20}^* - \left(\dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right)^* c_{\beta^20}^* \right) \left(\dfrac{{\rm d}}{{\rm d}\kern0.7pt x}\left(\dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right)\right)^* \nonumber\\ &\quad + \dfrac{1}{2 ({\rm d}\beta/{{\rm d}\kern0.7pt x})^*} \left(c_{N^21}^* - \left(\dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right)^* c_{\beta^21}^*\right). \end{align}The second derivative of the static pressure in the sonic section is finally given by
\begin{equation} \left(\dfrac{{\rm d}}{{\rm d}\kern0.7pt x}\left(\dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right)\right)^* = \dfrac{- c_{N^21}^* + \left(\dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right)^* c_{\beta^21}^*}{c_{N^20}^* - \left(\dfrac{1}{p}\dfrac{{\rm d}p}{{\rm d}\kern0.7pt x}\right)^* c_{\beta^20}^* - 2 \left(\dfrac{{\rm d}\beta}{{\rm d}\kern0.7pt x}\right)^*} , \end{equation}which requires (A25), (A28), (B33), (B34), (B43) and (B44). The second-order approximation of the static pressure gradient in the sonic section follows from the Taylor expansion in (3.2), for which the first and second derivatives are given by (A21), (A25), (B32) and (B42). A quasi-third order-approximation is obtained by computing the second derivatives on multiple grid points and applying finite differences to approximate the required third derivatives.
