2.1. Governing equations

We consider the isothermal and incompressible steady flow of a viscoelastic fluid which consists of a Newtonian solvent with constant shear viscosity $\eta _s^\ast $ and a polymeric material with constant shear viscosity $\eta _p^\ast $ and longest relaxation time ${\lambda ^\ast }$. The fluid has constant mass density denoted by ${\rho ^\ast }$. We use a cylindrical coordinate system $({z^\ast },{y^\ast },\theta )$ to describe the flow field, where z* is the main flow direction, y* is the radial direction and $\theta$ is the azimuthal angle; ${\boldsymbol{e}_\theta }$, ${\boldsymbol{e}_z}$ and ${\boldsymbol{e}_y}$ are the unit vectors in the $\theta$, z* and y* directions, respectively. The origin of the coordinate system is placed on the centre of the inlet cross-section with radius $h_0^\ast $ (see figure 1). The wall of the pipe is described by the shape function ${H^\ast } = {H^\ast }({z^\ast }) > 0$ for ${z^\ast } \in [0,{\ell ^\ast }]$, i.e. ${y^\ast } = {H^\ast }({z^\ast })$. The shape function ${H^\ast }$ is considered fixed and known in advance; for ${H^\ast }({z^\ast }) = h_0^\ast $ one gets a straight circular pipe. Also, we exclude variations of the flow field in the azimuthal angle which implies that the flow is axisymmetric and two-dimensional. The velocity vector in the flow domain is denoted by ${\boldsymbol{u}^\ast } = {V^\ast }(\kern0.08em {y^\ast },{z^\ast }){\boldsymbol{e}_y} + {U^\ast }(\kern0.08em {y^\ast },{z^\ast }){\boldsymbol{e}_z}$ and the total pressure is denoted by ${p^\ast } = {p^\ast }{(\kern0.08em {y^\ast },{z^\ast })^\ast }$. The latter is induced by the flow and the imposed constant flow rate Q* at the inlet,

(2.3)\begin{equation}{Q^\ast } = 2{\rm \pi} \int_0^{h_0^\ast } {{U^\ast }(\kern0.08em {y^\ast },0){y^\ast }\,\textrm{d}{y^\ast }} .\end{equation}

Using the unit tensor ${\boldsymbol{\mathsf{I}}}$, the rate of deformation tensor ${\dot{\boldsymbol{\gamma }}^\ast } = {\nabla ^\ast }{\boldsymbol{u}^\ast } + {({\nabla ^\ast }{\boldsymbol{u}^\ast })^\textrm{T}}$ where

(2.4)\begin{equation}{\nabla ^\ast } \equiv {\boldsymbol{e}_y}\frac{\partial }{{\partial {y^\ast }}} + {\boldsymbol{e}_z}\frac{\partial }{{\partial {z^\ast }}} + {\boldsymbol{e}_\theta }\frac{1}{{{y^\ast }}}\frac{\partial }{{\partial \theta }},\end{equation}

is the gradient operator, the viscoelastic extra-stress tensor ${\boldsymbol{\tau }^\ast }$, and the Reynolds stress tensor ${\boldsymbol{u}^\ast }{\boldsymbol{u}^\ast }$; we define the total momentum tensor as

(2.5)\begin{equation}{\boldsymbol{T}^\ast }: ={-} {P^\ast }{\boldsymbol{\mathsf{I}}} + \eta _s^\ast {\dot{\boldsymbol{\gamma }}^\ast } + {\boldsymbol{\tau }^\ast } - {\rho ^\ast }{\boldsymbol{u}^\ast }{\boldsymbol{u}^\ast }.\end{equation}

In the absence of any other external forces and torques, the conservation equations that govern the flow in the pipe are the mass and momentum balances, respectively,

(2.6a,b)\begin{equation}{\nabla ^\ast }\boldsymbol{\cdot }{\boldsymbol{u}^\ast } = 0,\quad {\nabla ^\ast }\boldsymbol{\cdot }{\boldsymbol{T}^\ast } = {\bf 0}.\end{equation}

Equations (2.6a,b) are accompanied with a constitutive equation that models the response of the polymeric material to the flow deformation, namely with an equation for ${\boldsymbol{\tau }^\ast }$. Specifically, we consider the fundamental nonlinear differential Oldroyd-B models, given by

(2.7)\begin{equation}{\boldsymbol{\tau }^\ast } + {\lambda ^\ast }\frac{{{\delta ^\ast }{\boldsymbol{\tau }^\ast }}}{{\delta {t^\ast }}} = \eta _p^\ast {\dot{\boldsymbol{\gamma }}^\ast },\end{equation}

where ${\delta ^\ast }{\boldsymbol{\tau }^\ast }/\delta {t^\ast }$ represents the upper convected derivative of ${\boldsymbol{\tau }^\ast }$ at steady state,

(2.8)\begin{equation}\frac{{{\delta ^\ast }{\boldsymbol{\tau }^\ast }}}{{\delta {t^\ast }}} = {\boldsymbol{u}^\ast }\boldsymbol{\cdot }{\nabla ^\ast }{\boldsymbol{\tau }^\ast } - {\boldsymbol{\tau }^\ast }\boldsymbol{\cdot }{\nabla ^\ast }{\boldsymbol{v}^\ast } - {({\nabla ^\ast }{\boldsymbol{v}^\ast })^\textrm{T}}\boldsymbol{\cdot }{\boldsymbol{\tau }^\ast }.\end{equation}

The domain of definition of (2.4)–(2.8) is

(2.9)\begin{equation}{\varOmega ^\ast } = \{ ({z^\ast },{y^\ast },\theta )|,\;0 < {z^\ast } < {\ell ^\ast },\;0 < {y^\ast } < {H^\ast }({z^\ast }),\;0 \le \theta < 2{\rm \pi} \} .\end{equation}

The governing equations are solved with no-slip and no-penetration boundary conditions along the wall of the pipe,

(2.10)\begin{equation}{U^\ast }({H^\ast }({z^\ast }),{z^\ast }) = {V^\ast }({H^\ast }({z^\ast }),{z^\ast }) = 0\quad \textrm{at}\;0 \le {z^\ast } \le {\ell ^\ast }.\end{equation}

All the equations are well-defined and all the dependent variables are finite in the entire axisymmetric domain satisfying the following conditions at ${y^\ast } = 0$:

(2.11a)\begin{gather}{V^\ast }(0,{z^\ast }) = {\left. {\frac{{\partial {U^\ast }}}{{\partial {y^\ast }}}} \right|_{{y^\ast } = 0}} = \tau _{yz}^\ast (0,{z^\ast }) = 0\quad \textrm{at}\;0 \le {z^\ast } \le {\ell ^\ast },\end{gather}

(2.11b)\begin{gather}\tau _{yy}^\ast (0,{z^\ast }) = \tau _{\theta \theta }^\ast (0,{z^\ast })\quad \textrm{at}\;0 \le {z^\ast } \le {\ell ^\ast }.\end{gather}

The integral constraint of mass at any distance from the inlet is also utilized,

(2.12)\begin{equation}{Q^\ast } = 2{\rm \pi} \int_0^{{H^\ast }({z^\ast })} {{U^\ast }(\kern0.08em {y^\ast },{z^\ast }){y^\ast }\,\textrm{d}{y^\ast }} = \textrm{constant}\quad \textrm{at}\;0 \le {z^\ast } \le {\ell ^\ast }.\end{equation}

Finally, a datum pressure, $p_{ref}^\ast $, is chosen at the wall of the outlet cross-section,

(2.13)\begin{equation}p_{ref}^\ast = {p^\ast }({H^\ast }({\ell ^\ast }),{\ell ^\ast }).\end{equation}

Note that when the Newtonian solvent is absent ($\eta _s^\ast = 0$) the governing equations reduce to the UCM model, or to the Oldroyd-B model when both the Newtonian solvent and the polymer molecules are present, $\eta _s^\ast ,\,\eta _p^\ast > 0$.

2.2. Scaling and non-dimensionalization

Dimensionless variables are introduced based on the lubrication theory using the transformation $X = {X^\ast }/X_c^\ast $ where ${X^\ast } = {z^\ast },\;{y^\ast },\;{H^\ast },\;{U^\ast },\;{V^\ast },\;{P^\ast },\;\tau _{zz}^\ast ,\;\tau _{yz}^\ast ,\;\tau _{yy}^\ast ,\;\tau _{\theta \theta }^\ast $ and $X_c^\ast $ is the relevant characteristic scale for ${X^\ast }$. The characteristic scales are defined with the aid of the aspect ratio of the pipe and one half of the cross-sectionally average velocity at the inlet of the pipe defined in (2.1a,b). The z*-coordinate is scaled by ${\ell ^\ast }$ and the ${y^\ast }$-coordinate and the shape function H* by $h_0^\ast $. The characteristic scale for the main velocity component ${U^\ast }$ is $u_c^\ast $, and from the continuity equation one finds that the characteristic scale for the vertical velocity component ${V^\ast }$ is $\varepsilon u_c^\ast $. With the aid of the momentum balance along the main flow direction, z*, the characteristic scale for the pressure difference ${P^\ast } - P_{ref}^\ast $ is $(\eta _s^\ast + \eta _p^\ast )u_c^\ast {\ell ^\ast }/h_0^{{\ast} 2}$ (Tavakol et al. Reference Tavakol, Froehlicher, Holmes and Stone2017; Housiadas & Tsagaris Reference Housiadas and Tsangaris2022, Reference Housiadas and Tsangaris2023; Housiadas & Beris Reference Housiadas and Beris2023, Reference Housiadas and Beris2024a,Reference Housiadas and Berisb,Reference Housiadas and Berisc). For the viscoelastic extra-stress components, $\tau _{zz}^\ast $, $\tau _{yz}^\ast $, $\tau _{yy}^\ast $ and $\tau _{\theta \theta }^\ast $, consistency of the governing equations for a negligible relaxation time of the polymer (${\lambda ^\ast } \to 0$) and a long pipe ($h_0^\ast\ll {\ell ^\ast }$) implies that different scales for the components of ${\boldsymbol{\tau }^\ast }$ must be used. Careful examination of the momentum balance and the individual components of the constitutive equation shows that the characteristic scales for $\tau _{yz}^\ast $, $\tau _{zz}^\ast $, $\tau _{yy}^\ast $ and $\tau _{\theta \theta }^\ast $ are $\eta _p^\ast u_c^\ast{/}h_0^\ast $, $\eta _p^\ast u_c^\ast {\ell ^\ast }/h_0^{{\ast} 2}$, $\eta _p^\ast u_c^\ast{/}{\ell ^\ast }$ and $\eta _p^\ast u_c^\ast{/}{\ell ^\ast }$, respectively. Based on these scales, the velocity vector and the Reynolds stress tensor are given as, respectively,

(2.14)\begin{gather}\boldsymbol{u} = U{\boldsymbol{e}_z} + \varepsilon V{\boldsymbol{e}_y},\end{gather}

(2.15)\begin{gather}\boldsymbol{uu} = {\boldsymbol{e}_z}{\boldsymbol{e}_z}{U^2} + {\boldsymbol{e}_z}{\boldsymbol{e}_y}\varepsilon UV + {\boldsymbol{e}_y}{\boldsymbol{e}_z}\varepsilon UV + {\boldsymbol{e}_y}{\boldsymbol{e}_y}{\varepsilon ^2}{V^2}.\end{gather}

Similarly, with the aid of the gradient operator

(2.16)\begin{equation}\nabla \equiv h_0^\ast {\nabla ^\ast } = {\boldsymbol{e}_y}\frac{\partial }{{\partial y}} + \varepsilon {\boldsymbol{e}_z}\frac{\partial }{{\partial z}} + {\boldsymbol{e}_\theta }\frac{1}{y}\frac{\partial }{{\partial \theta }},\end{equation}

we find that the rate of deformation tensor, scaled by $u_c^\ast{/}h_0^\ast $, is

(2.17)\begin{align}\dot{\boldsymbol{\gamma }} &= {\boldsymbol{e}_z}{\boldsymbol{e}_z}\left( {2\varepsilon \frac{{\partial U}}{{\partial z}}} \right) + {\boldsymbol{e}_z}{\boldsymbol{e}_y}\left( {\frac{{\partial U}}{{\partial y}} + {\varepsilon^2}\frac{{\partial V}}{{\partial z}}} \right) + {\boldsymbol{e}_y}{\boldsymbol{e}_z}\left( {\frac{{\partial U}}{{\partial y}} + {\varepsilon^2}\frac{{\partial V}}{{\partial z}}} \right) + {\boldsymbol{e}_y}{\boldsymbol{e}_y}\left( {2\varepsilon \frac{{\partial V}}{{\partial y}}} \right)\nonumber\\ &\quad + {\boldsymbol{e}_\theta }{\boldsymbol{e}_\theta }\left( {2\varepsilon \frac{V}{y}} \right).\end{align}

Also, the dimensionless polymer extra-stress tensor, scaled by $\eta _p^\ast u_c^\ast{/}h_0^\ast $, is

(2.18)\begin{equation}\boldsymbol{\tau } = {\boldsymbol{e}_z}{\boldsymbol{e}_z}{\tau _{zz}}/\varepsilon + {\boldsymbol{e}_z}{\boldsymbol{e}_y}{\tau _{yz}} + {\boldsymbol{e}_y}{\boldsymbol{e}_z}{\tau _{yz}} + {\boldsymbol{e}_y}{\boldsymbol{e}_y}\varepsilon {\tau _{yy}} + {\boldsymbol{e}_\theta }{\boldsymbol{e}_\theta }\varepsilon {\tau _{\theta \theta }}.\end{equation}

