2. Differential equations of the gas dynamics problem

Generalised Euler partial differential equations (PDEs) for a quasi-one-dimensional unsteady viscous diabatic compressible flow through nozzles and diffusers are given by (gravity is neglected) (Bermudez, Lopez & Vasquez-Cendon Reference Bermudez, Lopez and Vasquez-Cendon2017)

(2.1)\begin{equation}\left. {\begin{array}{c@{}} {\dfrac{{\partial (\rho A)}}{{\partial t}} + \dfrac{{\partial (\rho Au)}}{{\partial x}} = 0}\\ {\dfrac{{\partial (\rho uA)}}{{\partial t}} + \dfrac{{\partial (pA + \rho {\kern 1pt} {u^2}A)}}{{\partial x}} ={-} \mathrm{\pi }D{\tau_w} + p\dfrac{{\partial A}}{{\partial x}}}\\ {\dfrac{{\partial (\rho {e^0}A)}}{{\partial t}} + \dfrac{{\partial (A\rho u{h^0})}}{{\partial x}} = \mathrm{\pi }D{{\dot{q}}_f}} \end{array}} \right\},\end{equation}

where t is the time, x is the space coordinate, ρ, u and p are the cross-sectionally averaged 1-D density, velocity and pressure, respectively, A = A(x) is the variable pipe cross-section, D = D(x) is the variable diameter, e ⁰ is the stagnation internal energy (e ⁰ = e + u ²/2, where e is the internal energy), h ⁰ is the stagnation enthalpy (h ⁰ = h + u ²/2, where h is the enthalpy), τ_w is the wall friction shear stress and ${\dot{q}_f}$ is the heat flux with the walls. Equation (2.1) differs from the standard set of Euler equations for an isentropic flow because the friction and heat exchange are added on the right-hand side and evaluated as volume source terms (Hirsch Reference Hirsch2007; Maeda & Colonius Reference Maeda and Colonius2017) using Darcy–Weisbach's wall friction and convective heat transfer models. Equation (2.1) is a simplified approach compared to the Navier--Stokes equations, which use a multidimensional diffusive transport mechanism to model viscous stress and conductive heat. Wall friction stress can be modelled in (2.1) with (White Reference White2015)

(2.2)\begin{equation}{\tau _w} = \frac{f}{2}\rho {u^2},\end{equation}

where the friction coefficient f can be expressed as a function of the Reynolds number Re = ρuD/μ for a laminar flow, with μ being the dynamical viscosity of the fluid, and as a function of both Re and the relative roughness ε/D for a turbulent flow, with ε being the average roughness of the wall. Either Moody's diagram or both the 16/Re formula for a laminar flow (Re < 2000) and Colebrook–White's implicit equation for a turbulent flow (4000 < Re < 10⁸, 10⁻⁶ < ε/D < 10⁻²) can be used to evaluate f (Douglas et al. Reference Douglas, Gasiorek, Swaffield and Jack2005; White Reference White2015). An updated review of the most refined explicit formulas used in the literature to calculate f, as a function of Re and ε/D, is reported in Cheng (Reference Cheng2008). More accurate correlations, in which f generally depends on Re, ε/D and the Mach number are reported in Schwartz (Reference Schwartz1987), Asako et al. (Reference Asako, Pi, Turner and Faghri2003), Hong et al. (Reference Hong, Nakamura, Asako and Ueno2016).

According to a convective heat transfer model, the heat flux ${\dot{q}_f}$ in (2.1) is a function of the cross-sectionally averaged one-dimensional temperature (T) of the fluid, through a convective heat transfer coefficient that depends on Re and the Prandtl number (Bejan Reference Bejan2013). In other words, although a viscous friction and heat transfer simulation should require at least a 2-D model, to take into account the axial and radial variations of the velocity and temperature, the simplified model presented here uses a 1-D approach and then compensates for this shortcoming by applying empirical factors (f and a convective heat transfer coefficient) that need complex correlations.

Equation (2.1) is completed with the following state equations:

(2.3)\begin{equation}\left. {\begin{array}{c@{}} {p = \dfrac{{\gamma - 1}}{\gamma }\rho h,\quad h = {c_p}T}\\ {h = e + \dfrac{p}{\rho },\quad e = {c_v}T} \end{array}} \right\},\end{equation}

where γ = c_p/c_v is the ratio of the constant pressure specific heat (c_p) to the constant volume specific heat (c_v) and is here considered constant, as for a standard ideal gas.

When considering steady-state flows, the partial derivatives, with respect to time, vanish in (2.1), which reduces to a system of nonlinear ODEs with respect to x (Emmons Reference Emmons1958)

(2.4)\begin{equation}\left. {\begin{array}{c@{}} {\dfrac{{\textrm{d}\dot{m}}}{{\textrm{d}x}} = 0}\\ {A\dfrac{{\textrm{d}p}}{{\textrm{d}x}} + \dot{m}\dfrac{{\textrm{d}u}}{{\textrm{d}x}} ={-} \mathrm{\pi }D{\tau_w}}\\ {\dot{m}\dfrac{{\textrm{d}{h^0}}}{{\textrm{d}x}} = \mathrm{\pi }D{{\dot{q}}_f}} \end{array}} \right\},\end{equation}

where $\dot{m} = \rho Au$ is the constant mass flow rate through the pipe. These equations are applied to evaluate real flow evolutions in ducts, nozzles and diffusers with regular and symmetrical cross-section shapes and negligible curvature of the axis. The most accurate results for the simulation of real viscous diabatic (or adiabatic) compressible flows are obtained when $1/A|\textrm{d}A/\textrm{d}x|\mathrm{\Delta }x \ll 1$, where Δx is a characteristic distance, measured along the axis and related to the diameter size (Kamm & Shapiro Reference Kamm and Shapiro1979). When the radius of curvature is not large, compared with D or when $|\textrm{d}A/\textrm{d}x{|_{av}}\mathrm{\Delta }x \approx A$, a 2-D model is recommended for compressible flows.

3. Exact solution for compressible viscous flows in conical nozzles

