2. Set-up

The approach uses a control volume with a combination of heat transfer, friction and area change. A control volume that illustrates this combination is shown in figure 1. The surface of the control volume is divided into the inlet, the outlet and an impermeable wall. Additionally, there is no mass flow or velocity component normal to the wall. The flow properties are not assumed uniform across any control surface.

Figure 1. Control volume with simultaneous heat transfer, friction and area change.

The problem considers the $x$-direction momentum equation only. Further, the problem is assumed to be steady with only the $x$-component of velocity being non-zero at the inlet and exit locations. Additionally, the calorically perfect ideal gas assumption is used, which implies that the specific heats $c_P$ and $c_{\mathcal {V}}$ are constant, and the ideal gas law $P = \rho R T$ applies. The Mach number is given as $M = \sqrt {\gamma R T}$. The ratio of the stagnation enthalpy ($h_{o} \equiv h + {V^2}/{2}$) to the static value for a calorically perfect gas is

(2.1)\begin{equation} \frac{h_o}{h} = \frac{T_o}{T} = 1 + \frac{\gamma -1}{2}\,M^2 \equiv f(M),\end{equation}

where the function $f(M)$ is defined for convenience as it appears frequently in the analysis.

An important aspect of the analytical solution presented below is the ability to consider components in which the inflow or exit flow has non-uniform conditions. This non-uniformity could apply to any flow field variable, including velocity, density and temperature. Fluid variables at the inlet and exit of the control volume will be represented by appropriate average values. This will include conditions where more than one inflow stream is present, by considering a single inflow boundary with non-uniform conditions.

Two types of averaging will be considered here: the area average and the mass average. The area average of a quantity $\varphi$ is defined as

(2.2)\begin{equation} \bar{\varphi}^A \equiv \frac{1}{A} \int \varphi\, {\rm d} A. \end{equation}

The mass average of a quantity $\varphi$ is defined as

(2.3)\begin{equation} \bar{\varphi}^m \equiv \frac{1}{\dot{m}} \int \rho \varphi ( \boldsymbol{V} \boldsymbol{\cdot} \hat{\boldsymbol{n}})\, {\rm d} A. \end{equation}

Note that many other types of averages could be considered, as discussed in Cumpsty & Horlock (Reference Cumpsty and Horlock2005), and the choice of a particular averaging type is based on how it helps to simplify the analysis. Pressure, for example, most often appears in this analysis as an area average, $\bar {P}^A$, as this replaces an integral quantity that occurs naturally in the momentum equation (Cumpsty & Horlock Reference Cumpsty and Horlock2005). Similarly, the mass-averaged value of stagnation enthalpy, $\bar {h}^{m}_o$, replaces an integral that occurs in the control volume energy equation (Cumpsty & Horlock Reference Cumpsty and Horlock2005). The mass-averaged stagnation temperature for a calorically perfect ideal gas may then be calculated as

(2.4)\begin{equation} \bar{T}^m_o = \frac{\bar{h}^m_o}{c_P}. \end{equation}

For the present work, a combination of area averaging and mass averaging will be used most often. In cases where neither area or mass averaging is specified, a simple overbar will be used: $\overline {({\cdot })}$. Note that for uniform flows, the use of averaging is not required since all relevant quantities are constant along that particular control surface, and the overbar notation can be neglected, i.e. $\bar {\varphi }^A = \bar {\varphi }^m = \varphi$. Mass, area or any other averaging can be used for the Mach number, so it is left as $\bar {M}$. The solution will remain exact with this arbitrarily chosen average $\bar {M}$ by using ‘flux correction factors’, which are discussed later. Thus one should choose an averaging technique for $\bar {M}$ that provides convenience for the problem of interest.

The use of area-averaged static pressure, mass-averaged stagnation temperature, an average Mach number, and flux correction factors allows the problem to be solved for both uniform and non-uniform quasi-one-dimensional flow problems. The following quantities are defined so that integral forms of mass, momentum and energy conservation are visually similar to their uniform flow counterparts. The function $\tilde {f}$ is analogous to the stagnation-to-static enthalpy ratio and is defined as

(2.5)\begin{equation} \tilde{f} \equiv 1 + \frac{\gamma - 1}{2}\,\bar{M}^2. \end{equation}

Note that $\tilde {f}$ is not an average quantity itself; rather, it is defined using other average quantities. A representative static temperature can then be defined as

(2.6)\begin{equation} \tilde{T} \equiv \frac{\bar{T}^m_o}{\tilde{f}}. \end{equation}

A representative static density can then be defined as

(2.7)\begin{equation} \tilde{\rho} \equiv \frac{\bar{P}^A}{R\tilde{T}} = \frac{\bar{P}^A}{R}\, \frac{\tilde{f}}{\bar{T}^m_o}. \end{equation}

A representative velocity can then be defined as

(2.8)\begin{equation} \tilde{V} \equiv \bar{M} \sqrt{\gamma R \tilde{T}} = \bar{M} \sqrt{\gamma R\,\frac{\bar{T}^m_o}{\tilde{f}}}. \end{equation}

Similar to $\tilde {f}$, the quantities $\tilde {T}$, $\tilde {\rho }$ and $\tilde {V}$ are not average quantities themselves; rather, they are representative quantities defined in terms of other averages, i.e. $\bar {P}^{A}$, $\bar {T}^m_o$ and $\bar {M}$. They are analogous to their uniform-flow counterparts $f$, $T$, $\rho$ and $V$. The use of these definitions allows the formulation to be mathematically exact while maintaining a familiar form to the equations.

Conventionally, integral quantities in the control volume formulations with non-uniform profiles are represented using flux correction factors (also called non-uniformity coefficients). These are scalars that serve to modify values computed from spatially averaged quantities in order to preserve an integrated quantity of interest, e.g. momentum flux. Correction factors for compressible flow are discussed in Crocco (Reference Crocco and Emmons1958). For incompressible flow, these factors are described in more detail in White (Reference White2006) and Munson et al. (Reference Munson, Huebsch, Rothmayer, Gerhart, Gerhart and Hochstein2015).

