2.1. Governing equations

To non-dimensionalise the governing equations, dimensionless variables are introduced as follows:

(2.1a–e)\begin{equation} {\xi} = \frac{x}{D Re},\quad \eta = \frac{r}{D/2},\quad U= \frac{u}{u_b}, \quad V= \frac{v}{u_b /Re}, \quad P= \frac{p}{\rho u^2_b /2}, \end{equation}

where $\rho$, $D$ and ${u_b}$ represent the fluid density, the pipe diameter and the mean velocity, respectively. The Reynolds number, expressed as $Re= ( \rho D {u_b} )/\mu$, is defined using these reference values with the fluid viscosity $\mu$. All values follow conventional forms (Fargie & Martin Reference Fargie and Martin1971), and the reference pressure is assigned to the dynamic pressure consistent with the friction factor definition. The dimensionless continuity equation and full Navier–Stokes equations are transformed as follows:

(2.2)$$\begin{gather} \frac{\partial U}{\partial \xi} + \frac{2}{\eta} \frac{\partial (\eta V)}{\partial \eta} =0, \end{gather}$$

(2.3)$$\begin{gather}U \frac{\partial U}{\partial \xi} + 2V \frac{\partial U}{\partial \eta} =- \frac{1}{2} \frac{\partial P}{\partial \xi} + \frac{1}{Re^2} \frac{ {\partial}^2 U}{\partial {\xi}^2} + \frac{4}{\eta} \frac{\partial }{\partial \eta} \left(\eta\frac{\partial U}{\partial \eta} \right)\!, \end{gather}$$

(2.4)$$\begin{gather}\frac{1}{Re^2} \left( U \frac{\partial V}{\partial \xi} + 2V \frac{\partial V}{\partial \eta} \right)=- \frac{\partial P}{\partial \eta} + \frac{1}{Re^2} \left[ \frac{1}{Re^2} \frac{ {\partial}^2 V}{\partial {\xi}^2} + \frac{4}{\eta} \frac{\partial}{\partial \eta} \left( \eta \frac{\partial V}{\partial \eta} \right) - \frac{4V}{{\eta}^2} \right]\!. \end{gather}$$

The Reynolds number solely governs the flow. For $Re \gg$ 1, the axial gradient of linear dilatation in (2.3) can be ignored, while (2.4) shows the radial pressure variation to be negligible. Consequently, similar to the boundary layer equation, the axial momentum equation (2.3) simplifies to the parabolised Navier–Stokes equation as follows:

(2.5)\begin{equation} U \frac{\partial U}{\partial \xi} + 2V \frac{\partial U}{\partial \eta} =- \frac{1}{2} \frac{{\rm d} P}{{\rm d} \xi} + \frac{4}{\eta} \frac{\partial }{\partial \eta} \left( \eta \frac{\partial U}{\partial \eta} \right)\!. \end{equation}

The Reynolds number disappears and is grouped within the dimensionless axial coordinate and radial velocity, as shown in (2.1a–e), indicating the dimensionless coordinate system's self-similarity with respect to $Re$. In the entrance flow region, flows with different $Re$ values will exhibit similar velocity profile shapes at specific $\xi$ positions. However, as $Re$ decreases, the omitted terms in the full Navier–Stokes equations gain significance, undermining the self-similarity.

To integrate the axial momentum equation analytically, the flow field is divided into two regions, similar to the previous studies in category (ii); however, distinct assumptions are introduced for each region. Within the inertia-decaying core, the radial variation of axial velocity is markedly diminished, which is attributed to the restricted wall effect, rendering it insignificant in comparison with the axial variation; furthermore, the radial velocity is negligible. Thus, $U_c (\partial U_c / \partial \xi ) \gg 2 V_c ( \partial U_c / \partial \eta )$, where subscript ${\it c}$ indicates the inertia-decaying core. Ignoring radial variation, the inertia term can be approximated by the central axis velocity $U_0$, which is the usual assumption in the inviscid core (Fargie & Martin Reference Fargie and Martin1971; Gupta Reference Gupta1977; Mohanty & Asthana Reference Mohanty and Asthana1978; Smith & Cui Reference Smith and Cui2004), leading to

(2.6)\begin{equation} U_c \frac{\partial U_c}{\partial \xi} + 2 V_c \frac{\partial U_c}{\partial \eta} \approx \frac{1}{2} \frac{{\rm d} U_0^2}{{\rm d} \xi}. \end{equation}

Substituting (2.6) into the momentum equation (2.5), we obtain