Let us consider a steady-state compressible viscous adiabatic flow along a variable circular cross-section pipe with negligible curvature of the axis. The total energy conservation law in (2.4) reduces to $({\dot{q}_f} = 0)$

(3.1)\begin{equation}{c_p}T + {\textstyle{1 \over 2}}{u^2} = {h^0} = cost.\end{equation}

The momentum balance law in (2.4) can be rewritten as follows (the momentum balance is multiplied by u and the expression of τ_w in (2.2) is used):

(3.2)\begin{equation}Au\frac{{\textrm{d}p}}{{\textrm{d}x}} + \dot{m}\frac{\textrm{d}}{{\textrm{d}x}}\left( {\frac{{{u^2}}}{2}} \right) + \frac{{4f}}{D}\dot{m}\left( {\frac{{{u^2}}}{2}} \right) = 0,\end{equation}

where f is here considered uniform (it can be made equal to a space-averaged value over the pipe), in line with the original Fanno flow and with several approximated solutions of compressible flows with wall friction (Parker Reference Parker1989; Shiraishi et al. Reference Shiraishi, Murayama, Kawakami and Nakano2011; Cavazzuti, Corticelli & Karayiannis Reference Cavazzuti, Corticelli and Karayiannis2019).

The pressure can be expressed with the following formula, which was obtained using the state equations of ideal gases, i.e. (2.3), in conjunction with (3.1) and the flow-rate equation $\rho uA = \dot{m} = cost$:

(3.3)\begin{equation}p = p(u,A) = \frac{{\gamma - 1}}{\gamma }\frac{{\dot{m}}}{{Au}}h = \frac{{\gamma - 1}}{\gamma }\frac{{\dot{m}}}{{Au}}\left( {{h^0} - \frac{1}{2}{u^2}} \right).\end{equation}

The pressure gradient can thus be calculated by applying the chain rule. By using the relation ${h^0} = {c_p}{T^0} = \gamma /(\gamma - 1)R{T^0}$, one achieves

(3.4)\begin{equation}\frac{{\textrm{d}p}}{{\textrm{d}x}} ={-} \frac{{\dot{m}}}{{{A^2}u}}R{T^0}\frac{{\textrm{d}A}}{{\textrm{d}x}} + \frac{{\dot{m}}}{{{A^2}}}\frac{{\gamma - 1}}{{2\gamma }}u\frac{{\textrm{d}A}}{{\textrm{d}x}} - \frac{{\dot{m}}}{{{Au^2}}}R{T^0}\frac{{\textrm{d}u}}{{\textrm{d}x}} - \frac{{\dot{m}}}{{2A}}\frac{{\gamma - 1}}{\gamma }\frac{{\textrm{d}u}}{{\textrm{d}x}}.\end{equation}

By substituting (3.4) in (3.2), taking into account the relation

(3.5)\begin{equation}\frac{1}{A}\frac{{\textrm{d}A}}{{\textrm{d}x}} = \frac{2}{D}\frac{{\textrm{d}D}}{{\textrm{d}x}},\end{equation}

and performing some algebraic manipulations, the following nonlinear ODE is obtained with a single unknown function, that is, the kinetic energy per unit of mass ($y: = {u^2}/2$ for conciseness of notation):

(3.6)\begin{equation}\frac{{\gamma + 1}}{\gamma }y\frac{{\textrm{d}y}}{{\textrm{d}x}} - R{T^0}\frac{{\textrm{d}y}}{{\textrm{d}x}} = \frac{{4R{T^0}}}{D}\frac{{\textrm{d}D}}{{\textrm{d}x}}y - \left[ {\frac{{4(\gamma - 1)}}{{\gamma D}}\frac{{\textrm{d}D}}{{\textrm{d}x}} + \frac{{8f}}{D}} \right]{y^2}.\end{equation}

Equation (3.6) can be identified with the ${g_1}y\,\textrm{d}y/\textrm{d}x + {g_0}\,\textrm{d}y/\textrm{d}x = {j_0} + {j_1}y + {j_2}{y^2}$ ODE model, which is the standard form of an Abel equation of the second kind. In the case of (3.6), ${g_0} = -R{T^0}$ and ${g_1} = (\gamma + 1)/\gamma $ are constant coefficients, that is, they are independent of x, whereas coefficients ${j_1} = [4R{T^0}/D \cdot \textrm{d}D/\textrm{d}x]$ and ${j_2} ={-} [4(\gamma - 1)/(\gamma D) \cdot \textrm{d}D/\textrm{d}x + 8f/D]$ are variable with x (the free term j ₀ is null).

The variable is changed according to the following transformation (Bogouffa Reference Bogouffa2014):

(3.7a,b)\begin{align}\frac{1}{z}: = {g_1}y + {g_0} = \frac{{\gamma + 1}}{\gamma }y - R{T^0} \Rightarrow y = \frac{\gamma }{{\gamma + 1}}\frac{{1 + R{T^0}z}}{z},\quad \frac{{\textrm{d}y}}{{\textrm{d}x}} ={-} \gamma \frac{{\textrm{d}z/\textrm{d}x}}{{(\gamma + 1){z^2}}}.\end{align}

Equation (3.6) can be reduced, after some analytical steps, to the following Abel equation of the first kind in homogeneous form, that is, without the free term

(3.8)\begin{equation}\frac{{\textrm{d}z}}{{\textrm{d}x}} = {r_0}(x) + {r_1}(x)z + {r_2}(x){z^2} + {r_3}(x){z^3},\end{equation}

where

(3.9)\begin{align}\left. {\begin{array}{c@{}} {{r_1}(x) = \dfrac{4}{{D(x)}}\left[ {\dfrac{{\gamma - 1}}{{\gamma + 1}}\dfrac{{\textrm{d}D}}{{\textrm{d}x}} + \dfrac{{2\gamma f}}{{\gamma + 1}}} \right],\quad {r_2}(x) = \dfrac{{4R{T^o}}}{{D(x)}}\left[ {\dfrac{{\gamma - 3}}{{\gamma + 1}}\dfrac{{\textrm{d}D}}{{\textrm{d}x}} + \dfrac{{4\gamma f}}{{\gamma + 1}}} \right]}\\ {{r_3}(x) = \dfrac{{8{{(R{T^o})}^2}}}{{D(x)}}\left[ {\dfrac{{\gamma f}}{{\gamma + 1}} - \dfrac{1}{{\gamma + 1}}\dfrac{{\textrm{d}D}}{{\textrm{d}x}}} \right],\quad {r_0} = 0} \end{array}} \right\}.\end{align}