The mass flux correction factor is defined here as

(2.9)\begin{equation} \alpha \equiv \frac{\displaystyle \int \rho ( \boldsymbol{V} \boldsymbol{\cdot} \hat{\boldsymbol{n}})\, {\rm d} A}{\tilde{\rho}\tilde{V}A}. \end{equation}

This term is equal to unity for constant density flows when using an area-averaged velocity. However, when the density is not spatially uniform, it is possible for $\alpha$ to be greater than or less than unity. For example, if spatial variations in $\rho$ and $V$ are positively correlated over the inlet surface, then the mass flow rate is $\dot {m} > \tilde {\rho }\tilde {V}A$, hence $\alpha > 1$. Similarly, if density and velocity variations are negatively correlated over the area, then $\alpha < 1$. Inserting the definitions of $\tilde {\rho }$ and $\tilde {V}$ shows its dependence on normalized profiles of static pressure, stagnation temperature, Mach number and stagnation-to-static enthalpy:

(2.10)\begin{equation} \alpha = \frac{1}{A} \int \frac{P}{\bar{P}^A} \left(\frac{T_o}{\bar{T}^m_o} \right)^{-{1}/{2}} \frac{\ M\ }{\bar{M}} \sqrt{\frac{\ f\ }{\tilde{f}}}\, {\rm d} A.\end{equation}

This provides a useful method to estimate $\alpha$ values if there is knowledge regarding the shapes of the relevant dimensionless profiles.

The momentum flux correction factor is defined as

(2.11)\begin{equation} \beta \equiv \frac{\displaystyle \int \rho V ( \boldsymbol{V} \boldsymbol{\cdot} \hat{\boldsymbol{n}} )\, {\rm d} A}{\tilde{\rho}\tilde{V}^2 A}. \end{equation}

This follows from the classical definition given in White (Reference White2006) and Munson et al. (Reference Munson, Huebsch, Rothmayer, Gerhart, Gerhart and Hochstein2015). Substituting the definitions of $\tilde {\rho }$ and $\tilde {V}$ shows its dependence on the normalized profiles of static pressure and Mach number:

(2.12)\begin{equation} \beta = \frac{1}{A} \int \frac{P}{\bar{P}^A} \left(\frac{\ M }{\bar{M}}\right)^2 {\rm d} A. \end{equation}

Thus the momentum flux correction factor can be evaluated as the area average of a combination of the normalized profiles over the control surface.

These flux correction factors allow for any averaging method to be used for $\bar {M}$. It requires only consistent use throughout the equations. Careful selection of the averaging method can guarantee certain properties of $\beta$. For example, using Schwarz's inequality, it can be shown that if $\bar {M} = \bar {M}^A$, then $\beta \geq 1$.

3. Derivation of control volume solution

The conservation of mass for a steady, single-inlet, single-exit control volume is given as

(3.1)\begin{equation} \alpha_ 2 \tilde{\rho}_2 \tilde{V}_2 A_2 - \alpha_1 \tilde{\rho}_1 \tilde{V}_1 A_1 = 0, \end{equation}

where the subscript 2 represents exit conditions, and the subscript 1 indicates inlet conditions.

Substituting the ideal gas law and the definition of the Mach number into (3.1), and using (2.1), gives

(3.2)\begin{equation} \frac{\bar{P}_2^A A_2}{\bar{P}_1^A A_1}\,\frac{\alpha_2}{\alpha_1} \sqrt{\frac{\bar{T}_{o1}^m}{\bar{T}_{o2}^m}}\,\frac{\bar{M}_2}{\bar{M}_1} \left( \frac{\tilde{f}_2}{\tilde{f}_1} \right)^{1/2} - 1 = 0.\end{equation}

It should be noted that in (3.2), the static pressure is used rather than the stagnation pressure. Further, the stagnation temperature is used rather than the static temperature. The reason for this choice is that the static pressure and stagnation temperature occur in the momentum and energy equations, respectively. Additionally, the pressure occurs with an integer exponent. Thus those equations require no additional manipulation for all three governing equations to contain the same set of unknown quantities.

The steady $x$-momentum balance equation is

(3.3)\begin{equation} \beta_2 \tilde{\rho}_2 \tilde{V}_2^2 A_2 - \beta_1 \tilde{\rho}_1 \tilde{V}_1^2 A_1 = \bar{P}_1^A A_1 - \bar{P}_2^A A_2 + \bar{P}_3^A ( A_2 - A_1 ) - F_{fx}, \end{equation}

where $\bar {P}_1^A$ and $\bar {P}_2^A$ represent the area-averaged static pressure on inlet and outlet control surfaces, respectively. Here, $F_{fx}$ is the friction force acting on the control volume projected in the $x$-direction. The term $\bar {P}_3^A (A_2-A_1)$ represents the axial force on the fluid due to the wall pressure. Therefore, $\bar {P}_3^A$ represents an average over the projected area of the walls of the device.

The area-averaged wall pressure $\bar {P}_3^A$ will be of primary importance in the solution for any problems that involve area change. In many cases, this term will be the only unknown needed for closure of the equations, hence it must be specified in order to solve for the control volume exit conditions. Without loss of generality, $\bar {P}_3^A$ may be prescribed as a linear combination of the inlet static pressure $\bar {P}_1^A$, the exit static pressure $\bar {P}_2^A$, and an average inlet stagnation pressure $\bar {P}_{o1}$:

(3.4)\begin{equation} \frac{\bar{P}_3^A}{\bar{P}_1^A} = (1-\xi-\nu) + \xi\, \frac{\bar{P}_2^A}{\bar{P}_1^A} + \nu\, \frac{\bar{P}_{o1}}{\bar{P}_1^A},\end{equation}