Finally, the total momentum tensor, scaled by $(\eta _s^\ast + \eta _p^\ast )u_c^\ast {\ell ^\ast }/h_0^{{\ast} 2}$, is given by

(2.19)\begin{equation}\boldsymbol{T} ={-} P{\boldsymbol{\mathsf{I}}} + \varepsilon (1 - \eta )\dot{\boldsymbol{\gamma }} + \varepsilon \eta \boldsymbol{\tau } - R{e_m}\boldsymbol{uu}.\end{equation}

In (2.19), the modified Reynolds number, $R{e_m}$, and the polymer viscosity ratio, $\eta$, appear:

(2.20a,b)\begin{equation}R{e_m} \equiv \frac{{{\rho ^\ast }u_c^\ast h_0^{{\ast} 2}}}{{(\eta _s^\ast + \eta _p^\ast ){\ell ^\ast }}},\quad \eta \equiv \frac{{\eta _p^\ast }}{{\eta _s^\ast + \eta _p^\ast }}.\end{equation}

The dimensionless form of the balance equations in scalar form is

(2.21)\begin{gather}\frac{{\partial U}}{{\partial z}} + \frac{{\partial V}}{{\partial y}} + \frac{V}{y} = 0,\end{gather}

(2.22)\begin{gather}R{e_m}\frac{{\textrm{D}U}}{{\textrm{D}t}} ={-} \frac{{\partial P}}{{\partial z}} + (1 - \eta )\left( {\frac{{{\partial^2}U}}{{\partial {y^2}}} + \frac{1}{y}\frac{{\partial U}}{{\partial y}} + {\varepsilon^2}\frac{{{\partial^2}U}}{{\partial {z^2}}}} \right) + \eta \left( {\frac{{\partial {\tau_{yz}}}}{{\partial y}} + \frac{{{\tau_{yz}}}}{y} + \frac{{\partial {\tau_{zz}}}}{{\partial z}}} \right),\end{gather}

(2.23)\begin{gather}\begin{aligned}[b] {\varepsilon ^2}R{e_m}\frac{{\textrm{D}V}}{{\textrm{D}t}} &={-} \frac{{\partial P}}{{\partial y}} + (1 - \eta ){\varepsilon ^2}\left( {\frac{{{\partial^2}V}}{{\partial {y^2}}} + \frac{1}{y}\frac{{\partial V}}{{\partial y}} - \frac{V}{{{y^2}}} + {\varepsilon^2}\frac{{{\partial^2}V}}{{\partial {z^2}}}} \right)\\ &\quad + \eta {\varepsilon ^2}\left( {\frac{{\partial {\tau_{yy}}}}{{\partial y}} + \frac{{\partial {\tau_{yz}}}}{{\partial z}} + \frac{{{\tau_{yy}} - {\tau_{\theta \theta }}}}{y}} \right),\end{aligned}\end{gather}

where $\textrm{D}/\textrm{D}t$ is the material derivative defined at steady state as

(2.24)\begin{equation}\frac{\textrm{D}}{{\textrm{D}t}} \equiv U\frac{\partial }{{\partial z}} + V\frac{\partial }{{\partial y}}.\end{equation}

For $\eta = 0$, i.e. in absence of the polymeric molecules, (2.21)–(2.24) reduce to the well-known dimensionless Navier–Stokes equations for an incompressible Newtonian fluid. Notice that for $\eta = 1$ there is no solvent contribution, i.e. the fluid is a pure polymeric material, while for $0 < \eta < 1$ the fluid consists of a Newtonian solvent and a polymeric material.

The UCM/Oldroyd-B models in scalar form are given as

(2.25)\begin{gather}{\tau _{zz}} + De\left( {\frac{{\textrm{D}{\tau_{zz}}}}{{\textrm{D}t}} - 2{\tau_{zz}}\frac{{\partial U}}{{\partial z}} - 2{\tau_{yz}}\frac{{\partial U}}{{\partial y}}} \right) = 2{\varepsilon ^2}\frac{{\partial U}}{{\partial z}},\end{gather}

(2.26)\begin{gather}{\tau _{yy}} + De\left( {\frac{{\textrm{D}{\tau_{yy}}}}{{\textrm{D}t}} - 2{\tau_{yz}}\frac{{\partial V}}{{\partial z}} - 2{\tau_{yy}}\frac{{\partial V}}{{\partial y}}} \right) = 2\frac{{\partial V}}{{\partial y}},\end{gather}

(2.27)\begin{gather}{\tau _{\theta \theta }} + De\left( {\frac{{\textrm{D}{\tau_{\theta \theta }}}}{{\textrm{D}t}} - 2{\tau_{\theta \theta }}\frac{V}{y}} \right) = 2\frac{V}{y},\end{gather}

(2.28)\begin{gather}{\tau _{yz}} + De\left( {\frac{{\textrm{D}{\tau_{yz}}}}{{\textrm{D}t}} + {\tau_{yz}}\frac{V}{y} - {\tau_{zz}}\frac{{\partial V}}{{\partial z}} - {\tau_{yy}}\frac{{\partial U}}{{\partial y}}} \right) = \frac{{\partial U}}{{\partial y}} + {\varepsilon ^2}\frac{{\partial V}}{{\partial z}}.\end{gather}

The dimensionless domain of definition of (2.21)–(2.28) is

(2.29)\begin{equation}\varOmega = \{ (z,y,\theta )|0 < z < 1,\;0 < y < H(z),0 \le \theta < 2{\rm \pi} \} .\end{equation}

Various cases can be considered for the shape function but here we focus on the hyperbolic pipe,

(2.30)\begin{equation}H(z) = \frac{1}{{\sqrt {1 + ({\varLambda ^2} - 1)z} }}\quad 0 \le z \le 1,\end{equation}

where $\varLambda = h_0^\ast{/}h_f^\ast $. In the entrance region ($z \le 0$) $H(z) = 1$, while in the exit region ($z \ge 1$), $H(z) = 1/\varLambda$. Thus, $H = H(z)$ is a piecewise smooth function which is continuous in the whole domain but not differentiable at $z = 0$ and $z = 1$. For $\varLambda > 1$, (2.30) describes a converging pipe, while for $0 < \varLambda < 1$ it describes an expanding pipe; in the former case $\varLambda $ will be referred to as the contraction ratio. Note that for axisymmetric geometries, experimentalists relate the contraction ratio to the Hencky strain, ${\varepsilon _H}$, as ${\varepsilon _H} = \textrm{ln}(({\rm \pi} h_0^{{\ast} 2})/({\rm \pi} h_f^{{\ast} 2})) = \textrm{ln}({\varLambda ^2})$ which, of course, implies $\varLambda = \textrm{exp}({\varepsilon _H}/2)$. Typical contraction ratios used in experiments are $\varLambda \sim 3\hbox{--} 7$; for instance, James & Roos (Reference James and Roos2021) and James & Tripathi (Reference James and Tripathi2023) used $\varLambda = \sqrt {117/5.8} \approx 4.49$, Kamerkar & Edwards (Reference Kamerkar and Edwards2007) used Hencky strain values in the range $4 \le {\varepsilon _H} \le 7$ and Collier et al. (Reference Collier, Romanoschi and Petrovan1998) used ${\varepsilon _H}$ in the range $6 \le {\varepsilon _H} \le 7$.

The auxiliary dimensionless conditions (boundary conditions, total mass balance, constraints on the axis of symmetry and the datum pressure) are

(2.31a)\begin{gather}V = U = 0\quad \textrm{at}\;y = H(z),\;0 \le z \le 1,\end{gather}

(2.31b)\begin{gather}V = \frac{{\partial U}}{{\partial y}} = {\tau _{yz}} = 0\quad \textrm{at}\;y = 0,\;0 \le z \le 1,\end{gather}

(2.31c)\begin{gather}{\tau _{yy}} = {\tau _{\theta \theta }}\quad \textrm{at}\;y = 0,\;0 \le z \le 1,\end{gather}

(2.31d)\begin{gather}\int_0^{H(z)} {U(y,z)y\,\textrm{d}y} = 1,\quad 0 \le z \le 1,\end{gather}

(2.31e)\begin{gather}P(H(1),1) = 0.\end{gather}

Equations (2.31b,c) are needed so that (2.21)–(2.23) are well-defined at y = 0, while (2.31d) results from the integration of (2.21) along the y-direction, between the limits $y = 0$ and $y = H(z)$, applying the no-penetration boundary condition and using the dimensionless flow rate at the inlet of the pipe.

2.3. The conformation tensor

As mentioned above, the constitutive models for the polymeric molecules are frequently provided in terms of the conformation tensor ${\boldsymbol{\mathsf{C}}}$ (dimensionless). For the Oldroyd-B model, and its limiting form, the UCM model, the conformation tensor is related to the extra-stress tensor (Bird, Armstrong & Hassager Reference Bird, Armstrong and Hassager1987; Beris & Edwards Reference Beris and Edwards1994) as follows:

(2.32)\begin{equation}{\boldsymbol{\tau }^\ast } = \frac{{\eta _p^\ast }}{{{\lambda ^\ast }}}({\boldsymbol{\mathsf{C}}} - {\boldsymbol{\mathsf{I}}}).\end{equation}

For a Newtonian fluid, i.e. for ${\lambda ^\ast } \to 0$, the conformation tensor reduces to the identity tensor, ${\boldsymbol{\mathsf{C}}} \to {\boldsymbol{\mathsf{I}}}$, and the extra-stress tensor to the rate-of-deformation tensor, ${\boldsymbol{\tau }^\ast } \to {\dot{\boldsymbol{\gamma }}^\ast }$; for no-flow ${\boldsymbol{\mathsf{C}}} \to {\boldsymbol{\mathsf{I}}}$ and ${\boldsymbol{\tau }^\ast } \to {\bf 0}$. An important property of the real, second-order and symmetric tensor ${\boldsymbol{\mathsf{C}}}$ is that is positive definite, and therefore its eigenvalues and its determinant, are strictly positive.

Based on the dimensionless form of the extra-stress tensor, $\boldsymbol{\tau }$, the conformation tensor, ${\boldsymbol{\mathsf{C}}}$, is given through $\boldsymbol{\tau } = ({\boldsymbol{\mathsf{C}}} - {\boldsymbol{\mathsf{I}}})/Wi = \varepsilon ({\boldsymbol{\mathsf{C}}} - {\boldsymbol{\mathsf{I}}})/De$ where

(2.33a)\begin{equation}{\boldsymbol{\mathsf{C}}} = {\boldsymbol{e}_z}{\boldsymbol{e}_z}\underbrace{{({\varepsilon ^2}{C_{zz}})}}_{{{c_{zz}}}}/{\varepsilon ^2} + {\boldsymbol{e}_z}{\boldsymbol{e}_y}\underbrace{{(\varepsilon {C_{yz}})}}_{{{c_{yz}}}}/\varepsilon + {\boldsymbol{e}_y}{\boldsymbol{e}_z}\underbrace{{(\varepsilon {C_{yz}})}}_{{{c_{yz}}}}/\varepsilon + {\boldsymbol{e}_y}{\boldsymbol{e}_y}\underbrace{{{C_{yy}}}}_{{{c_{yy}}}} + {\boldsymbol{e}_\theta }{\boldsymbol{e}_\theta }\underbrace{{{C_{\theta \theta }}}}_{{{c_{\theta \theta }}}},\end{equation}

or

(2.33b)\begin{equation}{\boldsymbol{\mathsf{C}}} = {\boldsymbol{e}_z}{\boldsymbol{e}_z}{c_{zz}}/{\varepsilon ^2} + {\boldsymbol{e}_z}{\boldsymbol{e}_y}{c_{yz}}/\varepsilon + {\boldsymbol{e}_y}{\boldsymbol{e}_z}{c_{yz}}/\varepsilon + {\boldsymbol{e}_y}{\boldsymbol{e}_y}{c_{yy}} + {\boldsymbol{e}_\theta }{\boldsymbol{e}_\theta }{c_{\theta \theta }}.\end{equation}

From (2.18) and (2.33b), we find that the components of $\boldsymbol{\tau }$ and the rescaled components of ${\boldsymbol{\mathsf{C}}}$ are connected as

(2.34a–d)\begin{equation}{\tau _{zz}} = \frac{{{c_{zz}} - {\varepsilon ^2}}}{{De}},\quad {\tau _{yz}} = \frac{{{c_{yz}}}}{{De}},\quad {\tau _{yy}} = \frac{{{c_{yy}} - 1}}{{De}},\quad {\tau _{\theta \theta }} = \frac{{{c_{\theta \theta }} - 1}}{{De}}.\end{equation}

One is easy to confirm that the Oldroyd-B model can be equivalently and consistency expressed in terms of ${\tau _{zz}},\;{\tau _{yz}},\;{\tau _{yy}}$ and ${\tau _{\theta \theta }}$ or ${c_{zz}},\;{c_{yz}},\;{c_{yy}}$ and ${c_{\theta \theta }}$; the same expressions for ${\tau _{zz}}$, ${\tau _{yz}}$ and ${\tau _{yy}}$ have also been derived by Boyko & Stone (Reference Boyko and Stone2022), Ahmed & Biancofiore (Reference Ahmed and Biancofiore2021, Reference Ahmed and Biancofiore2023) and Housiadas & Beris (Reference Housiadas and Beris2023) for the viscoelastic lubrication flow in a hyperbolic channel. Note that an unstretched or equilibrium state for the polymer molecules implies that

(2.35a)\begin{equation}{c_{zz}} = {\varepsilon ^2},\quad {c_{yy}} = {c_{\theta \theta }} = 1,\quad {c_{yz}} = 0,\end{equation}