There is no general solution to Abel equations of the first or second kind, in terms of elementary functions. Many investigations have attempted to propose exact solutions for Abel equations under certain conditions. In general, the method to derive exact solutions involves applying appropriate transformations to reduce both kinds of Abel ODEs to equivalent canonical forms that can accept closed-form solutions, such as separable, linear, Bernoulli and Riccati ODEs, or that can admit implicit or parametrical solutions, provided their coefficients satisfy specific functional relations (Mak & Harko Reference Mak and Harko2002; Panayotounakos Reference Panayotounakos2005; Bogouffa Reference Bogouffa2010, Reference Bogouffa2014). A comprehensive collection of the cases of closed-form explicit, implicit and parametrical solutions of both kinds of Abel equations can be found in Polyanin & Zaitsev (Reference Polyanin and Zaitsev2003).

3.1. Exact solution of the homogeneous Abel equation of the first kind

Since the Abel equation of the first kind is homogeneous for the considered gas dynamics problem, the most appropriate transformation is the following one, as provided by Kamke (Reference Kamke1983):

(3.10a–c)\begin{equation}z(x) = t(x)q[w(x)],\quad t(x) = exp \left[ {\int {{r_1}(x)\,\textrm{d}x} } \right],\quad w(x) = \int {{r_2}(x)t(x)\,\textrm{d}x} .\end{equation}

Equation (3.8) therefore becomes (Kamke Reference Kamke1983)

(3.11)\begin{equation}\left. {\begin{array}{c@{}} {\dfrac{{\textrm{d}q}}{{\textrm{d}w}} = \dfrac{{{r_3}(x){\kern 1pt} {\kern 1pt} t(x)}}{{{r_2}(x)}}{q^3} + {q^2}}\\ {x = x(w)} \end{array}} \right\},\end{equation}

which is a canonical form of a parametric representation of Abel's equation, where w is the parameter. A sufficient condition for the derivation of a closed-form solution of (3.11), and thus of a homogenous Abel equation of the first kind, i.e. (3.8), is obtained from the following relation between the coefficients (Markakis Reference Markakis2009):

(3.12)\begin{equation}\frac{\textrm{d}}{{\textrm{d}x}}\left[ {\frac{{{r_3}(x){\kern 1pt} }}{{{r_2}(x)}}} \right] = \frac{c}{{{{(1 + c)}^2}}}{r_2}(x) - \frac{{{r_3}(x){\kern 1pt} {r_1}(x)}}{{{r_2}(x)}},\end{equation}

where c is an arbitrary constant.

If (3.12) holds and ${r_2}(x) \ne 0$ (which implies that $\textrm{d}D/\textrm{d}x \ne 4\gamma f/(3 - \gamma )$), the homogeneous Abel equation of the first kind admits an implicit solution of the following form (Markakis Reference Markakis2009):

(3.13)\begin{equation}\frac{{|(1 + c)v(x)z + 1{|^c}}}{{|(1 + c)v(x)z + c|}} - K|(1 + c)v(x)z{|^{c - 1}}exp \left[ {\frac{{c(1 - c)}}{{{{(1 + c)}^2}}}\int {\frac{{{r_2}(x)}}{{v(x)}}\,\textrm{d}x} } \right] = 0,\end{equation}

where $v(x) = {r_3}(x)/{r_2}(x)$ and K stands for the parameter of the family of solutions of the first-order ODEs.

If the expressions of r ₁, r ₂ and r ₃ in (3.9) are substituted in (3.12), the following relation is obtained:

(3.14)\begin{align}& \dfrac{1}{2}{(\gamma + 1)^2}\gamma fD\dfrac{{{\textrm{d}^2}D}}{{\textrm{d}{x^2}}} + \dfrac{c}{{{{(1 + c)}^2}}}{\left[ {(\gamma - 3)\dfrac{{\textrm{d}D}}{{\textrm{d}x}} + 4\gamma f} \right]^3}\nonumber\\ & \quad = 2\left[ {(\gamma - 1)\dfrac{{\textrm{d}D}}{{\textrm{d}x}} + 2\gamma f} \right]\left[ {(\gamma - 3)\dfrac{{\textrm{d}D}}{{\textrm{d}x}} + 4\gamma f} \right]\left[ {\gamma f - \dfrac{{\textrm{d}D}}{{\textrm{d}x}}} \right]. \end{align}

The objective is now to find those functions of D(x) that allow (3.14) to be satisfied in order to obtain an analytical solution from (3.13) for the viscous flow through a variable cross-section pipe. A linear function, D(x), can satisfy (3.14): in this case, dD/dx is constant and d²D/dx ² is null. If such a linear dependence of D on x is hypothesised, (3.14) reduces to

(3.15)\begin{align}&\frac{c}{{{{(1 + c)}^2}}}{\left[ {(\gamma - 3)\frac{{\textrm{d}D}}{{\textrm{d}x}} + 4\gamma f} \right]^3}\nonumber\\ &\quad = 2\left[ {(\gamma - 1)\frac{{\textrm{d}D}}{{\textrm{d}x}} + 2\gamma f} \right]\left[ {(\gamma - 3)\frac{{\textrm{d}D}}{{\textrm{d}x}} + 4\gamma f} \right]\left[ {\gamma f - \frac{{\textrm{d}D}}{{\textrm{d}x}}} \right].\end{align}

Equation (3.15) can be simplified to the following second-order algebraic equation, with respect to dD/dx;