where $\xi$ and $\nu$ represent wall-pressure scaling coefficients. Note that the wall pressure was prescribed in terms of the known inlet conditions and the unknown outlet static pressure, similar to Shapiro (Reference Shapiro1953) and John & Keith (Reference John and Keith2006). Therefore, the wall pressure is prescribed implicitly. This leads to a closed-form solution for the exit conditions rather than an iterative solution. Determining the appropriate wall-pressure scaling coefficients for various problems is discussed in later sections.

The net friction force per unit area is written as a friction coefficient multiplied by $\frac {1}{2} \rho V^2$. The net friction force is then calculated by integrating over the wall area and taking the $x$-component of the resulting force. The integration of $\frac {1}{2} \rho V^2$ cannot be performed directly, thus the integral is approximated using the mean of the inlet and outlet $\rho V^2$ values:

(3.5)\begin{equation} F_{fx} \equiv \int_1^2 \frac{c_f}{2}\,\rho V^2 \cos(\theta_3)\, {\rm d} A_3 \approx \frac{c_f}{4}\,\gamma ( \beta_1 \bar{P}_1^A \bar{M}_1^2 + \beta_2 \bar{P}_2^A \bar{M}_2^2 ) A_{3f}, \end{equation}

where $c_f$ is the friction coefficient, $\theta _3$ is the angle from the $x$-direction to the wall, and $A_{3f}$ is defined as

(3.6)\begin{equation} A_{3f} \equiv \int_1^2 \cos(\theta_3)\, {\rm d} A_3. \end{equation}

The momentum equation then becomes

(3.7)\begin{align} &1 + \beta_1 \gamma \left(1 - \frac{c_f}{4}\,\frac{A_{3f}}{A_1} \right) \bar{M}_1^2 + \left(1-\xi + \nu \left( \frac{\bar{P}_{o1}}{\bar{P}_1^A} - 1 \right) \right) \left( \frac{A_2}{A_1} - 1\right) \nonumber\\ &\quad =\frac{\bar{P}_2^A A_2}{\bar{P}_1^A A_1} \left( 1 + \beta_2 \gamma \left(1 + \frac{c_f}{4}\,\frac{A_{3f}}{A_2} \right) \bar{M}_2^2 + \xi \left( \frac{A_1}{A_2} - 1 \right)\right). \end{align}

The steady energy conservation equation without power transfer across the control surface, consistent with figure 1, is

(3.8)\begin{equation} \bar{h}_{o2}^m - \bar{h}_{o1}^m = \frac{\dot{Q}}{\dot{m}} \equiv q.\end{equation}

This can be written in terms of stagnation temperature as

(3.9)\begin{equation} \frac{\bar{T}_{o2}^m}{\bar{T}_{o1}^m} = \frac{q}{\bar{h}_{o1}^m} + 1.\end{equation}

Substitution of (3.9) into (3.2), solving for $(\bar {P}_2^A A_2 ) / (\bar {P}_1^A A_1)$, and substitution into (3.7) results in

(3.10)\begin{align} &\frac{ 1 + \beta_1 \gamma \left(1 - \dfrac{c_f}{4}\, \dfrac{A_{3f}}{A_1} \right) \bar{M}_1^2 + \left(1-\xi + \nu \left( \dfrac{\bar{P}_{o1}}{\bar{P}_1^A} - 1 \right) \right) \left( \dfrac{A_2}{A_1} - 1\right)}{\dfrac{\alpha_1}{\alpha_2}\,\bar{M}_1 \sqrt{\tilde{f}_{1} \left( \dfrac{q}{\bar{h}_{o1}^m} + 1 \right)}}\nonumber\\ &\quad =\frac{1 + \beta_2 \gamma \left(1 + \dfrac{c_f}{4}\, \dfrac{A_{3f}}{A_2} \right) \bar{M}_2^2 + \xi \left( \dfrac{A_1}{A_2} - 1 \right)}{\bar{M}_2 \sqrt{\tilde{f}_{2}}}. \end{align}

All conditions at the inlet are assumed to be known, as well as the friction coefficient, heat transfer and wall-pressure scaling coefficients for cases with area change. Equation (3.10) thus contains only a single unknown: $\bar {M}_2$. Squaring both sides of (3.10) and simplifying results in the quadratic equation

(3.11)\begin{equation} (\bar{M}_2^2)^2 + b (\bar{M}_2^2) + c =0.\end{equation}

The quadratic coefficients are

(3.12)\begin{equation} \left.\begin{array}{c} b \equiv \dfrac{\varLambda - 2\gamma \beta_2 \left[ 1 + \xi \left(\dfrac{A_1}{A_2} - 1\right)\right] \left( 1 + \dfrac{c_f}{4}\, \dfrac{A_{3f}}{A_2} \right)}{\varLambda\,\dfrac{\gamma - 1}{2} - \left[\gamma \beta_2 \left( 1 + \dfrac{c_f}{4}\,\dfrac{A_{3f}}{A_2} \right) \right]^2} \\ \text{and} \quad c \equiv{-} \dfrac{\left[ 1 + \xi \left(\dfrac{A_1}{A_2} - 1\right)\right]^2}{\varLambda\,\dfrac{\gamma - 1}{2} - \left[\gamma \beta_2 \left( 1 + \dfrac{c_f}{4}\,\dfrac{A_{3f}}{A_2} \right) \right]^2}, \end{array}\right\} \end{equation}

where the term $\varLambda$ was introduced only for brevity as