or, based on the components of the original conformation tensor as follows:

(2.35b)\begin{equation}{C_{zz}} = {C_{yy}} = {C_{\theta \theta }} = 1,\quad {C_{yz}} = 0.\end{equation}

Before proceeding with the derivation of some general results, we emphasize that with the method of solution used here, i.e. the lubrication approximation, the boundary conditions at the inlet and/or outlet of the hyperbolic section of the pipe require special attention. We will return to this issue with more details below, in §§ 4 and 7.

3.1. Pressure-drop

The average pressure-drop is defined as

(3.1)\begin{align}\mathrm{\Delta }{\varPi ^\ast }&:= 2{\rm \pi} \left( {\int_0^{{H^\ast }(0)} {{P^\ast }(\kern0.08em {y^\ast },0){y^\ast }\,\textrm{d}{y^\ast }} - \int_0^{{H^\ast }({\ell^\ast })} {{P^\ast }(\kern0.08em {y^\ast },1){y^\ast }\,\textrm{d}{y^\ast }} } \right) \nonumber\\ & ={-} 2{\rm \pi} \int_0^{{\ell ^\ast }} {\frac{\textrm{d}}{{\textrm{d}{z^\ast }}}\left( {\int_0^{{H^\ast }({z^\ast })} {{P^\ast }(\kern0.08em {y^\ast },z){y^\ast }\,\textrm{d}{y^\ast }} } \right)\textrm{d}{z^\ast }} ,\end{align}

where $\Delta {f^\ast }: = {f^\ast }({z^\ast } = 0) - f({z^\ast } = {\ell ^\ast })$ for any field variable ${f^\ast }$ (or $\Delta f: = f(z = 0) - f(z = 1)$ in dimensionless form). Using the characteristic scales reported in § 2.3, we find $\Delta \varPi \equiv \Delta {\varPi ^\ast }/\Delta \varPi _c^\ast $:

(3.2)\begin{align}\mathrm{\Delta }\varPi = 2\int_0^{H(0)} {P(y,0)y\,\textrm{d}y} - 2\int_0^{H(1)} {P(y,1)y\,\textrm{d}y} ={-} 2\int_0^1 {\frac{\textrm{d}}{{\textrm{d}z}}\left( {\int_0^{H(z)} {P(y,z)y\,\textrm{d}y} } \right)\textrm{d}z} ,\end{align}

where $\mathrm{\Delta }\varPi _c^\ast = (\eta _s^\ast + \eta _t^\ast ){\ell ^\ast }{Q^\ast }/(2h_0^{{\ast} 2})$. Applying Leibniz's rule for integrals in (3.2), and splitting the resulting expression, gives

(3.3)\begin{equation}\mathrm{\Delta }\varPi ={-} 2\int_0^1 {\left\{ {\int_0^H {\left( {\frac{{\partial P}}{{\partial z}}y{\kern 1pt} \,\textrm{d}y} \right)\textrm{d}z} } \right\}} - \int_0^1 {2H^{\prime}HP(H,z)\,\textrm{d}z} .\end{equation}

The second term on the right-hand side of (3.3) is simplified using $2H^{\prime}H = ({H^2} - 1)^{\prime}$ and performing integration by parts gives

(3.4)\begin{equation}\int_0^1 {2H^{\prime}HP(H,z)\,\textrm{d}z} = [({H^2}(z) - 1)P(H(z),z)]_0^1 - \int_0^1 {({H^2} - 1)\frac{\textrm{d}}{{\textrm{d}z}}P(H,z)\,\textrm{d}z} .\end{equation}

Since $H(0) = 1$ (due to dimensionalization) and $P(H(1),1) = 0$ (due to the reference pressure), the first term on the right-hand side of (3.4) vanishes. Thus, (3.3) and (3.4) yield

(3.5)\begin{equation}\mathrm{\Delta }\varPi ={-} 2\int_0^1 {\left( {\int_0^H {\frac{{\partial P}}{{\partial z}}y\,\textrm{d}y} } \right)\textrm{d}z} + \int_0^1 {({H^2} - 1)\frac{\textrm{d}}{{\textrm{d}z}}P(H,z)\,\textrm{d}z} .\end{equation}

Using the chain rule to find the derivative of $P(H(z),z)$ with respect to z, we find the general formula for the average pressure drop in the varying region of the pipe:

(3.6)\begin{equation}\mathrm{\Delta }\varPi ={-} 2\int_0^1 {\left( {\int_0^H {\left( {\frac{{\partial P}}{{\partial z}}y\,\textrm{d}y} \right)} } \right)\textrm{d}z} + \int_0^1 {({H^2} - 1)\left( {{{\left. {\frac{{\partial P}}{{\partial z}}} \right|}_{y = H}} + H^{\prime}{{\left. {\frac{{\partial P}}{{\partial y}}} \right|}_{y = H}}} \right)\textrm{d}z} .\end{equation}

Worth mentioning is that the second integral on the right-hand side of (3.6) results from the variation of the shape function with respect to the axial coordinate. All quantities in (3.6) can be obtained directly from the non-trivial components of the momentum balance, i.e. (2.22)–(2.23), which can also be used to reveal the individual contributions to the total average pressure drop.

3.2. First normal stress difference

The extensional character of the flow is investigated starting from the constraints at y = 0, given by (2.31b), along with continuity equation evaluated on the axis of symmetry. In this case, we find

(3.7)\begin{equation}\lim_{y \to 0} \frac{V}{y} = \lim_{y \to 0} \frac{{\partial V}}{{\partial y}} ={-} \frac{1}{2}{\left. {\frac{{\partial U}}{{\partial z}}} \right|_{y = 0}},\end{equation}

where the first equality is due to L'Hôpital's rule. Due to (3.7), the rate-of-deformation tensor is given by

(3.8)\begin{equation}{ {\dot{\boldsymbol{\gamma }}} |_{y = 0}} = \varepsilon {\left. {\frac{{\partial U}}{{\partial z}}} \right|_{y = 0}}(2{\boldsymbol{e}_z}{\boldsymbol{e}_z} - {\boldsymbol{e}_y}{\boldsymbol{e}_y} - {\boldsymbol{e}_\theta }{\boldsymbol{e}_\theta }).\end{equation}

Thus, the flow on the axis of symmetry is a pure uniaxial extensional flow and the quantity $\varepsilon { {\partial U/\partial z} |_{y = 0}}$ is a dimensionless rate of extension $\dot{E}$. We emphasize, however, that the flow in the hyperbolic pipe is a mixed flow with both extensional and shear characteristics. It is purely extensional in character along the axis of symmetry of the pipe only, while at the wall the shear characteristics dominate. When the rate-of-deformation tensor is given in the form of (3.8), the dimensionless extensional viscosity of the fluid can be determined merely from the first total normal stress difference ${N_1}: = ({T_{zz}} - {T_{yy}})/{\varepsilon ^2}$ which has been scaled by $\eta _0^\ast u_c^\ast {\ell ^\ast }/h_0^{{\ast} 2}$. Following the definition of the total stress tensor given in (2.19), the rate of deformation tensor given in (2.17) and the viscoelastic extra-stress tensor given in (2.18), we find

(3.9a)\begin{equation}{N_1} = 3(1 - \eta )\frac{{\partial U}}{{\partial z}} + \eta \left( {\frac{{{\tau_{zz}}}}{{{\varepsilon^2}}} - {\tau_{yy}}} \right)\quad \textrm{at}\;y = 0.\end{equation}

For a Newtonian fluid, i.e. for $De = 0$, (2.25) and (2.26) with the aid of (3.7) give ${\tau _{zz,N}} = 2{\varepsilon ^2}{ {\partial U/\partial z} |_{y = 0}}$ and ${\tau _{yy,N}} ={-} { {\partial U/\partial z} |_{y = 0}}$, and therefore (3.9a) reduces to ${N_{1,N}} = 3{ {\partial U/\partial z} |_{y = 0}}$. Equation (3.9a) can also be expressed in terms of the rescaled components of the conformation tensor (see (2.33b)) as follows:

(3.9b)\begin{equation}{N_1} = 3(1 - \eta )\frac{{\partial U}}{{\partial z}} + \frac{\eta }{{De}}\left( {\frac{{{c_{zz}}}}{{{\varepsilon^2}}} - {c_{yy}}} \right)\quad \textrm{at}\;y = 0,\end{equation}

or in terms of the original components of the conformation tensor,

(3.9c)\begin{equation}{N_1} = 3(1 - \eta )\frac{{\partial U}}{{\partial z}} + \frac{\eta }{{De}}({C_{zz}} - {C_{yy}})\quad \textrm{at}\;y = 0.\end{equation}

Considering the constraints on the axis of symmetry (see (2.31b,c)), the $zz$_-, $yy$- and $\theta \theta$-components of the Oldroyd-B model can be solved analytically. Denoting with $u(z) \equiv U(y = 0,z)$ and $u^{\prime}(z) \equiv { {\partial U/\partial z} |_{y = 0}}$, and taking into account (2.34) and (3.7), (2.10)–(2.12) are given in terms of the rescaled components of the conformation tensor (see (2.33a)),

(3.10)\begin{gather}{c_{zz}} + De(u{c^{\prime}_{zz}} - 2{c_{zz}}u^{\prime}) = {\varepsilon ^2},\end{gather}

(3.11)\begin{gather}{c_{yy}} + De(u{c^{\prime}_{yy}} + {c_{yy}}u^{\prime}) = 1,\end{gather}

(3.12)\begin{gather}{c_{\theta \theta }} + De(u{c^{\prime}_{\theta \theta }} + {c_{\theta \theta }}u^{\prime}) = 1.\end{gather}

The exact analytical solution of (3.10)–(3.12) is

(3.13)\begin{gather}{c_{zz}}(0,z) = {u^2}(z)\varphi (z)\left\{ {\frac{{{c_{zz}}(0,0)}}{{{u^2}(0)}} + \frac{{{\varepsilon^2}}}{{De}}\int_0^z {\frac{{\textrm{d}\kern 0.06em x}}{{\varphi (x){u^3}(x)}}} } \right\},\end{gather}

(3.14)\begin{gather}{c_{yy}}(0,z) = \frac{{\varphi (z)}}{{u(z)}}\left\{ {u(0){c_{yy}}(0,0) + \frac{1}{{De}}\int_0^z {\frac{{\textrm{d}\kern 0.06em x}}{{\varphi (x)}}} } \right\},\end{gather}

(3.15)\begin{gather}{c_{\theta \theta }}(0,z) = \frac{{\varphi (z)}}{{u(z)}}\left\{ {u(0){c_{\theta \theta }}(0,0) + \frac{1}{{De}}\int_0^z {\frac{{\textrm{d}\kern 0.06em x}}{{\varphi (x)}}} } \right\},\end{gather}

where

(3.16)\begin{equation}\varphi (z): = \exp \left( { - \frac{1}{{De}}\int_0^z {\frac{{\textrm{d}s}}{{u(s)}}} } \right).\end{equation}

From (2.34) and (3.13)–(3.16), the normal extra-stress components are trivially found as

(3.17a–c)\begin{align}{\tau _{zz}}(0,z) = \frac{{{c_{zz}}(0,z) - {\varepsilon ^2}}}{{De}},\quad {\tau _{yy}}(0,z) = \frac{{{c_{yy}}(0,z) - 1}}{{De}},\quad {\tau _{\theta \theta }}(0,z) = \frac{{{c_{\theta \theta }}(0,z) - 1}}{{De}}.\end{align}

Moreover, since ${C_{zz}} \equiv {c_{zz}}/{\varepsilon ^2}$, ${C_{yy}} \equiv {c_{yy}}$ and ${C_{\theta \theta }} \equiv {c_{\theta \theta }}$, the components of the original conformation tensor can easily be extracted too; worth mentioning is that the solution for ${C_{zz}}$, ${C_{yy}}$ and ${C_{\theta \theta }}$ at $y = 0$ depends on the aspect ratio only implicitly through the velocity profile $u = u(z)$. Equations (3.13)–(3.16) also show that the spatial evolution of ${C_{zz}}$, ${C_{yy}}$ and ${C_{\theta \theta }}$ at $y = 0$ depends on the conditions at the inlet of the pipe (as it should be expected) as well as on the fluid velocity at $y = 0$, $u = u(z)$; note that hereafter we will be referring to the inlet conditions as the ‘initial conditions’ (considered in a Lagrangian sense as such for the hyperbolic section of the flow). Noticing that $\int_0^z {(1/u(s))\,\textrm{d}s} > 0$, it is clear that $\varphi = \varphi (z)$ is a continuous and positive function with $\varphi (0) = 1$ which decreases exponentially to zero in the range $0 \le z \le 1$.

Finally, by plugging (3.13)–(3.17) in (3.9b), we find

(3.18a)\begin{align} {N_1} & = 3(1 - \eta )u^{\prime}(z) + \dfrac{\eta }{{D{e^2}}}\varphi (z)\left( {{u^2}(z)\int_0^z {\dfrac{{\textrm{d}\kern 0.06em x}}{{\varphi (x){u^3}(x)}}} - \dfrac{1}{{u(z)}}\int_0^z {\dfrac{{\textrm{d}\kern 0.06em x}}{{\varphi (x)}}} } \right)\nonumber\\ & \quad + \dfrac{\eta }{{De}}\varphi (z)\left( {\dfrac{{{u^2}(z){c_{zz}}(0,0)}}{{{u^2}(0){\varepsilon^2}}} - \dfrac{{u(0)}}{{u(z)}}{c_{yy}}(0,0)} \right). \end{align}