(3.16)\begin{align}& \left[ {\dfrac{c}{{{{(1 + c)}^2}}}{{(\gamma - 3)}^2} + 2(\gamma - 1)} \right]{\left( {\dfrac{{\textrm{d}D}}{{\textrm{d}x}}} \right)^2} + 2\gamma f(3 - \gamma ){\left( {\dfrac{{c - 1}}{{c + 1}}} \right)^2}\dfrac{{\textrm{d}D}}{{\textrm{d}x}}\nonumber\\ & \quad - 4{\left( {\dfrac{{c - 1}}{{c + 1}}} \right)^2}{\gamma ^2}{f^2} = 0.\end{align}

The solutions to (3.16), and thus to (3.14), expressed as functions of c, are provided by

(3.17a,b)\begin{equation}{\left. {\frac{{\textrm{d}D}}{{\textrm{d}x}}} \right|_1} = \frac{{2\gamma f(c - 1)}}{{2c + \gamma - 1}}\quad {\left. {\frac{{\textrm{d}D}}{{\textrm{d}x}}} \right|_2} = \frac{{2\gamma f(c - 1)}}{{c - \gamma c - 2}}.\end{equation}

Equation (3.16) always admits two real solutions, whatever the value of c, except for $c ={-} (\gamma - 1)/2$ and $c ={-} 2/(\gamma - 1)$. The expressions in (3.17) correspond to hyperboles, and the dD/dx slope can therefore take on any finite value (−∞ < dD/dx < +∞), depending on c, which is an arbitrary dimensionless constant in the procedure; however, c ≠ −1, otherwise r ₂(x) = 0, and (3.13) is no longer in force. Figure 1 qualitatively plots (dD/dx)₁ as a function of c, according to (3.17a): the corresponding graph for (dD/dx)₂ is similar and has therefore been omitted for conciseness. Hence, once the value of dD/dx has been selected for the considered conical pipe, c can be determined, by making use of either (3.17) or a graph like the one in figure 1. The corresponding exact solution for the velocity of the 1-D flow through the conical nozzle or diffusor (0≤x ≤ L) follows from (3.13), which is hereafter rewritten and customised for the specific problem, taking into account (3.7a,b) and (3.9) (u ² must be different from $2\gamma /[(\gamma + 1)R{T^0}]$, which corresponds to a choking condition, and dD/dx ≠ 0)

(3.18)\begin{align}\left.\begin{array}{c@{}}\dfrac{{|2(1 + c){\kern 1pt} \gamma v + (\gamma + 1){u^2} - 2\gamma R{T^o}{|^c}}}{{|2(1 + c){\kern 1pt} \gamma v + c[(\gamma + 1){u^2} - 2\gamma R{T^o}]|}} - K^{\prime\prime}{[D(x)]^a} = 0,\\ \hspace{-10pc} \textrm{where}\;D(x) = {D_1} + \dfrac{{\textrm{d}D}}{{\textrm{d}x}}x\\ v = 2R{T^o}\dfrac{{\gamma f - \dfrac{{\textrm{d}D}}{{\textrm{d}x}}}}{{4\gamma f + (\gamma - 3)\dfrac{{\textrm{d}D}}{{\textrm{d}x}}}},\quad a = \dfrac{{2c(1 - c)}}{{{{(1 + c)}^2}}}\dfrac{{{{\left[ {4\gamma f + (\gamma - 3)\dfrac{{\textrm{d}D}}{{\textrm{d}x}}} \right]}^2}}}{{({\gamma + 1} )\left[ {\gamma f\left( {\dfrac{{\textrm{d}D}}{{\textrm{d}x}}} \right) - {{\left( {\dfrac{{\textrm{d}D}}{{\textrm{d}x}}} \right)}^2}} \right]}},\\ \hspace{-17.5pc} K^{\prime\prime} = K|(1 + c)\gamma v{|^{c - 1}}\end{array}\right\}.\end{align}

In the considered problem, v and a are constant parameters, that is, they are independent of x, and D ₁ is the value of the diameter at the inlet of the conical pipe (x ₁ = 0).

Figure 1. Slope (dD/dx)₁ as a function of c (graph not to scale).

The provided boundary conditions for the steady-state problem allow parameter K′′ to be determined (K′′ is a constant because $|(1 + c){\kern 1pt} \gamma v{\kern 1pt} {|^{c - 1}}$ is here independent of x). In general, 3 boundary conditions are required to solve a steady-state 1-D gas dynamics problem: these three input data can be assigned at the x ₁ = 0 and x ₂ = L boundaries without any constraints (Urata Reference Urata2013).

Let us suppose that the stagnation temperature T ⁰, mass flow rate $\dot{m}$ and velocity u ₁ at x ₁ = 0 are assigned; although this choice is not mandatory, it simplifies the calculations. Parameter K′′ can in fact be evaluated from (3.18) as follows:

(3.19)\begin{equation}K^{\prime\prime} = \frac{1}{{D_1^a}}\frac{{|2(1 + c)\gamma v + (\gamma + 1)u_1^2 - 2\gamma R{T^o}{|^c}}}{{|2(1 + c){\kern 1pt} \gamma v + c[(\gamma + 1)u_1^2 - 2\gamma R{T^o}]|}}.\end{equation}

Once the distribution of u ²/2, with respect to x, has been fully determined by means of (3.18) and (3.19), the corresponding T(x) and ρ(x) can be determined from (3.1) and from the constant mass flow rate, respectively, while p(x) can be calculated with (3.3).

When dD/dx = 0, it results, from (3.17), that c = 1. Hence, (3.18) is always verified for K′′ = 1 (an identity is obtained) and the u(x) distribution cannot be calculated. Since the pipe features a constant diameter for c = 1, Fanno's model can therefore be used to obtain the exact solution for the viscous compressible flow.

3.2. Dimensionless representation of the exact solution