(3.13)\begin{equation} \varLambda \equiv \left[\dfrac{ 1 + \beta_1 \gamma\left(1 - \dfrac{c_f}{4}\, \dfrac{A_{3f}}{A_1} \right) \bar{M}_1^2 + \left(1-\xi + \nu \left( \dfrac{\bar{P}_{o1}}{\bar{P}_1^A} - 1 \right)\right)\left( \dfrac{A_2}{A_1} - 1\right) }{\dfrac{\alpha_1}{\alpha_2}\,\bar{M}_1 \sqrt{\tilde{f}_{1} \left( \dfrac{q}{\bar{h}_{o1}^m} + 1 \right)}}\right]^2.\end{equation}

Physically, $\varLambda$ is the squared sum of all known momenta non-dimensionalized by the square of the mass flow rate and the stagnation enthalpy. This can be shown to be equivalent to the non-dimensionalization presented by Denton (Reference Denton1993) and Greitzer et al. (Reference Greitzer, Tan and Graf2004) for the two-stream mixing problem.

4. Solution procedure

The unknown exit Mach number represents the full solution. This can now be found from the quadratic equation as

(4.1)\begin{equation} \bar{M}_2^2 ={-}\frac{b}{2} \pm \sqrt{\left(\frac{b}{2}\right)^2 - c}.\end{equation}

This represents a closed-form algebraic solution to the generalized internal flow problem that does not require a numerical solution or iterative methods.

The final issue to consider is that the solution was obtained from a bi-quadratic, and hence allows for four roots as potential solutions. Physically realizable solutions may be determined by considering two restrictions. First, the exit Mach number must be a real positive number. This requirement eliminates all roots for which $\bar {M}_2 < 0$ (i.e. reversed outlet flow), $\bar {M}_2^2 <0$ and $(b/2)^2 - c < 0$ (i.e. imaginary $\bar {M}_2$). A single solution is possible for any real $b$ value given that $|b/2| < \sqrt { ( b/2)^2 - c}$ . Two solutions are possible if $b<0$ and $|b/2| > \sqrt { ( b/2)^2 - c}$. Whether or not these possible solutions will be supersonic or subsonic will depend on $|{-b/2} - 1|$, $\sqrt {( b/2)^2-c}$, and the choice of ‘$\pm$’. Hereafter, the choice ‘$+$’ will be referred to as a ‘positive root’, and ‘$-$’ will be referred to as a ‘negative root’. The results are summarized in table 1. The final restriction requires eliminating non-physical solutions by satisfying the second law:

(4.2)\begin{equation} {\rm \Delta} s \ge \int_1^2 \frac{\delta q}{T}.\end{equation}

Table 1. Conditions for $\bar {M}_2$ solutions.

The stated restrictions reduce the number of allowable solutions to either one or two. In the event that two $\bar {M}_2$ solutions are obtained, both should be considered to be physically realizable for the boundary conditions provided, e.g. one subsonic and one supersonic. This will be demonstrated in the following sections, with examples of both canonical flow problems and several unique solutions that were made possible by the formulation given by (3.11).

Once the solution for $\bar {M}_2$ is known, all other quantities at the exit may be obtained. For example, the exit pressure $\bar {P}_2^A$ can be found from the conservation of mass, i.e. (3.1). Hence all exit conditions, including entropy change, can be found from the solution.

The control volume solution as derived can be used for any internal flow component, as long as the required inputs are known or can be approximated. However, there may be cases where there is an advantage to dividing the flow into multiple sub-elements in the streamwise direction. An example of this is very long geometries with friction, such as long pipe flow problems. In this case, the approximation to the friction force given by (3.5) is improved if the component is segmented into multiple short sections. The control volume solution still provides a closed-form solution, but a streamwise-direction marching approach is applied. This provides a closed-form solution to the discretized problem, which is defined as ‘semi-analytical’ (De Vidts & White Reference De Vidts and White1992). This can be shown to provide second-order convergence with increasing sub-element count.

A second example where multiple sub-elements may be useful is for flows with continuously varying cross-sectional area in which the wall-pressure scaling ($\xi$) is unknown. The wall-pressure scaling is set to $\xi =1/2$, $\nu =0$, which approximates the wall pressure as the mean of inlet and outlet static pressure. This is commonly used with the steady one-dimensional differential equations (Shapiro Reference Shapiro1953; John & Keith Reference John and Keith2006; Kundu et al. Reference Kundu, Cohen and Dowling2015). Effectively, this means that the wall pressure is being solved for along the length of overall control volume. Like the case of friction, this can be shown to provide second-order convergence with increasing sub-element count. Convergence to the solution occurs regardless of the value of $\xi$ as long as $\nu =0$, though the rate of convergence varies with the actual value of $\xi$. This result indicates that in the limit of infinitely many sub-elements (i.e. a differential size control volume), the solution presented in this work is consistent with the solution of the differential form of the governing equations for steady, compressible, quasi-one-dimensional flow.

Problems that include arbitrary area change may be solved with a single control volume provided that the values of $\xi$ and $\nu$ are given. For example, in a flow with gradual area change, it may be assumed that the wall pressure varies between the inlet and exit static pressure values. This would be independent of the inlet stagnation pressure. Hence $\nu =0$, and choosing $\xi$ in the range $0 \leq \xi \leq 1$ will set $\bar {P}_3^A$ between the inlet and exit static pressures. In problems where the inlet stagnation pressure is more appropriate, with a relevant constraint $1 - \xi - \nu = 0$, and choosing $0 \leq \nu \leq 1$, will set $\bar {P}_3^A$ between the exit static pressure and the inlet stagnation pressure. Two examples, which set $\xi$ and $\nu$ within these ranges, are the cases of sudden expansion and sudden contraction. These are discussed in § 6.

Alternatively, if the values of $\xi$ and $\nu$ are not readily apparent from the geometry, the values may be varied over some range and selected to achieve a specific outlet condition. For example, a known component efficiency or entropy generation may be used as a constraint.

Finally, special cases allow for closed-form solutions for the wall pressure. For example, a closed-form solution for isentropic flow can be found if $M_2$ is sufficiently small. Here, the overbar notation can be neglected since the flow is assumed uniform and the flux correction factors are unity. First, the momentum equation is rearranged to obtain an explicit expression for $\xi$ with $\nu =0$:

(4.3)\begin{equation} \xi = \frac{1}{\dfrac{P_2}{P_1}-1} \left( \frac{\dfrac{P_2}{P_1}\, \dfrac{A_2}{A_1}\,( 1 + \gamma M_2^2) - ( 1 + \gamma M_1^2 )}{\dfrac{A_2}{A_1}-1} - 1 \right). \end{equation}

Two new relationships are required to substitute for $P_2/P_1$ and $M_2$ such that a closed-form expression can be determined based on known inlet conditions. The assumption of isentropic flow can be used so that the static pressure ratio becomes equivalent to

(4.4)\begin{equation} \frac{P_2}{P_1} = \frac{P_2}{P_{o2}}\,\frac{P_{o1}}{P_1}. \end{equation}

The outlet-stagnation-to-static pressure ratio may be expressed as a Taylor series about $M_2^2=0$:

(4.5)\begin{equation} \frac{P_{o2}}{P_2} = \left( 1 + \frac{\gamma-1}{2}\, M_2^2\right)^{\gamma/({\gamma-1})} = 1 + \frac{\gamma}{2}\,M_2^2 + \frac{\gamma}{8}\,M_2^4 + \frac{\gamma}{48}\,(2-\gamma)\,M_2^6 + O(M_2^8). \end{equation}

Neglecting terms of order $M_2^4$ and higher allows a solution for $M_2$ to be found using mass and energy conservation. The roots of the resulting polynomial are the solutions to the two-term isentropic flow approximation

(4.6)\begin{equation} M_{2s} \approx \sqrt{\frac{M_1^2}{\left( \dfrac{A_2}{A_1} \right)^2 f_1^{({\gamma+1})/({\gamma-1})} - \gamma M_1^2 }},\end{equation}

where $M_{2s}$ is the isentropic exit Mach number. Substituting the first two terms of (4.5) and (4.6) into (4.3) yields the isentropic value of the wall-pressure scaling coefficient $\xi _s$:

(4.7)\begin{align} \xi_s &\approx \frac{1+ \dfrac{\dfrac{\gamma}{2}\, M_1^2}{\left(\dfrac{A_2}{A_1}\right)^2 f_1^{({\gamma+1})/({\gamma-1})}-\gamma M_1^2}}{f_1^{{\gamma}/({\gamma-1})}-\left(1 + \dfrac{\dfrac{\gamma}{2}\, M_1^2}{\left(\dfrac{A_2}{A_1}\right)^2 f_1^{({\gamma+1})/({\gamma-1})}-\gamma M_1^2}\right)} \nonumber\\ &\quad \times \left[ f_1^{{\gamma}/({\gamma-1})}\,\dfrac{\dfrac{A_2}{A_1}}{\dfrac{A_2}{A_1}-1} \left(\dfrac{ 1 + \dfrac{\gamma M_1^2}{\left(\dfrac{A_2}{A_1}\right)^2 f_1^{({\gamma+1})/({\gamma-1})}-\gamma M_1^2}}{1 + \dfrac{\dfrac{\gamma}{2}\, M_1^2}{\left(\dfrac{A_2}{A_1}\right)^2 f_1^{({\gamma+1})/({\gamma-1})}-\gamma M_1^2}}\right) - \dfrac{1+\gamma M_1^2}{\dfrac{A_2}{A_1}-1}-1\right]. \end{align}

This relationship for $\xi _s$, along with (4.1), provides a closed-form solution to isentropic flow with area change for arbitrary inlet Mach numbers and small exit Mach numbers. This result will be accurate as long as $P_{o2}/P_2$ is well approximated as $1 + (\gamma /2)M_2^2$. The leading-order error term of this approximation is $(\gamma /8)M_2^4$. This approximation is accurate to within approximately $0.1\,\%$ for $M_2<0.3$, and approximately $1\,\%$ for $M_2<0.5$ for $\gamma =7/5$. Retaining additional terms in (4.5) would improve the approximation for isentropic flow. In principle, closed-form solutions using more terms could be found as long as the roots of the resulting polynomial are able to be solved explicitly. An additional approximation for isentropic flow without the assumption of small $M_{2s}$ is given in Appendix A.

5. Comparison to classical solutions

The utility of the derived solution will be demonstrated with a series of example problems. First, several of the classical one-dimensional compressible flow problems will be shown. This is to show consistency with established theory, as well as to demonstrate the variety of boundary conditions that can be considered from the same solution, i.e. (4.1). A later section will compare the results of (4.1) with experimental measurements and other recent solutions from the available literature.

The first problem considered is isentropic flow with area change. The flow is assumed uniform for this analysis, so the overbar notation is neglected and the flux correction factors are unity. An analytical relationship between the area ratio $A_2/ A_1$, the inlet Mach number $M_1$, and the isentropic exit Mach number $M_{2s}$, can be found from (3.2) by setting the stagnation pressure, stagnation temperature and mass-flux correction factor ratios to unity, resulting in

(5.1)\begin{equation} \frac{A_2}{A_1} = \frac{M_1}{M_{2s}} \left( \frac{1+ \dfrac{\gamma-1}{2}\, M_{1}^2}{1+ \dfrac{\gamma-1}{2}\,M_{2s}^2} \right)^{{-(\gamma+1)}/{2(\gamma-1)}}. \end{equation}

The solution in the form of $M_{2s}$ for a given area ratio and $M_1$ value can be found numerically but has no closed-form solution.