Equation (3.18a) is the general formula for the evaluation of the dimensionless first normal stress difference based on the flow field on the axis of symmetry of the pipe. Restoring the original components of the conformation tensor, ${c_{zz}} \equiv {\varepsilon ^2}{C_{zz}}$ and ${c_{yy}} \equiv {C_{yy}}$, and using that the polymer molecules are unstretched at the origin of the coordinate system, i.e. at the centre of the inlet of the varying region of the pipe, namely for ${c_{zz}}(0,0) = {\varepsilon ^2}$ and ${c_{yy}}(0,0) = 1$ (see (2.35a)), or equivalently, ${C_{zz}}(0,0) = {C_{yy}}(0,0) = 1$ (see (2.35b)), gives

(3.18b)\begin{align}{N_1} &= 3(1 - \eta )u^{\prime}(z) \nonumber\\ &\quad + \frac{{\eta \varphi (z)}}{{De}}\left( {\frac{1}{{De}}\left( {{u^2}(z)\int_0^z {\frac{{\textrm{d}\kern 0.06em x}}{{\varphi (x){u^3}(x)}}} - \frac{1}{{u(z)}}\int_0^z {\frac{{\textrm{d}\kern 0.06em x}}{{\varphi (x)}}} } \right) + \frac{{{u^2}(z)}}{{{u^2}(0)}} - \frac{{u(0)}}{{u(z)}}} \right).\end{align}

Equations (3.18a) and (3.18b) are exact, namely no approximations whatsoever have been made for their derivation. Also, notice that (3.18b) does not depend directly on the aspect ratio; instead, there is an indirect dependence of ${N_1}$ on ${\varepsilon ^2}$ through u. We emphasize that the equilibrium conditions imposed at $y = z = 0$ are a consequence of the exact analytical solution at the entrance region (see below in § 4.1). Finally, worth noting is that the Trouton ratio is evaluated at a postprocessing step after the solution of the governing equations has been evaluated. It is calculated from (3.18b) based on the (exact, numerical or approximate) solution for the velocity profile of the fluid at the centreline.

3.3. Viscoelastic extra-stresses at the wall

Considering the boundary conditions for the velocity field at the wall eliminates the material derivative of the polymeric extra stress(es) along the wall, namely $\textrm{D}f/\textrm{D}t = 0$ for $f = {\tau _{zz}},\;{\tau _{yz}},\;{\tau _{yy}},\;{\tau _{\theta \theta }}$ at $y = H(z)$. One can also confirm that along the wall, all velocity gradients can be expressed in terms of the wall shear rate $\dot{\gamma } \equiv { {\partial U/\partial y} |_{y = H}}$ only. This can be shown starting from the no-slip and no-penetration conditions (see (2.31a)), differentiating with respect to z and using the continuity equation (see (2.21)),

(3.19)\begin{equation}\frac{{\partial U}}{{\partial z}} ={-} H^{\prime}\dot{\gamma },\quad \frac{{\partial V}}{{\partial y}} = H^{\prime}\dot{\gamma },\quad \frac{{\partial V}}{{\partial z}} ={-} H^{\prime 2}\dot{\gamma }\quad \textrm{at}\;y = H.\end{equation}

Using (3.19), the exact analytical solution of the constitutive model, (2.25)–(2.28), is

(3.20)\begin{equation}\left. {\begin{array}{@{}c@{}} {{\tau_{zz}}(H,z) = 2De\,{{\dot{\gamma }}^2} + 2{\varepsilon^2}\dot{\gamma }H^{\prime}( - 1 + De\,H^{\prime}\dot{\gamma }),}\\ {{\tau_{yy}}(H,z) = 2H^{\prime}\dot{\gamma }(1 + De\,H^{\prime}\dot{\gamma }) + 2{\varepsilon^2}\,De\,{{\dot{\gamma }}^2}{{H^{\prime}}^4},}\\ {{\tau_{\theta \theta }}(H,z) = 0,}\\ {{\tau_{yz}}(H,z) = \dot{\gamma }(1 + 2De\,H^{\prime}\dot{\gamma }) + {\varepsilon^2}\dot{\gamma }{{H^{\prime}}^2}( - 1 + 2De\,\dot{\gamma }H^{\prime}).} \end{array}} \right\}\end{equation}

Equation (3.20) is well-defined for any value of $\dot{\gamma }$ and $De \ge 0$. Moreover, based on the positive definiteness criterion of the confirmation tensor ${\boldsymbol{\mathsf{C}}}$, its diagonal elements and its determinant, $\textrm{det}({\boldsymbol{\mathsf{C}}}) = {C_{\theta \theta }}({C_{yy}}{C_{zz}} - C_{yz}^2) = {c_{\theta \theta }}({c_{yy}}{c_{zz}} - c_{yz}^2)/{\varepsilon ^2}$, must be positive. Indeed, by using (2.34) and (3.20) we have confirmed that the diagonal elements of ${\boldsymbol{\mathsf{C}}}$ and its determinant

(3.21)\begin{equation}\textrm{det}({\boldsymbol{\mathsf{C}}}) = \frac{{{\varepsilon ^2} + {{(De\,\dot{\gamma })}^2}{{(1 + {\varepsilon ^2}{{H^{\prime}}^2})}^2}}}{{{\varepsilon ^2}}} = 1 + {(Wi\,\dot{\gamma })^2}{(1 + {\varepsilon ^2}H^{\prime 2})^2},\end{equation}

are strictly positive for any value of $De$, or $Wi$, and any geometry.

4.1. Solution at the entrance region

The pipe at the entrance region has a constant radius, i.e. $H(z) = 1$ for $z \le 0$. In this case, (4.1)–(4.4a–d) can easily be solved analytically for all the field variables,

(4.5)\begin{equation}\left. {\begin{array}{@{}c@{}} {{U_{en}} = 4(1 - {y^2}),\quad {V_{en}} = 0,\quad {P^{\prime}_{en}} ={-} 16,}\\ {{\tau_{zz,en}} = 128De\,{y^2},\quad {\tau_{yz,en}} ={-} 8y,\quad {\tau_{yy,en}} = {\tau_{\theta \theta ,en}} = 0,} \end{array}} \right\}\end{equation}

where the subscript ‘en’ is used to denote the entrance region. It is interesting to mention here that on the axis of symmetry of the pipe, $y = 0$, all the components of the polymer extra-stress tensor vanish, i.e. the polymer molecules are upstretched; the same holds for the planar case too (see (56) in Housiadas & Beris (Reference Housiadas and Beris2023)).

4.2. Newtonian solution for creeping flow

For a Newtonian fluid (De = 0) and creeping flow ($R{e_m} = 0$) (4.1)–(4.3) can be solved analytically too (Housiadas & Tsangaris Reference Housiadas and Tsangaris2022; Sialmas & Housiadas Reference Sialmas and Housiadas2024). The solution is

(4.6a–c)\begin{equation}{U_N} = \frac{4}{{{H^2}}}\left( {1 - \frac{{{y^2}}}{{{H^2}}}} \right),\quad {V_N} = \frac{{H^{\prime}}}{H}y{U_N},\quad {P^{\prime}_N} ={-} \frac{{16}}{{{H^4}}},\end{equation}

where a subscript ‘N’ throughout the paper denotes the solution for the Newtonian fluid. Using (4.6), one can find the polymer extra-stresses through (4.4a–d):

(4.7)\begin{equation}\left. {\begin{array}{@{}c@{}} {{\tau_{zz,N}} = 0,\quad {\tau_{yz,N}} = \dfrac{{\partial {U_N}}}{{\partial y}} ={-} \dfrac{{8y}}{{{H^4}}},}\\[12pt] {{\tau_{yy,N}} = 2\dfrac{{\partial {V_N}}}{{\partial y}} = 8\dfrac{{H^{\prime}}}{{{H^3}}}\left( {1 - 3\dfrac{{{y^2}}}{{{H^2}}}} \right),\quad {\tau_{\theta \theta ,N}} = 2\dfrac{{{V_N}}}{y} = 8\dfrac{{H^{\prime}}}{{{H^3}}}\left( {1 - \dfrac{{{y^2}}}{{{H^2}}}} \right).} \end{array}} \right\}\end{equation}

Along the axis of symmetry of the pipe, i.e. at $y = 0$, (4.7) gives ${\tau _{zz,N}} = {\tau _{yz,N}} = 0$ and ${\tau _{yy,N}} = {\tau _{\theta \theta ,N}} = 8H^{\prime}/{H^3}$; thus for a hyperbolic shape function ${\tau _{yy,N}} = {\tau _{\theta \theta ,N}} = {-} 4({\varLambda ^2} - 1)$.

4.3. Extra-stresses discontinuity

As mentioned above for a straight pipe, the velocity, pressure and extra-stress profiles given in (4.5) satisfy the lubrication equations (4.1)–(4.4a–d). However, the entrance extra-stress profiles given in (4.5) are not compatible with the exact solution at the wall given by (3.20). Indeed, from (4.5) we find that ${\lim _{z \to {0^ - }}}\dot{\gamma }(z) ={-} 8$ and

(4.8a–c)\begin{equation}\lim_{z \to {0^ - }} {\tau _{zz}}(1,z) = 128De,\quad \lim_{z \to {0^ - }} {\tau _{yy}}(1,z) = \lim_{z \to {0^ - }} {\tau _{\theta \theta }}(1,z) = 0,\quad \lim_{z \to {0^ - }} {\tau _{yz}}(1,z) ={-} 8.\end{equation}

Since ${\lim _{z \to {0^ + }}}\dot{\gamma }(z) ={-} 8$ and ${\lim _{z \to {0^ + }}}H^{\prime}(z) \ne 0$, from (3.20) we also find:

(4.9)\begin{equation}\left. {\begin{array}{@{}c@{}} \displaystyle{\lim_{z \to {0^ + }} {\tau_{zz}}(1,z) = 128De + \underline {16{\varepsilon^2}H^{\prime}(0)(1 + 8De\,H^{\prime}(0))} ,}\\[12pt] \displaystyle{\lim_{z \to {0^ + }} {\tau_{yy}}(1,z) ={-} 16H^{\prime}(0) + 128De\,H^{\prime}{{(0)}^2} + \underline {128{\varepsilon^2}\,De\,H^{\prime}{{(0)}^4}} ,}\\[12pt] \displaystyle{\lim_{z \to {0^ + }} {\tau_{\theta \theta }}(1,z) = 0,}\\[12pt] \displaystyle{\lim_{z \to {0^ + }} {\tau_{yz}}(1,z) ={-} 8 + 128De\,H^{\prime}(0) + \underline {8{\varepsilon^2}H^{\prime}{{(0)}^2}(1 + 16De\,H^{\prime}(0))} ,} \end{array}} \right\}\end{equation}

where in (4.9) the underlying $O({\varepsilon ^2})$ terms are disregarded at the classic lubrication limit, i.e. at ${\varepsilon ^2} \to 0$, and have been included here for completeness only. Equations (4.8) and (4.9) show that the viscoelastic extra-stress ${\tau _{yz}}$ and ${\tau _{yy}}$ at the inlet are discontinuous, and therefore a substantial rearrangement of the stresses is expected to take place when the fluid enters the varying (hyperbolic) section of the pipe.

4.4. Average pressure drop

The general formula for the average pressure-drop required to drive the flow, at the classic lubrication limit, is found from (3.6) by taking into account that the leading-order pressure field is independent of the transverse direction y:

(4.10)\begin{equation}\mathrm{\Delta }\varPi ={-} \int_0^1 {P^{\prime}(z){H^2}(z)\,\textrm{d}z} + \int_0^1 {({H^2}(z) - 1)P^{\prime}(z)\,\textrm{d}z} ={-} \int_0^1 {P^{\prime}(z)\,\textrm{d}z} .\end{equation}

Although (4.10) may seem trivial, it is a consequence of the classic lubrication approximation according to which the pressure gradient is independent of the radial coordinate. For Re_m = 0, we proceed by multiplying (4.2) with y, integrating with respect to y from $y = 0$ to $y = H(z)$ and applying the conditions on the axis of symmetry. Simplifying the result, gives

(4.11)\begin{equation}P^{\prime}(z) = \frac{2}{H}{\left( {(1 - \eta )\frac{{\partial U}}{{\partial y}} + \eta ({\tau_{yz}} - H^{\prime}{\tau_{zz}})} \right)_{y = H}} + \frac{{2\eta }}{{{H^2}}}\frac{\textrm{d}}{{\textrm{d}z}}\underbrace{{\left( {\int_0^H {{\tau_{zz}}y\,\textrm{d}y} } \right)}}_{{I(z)}},\end{equation}

where the quantity $I(z): = \int_0^{H(z)} {{\tau _{zz}}(y,z)y\,\textrm{d}y}$ has been defined for convenience. Equation (4.11) can be simplified further by considering (3.20) at the limit ${\varepsilon ^2} \to 0$, i.e. by omitting the O(ε ²) and O(ε ⁴) terms. In this case, ${({\tau _{yz}} - H^{\prime}{\tau _{zz}})_{y = H}} = \dot{\gamma }$, where $\dot{\gamma } = { {\partial U/\partial y} |_{y = H}}$ and thus, (4.10) and (4.11) give

(4.12)\begin{equation}\mathrm{\Delta }\varPi = \underbrace{{ - 2\int_0^1 {\frac{{\dot{\gamma }(z)}}{{H(z)}}\,\textrm{d}z} }}_{{{{\dot{\gamma }}_w}}} + \underbrace{{2\eta \left( {I(0) - {\varLambda ^2}I(1) - \int_0^1 {\frac{{2H^{\prime}(z)}}{{{H^3}(z)}}I(z)\,\textrm{d}z} } \right)}}_{\tau }.\end{equation}