Equations (3.18) and (3.19) and the geometrical domain include 8 independent dimensional sizes, that is, x, T ⁰, D, D ₁, L, u, u ₁, R, and 2 independent dimensionless parameters, that is, γ and f . Since both mechanical and thermal quantities are involved in the problem, the rank of the dimensional problem is equal to 4, i.e. m, s, kg and K are required as measurement units to express all the quantities that appear in the problem. According to Buckingham's theorem (Yarin Reference Yarin2012), it is possible to express the solution and boundary conditions in terms of (8 + 2) − 4 = 6 dimensionless groups or parameters. The dimensionless representation of the solution can be useful to identify the physical factors that drive the analysed phenomenon. If the expression of $K^{\prime\prime}$ given by (3.19) is substituted in (3.18), and the thus resulting equation is divided by ${(\gamma R{T^0})^{c - 1}}$, one obtains

(3.20) \begin{equation}\left.\begin{gathered} \hspace{-4pc}\dfrac{{{{\left|{{\kern 1pt} (1 + c){\kern 1pt} {v^\ast } + \dfrac{{(\gamma + 1)}}{2}\dfrac{{{u^2}}}{{\gamma R{T^o}}} - 1} \right|}^c}}}{{\left|{(1 + c){v^\ast } + c\left[ {\dfrac{{(\gamma + 1)}}{2}\dfrac{{{u^2}}}{{\gamma R{T^o}}} - 1} \right]} \right|}}\\ \hspace{4pc} - \dfrac{{{{\left|{(1 + c){v^\ast } + \dfrac{{(\gamma + 1)}}{2}\dfrac{{u_1^2}}{{\gamma R{T^o}}} - 1} \right|}^c}}}{{\left|{(1 + c){v^\ast } + c\left[ {\dfrac{{(\gamma + 1)}}{2}\dfrac{{u_1^2}}{{\gamma R{T^o}}} - 1} \right]{\kern 1pt} {\kern 1pt} } \right|}}{\left[ {\dfrac{{D(x)}}{{{D_1}}}} \right]^{^a}} = 0\\ {{v^\ast } = 2\frac{{\gamma f - \frac{{\textrm{d}D}}{{\textrm{d}x}}}}{{4\gamma f + (\gamma - 3)\frac{{\textrm{d}D}}{{\textrm{d}x}}}},\quad a = \frac{{2c(1 - c)}}{{{{(1 + c)}^2}}}\frac{{{{\left[ {4\gamma f + (\gamma - 3)\frac{{\textrm{d}D}}{{\textrm{d}x}}} \right]}^2}}}{{({\gamma + 1} )\left[ {\gamma f\left( {\frac{{\textrm{d}D}}{{\textrm{d}x}}} \right) - {{\left( {\frac{{\textrm{d}D}}{{\textrm{d}x}}} \right)}^2}} \right]}}} \end{gathered}\right\}.\end{equation}

Let us now divide both the numerator and denominator of v* by f, both the numerator and denominator of a by f ² and replace T ⁰ by either $T + {u^2}/(2{c_p})$ or ${T_1} + u_1^2/(2{c_p})$, according to (3.1). If all the derivatives of D, with respect to x, are expressed as $\textrm{d}D/\textrm{d}x = {L^{ - 1}}\textrm{d}D/\textrm{d}{x^\ast }$, where ${x^\ast }: = x/L$ is a dimensionless abscissa, the following equation can be obtained from (3.20) and (3.17) after some algebraic manipulation:

(3.21) \begin{equation}\left.\begin{gathered} \hspace{-4pc}\dfrac{{{{\left|{(1 + c){v^\ast } + \dfrac{{(\gamma + 1)}}{2}\dfrac{{{u^2}}}{{\gamma RT + \dfrac{{\gamma - 1}}{2}{u^2}}} - 1} \right|}^c}}}{{\left|{{\kern 1pt} (1 + c){v^\ast } + c\left[ {\dfrac{{(\gamma + 1)}}{2}\dfrac{{{u^2}}}{{\gamma RT + \dfrac{{\gamma - 1}}{2}{u^2}}} - 1} \right]{\kern 1pt} } \right|}}\\ \hspace{4pc} - \dfrac{{{{\left|{{\kern 1pt} (1 + c){v^\ast } + \dfrac{{(\gamma + 1)}}{2}\dfrac{{u_1^2}}{{\gamma R{T_1} + \dfrac{{\gamma - 1}}{2}u_1^2}} - 1} \right|}^c}}}{{\left|{(1 + c){v^\ast } + c\left[ {\dfrac{{(\gamma + 1)}}{2}\dfrac{{u_1^2}}{{\gamma R{T_1} + \dfrac{{\gamma - 1}}{2}u_1^2{\kern 1pt} }} - 1} \right]{\kern 1pt} {\kern 1pt} } \right|}}{\left[ {\dfrac{{D(x)}}{{{D_1}}}} \right]^a} = 0\\ {v^\ast } = 2\dfrac{{\gamma - \dfrac{\textrm{d}}{{\textrm{d}{x^\ast }}}\left( {\dfrac{D}{{f\,L}}} \right)}}{{4\gamma + (\gamma - 3)\dfrac{\textrm{d}}{{\textrm{d}{x^\ast }}}\left( {\dfrac{D}{{fL}}} \right)}},\\ a = \dfrac{{2c(1 - c)}}{{{{(1 + c)}^2}}}\dfrac{{{{\left[ {4\gamma + (\gamma - 3)\dfrac{\textrm{d}}{{\textrm{d}{x^\ast }}}\left( {\dfrac{D}{{fL}}} \right)} \right]}^2}}}{{(\gamma + 1)\left\{ {\gamma \dfrac{\textrm{d}}{{\textrm{d}{x^\ast }}}\left( {\dfrac{D}{{fL}}} \right) - {{\left[ {\dfrac{\textrm{d}}{{\textrm{d}{x^\ast }}}\left( {\dfrac{D}{{fL}}} \right)} \right]}^2}} \right\}}},\\ c = \dfrac{{(\gamma - 1)\dfrac{\textrm{d}}{{\textrm{d}{x^\ast }}}\left( {\dfrac{D}{{fL}}} \right) + 2}}{{2\left[ {1 - \dfrac{\textrm{d}}{{\textrm{d}{x^\ast }}}\left( {\dfrac{D}{{fL}}} \right)} \right]}}\end{gathered}\right\}.\end{equation}