(2.7)\begin{equation} \frac{1}{2} \frac{{\rm d} U_0^2}{{\rm d} \xi} =- \frac{1}{2} \frac{{\rm d} P}{{\rm d} \xi} + \frac{4}{\eta} \frac{\partial }{\partial \eta} \left( \eta \frac{\partial U}{\partial \eta} \right)\!. \end{equation}

In contrast to the case of the inviscid core, the shear stress derivative remains, prevailing even at the symmetric central axis where the shear stress is zero. This is confirmed in the fully developed flow, where the shear stress derivative is constant across the entire pipe section. The velocity profile changes from concave to convex depending on the relative magnitudes of the pressure gradient and the inertia term.

The shear stress prevails in the force balance near the pipe wall owing to the no-slip condition, and the inertia terms are disregarded different from the previous studies in category (ii). A sudden change of the uniform velocity profile at the pipe inlet edge may cause a significant momentum difference. This effect is assumed to be confined locally at the inlet edge, and it may cause inaccuracies of present analytical solution in this region. The momentum equation for the wall shear layer is simplified to

(2.8)\begin{equation} 0 =- \frac{1}{2} \frac{{\rm d} P}{{\rm d} \xi} + \frac{4}{\eta} \frac{\partial }{\partial \eta} \left( \eta \frac{\partial U_w}{\partial \eta} \right)\!, \end{equation}

where subscript ${\it w}$ indicates the wall shear layer. An inertia-decaying core, where the flow progresses with increasing velocity and decreasing pressure gradient, is situated adjacent to the wall shear layer. Velocities at the interface of the two flow regions overlap, and the same pressure gradient governs both flows. Under these conditions, the flow in the wall shear layer can be regarded as a Couette flow with varying boundary velocity and pressure gradient. The forward pressure gradient results in a convex velocity profile, and a velocity overshoot arises from the prevailing effect of the pressure gradient over the interface velocity.

2.2. Velocity solutions

The interface radius between the two flow regions, denoted as ${\eta }_{\delta } ( \xi )$, is a new unknown variable, which separates the flow velocity as follows:

(2.9)\begin{equation} U (\xi , \eta ) = \left\{\begin{array}{@{}ll@{}} U_c ( \xi , \eta ), & 0 \le \eta \le {\eta}_{\delta} \\ U_w ( \xi , \eta ), & {\eta}_{\delta} \le \eta \le 1 \end{array}\right.\!. \end{equation}

As the wall shear layer extends to the central axis, the inertia-decaying core shrinks. When the inertia term in the inertia-decaying core sufficiently diminishes, the momentum equations for each flow region are approximately equivalent, and the velocity profiles overlap significantly. Beyond this point, the interface radius no longer changes, and the flow becomes fully developed as the inertia term vanishes completely.

The momentum equation (2.5) is approximated for each flow region, and the resultant equations (2.7) and (2.8) can be integrated analytically. Two boundary conditions are available – one for the inertia-decaying core at the symmetric central axis and the other for the wall shear layer at the no-slip wall – expressed as follows:

(2.10a,b)\begin{equation} \frac{\partial U_c ( \xi ,0)}{\partial \eta} =0, \quad U_w ( \xi ,1) =0. \end{equation}

At the interface radius, the velocities and their slopes (shear stresses) must be equal, ensuring a continuous function of velocity across the two regions

(2.11a,b)\begin{equation} U_c ( \xi , {\eta}_{\delta}) = U_w ( \xi , {\eta}_{\delta}), \quad \frac{\partial U_c ( \xi , {\eta}_{\delta})}{\partial \eta} = \frac{\partial U_w ( \xi , {\eta}_{\delta})}{\partial \eta}. \end{equation}

Applying the boundary conditions (2.10a,b) and interface-matching conditions (2.10a,b), four integral coefficients are determined, and the velocities are yielded as

(2.12)$$\begin{gather} U_c ( \xi ,\eta )=- \frac{D_{\! P}}{32} (1- {\eta}^2) + \frac{D_U}{32} {\eta}^2 - \frac{D_U}{32} {\eta}_{\delta}^2 (1- \ln {\eta}_{\delta}^2), \end{gather}$$

(2.13)$$\begin{gather}U_w ( \xi ,\eta )=- \frac{D_{\! P}}{32} (1- {\eta}^2) + \frac{D_U}{32} {\eta}_{\delta}^2 \ln {\eta}^2, \end{gather}$$

where $D_{\! P} ( \xi )$ and $D_U ( \xi )$ are introduced to reduce the complexity of the expression and are defined as follows:

(2.14a,b)\begin{equation} D_{\! P} ( \xi ) \equiv \frac{{\rm d} P}{{\rm d} \xi},\quad D_U ( \xi ) \equiv \frac{{\rm d} U_0^2}{{\rm d} \xi}. \end{equation}

The central axis velocity can be deduced from (2.12) at $\eta = 0$, from which the pressure gradient $D_P ( \xi )$ is expressed as

(2.15)\begin{equation} D_{\! P} ( \xi ) =-32 U_0 - D_U {\eta}_{\delta}^2 (1- \ln {\eta}_{\delta}^2). \end{equation}

The pressure gradient is removed from (2.12) and (2.13) to obtain

(2.16)$$\begin{gather} U_c ( \xi ,\eta )= U_0 ( 1- {\eta}^2 ) + \frac{D_U}{32} [ 1- {\eta}_{\delta}^2 ( 1- \ln {\eta}_{\delta}^2 )] {\eta}^2, \end{gather}$$

(2.17)$$\begin{gather}U_w ( \xi ,\eta )= U_c ( \xi ,\eta ) + \frac{D_U}{32} [{\eta}_{\delta}^2 (1- \ln {\eta}_{\delta}^2)- {\eta}^2 + {\eta}_{\delta}^2 \ln {\eta}^2]. \end{gather}$$

The velocity at an axial position $\xi$ is represented by a function of $\eta$, with two axially varying parameters $U_0 ( \xi )$ and ${\eta }_{\delta } ( \xi )$. The remaining $D_U ( \xi )$ is calculated from $U_0 ( \xi )$. Two additional equations are needed to determine these unknown parameters.

The continuity and the momentum equations integrated over the pipe cross-section are referred to as global equations. They are useful in analysing the internal flow, and the equations necessary to solve $U_0 ( \xi )$ and ${\eta }_{\delta } ( \xi )$ are derived from them. The continuity equation (2.2) is integrated over the pipe cross-section

(2.18)\begin{equation} \frac{{\rm d}}{{\rm d} \xi} \int_0^1 U( \xi,\eta ) \eta \,{\rm d} \eta + 2 \int_0^1 \frac{\partial ( \eta V)}{\partial \eta} \, {\rm d} \eta =0. \end{equation}

Due to the boundary conditions, the second term of the left-hand side of (2.18) becomes zero. Integrating equation (2.18) over $\xi$ results in

(2.19)\begin{equation} \int_0^1 U( \xi ,\eta ) \eta \,{\rm d} \eta = {\rm constant}. \end{equation}

Because the left-hand side of (2.19) represents the total mass flow rate at any $\xi$, by applying the uniform inlet velocity condition, the constant on the right-hand side is 1/2. Eventually, the global continuity equation can be derived by substituting the velocity equations (2.16) and (2.17) and integrating analytically for each region

(2.20)\begin{equation} \int_0^{{\eta}_{\delta}} U_c( \xi ,\eta ) \eta \, {\rm d} \eta + \int_{{\eta}_{\delta}}^1 U_w( \xi ,\eta ) \eta \, {\rm d} \eta = \tfrac{1}{2}. \end{equation}

Following mathematical modification, the global continuity equation can be expressed as

(2.21)\begin{equation} \frac{{\rm d} U_0^2 ( \xi )}{{\rm d} \xi} = \frac{32 ( U_0 -2) }{{\eta}_{\delta}^2 ( 1- {\eta}_{\delta}^2 + \ln {\eta}_{\delta}^2)}. \end{equation}

In addition, the momentum equation (2.5) is integrated over the pipe cross-section

(2.22)\begin{equation} \frac{{\rm d}} {{\rm d} \xi} \int_0^1 U^2 \eta \, {\rm d} \eta +2 \int_0^1 \frac{\partial ( \eta UV)}{\partial \eta} {\rm d} \eta =- \frac{D_{\! P}}{2} \int_0^1 \eta \, {\rm d} \eta + 4 \int_0^1 \frac{\partial}{\partial \eta} \left( \eta \frac{\partial U}{\partial \eta} \right) {\rm d} \eta. \end{equation}

Here, the left-hand side of (2.5) is switched to conservative form using the continuity equation. The second term of the left-hand side is removed owing to the boundary conditions, and the global momentum equation is obtained

(2.23)\begin{equation} \frac{{\rm d} \varTheta ( \xi )}{{\rm d} \xi} =-\frac{D_{\! P}}{4} + 4 \left( \eta \frac{\partial U_w}{\partial \eta} \right)_{\eta = 1}, \end{equation}

where the global momentum $\varTheta ( \xi )$ is defined as

(2.24)\begin{equation} \varTheta ( \xi ) = \int_0^1 U^2 \eta \, {\rm d} \eta. \end{equation}

Substituting (2.15), (2.17) and (2.21) into (2.23), the global momentum equation becomes

(2.25)\begin{equation} \frac{{\rm d} \varTheta (\xi )}{{\rm d} \xi} = \frac{8 ( U_0 -2)}{1- {\eta}_{\delta}^2 + \ln {\eta}_{\delta}^2}. \end{equation}

The global momentum in (2.24) is partially integrated by applying (2.16) and (2.17), resulting in

(2.26)\begin{align} \varTheta (\xi ) = \frac{2}{3} - \frac{( U_0 -2) ( 1- {\eta}_{\delta}^2)^2 }{6 ( 1- {\eta}_{\delta}^2 + \ln {\eta}_{\delta}^2)}+ \frac{( U_0 -2)^2 ( 2+3 {\eta}_{\delta}^2 -6 {\eta}_{\delta}^4 + {\eta}_{\delta}^6 +6 {\eta}_{\delta}^2 \ln {\eta}_{\delta}^2) }{12 ( 1- {\eta}_{\delta}^2 + \ln {\eta}_{\delta}^2)^2}. \end{align}

There are three unknowns ($U_0 ( \xi )$, ${\eta }_{\delta } ( \xi )$ and ${\it \varTheta } ( \xi )$) to solve, and we have two first-order ordinary differential equations (2.21) and (2.25) and one algebraic equation (2.26). To create a closed form for the equation set, three inlet conditions are required to integrate the differential equations. The central axis velocity at the inlet is $U_0 (0)= 1$ due to the uniform inlet velocity condition, and applying it to (2.24) results in $\varTheta ( 0 )= 1/2$. Because the wall shear layer takes place from the pipe wall, the interface radius at the inlet is ${\eta }_{\delta } ( 0 )= 1$. The ordinary differential equations (2.21) and (2.25) were integrated using the fourth-order Runge–Kutta method. At each prediction and correction step of the Runge–Kutta integration, ${\eta }_{\delta } ( \xi )$ was calculated from (2.26) using the bisection method.

2.3. Supplementary equations

The total pressure drop from the inlet can be calculated by integrating equation (2.23) with respect to $\xi$, which is expressed as follows:

(2.27)\begin{equation} {\rm \Delta} P ( \xi ) = 4 \left[ {\it \varTheta} ( \xi ) - \tfrac{1}{2} \right] + \int_0^{\xi} (32 U_0 - D_U {\eta}_{\delta}^2 \ln {\eta}_{\delta}^2)\,{\rm d} \hat{\xi}. \end{equation}

The term within the square brackets represents the net change in the inertia force, and the integral is the shear force at the pipe wall. The fourth-order Runge–Kutta method was adopted for the integration.

The velocity functions are continuous across two flow regions, and the streamlines can be obtained using the streamfunction defined as

(2.28)\begin{equation} \varPsi ( \xi, \eta ) = \tfrac{1}{2} \int_0^{\eta} U( \xi , \eta ) \hat{\eta} \, {\rm d} \hat{\eta}. \end{equation}

Substituting (2.16) and (2.17) into (2.28) yields the following streamfunctions:

(2.29)$$\begin{gather} \varPsi_c ( \xi, \eta ) = \frac{U_0}{8} ( 2- \eta^2) \eta^2 + \frac{D_U}{256} [(1- {\eta}_{\delta}^2)^2 + {\eta}_{\delta}^2 ( 1- {\eta}_{\delta}^2 + \ln {\eta}_{\delta}^2)] {\eta}^4, \end{gather}$$

(2.30)$$\begin{gather}\varPsi_w ( \xi , \eta ) = {\it \varPsi_c} ( \xi , \eta ) + \frac{D_U}{256} [ {\eta}_{\delta}^4 +2 {\eta}_{\delta}^2 (1- \ln {\eta}_{\delta}^2) {\eta}^2 - {\eta}^4 - 2 {\eta}_{\delta}^2 (1 - \ln {\eta}^2) {\eta}^2]. \end{gather}$$

The radial velocity can be determined by differentiating the streamfunction

(2.31)\begin{equation} V ( \xi, \eta ) =- \frac{1}{\eta} \frac{\partial \varPsi ( \xi , \eta ) }{\partial \xi}. \end{equation}

The central differencing method was used for this differentiation.