Due to the definitions shown in (4.12), $\Delta \varPi = {\dot{\gamma }_w} + \tau$, namely the total average pressure drop along the pipe is a consequence of the tangential viscous stresses along the wall, ${\dot{\gamma }_w}$, plus an additional term, $\tau$, that comes from the axial viscoelastic extra-stress along the main flow direction. Furthermore, using (2.30) gives $2H^{\prime}(z)/{H^3}(z) ={-} ({\varLambda ^2} - 1)$ and thus (4.12) reduces to

(4.13)\begin{equation}\mathrm{\Delta }\varPi ={-} 2\int_0^1 {\frac{{\dot{\gamma }(z)}}{{H(z)}}\,\textrm{d}z} + 2\eta \left( {I(0) - {\varLambda ^2}I(1) + ({\varLambda ^2} - 1)\int_0^1 {I(z)\,\textrm{d}z} } \right).\end{equation}

For a Newtonian fluid, ${\tau _{zz}} = 0$ (see (4.4a)) and (4.12), or (4.13), gives

(4.14a)\begin{equation}\mathrm{\Delta }{\varPi _N} = 16\int_0^1 {\frac{{\textrm{d}z}}{{{H^4}(z)}}} .\end{equation}

For the hyperbolic pipe (4.14a) reduces to

(4.14b)\begin{equation}\mathrm{\Delta }{\varPi _N} = 16(1 + {\varLambda ^2} + {\varLambda ^4})/3.\end{equation}

It is trivial to confirm that for a straight pipe $\varLambda = 1$, or $H(z) = 1$, and (4.14a,b) give $\mathrm{\Delta }{\varPi _N} = 16$. Finally, we mention that (4.13) can also be derived from the total force balance on the entire hyperbolic section of the flow domain at the classic lubrication limit.

4.5. Mechanical energy

Here, we present the analysis for the total mechanical energy of the system at the classic lubrication limit, in terms of its individual contributions. For a planar hyperbolic channel, the same analysis has also been performed by Housiadas & Beris (Reference Housiadas and Beris2024c), including fluid slip along the walls of the channel. We start from the momentum equation for creeping flow, i.e. (4.2) with $R{e_m} = 0$, which we multiply with $yU$. Then, we integrate along the y-direction and perform integration by parts on the viscous and viscoelastic terms applying the no-slip and no-penetration conditions at the walls. We also use the constraints along the axis of symmetry, integrate the resulting equation with respect to the axial coordinate, and use Leibniz's rule for integrals in the viscoelastic terms, to find

(4.15)\begin{equation}\mathrm{\Delta }\varPi = {\varPhi _V} + {\varPhi _{el}} + {W_{el}},\end{equation}

where ${\varPhi _V}$ is the viscous dissipation, ${\varPhi _{el}}$ is the elastic bulk contribution, and ${W_{el}}$ is the work done by the elastic forces. Reiterating that $\Delta f = f(z = 0) - f(z = 1)$, these quantities are given as

(4.16a)\begin{gather}{\varPhi _V} = (1 - \eta )\int_0^1 {\int_0^H {{{\left( {\frac{{\partial U}}{{\partial y}}} \right)}^2}y\,\textrm{d}y\,\textrm{d}z} } ,\end{gather}

(4.16b)\begin{gather}{\varPhi _{el}} = \eta \int_0^1 {\int_0^H {\left( {{\tau_{zz}}\frac{{\partial U}}{{\partial z}} + {\tau_{yz}}\frac{{\partial U}}{{\partial y}}} \right)y\,\textrm{d}y\,\textrm{d}z} } ,\end{gather}

(4.16c)\begin{gather}{W_{el}} = \eta \mathrm{\Delta }\left( {\int_0^H {{\tau_{zz}}Uy\,\textrm{d}y} } \right).\end{gather}

Obviously, ${\varPhi _V}$ is strictly positive, while ${\varPhi _{el}}$ can theoretically take any value. It is positive, namely purely dissipative, only at the limit $De \to 0$ in which case $(1 - \eta ){\varPhi _{el}} = \eta {\varPhi _v}$. However, we need to emphasize that in the parameters range investigated here, ${\varPhi _{el}}$ is indeed dissipative (see also below in § 6.4).

Equation (4.15) is an alternative expression for $\mathrm{\Delta }\varPi $ to that given by (4.13). For a Newtonian fluid in a hyperbolic axisymmetric pipe, (4.16a–c) reduce to

(4.17a,b)\begin{equation}\frac{{{\varPhi _{v,N}}}}{{1 - \eta }} = \frac{{{\varPhi _{el,N}}}}{\eta } = \mathrm{\Delta }{\varPi _N},\quad {W_{el,N}} = 0.\end{equation}

Although in the Newtonian limit ${\tau _{zz,N}} = 0$ everywhere, as seen in (4.4a), resulting in ${W_{el,N}} = 0$, in the viscoelastic case, ${\tau _{zz}}$ is non-zero in the bulk except on the axis of symmetry of the pipe (${\tau _{zz}}(0,z) = 0$). Also note that ${\tau _{yz,N}} = \partial {U_N}/\partial y$ (see (4.4d)). Taken together the above explain the viscous dissipation and elastic contributions as given in (4.17a).

6.1. Streamlines

The stream function is given with of the aid of $\bar{y} \equiv y/H(z)$ and using the hyperbolic shape function for H, i.e. (2.30), and its derivatives $H^{\prime} ={-} ({\varLambda ^2} - 1){H^3}/2$ and $H^{\prime\prime} = 3{({\varLambda ^2} - 1)^2}{H^5}/4$. Thus, (6.8a) is simplified substantially because functions $\hat{\varPsi }_2^{(2j)}(z)$, $j = 1,2,3,4$, become constants, namely they turn out to be independent of z (see the Appendix A):

(6.9)\begin{equation}\varPsi = 2{\bar{y}^2} - {\bar{y}^4} + \eta De_m^2( - {\textstyle{1 \over 3}}{\bar{y}^2} + {\textstyle{1 \over 2}}{\bar{y}^4} - {\textstyle{1 \over 6}}{\bar{y}^8}) + O(De_m^3).\end{equation}

Recall that $\varPsi (\kern0.08em \bar{y} = 0) = 0$ (along the axis of symmetry) and $\varPsi (\kern0.08em \bar{y} = 1) = 1$ (at the wall) by construction. The shape of the streamlines, ${y_s} = {y_s}(z)$ or ${\bar{y}_s} \equiv {y_s}(z)/H(z)$, can be determined from (6.9) by considering the new transformed variable ${\bar{Y}_s}: = \bar{y}_s^2 = y_s^2(z)/{H^2}(z)$ and solving the equation

(6.10)\begin{equation}\varPsi \approx 2{\bar{Y}_s} - \bar{Y}_s^2 + \eta \,De_m^2( - {\textstyle{1 \over 3}}{\bar{Y}_s} + {\textstyle{1 \over 2}}\bar{Y}_s^2 - {\textstyle{1 \over 6}}\bar{Y}_s^4) = c\quad 0 \le c \le 1,\end{equation}

from which we can find ${\bar{Y}_s} = {\bar{Y}_s}(c,\eta \,De_m^2) \ge 0$. Notice that since z or $H(z)$ do not appear in (6.10), the z-dependence of the streamlines is hidden in ${\bar{Y}_s}$. Equation (6.10) has four roots, two of which are complex conjugate, the third is negative and the fourth is positive. Thus, the positive root is the relevant one for this problem, and although available analytically, is too long to be printed here. A very good approximation, however, is given by the following expression:

(6.11)\begin{equation}{\bar{Y}_s} \approx 1 - \sqrt {1 - c} + \frac{{\eta \,De_m^2}}{{12}}(4\sqrt {1 - c} - 4 + (4 - \sqrt {1 - c} )c).\end{equation}

Therefore, the streamlines are described by

(6.12)\begin{equation}{y_s}(z) \approx H(z)\sqrt {{{\bar{Y}}_s}} ,\end{equation}

Equation (6.12) shows that at the classic lubrication limit the shape of the streamlines is not affected by viscoelasticity. The streamlines are determined exclusively by the hyperbolic shape of the wall, while viscoelasticity influences only their position. In figure 3, we plot $\sqrt {{{\bar{Y}}_s}} \equiv {\bar{y}_s}$ as a function of c, for $\eta \,De_m^2 = 0$, 0.5 and 1. Both the full analytical solution (solid lines) and its approximation given by (6.11)–(6.12) (dashed lines) are shown in figure 3(a). One can hardly see the differences, and thus (6.11)–(6.12) are an excellent approximation of the exact solution, at least in the range $0 \le \eta \,De_m^2 \le 1$. Using the exact solution, or the approximate one, we can draw the shape of the streamlines, ${y_s}$, as function of the distance from the inlet, z, for various values of the stream function, c, as illustrated in figure 3(b).

Figure 3.

(a) The solution of (6.10) for a Newtonian fluid (De_m = 0), and the two viscoelastic cases ($\eta \,De_m^2 = 0.5\;\textrm{and}\;1$) as function of c. Both the exact (solid lines) and approximate ((6.11); dashed lines) solutions are shown. (b) The shape of the streamlines for a Newtonian fluid (solid lines) and a highly viscoelastic fluid with $\eta \,De_m^2 = 1$ (dashed lines) as function of z; the arrow shows in the direction of increasing c.

The results for the Newtonian ($\eta \,De_m^2 = 0$) and a highly viscoelastic fluid ($\eta \,De_m^2 = 1$) are shown for values of the steam function from zero to one in increments of 0.25, i.e. for $c = 0$, 0.25, 0.50, 0.75 and 1. As expected, the differences between the Newtonian and viscoelastic fluids are minute. The streamlines are slightly pushed towards the wall, but the displacement is inconsequential.

6.2. Average pressure-drop

Based on the series expansion in terms of De, the solution for the pressure gradient is

(6.13)\begin{equation}P^{\prime}(z) \approx {P^{\prime}_N}(z) + \eta \sum_{k = 1}^8 {{{\hat{P^{\prime}}}_k}(z)\,D{e^k}} ,\end{equation}

where $\eta {\hat{P^{\prime}}_k}(z) \equiv {P^{\prime}_k}(z)$, k = 1,2, …,8. The individual components are provided below up to $O(D{e^4})$, while the higher-order components are too long to be printed here,

(6.14a)\begin{gather}{P^{\prime}_N} ={-} \frac{{16}}{{{H^4}}},\end{gather}

(6.14b)\begin{gather}{\hat{P^{\prime}_1}} ={-} 256\frac{{H^{\prime}}}{{{H^7}}} = \frac{{128}}{3}{\left( {\frac{1}{{{H^6}}}} \right)^\prime },\end{gather}

(6.14c)\begin{gather}{\hat{P^{\prime}_2}} = 768\left( {\frac{{H^{\prime\prime}}}{{{H^9}}} - \frac{{19}}{3}\frac{{{{H^{\prime}}^2}}}{{{H^{10}}}}} \right),\end{gather}

(6.14d)\begin{gather}{\hat{P^{\prime}_3}} = 3072\left( { - \frac{{152}}{3}\frac{{{{H^{\prime}}^3}}}{{{H^{13}}}} + \frac{{296}}{{15}}\frac{{H^{\prime}H^{\prime\prime}}}{{{H^{12}}}} - \frac{{16}}{{15}}\frac{{16H^{\prime\prime\prime}}}{{15{H^{11}}}} + \frac{\eta }{5}\left( {52\frac{{{{H^{\prime}}^3}}}{{{H^{13}}}} - \frac{{59}}{3}\frac{{H^{\prime}H^{\prime\prime}}}{{{H^{12}}}} + \frac{{H^{\prime\prime\prime}}}{{{H^{11}}}}} \right)} \right),\end{gather}

(6.14e)\begin{gather} {{\hat{P^{\prime}_4}}} ={-} \dfrac{{840{{H^{\prime}}^4}}}{{{H^{16}}}} + \dfrac{{1586{{H^{\prime}}^2}H^{\prime\prime}}}{{3{H^{15}}}} - \dfrac{{85{{H^{\prime\prime}}^2}}}{{3{H^{14}}}} - \dfrac{{46H^{\prime}H^{\prime\prime\prime}}}{{{H^{14}}}} + \dfrac{{5{H^{(4 )}}}}{{3{H^{13}}}}\nonumber\\ \ \ \ \qquad\qquad\qquad + \,\eta \left( {\dfrac{{9221{{H^{\prime}}^4}}}{{20{H^{16}}}} - \dfrac{{2729{{H^{\prime}}^2}H^{\prime\prime}}}{{10{H^{15}}}} + } \right.\left. {\dfrac{{109{{H^{\prime\prime}}^2}}}{{8{H^{14}}}} + \dfrac{{133H^{\prime}H^{\prime\prime\prime}}}{{6{H^{14}}}} - \dfrac{{89{H^{(4 )}}}}{{120{H^{13}}}}} \right)\nonumber\\ \ \ \ \quad\qquad\qquad\qquad + \,{\eta ^2}\left( { - \dfrac{{136{{H^{\prime}}^4}}}{{{H^{16}}}} + \dfrac{{1178{{H^{\prime}}^2}H^{\prime\prime}}}{{15{H^{15}}}} - \dfrac{{11{{H^{\prime\prime}}^2}}}{{3{H^{14}}}} - \dfrac{{94H^{\prime}H^{\prime\prime\prime}}}{{15{H^{14}}}} + \dfrac{{{H^{(4 )}}}}{{5{H^{13}}}}} \right), \end{gather}