By introducing the Mach numbers $Ma = u/\sqrt {\gamma RT} $ and $M{a_1} = {u_1}/\sqrt {\gamma R{T_1}} $, (3.21) can be presented in the following dimensionless form (0 ≤ x* ≤ 1):

(3.22) \begin{gather}\left.\begin{gathered} \hspace{-5pc}\dfrac{{{{\left|{(1 + c){v^\ast } + \dfrac{{(\gamma + 1)}}{2}\dfrac{{M{a^2}}}{{1 + \dfrac{{\gamma - 1}}{2}M{a^2}}} - 1{\kern 1pt} {\kern 1pt} } \right|}^c}}}{{\left|{(1 + c){v^\ast } + c\left[ {\dfrac{{(\gamma + 1)}}{2}\dfrac{{M{a^2}}}{{1 + \dfrac{{\gamma - 1}}{2}M{a^2}}} - 1} \right]{\kern 1pt} } \right|}}\\ \hspace{5pc} - \dfrac{{{{\left|{(1 + c){v^\ast } + \dfrac{{(\gamma + 1)}}{2}\dfrac{{Ma_1^2}}{{1 + \dfrac{{\gamma - 1}}{2}Ma_1^2}} - 1{\kern 1pt} } \right|}^c}}}{{\left|{(1 + c){v^\ast } + c\left[ {\dfrac{{(\gamma + 1)}}{2}\dfrac{{Ma_1^2}}{{1 + \dfrac{{\gamma - 1}}{2}Ma_1^2}} - 1} \right]{\kern 1pt} } \right|}}{\left\{ {\dfrac{{\left( {\dfrac{{fL}}{{{D_1}}}} \right)}}{{\left[ {\dfrac{{fL}}{{D({x^\ast })}}} \right]}}} \right\}^a} = 0\\ {v^\ast } = 2\dfrac{{\gamma - \dfrac{\textrm{d}}{{\textrm{d}{x^\ast }}}\left( {\dfrac{D}{{fL}}} \right)}}{{4\gamma + (\gamma - 3)\dfrac{\textrm{d}}{{\textrm{d}{x^\ast }}}\left( {\dfrac{D}{{fL}}} \right)}},\\ a = \dfrac{{2c(1 - c)}}{{{{(1 + c)}^2}}}\dfrac{{{{\left[ {4\gamma + (\gamma - 3)\dfrac{\textrm{d}}{{\textrm{d}{x^\ast }}}\left( {\dfrac{D}{{fL}}} \right)} \right]}^2}}}{{({\gamma + 1} )\left\{ {\gamma \dfrac{\textrm{d}}{{\textrm{d}{x^\ast }}}\left( {\dfrac{D}{{fL}}} \right) - {{\left[ {\dfrac{\textrm{d}}{{\textrm{d}{x^\ast }}}\left( {\dfrac{D}{{fL}}} \right)} \right]}^2}} \right\}}},\\ c = \dfrac{{(\gamma - 1)\dfrac{\textrm{d}}{{\textrm{d}{x^\ast }}}\left( {\dfrac{D}{{fL}}} \right) + 2}}{{2\left[ {1 - \dfrac{\textrm{d}}{{\textrm{d}{x^\ast }}}\left( {\dfrac{D}{{fL}}} \right)} \right]}}\end{gathered} \right\}.\end{gather}

Six dimensionless factors appear in (3.22): fL/D ₁, fL/D, Ma, γ, Ma ₁ and x*. The Ma group, as expected, plays a fundamental role in the examined phenomenon. The steady-state evolution of the flow in the pipe is driven by both wall friction and a variable flow area, and the Mach number is an essential factor in both Fanno flow and nozzle flow problems. Moreover, the f L/D quantity is assembled and introduced into the dimensionless equation, since it is also an important group for Fanno flows (Shapiro Reference Shapiro1953).

3.3. Exact solution for the $\textrm{d}D/\textrm{d}x = 4\gamma f/(3 - \gamma )$ case

As already mentioned, when $\textrm{d}D/\textrm{d}x = 4\gamma f/(3 - \gamma )$, r ₂ = 0 and (3.18) cannot be applied. Therefore, an exact solution for this particular positive value of the diameter slope (for ideal gases γ < 3) should be derived to complete the model of the conical nozzle with wall friction. It is possible to prove that if the coefficients of an Abel equation of the second kind satisfy the subsequent relation, where ${g^{\prime}_0}$ and ${g^{\prime}_1}$ are the synthetic notations for the derivatives of g ₀ and g ₁, respectively,

(3.23)\begin{equation}{g_0}(x)[2{j_2}(x) + {g^{\prime}_1}(x)] = {g_1}(x)[{j_1}(x) + {g^{\prime}_0}(x)],\end{equation}

then the following Julia solution of Abel's ${g_1}y\,\textrm{d}y/\textrm{d}x + {g_0}\,\textrm{d}y/\textrm{d}x = {j_0} + {j_1}y + {j_2}{y^2}$ equation holds (Bogouffa Reference Bogouffa2014)

(3.24)\begin{equation}\frac{{{g_1}(x){y^2} + 2{g_0}(x)y}}{{{g_1}\varGamma }} = 2\int {\frac{{{j_0}(x)}}{{{g_1}\varGamma }}\,\textrm{d}x + K,}\end{equation}

where K is an integration constant and Γ is equal to

(3.25)\begin{equation}\varGamma = exp \left( {\int {\frac{{2{j_2}(x)}}{{{g_1}}}\,\textrm{d}x} } \right).\end{equation}

As a result of the identification of the various coefficients made in § 3, (3.23) leads to the following condition on the diameter slope:

(3.26)\begin{equation}\frac{{\textrm{d}D}}{{\textrm{d}x}} = \frac{{4\gamma f}}{{3 - \gamma }},\end{equation}

and, recalling that j ₀ = 0 in the considered gas dynamics problem, (3.24) and (3.25) reduce to