The isentropic flow solution of (4.1) requires incorporating a constraint ${\rm \Delta} s=0$. This additional constraint effectively sets the wall-pressure coefficients to their isentropic values for the given boundary conditions. There are three methods to incorporate this new constraint into the solution: approximate the solution using isentropic relations, multiple sub-elements or iteration. The isentropic value of the wall pressure is independent of the inlet stagnation pressure, hence $\nu =0$ for each method. The first method used (4.7) to set $\xi$ and was found to be accurate for small $M_2$ values. The second method implemented multiple sub-elements using $\xi =0.5$. The last method used a single control volume to evaluate the solution as a function of $\xi$. This approach used iteration to determine the value of $\xi$ such that ${\rm \Delta} s=0$.

The results of iteration using a single control volume for isentropic flow are shown in figure 2. The abscissa is the inlet Mach number, and the ordinate is the isentropic exit Mach number. The ‘$\times$’ symbols refer to the negative root, and the ‘$\circ$’ symbols refer to the positive root for the isentropic iterative results from (4.1). These results are compared with (5.1) for area ratio in figure 2(a). The corresponding values of $\xi$ are shown in figure 2(b) and are compared against (4.3). Both the area ratio and $\xi$ values show the iterative results agree with the analytical expressions to numerical precision.

Figure 2. Isentropic flow for $\gamma =7/5$ using a single control volume. (a) Contours of area ratio; lines are from (5.1), and symbols are from (4.1). (b) Contours of $\xi$ with $\nu =0$; lines are from (4.3), and symbols are from (4.1).

The next problem considered is the flow across a normal shock. Again, the overbar notation can be neglected since the flow is assumed uniform and the flux correction factors are unity. Equation (4.1) was solved using $A_2/A_1 = 1$ and $q/\bar {h}_{o1}^m=0$. The entropy change across the shock is used to compare the solution from this work against the literature. This is given by John & Keith (Reference John and Keith2006) as

(5.2)\begin{equation} \frac{s_2-s_1}{c_P} ={-}\frac{\gamma-1}{\gamma} \ln \left[ \left(\frac{\dfrac{\gamma+1}{2}\,M_1^2}{1 + \dfrac{\gamma-1}{2}\,M_1^2} \right)^{{\gamma}/({\gamma-1})} \left( \frac{2\gamma}{\gamma+1}\,M_1^2 - \frac{\gamma-1}{\gamma+1} \right)^{{-1}/({\gamma-1})} \right].\end{equation}

The entropy change is shown over a range of $M_1$ values in figure 3 as a dashed line. Note that (5.2) admits both trivial solutions (${\rm \Delta} s =0$) and non-physical shocks with $M_1 < 1$. The solutions from (4.1) are shown with symbols; the positive and negative roots are shown separately. Interestingly, the negative root gives the non-trivial solution for $M_1>1$ and trivial solution for $M_1 < 1$. The positive root results in a subsonic shock with negative entropy change or the trivial supersonic solution. The results from (4.1) agree with (5.2) to numerical precision.

Figure 3. Entropy generation across a normal shock; the dashed line is from (5.2), and symbols are from (4.1).

Fanno and Rayleigh flows are now considered. Again, the overbar notation can be neglected since the flow is assumed uniform and the flux correction factors are unity. Entropy change is again used as the dependent variable of interest, as it is easy to visualize and compare solutions. The entropy change for Fanno flow is given by Shapiro (Reference Shapiro1953) as

(5.3)\begin{equation} \frac{s-s^*}{c_P} = \ln \left[ M^2 \left(\frac{\dfrac{\gamma+1}{2}}{M^2 \left(1 + \dfrac{\gamma-1}{2}\, M^2\right)}\right)^{({\gamma+1})/{2\gamma}}\right], \end{equation}

where $s^*$ refers to the sonic condition. The entropy change for Rayleigh flow is given in Shapiro (Reference Shapiro1953) and John & Keith (Reference John and Keith2006) as

(5.4)\begin{equation} \frac{s-s^*}{c_P} = \ln \left[ M^2 \left( \frac{\gamma+1}{1+\gamma M^2}\right)^{({\gamma+1})/{\gamma}}\right]. \end{equation}

These results compared to (4.1) are shown in figure 4. The solutions from Shapiro (Reference Shapiro1953) and John & Keith (Reference John and Keith2006) are shown as dashed lines, with the positive and negative roots from (4.1) shown with symbols. The results show the consistency of (4.1) with Rayleigh flows to within numerical precision of the solution. The results also show that the solution follows the correct non-dimensional entropy change versus Mach number for Fanno flow.

Figure 4. Entropy relative to the reference condition ($M=1$) for Fanno flow, shown as a dotted line using (5.3), and Rayleigh flow, shown as a dashed line using (5.4); symbols are from (4.1).

6. Additional example problems

Several distinctive internal flows are now considered that can be considered with the present formulation. This includes sudden expansion, sudden contraction, the two-stream mixing layer, and combinations of friction and heat transfer with area change.

The problem of sudden expansion with compressible flow provides a useful demonstration of the solution presented in this paper. The control volume and problem set-up are shown in figure 5. There are three control surfaces: inlet, outlet and an impermeable wall. There is no heat addition through the walls, and there is no friction present at the walls. The flow is assumed to be uniform at both the inlet and exit of the control volume. Note that this requires the outlet of the control volume to be sufficiently downstream of the inlet so that the uniform assumption is justified. Further, the overbar notation can be neglected and the flux correction factors are set to unity. The pressure on the wall is specified based on the assumption that the streamlines separate horizontally at the expansion point. Hence it can be assumed that $\bar {P}_3^A = \bar {P}_1^A$.

Figure 5. Sudden expansion control volume.

The assumptions noted lead to the parameters given by $q=0$, $\nu =0$ and $\xi = 0$. The quantities $M_1$ and $A_2/A_1$ were considered to be independent variables so that $M_2$ may be determined. Figure 6 shows the solution to (4.1) along with experimental data from Hall & Orme (Reference Hall and Orme1955). The abscissa and ordinate of the figure represent the inlet and exit Mach numbers, respectively. The different lines correspond to varied area ratios.