where the explicit dependence on the axial coordinate is omitted for simplicity. For a hyperbolic pipe, i.e. when the shape function is given by (2.30), (4.2) gives up to eighth-order

(6.15)\begin{align} \dfrac{{\mathrm{\Delta }\varPi }}{{\mathrm{\Delta }{\varPi _N}}} & = 1 - 2\eta \,D{e_m} + \dfrac{{5\eta }}{2}\,De_m^2 + \eta \left( { - \dfrac{{14}}{5} + \dfrac{{3\eta }}{5}} \right)De_m^3 + \eta \left( {3 - \dfrac{{49}}{{20}}\eta + \dfrac{4}{5}{\eta^2}} \right)De_m^4\nonumber\\ & \quad + \sum_{k = 5}^8 {\eta {\delta _k}(\eta )\,De_m^k} , \end{align}

where $\Delta {\varPi _N}$ is given in (4.14b), and ${\delta _k} = {\delta _k}(\eta )$, $k = 5,6,7,8$, are given by

(6.16a)\begin{gather}{\delta _5}(\eta ) ={-} \frac{{22}}{7} + \frac{{653}}{{105}}\eta - \frac{{506}}{{105}}{\eta ^2} + \frac{8}{7}{\eta ^3},\end{gather}

(6.16b)\begin{gather}{\delta _6}(\eta ) = \frac{{13}}{4} - \frac{{4399}}{{350}}\eta + \frac{{9049}}{{525}}{\eta ^2} - \frac{{4849}}{{525}}{\eta ^3} + \frac{{12}}{7}{\eta ^4},\end{gather}

(6.16c) \begin{gather}{\delta _7}(\eta ) ={-} \frac{{10}}{3} + \frac{{4649}}{{210}}\eta - \frac{{16\,609}}{{350}}{\eta ^2} + \frac{{270\,569}}{{6300}}{\eta ^3} - \frac{{5512}}{{315}}{\eta ^4} + \frac{8}{3}{\eta ^5},\end{gather}

(6.16d) \begin{gather}{\delta _8}(\eta ) = \frac{{17}}{5} - \frac{{6397}}{{180}}\eta \!+ \frac{{2\,095\,259}}{{18\,900}}{\eta ^2} - \frac{{56\,720\,593}}{{378\,000}}{\eta ^3} \!+ \frac{{791\,519}}{{7875}}{\eta ^4} - \frac{{259\,843}}{{7875}}{\eta ^5} \!+ \frac{{64}}{{15}}{\eta ^6}.\end{gather}

It is interesting, and somehow unexpected, that the final series for the reduced average pressure drop can be recast in terms of $D{e_m} \equiv 4({\varLambda ^2} - 1)De$, and the polymer viscosity ratio, η, only. A similar formula for $\Delta \varPi /\Delta {\varPi _N}$ holds for the planar case, studied theoretically previously by Boyko & Stone (Reference Boyko and Stone2022) and Housiadas & Beris (Reference Housiadas and Beris2023), albeit these authors did not report the corresponding formula; the latter can be found in Housiadas & Beris (Reference Housiadas and Beris2024a) though.

To get more insight about the results, we consider η = 4/10 and η = 1 and we calculate the coefficients of the $O(De_m^k)$ terms, $k = 0,1,2, \ldots ,8$, in (6.15); the results are reported in table 1. We notice that the numerical coefficients for η = 0.4 are smaller than those for η = 1, as well as that the magnitude of the coefficients decreases very slowly. Thus, and to increase the accuracy of our eighth-order formula, (6.15), we proceed by applying a technique that rearranges nonlinearly the terms of a truncated series. Specifically, we apply the diagonal [M/M] Padé approximants (Padé Reference Padé1892) on (6.15), where M = 1, 2, 3 and 4 and we derive new, transformed, analytical solutions. We remind the reader that the diagonal [M/M] Padé approximant of a function $f = f(\delta )$ agrees with the corresponding Taylor series of f about the point $\delta = 0$ up to $O({\delta ^{2M}})$. The successive approximants for M = 1, 2, 3 and 4 can be used to check the convergence of the transformed solutions.

Table 1.

The coefficients of the $O(De_m^k)$ terms, $k = 0,1,2,. \ldots ,8$, in (6.15).

In figure 4(a), we present the reduced pressure drop, $\Delta \varPi /\Delta {\varPi _N}$, for $\eta = 4/10$ up to second-, fourth-, sixth and eighth-order perturbation solutions as functions of $D{e_m}$. Convergence of the results is clearly observed gradually as more terms are included in the series, and we safely claim that the radius of convergence of the series is at least 0.85. The results show a decrease of $\Delta \varPi /\Delta {\varPi _N}$, as previously predicted for a hyperbolic symmetric channel (Boyko & Stone Reference Boyko and Stone2022; Housiadas & Beris Reference Housiadas and Beris2023, Reference Housiadas and Beris2024a). Furthermore, the accelerated (transformed) solutions which are depicted in figure 4(b), clearly show the convergence of the perturbation results. The convergence is achieved when at least five terms in the series are taken into account for the construction of the Padé [M/M] diagonal approximant (i.e. for M = 2). Indeed, the curves with M = 2, 3 and 4, which are constructed using the first five, seven and nine terms in the series (6.15), respectively, are indistinguishable.

Figure 4.

Reduced pressure drop versus De_m for the Oldroyd-B model with η = 0.4. (a) Perturbation solutions: solid (black) line, second-order solution; dashed (red) line, fourth-order solution; dotted (blue) line, sixth-order solution; dot–dashed (magenta) line, eighth-order solution. (b) Accelerated solutions: solid (black) line, up to second order; dashed (red) line, up to fourth order; dotted (blue) line, up to sixth order; dot–dashed (magenta) line, up to eighth order. The accelerated solutions up to fourth, sixth and eighth orders are indistinguishable.

In figure 5, we present the individual contributions to the average pressure drop, i.e. ${\dot{\gamma }_w}$ and $\tau$ (see (4.12)) normalized by the Newtonian value $\Delta {\varPi _N}$, as function of the modified Deborah number for a polymer viscosity ratio η = 4/10. First, we notice the clear convergence of both contributions to the average pressure drop. We also see that with increasing De_m, an increase of the magnitude of the positive viscous contribution at the wall (${\dot{\gamma }_w} > 0$) is predicted, as well as the increase of the magnitude of the negative viscoelastic contribution ($\tau < 0$). However, it appears that the negative viscoelastic contribution overwhelms the positive viscous contribution leading to a decrease in the average pressure drop. Similar results were observed for the planar case too (Housiadas & Beris, Reference Housiadas and Beris2023, Reference Housiadas and Beris2024a). The viscous and viscoelastic contributions are also given up to fourth order in De_m,

(6.17)\begin{align}\frac{{{{\dot{\gamma }}_w}}}{{\mathrm{\Delta }{\varPi _N}}} = 1 + \frac{\eta }{2}\,De_m^2 + \frac{\eta }{5}(3\eta - 4)\,De_m^3 + \eta \left( {1 + \frac{\eta }{{20}}(16\eta - 37)} \right)De_m^4 + O(De_m^5),\end{align}

and

(6.18)\begin{equation}\frac{\tau }{{\mathrm{\Delta }{\varPi _N}}} = 2\eta \,D{e_m}\left( { - 1 + D{e_m} - De_m^2 + \left( {1 - \frac{{3\eta }}{{10}}} \right)De_m^3} \right) + O(De_m^5).\end{equation}

Figure 5.

Decomposition of the average pressed drop, ΔΠ, based on (4.12), or (4.13), in (a) wall shear stress contribution, γ_w, (b) viscoelastic extra-stress contribution, τ, and (c) surface (S_el) and volume (V_el) viscoelastic contributions. All contributions are normalized by ΔΠ_N (4.14b) and the viscosity ratio is η = 4/10. Solid (black), dashed (red) and dotted (blue) lines are acceleration formulae up to O(De ⁴), O(De ⁶) and O(De ⁸), respectively.

In order to identify the terms responsible for the decrease of $\Delta \varPi$, it is helpful to further analyse the viscoelastic component, $\tau$, as defined in (4.12)–(4.13), in terms of a net surface contribution between the inlet and the outlet, ${S_{el}} \equiv I(0) - {\varLambda ^2}I(1)$, and a volume contribution, ${V_{el}} \equiv ({\varLambda ^2} - 1)\int_0^1 {I(z)\,\textrm{d}z}$, resulting in

(6.19)\begin{equation}\tau = {S_{el}} + {V_{el}}.\end{equation}

Based on the high-order analytical solution, we report these quantities up to $O(De_m^4)$,

(6.20)\begin{equation}\frac{{{S_{el}}}}{{\mathrm{\Delta }{\varPi _N}}} ={-} 3\eta \,D{e_m}\left( {1 - D{e_m} + De_m^2 - \left( {1 - \frac{{3\eta }}{{10}}} \right)De_m^3} \right) + O(De_m^5),\end{equation}

and

(6.21)\begin{equation}\frac{{{V_{el}}}}{{\mathrm{\Delta }{\varPi _N}}} = \eta \,D{e_m}\left( {1 - D{e_m} + De_m^2 - \left( {1 - \frac{{3\eta }}{{10}}} \right)De_m^3} \right) + O(De_m^5).\end{equation}

Note also, in reference to (4.13), that the analytical solution at the entrance region gives $I({0^ - }) = 32De$, while the high-order analytical solution in the hyperbolic section gives $I({0^ + }) = 32De + O(D{e^2})$.

The results for ${V_{el}}/\Delta {\varPi _N}$ and ${S_{el}}/\Delta {\varPi _N}$ are shown in figure 5(c) as function of $D{e_m}$. An interesting finding from our analysis that is worth mentioning is that the surface contribution is minus three times the positive volume contribution:

(6.22)\begin{equation}{S_{el}} ={-} 3{V_{el}}.\end{equation}

By verifying that (6.22) holds up to $O(De_m^8)$, we conjecture that is exact, although we can neither explain its origin, nor we have proved it formally. Equations (6.17), (6.20) and (6.21) reveal that both the pure viscous contribution along the wall of the tube and the volume viscoelastic contribution result to an increase of the pressure drop required to drive the flow, while the net surface contribution due to viscoelastic has the opposite effect, namely it facilitates the transport of the fluid through the tube.

6.3. First normal stress difference and Trouton ratio

As far as the exact solution for the Trouton ratio, given by (3.18b), is concerned we need to emphasize that the regular perturbation scheme in terms of De employed above cannot be applied directly on (3.18b). The reason is that the function $\varphi (z) = \exp \left( { - D{e^{ - 1}}\int_0^z {{u^{ - 1}}(x)\,\textrm{d}\kern 0.06em x} } \right)$ which appears in (3.18a,b) goes to zero faster than any power of the Deborah number, namely $\varphi = o(D{e^n})$ for any $n = 0,1,2,3, \ldots$ as $De \to {0^ + }$. Provided that the fluid velocity on the axis of symmetry, $u = u(z)$, is a continuous and strictly positive function for any $z \in [0,1]$, $\varphi$ is smooth and zero over an infinitely long interval, i.e. for $0 < De < \infty$, and yet non-zero, but it cannot be described by a Taylor series because it is not holomorphic.

However, and even if we are aware that the exact solution cannot be approximated by a simple power series in terms of De, we ignore this information and proceed by going one step back, i.e. we start from (3.9b) and apply the lubrication expansion,

(6.23)\begin{equation}{N_1} = 3(1 - \eta ){u^{\prime}_{(0)}} + \frac{\eta }{{De}}\left( {\frac{{{c_{zz(0)}}}}{{{\varepsilon^2}}} + {c_{zz(2)}} - {c_{yy(0)}}} \right) + O({\varepsilon ^2})\quad \textrm{at}\;y = 0,\end{equation}

where we retain the lubrication subscript to avoid confusion, ${u^{\prime}_{(0)}} \equiv { {\partial {U_{(0)}}/\partial z} |_{y = 0}}$ and the equations that govern ${c_{zz(0)}}(0,z)$, ${c_{zz(2)}}(0,z)$ and ${c_{yy(0)}}(0,z)$ are found from (3.10) and (3.11). For completeness, we give the leading-order velocity profile of the fluid on the axis of symmetry, ${u_{(0)}}$, which, interestingly, can be recast as follows:

(6.24a,b)\begin{equation}{u_{(0)}} = {u_{(0)N}}(z){q_{[n]}}(\eta ,D{e_m}),\quad {q_{[n]}} = 1 + {\hat{q}_2}(\eta )\,De_m^2 + \cdots + {\hat{q}_n}(\eta )\,De_m^n,\end{equation}

where ${u_{(0)N}} = 4/{H^2} = 4(1 + ({\varLambda ^2} - 1)z)$, and therefore ${u^{\prime}_{(0)}} = 4({\varLambda ^2} - 1){q_{[n]}}(\eta ,D{e_m})$; recall that a prime denotes differentiation with respect to the axial coordinate. The first three functions ${\hat{q}_k}$, $k = 2,3,4$ are

(6.25a–c)\begin{equation}{\hat{q}_2} ={-} \frac{\eta }{6},\quad {\hat{q}_3} = \frac{\eta }{{30}}(11 - 7\eta ),\quad {\hat{q}_4} ={-} \frac{\eta }{{300}}(170 - 283\eta + 104{\eta ^2}).\end{equation}

First, the initial value problem that determines ${c_{zz(0)}}(0,z)$ is