(3.27a,b)\begin{equation}{y^2} - \frac{{2\gamma }}{{\gamma + 1}}R{T^0}y = K{[D(x)]^b}\quad \textrm{with } b ={-} 4,\end{equation}

where the boundary data allow K to be determined.

Equations (3.26) and (3.27) can be rewritten in dimensionless form as follows (0 ≤ x* ≤ 1):

(3.28) \begin{align} &\dfrac{{M{a^4}}}{{{{\left( {1 + \dfrac{{\gamma - 1}}{2}M{a^2}} \right)}^2}}} - \dfrac{4}{{\gamma + 1}}\dfrac{{M{a^2}}}{{\left( {1 + \dfrac{{\gamma - 1}}{2}M{a^2}} \right)}} = \left[ {\dfrac{{Ma_1^4}}{{{{\left( {1 + \dfrac{{\gamma - 1}}{2}Ma_1^2} \right)}^2}}}}\right.\nonumber\\ & \quad \left. { - \dfrac{4}{{\gamma + 1}}\dfrac{{Ma_1^2}}{{\left( {1 + \dfrac{{\gamma - 1}}{2}Ma_1^2} \right)}}} \right]{\left\{ {\dfrac{{\left[ {\dfrac{{fL}}{{D({x^\ast })}}} \right]}}{{\left( {\dfrac{{fL}}{{{D_1}}}} \right)}}} \right\}^4},\quad \dfrac{\textrm{d}}{{\textrm{d}{x^\ast }}}\left( {\dfrac{D}{{fL}}} \right) = \dfrac{{4\gamma }}{{3 - \gamma }}.\end{align}

4. The fluid dynamics of compressible viscous adiabatic flows in conical pipes

Let us again refer to (2.4), which governs the general steady-state compressible quasi-one-dimensional flow along a pipe with a variable cross-section and wall friction, under the hypothesis ${\dot{q}_f} \approx 0$. The objective of the present analysis is to design suitable nozzle cross-section variations in order to obtain the desired evolution of the compressible viscous adiabatic flow or to predict the compressible viscous adiabatic flow transformation for an assigned geometry of the nozzle. In both cases, it is necessary to define the relationships between the variations along x of the 1-D flow properties and of the nozzle cross-section.

By adding (2.3), in the form p = ρRT, to (2.4) and by taking into account (3.1), it is possible to derive the following differential equation after some analytical steps:

(4.1)\begin{equation}(M{a^2} - 1)\frac{\textrm{d}}{{\textrm{d}x}}(\ln u) = \frac{\textrm{d}}{{\textrm{d}x}}(\ln A) - \frac{{2f}}{D}\gamma M{a^2}.\end{equation}

This equation is of fundamental importance because it relates the driving factors of the flow evolution, that is, the nozzle cross-section variation and the friction term, to the quasi-one-dimensional flow properties and their variations along the nozzle. When the friction is negligible (f ≈ 0), (4.1) reduces to Hugoniot's differential equation for isentropic nozzles, whereas if a straight duct is considered for which d(ln A)/dx = 0, (4.1) coincides with the Fanno flow relation that links the velocity and space coordinate by means of the Mach number (Shapiro Reference Shapiro1953).

Applying the logarithmic derivative rule to A(x) in (4.1) and using (3.5), one obtains

(4.2)\begin{equation}(M{a^2} - 1)\frac{\textrm{d}}{{\textrm{d}x}}(\ln u) = \frac{2}{D}\left( {\frac{{\textrm{d}D}}{{\textrm{d}x}} - f\gamma M{a^2}} \right).\end{equation}

Equation (4.2) generally guides the design of nozzles and diffusers:

• • in order to accelerate a steady-state viscous flow (du/dx > 0) along a convergent nozzle (dD/dx < 0), the flow must be subsonic (Ma < 1 is a necessary condition). In fact, if the viscous flow were supersonic, it could only be decelerated along a convergent pipe, which therefore would act as a diffuser;

• • a divergent nozzle (dD/dx > 0) is necessary to accelerate a steady-state supersonic viscous flow (du/dx > 0, Ma > 1), but it is also required that $\textrm{d}D/\textrm{d}x \gt f\gamma M{a^2}\,\forall x$;

• • a subsonic viscous flow (Ma < 1) can also be accelerated (du/dx > 0) along a divergent nozzle, provided that $0 \lt \textrm{d}D/\textrm{d}x \lt f\gamma M{a^2}\forall x$;

• • a critical flow (Ma = 1) can only be reached along a divergent pipe, at a section where $\textrm{d}D/\textrm{d}x = f\gamma$;

• • it is possible to decelerate supersonic and subsonic viscous flows by means of divergent diffusers that satisfy the $0 \lt \textrm{d}D/\textrm{d}x \lt f\gamma M{a^2}$ ∀x and $\textrm{d}D/\textrm{d}x \gt f\gamma M{a^2}\,\forall x$ conditions, respectively;

• • a local maximum point or a local minimum one in the velocity distribution (du/dx = 0) can only occur when the $\textrm{d}D/\textrm{d}x = f\gamma M{a^2}$ condition is satisfied at a certain x for a nozzle or a diffuser, respectively;

• • the flow is subsonic in correspondence to the throat ($\textrm{d}D/\textrm{d}x = 0$), and is accelerated (du/dx > 0) for a nozzle, whereas it is supersonic and decelerated for a diffuser.

Let us now limit the discussion to conical pipes for which dD/dx is a constant term, and let us suppose that a viscous flow moves through a conical nozzle of assigned geometry. The flow at the inlet of a conical convergent nozzle ($\textrm{d}D/\textrm{d}x = \textrm{const}. \lt 0$) must be subsonic (otherwise the convergent pipe would be a diffuser) and can only be accelerated up to $Ma \to {1^ - }$, provided the nozzle is sufficiently long. In fact, for $Ma \to 1$, $\textrm{d}u/\textrm{d}x \to + \infty$ in (4.2), and this confirms that a critical viscous flow cannot be reached in a conical convergent nozzle as a choked flow occurs.