Figure 6. Sudden expansion solution using (4.1), shown as lines, compared to experimental data, shown as symbols, from Hall & Orme (Reference Hall and Orme1955) for $\xi =0$, $\nu =0$ and $\gamma =7/5$.

The exit Mach number increases as $M_1$ is increased and as $A_2 / A_1$ decreases, as expected. Generally, the agreement between the analytical solution and the experimental data is excellent over the full range of inlet Mach numbers and area ratios. It is noted that the lowest area ratio values do show a small bias in the prediction. This is likely due to the assumption $\bar {P}_3^A = \bar {P}_1^A$, which is likely invalid for smaller area ratios, since the flow reattachment location will be closer to the expansion location, thereby altering the streamline pattern. A more accurate estimate of the wall pressure, using theoretical values or empirical relations (Marjanovic & Djordjevic Reference Marjanovic and Djordjevic1994; Morimune, Hirayama & Maeda Reference Morimune, Hirayama and Maeda1980), would improve the predictions for area ratios closer to unity.

The compressible sudden contraction, like sudden expansion, provides another useful demonstration of the solution presented in this paper. The control volume and problem set-up are shown in figure 7. The assumptions are identical to the sudden expansion just discussed. This includes the assumption that the measurements used are sufficiently downstream so that the uniform flow assumption is justified. Similar to the case of sudden expansion, the overbar notation is also dropped and the flux correction factors are set to unity. The wall pressure in this case, however, can be assumed to be close to the inlet stagnation pressure (Marjanovic & Djordjevic Reference Marjanovic and Djordjevic1994). That is, $\bar {P}_3^A = \bar {P}_{o1}$. Further, the outlet control surface is placed sufficiently downstream of the vena contracta to justify the uniform flow assumption. The problem set-up leads to $q=0$, $\nu =1$ and $\xi = 0$.

Figure 7. Sudden contraction control volume.

Figure 8 shows the theoretical solution along with experimental data from Benedict et al. (Reference Benedict, Carlucci and Swetz1966). These measurements were also used by Marjanovic & Djordjevic (Reference Marjanovic and Djordjevic1994). The abscissa represents the exit static-to-stagnation pressure ratio ($P_2/P_{o2}$), which is related directly to $M_2$. The ordinate is the exit-to-inlet stagnation pressure ratio ($P_{o2}/P_{o1}$), which can be related directly to the net entropy change across the control volume. Results from four different area ratio values are shown.

Figure 8. Sudden contraction solution using (4.1), shown as lines, compared to experimental data, shown as symbols, from Benedict et al. (Reference Benedict, Carlucci and Swetz1966) for $\xi =0$, $\nu =1$ and $\gamma =7/5$.

The exit-to-inlet stagnation pressure ratio is observed to decrease with decreasing exit static-to-stagnation pressure ratio (increasing exit Mach number). The theoretical solution agrees well with the experimental data in magnitude and trend, but underpredicts the exit-to-inlet stagnation pressure ratio by as much as $\approx 0.03$. This implies that using the inlet stagnation pressure value for the wall pressure is a useful upper bound, and provides a lower bound to the stagnation pressure ratio through a sudden contraction. It is worth noting that this prediction did not require any information regarding the vena contracta in the outlet pipe section. Similar to the sudden expansion, a more accurate estimate of the wall pressure, using theoretical values or empirical relations (Marjanovic & Djordjevic Reference Marjanovic and Djordjevic1994), would improve the predictions.

The two-stream mixing layer, like the previous examples, provides another useful demonstration of the solution presented in this paper. The control volume and problem set-up are shown in figure 9, applied to the case of two supersonic streams mixing. There are two control surfaces: one inlet and one outlet. Note that the inlet was considered using a non-uniform flow, while the outlet was uniform, such that all outlet flux correction factors were assumed to be unity. Further, the inlet non-uniform profile was considered as a combination of two uniform flow streams, as shown. The ‘primary stream’ is given as the lower inlet stream, and the ‘secondary stream’ is the upper inlet stream, as shown in figure 9. The primary and secondary streams together compose the overall inlet control surface. There was no heat addition through the walls and the area ratio is unity (i.e. $A_2/A_1=1$). Friction was assumed to be significant due to the high flow speed. The friction coefficient was assumed to be $c_f\approx 0.005$, which was based on an estimate by Chow & Addy (Reference Chow and Addy1964). The inlet static pressure was not assumed to be uniform. The inlet mass-averaged Mach number was used as the normalizing Mach number. The presence of a shock structure inside the control volume was handled implicitly by the selection of subsonic or supersonic $\bar {M}_2$, i.e. selection of the positive or negative root in (4.1). This particular case required the subsonic $\bar {M}_2$.

Figure 9. Non-uniform inlet flow with constant area mixing.

Figure 10 shows the theoretical solution along with experimental data from Dutton et al. (Reference Dutton, Mikkelsen and Addy1982). The abscissa represents the primary-stagnation-to-secondary-stream static pressure ratio ($P_{op}/P_{s}$), and the ordinate is the exit-to-secondary-stream static pressure ratio ($P_2/P_{s}$). Results from four different secondary stream Mach numbers and two different primary-to-secondary-stream area ratios are shown. The primary stream Mach number was set to $2.5$. Solid symbols and open symbols of the same shape denote two different Reynolds numbers for the same area ratio and inlet Mach numbers.

Figure 10. Supersonic constant-area two-stream mixing with friction solution using (4.1), shown as lines, compared to experimental data, shown as symbols, from Dutton et al. (Reference Dutton, Mikkelsen and Addy1982) for $M_p=2.50$, $c_f=0.005$, $A_3/A_1=40$ and $\gamma =7/5$.