(6.26)\begin{equation}\left. {\begin{array}{@{}c@{}} {{c_{zz(0)}} + De({u_{(0)}}{{c^{\prime}}_{zz(0)}} - 2{c_{zz(0)}}{{u^{\prime}}_{(0)}}) = 0\quad 0 < z \le 1,}\\ {{c_{zz(0)}}(0,z = 0) = 0,} \end{array}} \right\}\end{equation}

where an unstretched inlet condition for the polymer molecules is used in accordance with the analytical solution at the entrance region, (4.5). The perturbation solution of (6.26) in terms of De (or De_m), up to any order, is

(6.27)\begin{equation}{c_{zz(0)}}(0,z) = 0\quad 0 \le z \le 1.\end{equation}

Using (6.27), we find that the initial value problem that determines ${c_{zz(2)}}(0,z)$ is

(6.28) \begin{equation}\left. {\begin{array}{@{}c@{}} {{c_{zz(2)}} + De({u_{(0)}}{c^{\prime}_{zz(2)}} - 2{c_{zz(2)}}{u^{\prime}_{(0)}}) = 1\quad 0 < z \le 1,}\\ {{c_{zz(2)}}(0,z = 0) = 1.} \end{array}} \right\}\end{equation}

Unexpectedly, we see that the solution for ${c_{zz(2)}}(0,z)$ can be found using only ${u_{(0)}}$. The solution of (6.28) up to fourth order in $D{e_m}$ is

(6.29)\begin{align}{c_{zz(2)}}(0,z) &= 1 + 2D{e_m} + 4De_m^2 + \left( {8 - \frac{\eta }{3}} \right)De_m^3 \nonumber\\ &\quad + \left( {16 - \frac{1}{{15}}\eta (9 + 7\eta )} \right)De_m^4 + O(De_m^5).\end{align}

Similarly, the initial value problem that determines ${c_{yy(0)}}(0,z)$ is

(6.30) \begin{equation}\left. {\begin{array}{@{}c@{}} {{c_{yy(0)}} + De({u_{(0)}}{c^{\prime}_{yy(0)}} + {c_{yy(0)}}{u^{\prime}_{(0)}}) = 1\quad 0 < z \le 1,}\\ {{c_{yy(0)}}(0,z = 0) = 1} \end{array}} \right\}. \end{equation}

The solution for ${c_{yy(0)}}(0,z)$ up to fourth order in $D{e_m}$ is

(6.31)\begin{equation}{c_{yy(0)}}(0,z) = 1 + De_m^2 + {\textstyle{1 \over 6}}(\eta - 6)\,De_m^3 + {\textstyle{1 \over {30}}}(30 + 7(\eta - 3)\eta )\,De_m^4 + O(De_m^5).\end{equation}

We reiterate, however, that ${c_{zz(0)}}$, ${c_{zz(2)}}$, ${c_{yy(0)}}$ and ${u^{\prime}_{(0)N}}$ are independent of z. Based on (6.23), the leading-order first normal stress difference is

(6.32)\begin{equation}{N_{1(0)}} = 3(1 - \eta ){u^{\prime}_{(0)}} + \frac{\eta }{{De}}({c_{zz(2)}} - {c_{yy(0)}}).\end{equation}

By substituting all known quantities (${c_{zz(2)}}$, ${c_{yy(0)}}$ and ${u^{\prime}_{(0)}}$) into (6.32) and simplifying the result, we find the Trouton ratio in the lubrication limit, $\textrm{Tr} \equiv {N_{1(0)}}/{u^{\prime}_{(0)}}$:

(6.33a)\begin{align}\textrm{Tr} &= 3\left( {1 + \eta \,D{e_m} + 3\eta \,De_m^2 + \frac{\eta }{6}(30 - \eta )\,De_m^3 + \frac{\eta }{{30}}(330 - 19\eta - 7{\eta^2})\,De_m^4} \right)\nonumber\\ &\quad + O(De_m^5).\end{align}

Instead, if the Newtonian strain rate on the symmetry axis is used, we form the quantity $\textrm{T}{\textrm{r}_\ast } \equiv {N_{1(0)}}/{u^{\prime}_{(0)N}}$ and find

(6.33b)\begin{align}\textrm{T}{\textrm{r}_\ast } &= 3\left( 1 + \eta \,D{e_m} + \frac{{17}}{6}\eta \,De_m^2 + \frac{\eta }{{30}}(161 - 17\eta )\,De_m^3\right.\nonumber\\ &\quad \left. +\, \frac{\eta }{{300}}(3130 + 53\eta - 244{\eta^2})\,De_m^4 \right) + O(De_m^5).\end{align}

As previously observed for $\Delta \varPi$ given by (6.15), the perturbation solutions for the Trouton ratio are expressed in terms of $D{e_m}$ and $\eta$ only. For the UCM model, i.e. for η = 1, we find up to $O(De_m^8)$,

(6.34a) \begin{align}\textrm{Tr} & \approx 3\left( {1 + D{e_m} + 3De_m^2 + \dfrac{{29}}{6}\,De_m^3 + \dfrac{{152}}{{15}}\,De_m^4 + \dfrac{{1933}}{{100}}\,De_m^5 + \dfrac{{34\,019}}{{900}}\,De_m^6} \right.\nonumber\\ & \quad \left. +\, \dfrac{{3\,235\,319}}{{44\,100}}\,De_m^7 + \dfrac{{200\,133}}{{1400}}\,De_m^8 \right), \end{align}

and

(6.34b) \begin{align} \textrm{T}{\textrm{r}_\ast } & \approx 3\left( {1 + D{e_m} + \dfrac{{17}}{6}\,De_m^2 + \dfrac{{24}}{5}\,De_m^3 + \dfrac{{2939}}{{300}}\,De_m^4 + \dfrac{{5647}}{{300}}\,De_m^5 + \dfrac{{1\,621\,259}}{{44\,100}}\,De_m^6} \right.\nonumber\\ & \quad + \left. {\dfrac{{6\,297\,799}}{{88\,200}}\,De_m^7 + \dfrac{{1\,103\,656\,411}}{{7\,938\,000}}De_m^8} \right). \end{align}

The first two terms in (6.34a) and (6.34b) are the same but at higher orders in $D{e_m}$ the coefficients are slightly different (as it should be expected). Also, in contrast to (6.15), all terms in (6.33a,b) and (6.34a,b) are positive which is indicative of the existence of a function irregularity at a critical value of $D{e_m}$. Unlike the coefficients of (6.15), the magnitude of the coefficients of $D{e_m}$ in (6.33a,b) and (6.34a,b) increases. Applying the lowest-order diagonal [M/M] Padé approximant (the case with M = 1) in the truncated series (6.33a) and (6.33b) shows that the approximants become singular at the critical value $D{e_{m,c}} = 1/3$ and $D{e_{m,c}} = 6/17$, respectively. Increasing the parameter M does not eliminate the singular point but pushes it towards 1/2 (for M = 2, one finds $D{e_{m,c}} \approx 0.507$ for both (6.33a) and (6.33b)). This shows the consistency of the approximants and implies that accurate predictions as the $D{e_m} = 1/2$ is approached cannot be expected.

In figure 6, we present the Trouton ratio as predicted by (6.33a) and increasing gradually the number of terms that are included in the series (shown with dashed lines), the corresponding Padé approximants with M = 2, 3 and 4 (shown with solid lines), along with the solution for the pure homogeneous steady uniaxial extension, $\textrm{T}{\textrm{r}_h}$, (dotted red line) as indicated in (5.4b). First, we observe that the Padé approximants are indistinguishable, as well as that they diverge near $D{e_m} \approx 1/2$. Second, we see that the homogeneous solution is practically the same as the transformed/accelerated solutions (i.e. with the Padé approximants for M = 2, 3 and 4). Third, we see that the perturbation solutions converge slowly and only for $D{e_m} < 1/2$. This information in conjunction with the convergence of the successive Padé approximants implies that the radius of convergence of the perturbation series for the Trouton ratio, (6.33a), is less than one half. Moreover, in figure 6(b) we compare the accelerated solution up to $O(De_m^8)$ with the exact solution based on the Newtonian velocity ${u_{(0)N}} \equiv {U_N}(0,z)$, i.e. (5.3a) or (5.3b), for different values of the contraction ratio $\varLambda $. It is seen that for $D{e_m} > 0.3$ approximately, the accelerated solution cannot follow closely the exact solution as it should be expected due to its divergence near $D{e_m} = 1/2$; recall that when $\textrm{Tr}$ is expressed in terms of $D{e_m}$, Λ does not appear explicitly in (6.33a), or in (6.33b).

Figure 6.

Trouton ratio along the cenreline at the exit of the pipe as function of De_m = 4($\varLambda $² − 1)De for η = 4/10. (a) The perturbation solutions are calculated up to O(De ²), O(De ⁴), O(De ⁶) and O(De ⁸) (dashed lines), while the accelerated solutions are constructed based on the perturbation solutions up to O(De ⁴) and O(De ⁸). (b) The Trouton ratio evaluated from the exact solution (3.18b), is shown with the velocity approximated by the Newtonian lubrication solution ((4.6) at y = 0), at the exit of the hyperbolic section of the pipe (z = 1) and for $\varLambda $ = 2 (solid), 4 (dashed) and 8 (dot–dashed). For comparison, the accelerated evaluation of the eighth-order perturbation solution is also shown and indicated with a solid (blue) line.

In figure 7, we investigate the sensitivity of the Trouton ratio, $\textrm{Tr} \equiv {N_1}/u^{\prime}(z)$ where ${N_1}$ is given by (3.18b), on the velocity profile on the axis of symmetry. Reiterating that (3.18b) is exact, we consider two cases, both of which have the same base velocity, i.e. the Newtonian velocity profile at the classic lubrication limit. In the first case, we consider higher-order corrections due to viscoelasticity, namely in terms of the Deborah number, and in the second we consider higher-order corrections due to variations in geometry, namely in terms of the square of the aspect ratio of the tube.

Figure 7.

The Trouton ratio at the exit of the pipe as function of De_m = 4($\varLambda $² − 1)De for η = 4/10, using (a) the Newtonian velocity profile ((4.6) at y = 0), i.e. (5.4b), solid (black) lines; the viscoelastic velocity profile up to O(De ⁸) with acceleration, dashed (red) lines, (b) the Newtonian velocity profile ((4.6) at y = 0), solid (black) line; the Newtonian velocity profile at y = 0 up to O(ε ⁴) ((6.39a)), dashed (red) line; the Newtonian velocity profile at y = 0 up to O(ε ⁴) with Padé acceleration (6.39b), dotted (blue) line.

In the first case, the starting point of the analysis is (6.24a), ${u_{(0)}} \approx {u_{(0)N}}(z){q_{[2M]}}(\eta ,D{e_m})$, and its corresponding transformed formula based on the Padé [M/M] approximant, ${u_{(0)}} \approx {u_{(0)N}}(z){q_{[M/M]}}(\eta ,D{e_m})$ where $M = 1,2,3,4$; recall that ${q_{[M/M]}}(\eta ,D{e_m})$ agrees with ${q_{[2M]}}(\eta ,D{e_m})$ up to $O(De_m^{2M})$. When these two approximate velocity profiles are substituted into (3.18b), gives

(6.35a)\begin{equation}\textrm{Tr} = 3(1 - \eta ) + \frac{\eta }{a}\left( {\frac{{1 - 2a{{(1 + ({\varLambda ^2} - 1)z)}^{2 - 1/.a}}}}{{1 - 2a}} - \frac{{1 + a{{(1 + ({\varLambda ^2} - 1)z)}^{ - 1 - 1/a}}}}{{1 + a}}} \right),\end{equation}

or, after some rearrangement to show there is no singularity as $a \to {0^ + }$,

(6.35b)\begin{align}\textrm{Tr} &= 3\left( {1 + \eta \frac{{(1 + 2a)a}}{{(1 - 2a)(1 + a)}}} \right) \nonumber\\ &\quad - \eta \left( {\frac{{{{(1 + ({\varLambda ^2} - 1)z)}^{2 - 1/a}}}}{{1 - 2a}} - \frac{{{{(1 + ({\varLambda ^2} - 1)z)}^{ - 1 - 1/.a}}}}{{1 + a}}} \right),\end{align}

where

(6.36a,b)\begin{equation}a \equiv D{e_m}{q_{[2M]}}(\eta ,D{e_m})\quad \textrm{or}\quad a \equiv D{e_m}{q_{[M/M]}}(\eta ,D{e_m}).\end{equation}

At the exit of the pipe, (6.35b) reduces to

(6.37)\begin{equation}\textrm{T}{\textrm{r}_{ex}} = 3\left( {1 + \eta \frac{{(1 + 2a)a}}{{(1 - 2a)(1 + a)}}} \right) - \eta \left( {\frac{{{\varLambda ^{2(2 - 1/a)}}}}{{1 - 2a}} - \frac{{{\varLambda ^{ - 2(1 + 1/a)}}}}{{1 + a}}} \right).\end{equation}

The first expression for a, (6.36a), is based on the truncated perturbation series up to $O(De_m^{2M})$, while the second expression, (6.36b), is based on the corresponding Padé [M/M] approximant. Notice that (6.35a,b) are identical with (5.3a,b), respectively, and (6.37) is identical with (5.4b), with the only difference that a is in place of $D{e_m}$. Obviously, a is a new dimensionless group which combines both viscoelastic (through $De$ and $\eta$) and geometric effects (through $\varLambda$) and therefore is a quantity that characterizes better and more accurately, compared with $D{e_m}$, the steady viscoelastic flow in the hyperbolic pipe. Notice, however, that due to (6.24a,b) and (6.25),

(6.38)\begin{equation}a = D{e_m} - \eta \,De_m^3/6 + O(\eta \,De_m^4),\end{equation}