If the viscous flow is supersonic at the inlet of the conical divergent nozzle $\textrm{(d}D/\textrm{d}x = \textrm{const}\textrm{.} \gt 0)$, we have $1 \lt M{a_1} \lt \sqrt {(\textrm{d}D/\textrm{d}x)/(f\gamma )}$ at its inlet (otherwise the divergent pipe would be a diffuser) and the viscous flow can then be accelerated up to a maximum Mach number according to the $M{a_{max}} = \sqrt {(\textrm{d}D/\textrm{d}x)/(f\gamma )}$ formula. As already mentioned, one obtains du/dx = 0 for this Mach number.

If the viscous flow at the inlet of the divergent pipe is subsonic, $\sqrt {(\textrm{d}D/\textrm{d}x)/(f\gamma )} \lt M{a_1} \lt 1$ is a necessary condition for the flow to be accelerated up to $Ma = 1$, where $\textrm{d}u/\textrm{d}x \to \infty$. This means that a divergent conical nozzle cannot undergo transition from a subsonic viscous flow to a supersonic one. The $\textrm{d}D/\textrm{d}x \lt f\gamma M{a^2}$ with Ma < 1, and $\textrm{d}D/\textrm{d}x \gt f\gamma M{a^2}$ relations with Ma > 1 should be satisfied simultaneously for the same dD/dx constant value, which is clearly inconsistent. A similar situation occurs for Fanno's flow, for which choking makes it impossible for a subsonic flow to undergo a steady-state transition to a supersonic flow along a pipe of constant diameter. Hence, a divergent nozzle, with variable dD/dx with respect to x, is necessary to change a viscous flow from subsonic to supersonic conditions, whereas neither dD/dx = 0 nor dD/dx = const. > 0 are sufficient for this purpose.

4.1. Analytical solution in the presence of shocks

Pressure p_E in the environment downstream of the conical pipe is an important parameter for fluid dynamics analysis. The pressure value at the final section of a conical pipe, calculated according to the exact solution developed in § 3, is referred to as p ₂.

When the flow is subsonic in the final section of a conical pipe, p ₂ must be equal to p_E. This constraint is satisfied by the proposed exact solution, via the boundary data, which must be consistent with the p_E value.

When the flow is supersonic at the exit of a conical nozzle (dD/dx > 0), p ₂ may be different from p_E, but the exact solution developed in § 3 may become unacceptable. If p_E < p ₂, the flow is under-expanded for assigned T ⁰, $p_1^0$ and Ma ₁ > 1: the continuous exact solution developed in § 3 is still acceptable, and the transition from p ₂ to p_E occurs downstream of the nozzle, as a result of 2-D Prandtl–Meyer like expansion waves. Instead, if p_E > p ₂, either a shock within the conical nozzle or 2-D compression waves outside the nozzle can take place, depending on the value of p_E. In the case of the 2-D compression waves outside the nozzle, the continuous exact solution developed in § 3 is still acceptable within the divergent nozzle, while this solution is not acceptable in the case of a shock.

Since a 2-D approach is required to model an oblique shock, only a normal shock can be treated with the here presented quasi-one-dimensional approach. When a normal shock occurs in a divergent nozzle, the continuous exact solution calculated for the assigned T ⁰, $p_1^0$ and Ma ₁ holds until the location of the normal shock. The flow properties then change abruptly across the normal shock, based on the following Rankine–Hugoniot jump conditions (Anderson Reference Anderson1995; Toro Reference Toro2009)

(4.3) \begin{equation} \left.\begin{gathered} {\dfrac{{{p_{do,sh}}}}{{{p_{up,sh}}}} = \dfrac{{2\gamma }}{{\gamma + 1}}Ma_{up,sh}^2 - \dfrac{{\gamma - 1{\kern 1pt} }}{{\gamma + 1}}}\\ {Ma_{do,sh}^2 = \dfrac{{Ma_{up,sh}^2 + \dfrac{2}{{\gamma - 1}}}}{{\dfrac{{2\gamma {\kern 1pt} }}{{\gamma - 1}}Ma_{up,sh}^2 - 1}}}\\ {\dfrac{{{u_{do,sh}}}}{{{u_{up,sh}}}} = \dfrac{{2{\kern 1pt} }}{{\gamma + 1}}\dfrac{1}{{Ma_{up,sh}^2}} + \dfrac{{\gamma - 1{\kern 1pt} }}{{\gamma + 1}}}\end{gathered} \right\},\end{equation}

where p_up,sh, u_up,sh and Ma_up,sh are the flow properties upstream of the normal shock, and p_do,sh, u_do,sh and Ma_do,sh are the corresponding properties downstream of the normal shock. The flow property distributions in the diffuser downstream of the normal shock can be calculated with the solution developed in § 3, with T ⁰, p_do,sh and Ma_do,sh < 1 as the assigned boundary data.

The piecewise analytical solution that includes the normal shock requires that p_E satisfies the following condition:

(4.4)\begin{equation}{p_E} \ge {p_{sh,2}} = \left( {\frac{{2{\kern 1pt} \gamma }}{{\gamma + 1}}Ma_2^2 - \frac{{\gamma - 1{\kern 1pt} }}{{\gamma + 1}}} \right){p_2},\end{equation}

where p_{sh, 2} is the value of the environmental pressure that corresponds to the exact solution with the normal shock at the final section of the nozzle (p ₂ in that case coincides with the pressure value upstream of the shock) and Ma ₂ > 1 is the Mach number value in the final section of the divergent pipe, according to the continuous exact solution determined for the assigned T ⁰, $p_1^0$ and Ma ₁ > 1.

The higher the value of p_E (> p_{sh ,2}) is, the more advanced the position of the normal shock along the divergent nozzle for fixed T ⁰, $p_1^0$ and Ma ₁ > 1. As soon as the normal shock reaches the inlet section, a choking condition is reached for the conical divergent pipe. If the value of p_E is further increased, the boundary conditions at the inlet section of the pipe must be changed in order to obtain a steady-state solution.

If (4.4) is not satisfied, and thus ${p_2} \lt {p_E} \lt {p_{sh,2}}$, oblique shocks or 2-D compression waves occur.