The exit-to-secondary-stream static pressure ratio values were observed to increase linearly with primary-stagnation-to-secondary-static pressure ratio. The predictions show the same qualitative behaviour and are within $\approx 8\,\%$ error relative to the data. The larger errors are likely due to an increase in the friction coefficient as the Reynolds number is reduced. The direction of the error is also consistent with an increasing friction coefficient, meaning that a larger $c_f$ would lower the predicted exit-to-secondary-stream static pressure ratio. An additional test was reported for $M_p=4.0$ and $M_s=2.0$, which also had maximum error of approximately $8\,\%$. However, this was not included due to the large ordinate and abscissa scale differences relative to the other data.

The combination of simultaneous heat transfer and friction in a constant-area duct provides a useful example of the solution when using multiple sub-elements. A schematic of this problem is shown in figure 11. This problem was recently reduced by Ferrari (Reference Ferrari2021a) to a single implicit algebraic equation, which can then be solved iteratively. Ferrari's solution for simultaneous heat transfer and friction provides a useful comparison for the multiple sub-element solution method presented here.

Figure 11. Schematic of simultaneous friction and heat transfer.

Two separate cases are examined: one with a subsonic inlet, and the other with a supersonic inlet. This comparison is shown in figure 12. The subsonic inlet case is shown in figure 12(a), and the supersonic inlet case is shown in figure 12(b). The abscissa is non-dimensionalized length along the duct, $c_f x/D$. The ordinate is Mach number. The overbar notation can be neglected since the flow is assumed uniform and the flux correction factors are unity. There are multiple curves for various values of a non-dimensionalized heat-transfer-to-friction ratio. The heat-transfer-to-friction ratio is defined as

(6.1)\begin{equation} \varGamma \equiv \frac{\dfrac{q}{\bar{h}_{o1}^m} }{c_f \dfrac{L}{D}}, \end{equation}

where the overall length of the duct is $L$, and the diameter is $D$. Note that the $\varGamma =0$ case corresponds to Fanno flow. The open symbols in figure 12 are predictions using the solution as described above, but with multiple sub-elements. Both the subsonic and supersonic cases show a reduced non-dimensional length prior to choking with heat transfer into the fluid. Further, increasing heat transfer reduces the non-dimensional length at which choking occurs. This is consistent with the direction of increasing entropy.

Figure 12. Comparison of the multiple sub-element solution using (4.1), shown as symbols, and Ferrari's solution, shown using lines for various $\varGamma$ values, for the case of simultaneous friction and heat transfer (Ferrari Reference Ferrari2021a). (a) Subsonic flow with heat transfer and friction. (b) Supersonic flow with heat transfer and friction.

The prediction using multiple sub-elements is consistent both qualitatively and quantitatively with Ferrari's exact solution for both subsonic and supersonic flows. The use of multiple sub-elements to reproduce Ferrari's solution is required due to the friction force approximation. The subsonic and supersonic inlet cases both used $1000$ sub-elements. This shows that in the limit of many sub-elements, the result will converge to the appropriate solution for the case of simultaneous friction and heat transfer. This result is consistent with previous statements that the solution of the differential form is recovered in the limit of infinitely many control volumes, so that any one volume is of infinitesimal size.

The combination of simultaneous linear diameter change and friction provides another useful example of the solution when using multiple sub-elements. A schematic of this problem is shown in figure 13. This problem was recently reduced by Ferrari (Reference Ferrari2021b) to a single, implicit, algebraic equation, which can then be solved iteratively. Obtaining this algebraic equation required restricting the diameter to vary linearly over the length of the nozzle or diffuser. This restriction means that it is valid only for conical geometries. Ferrari's solution provides a benchmark with which to compare the solution presented here using multiple sub-elements.

Figure 13. Schematic of simultaneous area change, with linear variation in diameter, and friction.

Two separate cases are examined: one with a subsonic inlet, and the other with a supersonic inlet. This comparison is shown in figure 14. The subsonic inlet case for a linearly convergent diameter subsonic nozzle is shown in figure 14(a); the supersonic inlet case for a linearly divergent diameter supersonic nozzle is shown in figure 14(b). The abscissa is a non-dimensional length, $c_f L/D(x)$. The ordinate is Mach number. Again, the overbar notation can be neglected since the flow is assumed uniform and the flux correction factors are unity. Note that for the subsonic inlet case, the diameter is decreasing along the length of the subsonic nozzle, so the flow direction is in the direction of increasing $c_f L/D(x)$. However, for the supersonic inlet case, the diameter is increasing along the length of the supersonic nozzle, so the flow direction is in the direction of decreasing $c_f L/D(x)$. The subsonic inlet case shows an increasing Mach number in the flow direction, which is consistent with the decreasing area. The subsonic nozzle stops just short of a Mach number of unity. The supersonic case shows an increasing Mach number along the flow direction, which is consistent with an increasing area for supersonic flow.

Figure 14. Comparison of the multiple sub-element solution using (4.1), shown as symbols, and Ferrari's solution, shown using a dashed line, for the case of simultaneous area change and friction (Ferrari Reference Ferrari2021b). (a) Subsonic flow with area change and friction. (b) Supersonic flow with area change and friction.

The prediction using multiple sub-elements is consistent both qualitatively and quantitatively with Ferrari's exact solution for both the subsonic and supersonic inlet cases. The use of multiple sub-elements to reproduce Ferrari's solution is required due to both the friction force approximation and the unknown wall pressure along the duct length. The subsonic inlet case used $900$ sub-elements, while the supersonic inlet case used $1200$ sub-elements. These sub-element counts are identical to the number of elements used by Ferrari when comparing the exact solution to numerical integration results. This shows that in the limit of many sub-elements, the result will converge to the appropriate solution in the case of simultaneous area change and friction. This is true regardless of the value of $\xi$ as long as $\nu =0$, though the rate of convergence varies with $\xi$. This result is consistent with the previous result. More specifically, the solution of the differential form is recovered in the limit of infinitely many control volumes, so that any one volume is of infinitesimal size.