namely the higher-order corrections are very small (for the typical parameters $D{e_m} = 1/2$ and $\eta = 4/10$, $a \approx 1/2 - 1/120 \approx 0.4916$). The Trouton ratio as function of $D{e_m}$ is depicted in figure 7(a) for η = 4/10 and $\varLambda $ = 2 and 3. Clearly, the differences between the results are inconsequential.

In the second case, we use the Newtonian velocity profile on the axis of symmetry, ${u_N}(z) \equiv {U_N}(y = 0,z)$, up to $O({\varepsilon ^4})$ as previously found by Housiadas & Tsangaris (Reference Housiadas and Tsangaris2022) and further studied recently up to $O({\varepsilon ^{20}})$ by Sialmas & Housiadas (Reference Sialmas and Housiadas2024), as well as the corresponding Padé [2/2] approximant:

(6.39a,b)\begin{equation}{u_N} = {u_{(0)}} + {\varepsilon ^2}{u_{(2)}} + {\varepsilon ^4}{u_{(4)}} + O({\varepsilon ^6}),\quad {u_{N[2 \times 2]}} = {u_{(0)}} + \frac{{{\varepsilon ^2}u_{(2)}^2}}{{{u_{(2)}} - {\varepsilon ^2}{u_{(4)}}}}.\end{equation}

For the hyperbolic pipe, the velocity components in (6.39a,b) are

(6.40a–c)\begin{equation}{u_{(0)}} = \frac{4}{{{H^2}}},\quad {u_{(2)}} = \frac{1}{3}{H^4}{({\varLambda ^2} - 1)^2},\quad {u_{(4)}} = \frac{1}{9}{H^{10}}{({\varLambda ^2} - 1)^4}.\end{equation}

Using the velocity profiles shown in (6.39a,b)–(6.40), and η = 4/10, $\varLambda $ = 3, we present the results for the Trouton ratio as function of the modified Deborah number in figure 7(b). The results are shown at the classic lubrication (${\varepsilon ^2} \to 0$) and for aspect ratio $\varepsilon = 0.15$; note that for $\varepsilon \le 0.1$ no observable differences are predicted. It is seen that the as the aspect ratio increases, the Trouton ratio drops, as well as that the differences between ${u_N}$ and the corresponding Padé approximant cause very small differences on the Trouton ratio too.

Based on the sensitivity analysis conducted here and the results presented in figures 6 and 7, we conclude that the exact analytical solution for the Trouton ratio, which is extracted with the aid of (3.18b) and is based on the velocity profile on the axis of symmetry of the pipe, is insensitive to small changes from the Newtonian velocity field. Therefore, the approximate analytical solution(s) at the classic lubrication limit, (5.3a,b) and (5.4a,b), can safely be used for viscoelastic fluids in axisymmetric hyperbolic pipes with $\varepsilon \le 0.1$. More comments and discussion on the theoretical features of the perturbation solution can be found in § 8.

6.4. Mechanical energy decomposition

In figure 8, we present the individual contributions on the pressure-drop resulting from the mechanical energy decomposition of the flow system as given in (4.15), i.e. the viscous dissipation, ${\varPhi _v}$, elastic dissipation, ${\varPhi _{el}}$, and the work done by the elastic forces, ${W_{el}}$, as defined in (4.16a,b,c), respectively. Based on the high-order perturbation scheme, ${\varPhi _v}$, ${\varPhi _{el}}$ and ${\varPhi _v}$ are

(6.41a)\begin{gather}\frac{{{\varPhi _V}}}{{(1 - \eta )\mathrm{\Delta }{\varPi _N}}} = 1 + \frac{{{\eta ^2}}}{{12}}\,De_m^4 + O(De_m^5),\end{gather}

(6.41b)\begin{gather}\frac{{{\varPhi _{el}}}}{{\eta \,\Delta {\varPi _N}}} = 1 - \frac{1}{2}\,De_m^2 + \frac{4}{5}\,De_m^3 + \left( { - 1 + \frac{{4\eta }}{{15}} + \frac{{{\eta^2}}}{{12}}} \right)De_m^4 + O(De_m^5),\end{gather}

(6.41c)\begin{gather}\frac{{{W_{el}}}}{{\eta \mathrm{\Delta }{\varPi _N}}} ={-} 2D{e_m} + 3De_m^2 + \frac{1}{5}( - 18 + 3\eta )\,\textrm{De}_m^3 + \left( {4 - \frac{{14\eta }}{5} + \frac{4}{5}{\eta^2}} \right)De_m^4 + O(De_m^5).\end{gather}

Figure 8.

Contributions to the pressure drop based on the total mechanical energy of the system. The results are shown as function of De_m for the viscous dissipation, Φ_v, elastic dissipation, Φ_el, and inlet/outlet work done by the elastic forces, W_el. Solid lines, η = 8/10; dashed lines, η = 4/10.

Equation (6.41a) shows that the energy lost due to viscous dissipation of the flow is affected at fourth order in $D{e_m}$, and thus ${\varPhi _v}$ remains almost constant with respect to $D{e_m}$ (for $\eta = 4/10$ ${\varPhi _v}$ increases from 1 at the Newtonian limit to 1.0133 at $D{e_m} = 1$), (6.41b) shows that the leading-order contribution to ${\varPhi _{el}}$ is positive with the first correcting term being $O(De_m^2)$ and negative, and (6.41c) shows that the leading-order term is negative. Note also that the asymptotic series for ${\varPhi _{el}}$ is alternating for all values of the polymer viscosity ratio, $0 < \eta \le 1$.

The high-order asymptotic solutions up to $O(De_m^8)$ for ${\varPhi _v}$, ${\varPhi _{el}}$ and ${\varPhi _v}$ have been processed with the Padé [M/M] approximant with M = 2, 3 and 4 to achieve convergence, and the results are shown as function of the modified Deborah number for polymer viscosity ratios 4/10 and 8/10. As previously revealed by (6.41a–c), both ${\varPhi _v}$ and ${\varPhi _{el}}$ are dissipative, with ${\varPhi _v}$ being practically constant in the range of $D{e_m}$ covered in the figure, and ${\varPhi _{el}}$ decreasing slightly with increasing $D{e_m}$. On the contrary, ${W_{el}}$ is negative and decreases substantially with the increase of the modified Deborah number leading to the decrease of the total pressure drop required to maintain the constant flowrate through the pipe.

Regarding our conjecture that for this type of flow ${\varPhi _{el}}$ is strictly positive, i.e. purely dissipative, we mention that the Padé [M/M] approximants for ${\varPhi _{el}}$ reveal that indeed, ${\varPhi _{el}} > 0$ for any value of $D{e_m}$ and $\eta$. For instance, for M = 2, we find

(6.42)\begin{equation}{\left. {\frac{{{\varPhi _{el}}}}{{\eta \mathrm{\Delta }{\varPi _N}}}} \right|_{[2 \times 2]}} = \frac{{1 + {\textstyle{8 \over 5}}\,D{e_m} + {\textstyle{1 \over {150}}}(9 + 80\eta + 25{\eta ^2})\,De_m^2}}{{1 + {\textstyle{8 \over 5}}\,D{e_m} + {\textstyle{1 \over {150}}}(84 + 80\eta + 25{\eta ^2})\,De_m^2}}.\end{equation}

Similarly, for M = 3 and 4 the approximants, which are too lengthy to print here, are strictly positive. Therefore, although a formal proof for the dissipative character of the elastic contribution to the system's total mechanical energy cannot be provided, the analytical accelerated formulae clearly indicate this.

6.5. A note on the purely elastic instabilities

It is known that flow instabilities occur in the motion of non-Newtonian polymeric liquids even in the absence of inertia ($Re \to 0$) (Shaqfeh Reference Shaqfeh1996). Pakdel & McKinley (Reference Pakdel and McKinley1996) developed a dimensionless criterion that characterizes the critical conditions for onset of elastic instabilities in two-dimensional, single-phase isothermal viscoelastic flows. The instabilities arise due to perturbations which cause the polymer molecules to stretch non-uniformly and deviate from their base state motion along a curvilinear streamline. Pakdel & McKinley (Reference Pakdel and McKinley1996) investigated in detail the lid-driven cavity flow problem of two different ideal elastic fluids, while McKinley et al. (Reference Mckinley, Pakdel and Oztekin1996) extended the analysis of Pakdel & McKinley (Reference Pakdel and McKinley1996) in order to describe the geometric and rheological sensitivity that has been observed in a variety of purely elastic instabilities.

Using the notation of the present work, the proposed dimensionless criterion of McKinley et al. (Reference Mckinley, Pakdel and Oztekin1996) is represented in the general form,

(6.43) \begin{equation}{M_{cr}} = \sqrt {\frac{{{\lambda ^\ast }u_c^\ast }}{{{\mathrm{\Re }^\mathrm{\ast }}}}\frac{{\tau _{tt}^\ast }}{{(\eta _s^\ast + \eta _p^\ast )(u_c^\ast{/}h_0^\ast )}}} \Leftrightarrow {M_{cr}} = \sqrt {\frac{1}{{\textrm{(}{\mathrm{\Re }^\mathrm{\ast }}/h_0^\ast )}}\frac{{{\lambda ^\ast }\tau _{tt}^\ast }}{{(\eta _s^\ast + \eta _p^\ast )}}} ,\end{equation}

where ${\mathrm{\Re }^\mathrm{\ast }}$ is a characteristic radius of curvature of the streamlines, and we have used that the characteristic value for the local deformation of the flow is $u_c^\ast{/}h_0^\ast $ . Also, $\tau _{tt}^\ast $ is the tensile extra-stress tangentially to the streamlines, i.e. $\tau _{tt}^\ast \equiv ({\boldsymbol{\tau }^\ast }\boldsymbol{\cdot }{\boldsymbol{t}_s})\boldsymbol{\cdot }{\boldsymbol{t}_s}$ with ${\boldsymbol{t}_s}$ being the tangential unit vector along a streamline. Both ${\Re ^ \ast }$ and ${\boldsymbol{t}_s}$ can be found provided that the shape of the streamlines is known. Indeed, following the analysis in § 6.1, and in particular (6.12) in dimensional form, i.e. $y_s^\ast (z) \approx {H^\ast }({z^\ast })\bar{Y}_s^{1/2}$ where ${\bar{Y}_s}$ is given by (6.11), we find

(6.44a,b)\begin{equation}{\boldsymbol{t}_s} \approx \frac{{{\boldsymbol{e}_z} + {H^\ast }^\prime \bar{Y}_s^{1/2}{\boldsymbol{e}_y}}}{{1 + {{({H^\ast }^\prime )}^2}{{\bar{Y}}_s}}},\quad {\Re ^ \ast } \approx \frac{{{{(1 + {{({H^\ast }^\prime )}^2}{{\bar{Y}}_s})}^{3/2}}}}{{H^{\ast \prime \prime }\bar{Y}_s^{1/2}}}\end{equation}

or

(6.45a,b)\begin{equation}{\boldsymbol{t}_s} \approx \frac{{{\boldsymbol{e}_z} + \varepsilon {H^\prime }\bar{Y}_s^{1/2}{\boldsymbol{e}_y}}}{{1 + {\varepsilon ^2}{H^\prime }^2{{\bar{Y}}_s}}},\quad \frac{{{\Re ^ \ast }}}{{h_0^\ast }} \approx \frac{{{{(1 + {\varepsilon ^2}{H^\prime }^2{{\bar{Y}}_s})}^{3/2}}}}{{{\varepsilon ^2}{H^{\prime \prime }}\bar{Y}_s^{1/2}}}.\end{equation}

Using the definition of $\tau _{tt}^\ast $, as well as (6.45a,b), and after some algebraic manipulations, we find (6.43) at the lubrication limit,

(6.46)\begin{equation}{M_{cr}} = \sqrt {\eta \,De\,\bar{Y}_s^{1/2}H^{\prime\prime}{\tau _{zz}}} ,\end{equation}

where the zero subscript has been dropped from ${\tau _{zz}}$ in accordance with the notation in § 4 and the present § 6. However, due to the hyperbolic shape function $H^{\prime\prime} = 3{({\varLambda ^2} - 1)^2}{H^5}/4$, while the leading-order term of ${\tau _{zz}}$ is given in (6.7) as ${\tau _{zz}} \approx 128\,De\,{\bar{y}^2}/{H^6}$. Substituting $H^{\prime\prime}$ and ${\tau _{zz}}$ in (6.46), using ${\bar{y}^2} \equiv {\bar{Y}_s}$ along a streamline and taking into account that $H \le 1$ gives

(6.47)\begin{equation}{M_{cr}} = D{e_m}\sqrt {\frac{{6\eta \bar{Y}_s^{3/2}}}{H}} \le \sqrt 6 \,\bar{Y}_s^{3/4}\sqrt {\eta \,De_m^2} .\end{equation}

Recall that ${\bar{Y}_s}$ varies in the range [0, 1], as well as that ${\bar{Y}_s} = {\bar{Y}_s}(\eta \,De_m^2)$. We emphasize that it is not possible to know in advance ${M_{cr}}$ for the steady-state flow studied here. This requires a linear stability analysis to be determined. However, if we consider Pakdel & McKinley's (Reference Pakdel and McKinley1996) heuristic criterion, we find that the elastic instability may appear for $\eta \,De_m^2 > \textrm{min}(M_{cr}^2/(6\bar{Y}_s^{3/2})) \approx 74.7$ where ${M_{cr}} \approx 21.17$; of course this limit is much beyond the range of parameters of interest which is covered in the present analysis ($0 \le De_m^2 \le 1$ and $0 \le \eta \le 1$).