Pressure-driven viscoelastic flow in axisymmetric geometries with application to the hyperbolic pipe

Abstract

We investigate theoretically the steady incompressible viscoelastic flow in a rigid axisymmetric cylindrical pipe with varying cross-section. We use the Oldroyd-B viscoelastic constitutive equation to model the fluid viscoelasticity. First, we derive exact general formulae: for the total average pressure-drop as a function of the wall shear rate and the viscoelastic axial normal extra-stress; for the viscoelastic extra-stress tensor and the Trouton ratio as functions of the fluid velocity on the axis of symmetry; and for the viscoelastic extra-stress tensor along the wall in terms of the shear rate at the wall. Then we exploit the classic lubrication approximation, valid for small values of the square of the aspect ratio of the pipe, to simplify the original governing equations. The final equations are solved analytically using a regular perturbation scheme in terms of the Deborah number, De, up to eighth order in De. For a hyperbolically shaped pipe, we reveal that the reduced pressure-drop and the Trouton ratio can be recast in terms of a modified Deborah number, De_m, and the polymer viscosity ratio, η, only. Furthermore, we enhance the convergence and accuracy of the eighth-order solutions by deriving transformed analytical formulae using Padé diagonal approximants. The results show the decrease of the pressure drop and the enhancement of the Trouton ratio with increasing De_m and/or increasing η. Comparison of the transformed solutions with numerical simulations of the lubrication equations using pseudospectral methods shows excellent agreement between the results, even for high values of De_m and all values of η, revealing the robustness, validity and efficiency of the theoretical methods and techniques developed in this work. Last, it is shown that the exact solution for the Trouton ratio gives a well-defined and finite solution for any value of De_m and reveals the reason for the failure of the corresponding high-order perturbation series for De_m > 1/2.

1. Introduction

Pressure-driven flows of viscoelastic fluids in narrow and long tubes (channel or pipes) are widely encountered in industrial processes, such as extrusion (Pearson 1985; Tadmor & Gogos 2013), in applications such as microfluidic extensional rheometers (Ober et al. 2013), in devices for subcutaneous drug administration (Allmendinger et al. 2014; Fischer et al. 2015) and many others. The complex rheological behaviour of viscoelastic fluids affects various features and properties of these flows, such as the average pressure drop along the tube, $\Delta {\varPi ^\ast }$, as function of the flow rate, ${Q^\ast }$, and the effective elongational viscosity of the fluid, $\eta _{el}^\ast $, as function of the extensional rate, ${\dot{E}^\ast }$; throughout the paper a star superscript denotes a dimensional quantity.

Contraction flows in internal confined geometries, and especially flows through hyperbolic geometries, have attracted much attention because of the almost constant extensional rate developed on the midplane for planar geometries or along the axis of symmetry for axisymmetric ones. These flows have been used to measure $\eta _{el}^\ast $ by relating $\Delta {\varPi ^\ast }$ to ${Q^\ast }$ (see e.g. James, Chandler & Armour 1990; Collier, Romanoschi & Petrovan 1998; Gotsis & Adriozola 1998; Koppol et al. 2009; Campo-Deaño et al. 2011; Wang & James 2011; Ober et al. 2013; Nyström et al. 2016; Kim et al. 2018). It should be emphasized that both the measurement of extensional viscosity (Binding & Walters 1988; Binding & Jones 1989; James & Walters 1993; Oliveira et al. 2007; Ober et al. 2013; Keshavarz & McKinley 2016) and its theoretical prediction, as for instance presented by the preliminary works of Cogswell (1972, 1978) and further improved by Binding (1988) and later by James (2016), are challenging because a steady and spatially uniform (i.e. homogeneous) flow from a Lagrangian point of view cannot be easily achieved (Petrie 2006). Other published works on the subject also include the experimental work of Lee & Muller (2017), the simulations of Nyström et al. (2012, 2016, 2017) and Feigl et al. (2003), and the optimization numerical work with respect to the geometry of Zografos et al. (2020), with goal to produce a constant strain-rate along the midplane of a symmetric channel. However, note that some authors have claimed that given that the flow in this geometry is not purely extensional, it is difficult to estimate the resistance to extensional motion directly from the $\Delta {\varPi ^\ast } - {Q^\ast }$ experimental data (James 2016; Nyström et al. 2016; Hsiao et al. 2017; James & Roos 2021).

Most of the available experimental data in the literature for viscoelastic fluids in contraction flows show an increase of the total pressure drop with increasing the fluid viscoelasticity. However, these data are mainly concerned with flows at high Weissenberg number, or in geometric configurations with singular points (sharp corners) such as the sudden contraction–expansion planar two-dimensional channels (Nigen & Walters 2002) or three-dimensional channels with square cross-section (Sousa et al. 2009), and the contraction–expansion circular axisymmetric tubes (Rothstein & McKinley 1999, 2001; Nigen & Walters 2002). In these cases, secondary motion, vortices and/or elastic instabilities arise which cause the increase of the average pressure drop in the tube although in the experiments conducted by Rothstein & McKinley (1999, 2001) in the sudden contraction–expansion geometry; it was mentioned that the increase of the pressure drop should not be directly connected with the onset of elastic instabilities. We also need to emphasize that those experimental data are in contrast with the numerical results obtained using macroscopic constitutive models, for instance those of Keiller (1993) using the Oldroyd-B and FENE models in planar and axisymmetric contractions; the work of Alves, Oliveira & Pinho (2003) using the Oldroyd-B and Phan-Thien–Tanner (PTT) models in planar contraction flow; the work of Aguayo, Tamaddon-Jahromi & Webster (2008) using the Oldroyd-B model in contraction and contraction/expansion geometries; and the work of Aboubacar, Matallah & Webster (2002) using the PTT model in planar contraction geometry. Boyko & Stone (2022) attributed the discrepancy between the experimental data and the numerical results to the inability of the macroscopic models to capture the complicated flow dynamics of viscoelastic fluids. These authors also speculated that a possible answer to this problem is the development of multiscale models, such as the mesoscopic numerical simulations with the bead–rod and bead–spring models performed by Koppol et al. (2009) in an abrupt contraction.

The abovementioned flows and geometries are different from the pure hyperbolic case in which no corners exist. In this geometry, the fluid enters the hyperbolic section of the axisymmetric tube smoothly in a fully developed state without vorticial or secondary motion. Two recent experimental papers by James & Roos (2021) and James & Tripathi (2023) provide, among other quantities, measurements of the pressure-drop versus the flow rate in this geometry; the former deals with the Boger-type, i.e. constant viscosity viscoelastic fluids, while the later deals with viscoelastic fluids with power-law viscous behaviour. The data showed that the pressure drop for the viscoelastic fluids remains practically the same with the corresponding Newtonian value, whereas the pressure drop increases with the flow rate for the power-law viscoelastic fluid.

For a very long time now, the theoretical study of viscoelastic flows has been proven to be a very demanding and complicated task. This is mainly because the relevant governing equations which predict the spatial and/or temporal evolution of the viscoelastic extra-stresses due to flow deformation are highly nonlinear. This holds true even when the most basic differential constitutive equation, i.e. the upper convected Maxwell (UCM) model, is used, and under the simplest possible flow conditions (steady, laminar and creeping conditions). Thus, the development and application of theories and techniques that further simplify the original governing equations is crucial. One such theory is the classic lubrication theory, a simple and efficient asymptotic technique widely used for modelling fluid lubricants (Szeri 2005; Tichy 2012) and thin fluid films (Ockendon & Ockendon 1995; Leal 2007; Langlois & Deville 2014), the motion of particles near surfaces (Goldman et al. 1967; Stone 2005), free surface flows (Ro & Homsy 1995) and flow in microchannels with known geometry (Plouraboué et al. 2004; Stone et al. 2004; Amyot & Plouraboué 2007). It is also commonly applied in the theoretical study of slow flows in narrow and confined geometries with small curvature. In these confined geometries, advances on the classic lubrication theory, referred to as high-order or extended lubrication theory, have been made recently by Tavakol et al. (2017), Housiadas & Tsangaris (2022, 2023) and Sialmas & Housiadas (2024).

The application of the lubrication approximation to study the flow of viscoelastic lubricants in confined tubes with solid walls was initiated and developed in the field of tribology by Tichy (1996), who was the first to derive the lubrication equations for the UCM model. Tichy's work was subsequently followed by others such as Zhang, Matar & Craster (2002) who investigated the dynamic spreading of a surfactant on a thin viscoelastic film using the same constitutive model; Li (2014) who investigated the effect of viscoelasticity on the lubrication performance in thin film flows using the PTT model; Gamaniel, Dini & Biancofiore (2021) who studied thin-film sliding lubricant flow in the presence of cavitation using the Oldroyd-B model; Ahmed & Biancofiore (2021, 2023) who studied lubricants viscoelastic flows using the UCM/Oldroyd-B and FENE type models; and Sari et al. (2024) who studied the effect of fluid viscoelasticity in soft lubrication using the Oldroyd-B model. Also, theoretical investigation of the effect of viscoelasticity, augmented by numerical simulations, on contracting channels with hyperbolic geometry was carried out by Pérez-Salas et al. (2019) using the PTT model; Boyko & Stone (2022) using the UCM/Oldroyd-B models; and Housiadas & Beris (2023, 2024a) who advanced the work of Boyko & Stone (2022) by deriving high-order asymptotic solutions (up to eighth order in the perturbation parameter) using the UCM/Oldroyd-B, PTT, Giesekus and FENE-P models. An extended list of the literature works for viscoelastic flows in internal and confined geometries (both planar and axisymmetric) has been compiled and can be found in the paper of Boyko & Stone (2022).

Although plenty of computational and experimental works are available in the literature for axisymmetric geometries (Rothstein & McKinley 1999, 2001; Nigen & Walters 2002; Oliveira et al. 2007; Sousa et al. 2009; James & Roos 2021; James & Tripathi 2023), a formal theoretical analysis for the axisymmetric hyperbolic case is still lacking. Additionally, while the axisymmetric geometry is similar to the planar symmetric case, experiments have shown differences (Nigen & Walters 2002; Rodd et al. 2005; Binding et al. 2006; James & Roos 2021). Moreover, the singular mathematical nature of the governing equations along the axis of symmetry of the tube requires special attention and can reveal interesting features of the solution. This is the motivation for the present work, in which we undertake a detailed theoretical analysis of the viscoelastic flow in an axisymmetric pipe with a variable (non-uniform) cross-section. The pipe radius is described with an appropriate smooth, continuous and adequately differentiable function which we call the ‘shape function’. Although the general analysis is performed in terms of this shape function, we focus on the hyperbolic contracting pipe given its importance to applications and experiments. An entrance and an exit region are also considered.

A major goal of our work is to develop a general theoretical framework for the evaluation of the dimensionless elongational viscosity (or Trouton ratio) of the fluid which is not linked to the relationship between $\Delta {\varPi ^\ast }$ and ${Q^\ast }$. This is achieved first by showing that the velocity field on the axis of symmetry of the pipe corresponds to pure uniaxial extension, and then by solving exactly the constitutive model on the axis of symmetry. The analytical solution allows for the development of a general formula for the Trouton ratio in terms of the velocity field only.

Furthermore, we aim to derive approximate analytical solutions for the flow field in this geometry. To this end, we exploit the classic lubrication theory to derive a simplified set of governing equations, valid at the limit of a vanishing small aspect ratio of the pipe; the latter is defined as the ratio of the radius of the cross-section at the inlet to the length of the pipe. Furthermore, we investigate the validity and accuracy of the approximate solutions, and we study the effect of the rheological parameters entering in the constitutive model. All these are of fundamental importance for comparison with experimental data and for building reliable mathematical model(s) with good predictive capabilities.

In this work, we restrict the analysis to the UCM/Oldroyd-B models, since previous works for the planar hyperbolic geometry (Housiadas & Beris 2023, 2024a,b,c) showed that for weakly viscoelastic fluids the results are insensitive to the nonlinear extensions of the UCM/Oldroyd-B models involved in more realistic macroscopic constitutive models (such as the Giesekus, PTT and FENE-P models). This is attributed mainly to the lubrication approximation within which most of the nonlinear terms of the polymer extra-stress tensor in the constitutive equations are eliminated. Moreover, additional terms are eliminated at the midplane of the channel or the axis of symmetry of the pipe because of the symmetries and extensional character of the flow, reducing the equations for all models there to those obtained for the UCM/Oldroyd models. For more details the interested reader is referred to Housiadas & Beris (2023, 2024a,b,c).

The range of validity and accuracy of the high-order asymptotic expansion in terms of the Deborah number is investigated in three ways (see § 2 for the definition of the Deborah number). First, by employing a method that accelerates the convergence of series such as the diagonal Padé approximants (Padé 1892); this method has been developed and extensively studied in the literature (see, for instance, Baker & Graves-Morris (1981, 1996)). The performance of Padé approximants, as well as a few other methods that accelerate the convergence of series, for pure shear and pure elongational viscoelastic flows has been investigated systematically and presented by Housiadas (2017, 2023). These methods have also been applied successfully in a variety of more complex viscoelastic problems (Housiadas 2021, 2023; Housiadas & Beris 2023, 2024c). Second, by developing and using highly accurate numerical methods for the solution of the lubrication equations, and third by comparing the high-order asymptotic results for the Trouton ratio, and the corresponding transformed formulae, with the exact analytical solution which is valid beyond the classic lubrication limit.

The rest of the paper is organized as follows. In § 2, we describe the flow geometry, the main dimensionless quantities, the governing equations and accompanied auxiliary conditions for the steady flow of an incompressible viscoelastic fluid in an axisymmetric non-uniform pipe. In § 3, we develop the general theoretical framework leading to general analytical formulae for the average pressure drop and the dimensionless elongational viscosity of the fluid (usually known as the Trouton ratio). In § 4 we derive a simplified set of governing equations following the lubrication approximation. This is followed in § 5 by the development of analytical expressions for the first normal stress difference and the Trouton ratio at the centreline, using an exact expression that involves only the flow velocity and its approximation by the closed-form analytical expression for the Newtonian lubrication solution. The method of developing high-order regular perturbation solutions of the classic lubrication equations in terms of the Deborah number is described in § 6. There, we report the most important new asymptotic formulae along with suitable convergence acceleration results for the stream function, the total average pressure-drop, the Trouton ratio, the total force balance and its individual contributions, and the mechanical energy decomposition of the flow system. In the same section we also demonstrate through detailed comparisons the accuracy of the Newtonian velocity approximation entering the evaluation of the Trouton ratio based in the exact formula developed before, as well as its extended validity with respect to the Deborah number in comparison with more standard asymptotic evaluations. In § 7, we present the numerical methodology developed to solve the final lubrication equations, along with a comparison of the numerical results with the formulae derived in § 6. In § 8 we discuss an important finding of the analysis, namely a non-algebraic dependence of the solution for the Trouton ratio with respect to the Deborah number, and finally, in § 9 we provide our conclusions.

2. Problem formulation

The geometry and its main features are shown in figure 1. A varying axisymmetric pipe which consists of three segments is depicted; an entrance region with constant radius $h_0^\ast $, a varying region with radius $h_0^\ast $ at the inlet and radius $h_f^\ast $ at the outlet and an exit region with constant radius $h_f^\ast $; the length of varying section of the pipe is ${\ell ^\ast }$. We also assume a constant volumetric flow rate, ${Q^\ast }$. The steady state flow in the entrance region of the tube is fully developed. Given the smooth transition to the hyperbolic region and the absence of inertia, the flow in the entrance region is unaltered from the steady state uniform pipe flow (Poiseuille). This is very different from, for instance, the inlet conditions of the sudden contraction–expansion geometry of Rothstein & McKinley (1999, 2001) according to which major vortices appear on the corner region right before the contraction section of the axisymmetric tube, or the work of Lubanssky et al. (2007) for the viscoelastic flow in an abrupt axisymmetric contraction.

Figure 1. Geometrical configuration and cylindrical coordinate system (y*, z*) for an axisymmetric hyperbolic pipe.

For later convenience, and to facilitate the discussion and analysis, we define the aspect ratio of the varying region of the pipe, as well as one half of the average inlet velocity,

(2.1a,b)\begin{equation}\varepsilon \equiv \frac{{h_0^\ast }}{{{\ell ^\ast }}},\quad u_c^\ast \equiv \frac{{{Q^\ast }}}{{2{\rm \pi} h_0^{{\ast} 2}}}.\end{equation}

For narrow and confined geometries such as that shown in figure 1 the aspect ratio is small, $\varepsilon < 1$, a feature that will be exploited in § 4 to simplify the original governing equations. We consider a complex viscoelastic fluid (a polymeric material into a Newtonian solvent) with longest relaxation time ${\lambda ^\ast }$, and we recognize two distinct time scales associated with the flow conditions for the process. The first is an inverse shear-rate $h_0^\ast{/}u_c^\ast $, and the second is an average residence time of the fluid in the pipe ${\ell ^\ast }/u_c^\ast $. Based on these three characteristic times, we define two dimensionless groups of importance: the Weissenberg and Deborah numbers,

(2.2a,b)\begin{equation}Wi \equiv \frac{{{\lambda ^\ast }}}{{h_0^\ast{/}u_c^\ast }} = \frac{{{\lambda ^\ast }u_c^\ast }}{{h_0^\ast }},\quad De \equiv \frac{{{\lambda ^\ast }}}{{{\ell ^\ast }/u_c^\ast }} = \frac{{{\lambda ^\ast }u_c^\ast }}{{{\ell ^\ast }}},\end{equation}

from where we observe that $De = \varepsilon \,Wi$. Since the analysis performed here assumes a confined and narrow geometry, $0 < \varepsilon < 1$, and in order to be able to observe the effect of viscoelasticity, we demand that the Deborah number is finite when the aspect ratio goes to zero, which in turn implies that the Weissenberg number is large, i.e. $De = O(1)$ and $Wi = O(1/\varepsilon )$ as $\varepsilon \to 0$. Note that if we had assumed $Wi = O(1)$ as $\varepsilon \to 0$, then $De \to 0$ and either the relaxation time would be very small or the residence time would be very large. In the first case, the effect of viscoelasticity would be negligible, while in the second its effect would be observed in a very small region near the entrance of the pipe. Hence, we proceed with $De = O(1)$ and $Wi = O(1/\varepsilon )$ as $\varepsilon \to 0$, and we keep in mind that a small value of $De$ indicates that the relaxation time of the polymeric material is small compared with the residence time (i.e. the flow is ‘slow’), while a large value of $De$ suggests that the polymer molecules do not have enough time to relax before their exit from the tube (i.e. the flow is ‘fast’). Note that very large values of the Deborah number are not permitted, as they correspond to an elastic solid where flow is not possible, and the Oldroyd-B model is not valid. In this work, we are primarily interested in creeping flow and small values of the Deborah number.

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. 2017; Housiadas & Tsagaris 2022, 2023; Housiadas & Beris 2023, 2024a,b,c). 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 (2021) and James & Tripathi (2023) used $\varLambda = \sqrt {117/5.8} \approx 4.49$, Kamerkar & Edwards (2007) used Hencky strain values in the range $4 \le {\varepsilon _H} \le 7$ and Collier et al. (1998) 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 1987; Beris & Edwards 1994) 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 (2022), Ahmed & Biancofiore (2021, 2023) and Housiadas & Beris (2023) 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. General analytical results

Main quantities of interest for this type of flow are the average pressure-drop required to maintain the constant flowrate through the pipe, the first normal stress difference and the Trouton ratio (or dimensionless elongational viscosity) of the fluid, as well as the viscoelastic extra-stresses at the wall.

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. Lubrication approximation

The evaluation of the total average pressure drop from (3.6), the first normal stress difference from (3.18b), and the viscoelastic extra-stresses at the wall from (3.20), require the velocity field in the flow domain $\varOmega$. We will not proceed by solving the full nonlinear governing equations reported in § 2, but we will find good approximations of the field variables. This can be achieved by considering pipes with small aspect ratio, i.e. $\varepsilon < 1$. From an asymptotic point of view and based on the magnitude of ε, all terms in (2.21)–(2.23) and the constitutive equation, (2.25)–(2.28), multiplied with ${\varepsilon ^2}$ or ${\varepsilon ^4}$ are much smaller compared with the other terms and can be ignored as a first approximation. Therefore, the solution for each dependent field variable X can be expressed formally as a perturbation power series in terms of ${\varepsilon ^2}$, $X(y,z) \approx \sum\nolimits_{i \ge 0} {{X_{(2i)}}(y,z){\varepsilon ^{2i}}}$ as ${\varepsilon ^2} \to 0$. Dropping the zero subscript in parenthesis throughout the paper for simplicity, unless otherwise noted, the leading-order balance equations are

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

(4.2)\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}}} \right) + \eta \left( {\frac{{\partial {\tau_{zz}}}}{{\partial z}} + \frac{{\partial {\tau_{yz}}}}{{\partial y}} + \frac{{{\tau_{yz}}}}{y}} \right),\end{gather}

(4.3)\begin{gather}0 ={-} \frac{{\partial P}}{{\partial y}},\end{gather}

(4.4a)\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) = 0,\end{gather}

(4.4b)\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}

(4.4c)\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}

(4.4d)\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}},\end{gather}

where the convective derivative $\textrm{D}/\textrm{D}t$ is given by (2.24). The same type of lubrication analysis has been performed before for the planar case by Boyko & Stone (2022), Ahmed & Biancofiore (2021, 2023) and Housiadas & Beris (2023, 2024a,b,c). However, as it will be become clear later, the analysis for the cylindrical pipe requires special attention because from the mathematical point of view, the governing equations at the axis of symmetry of the pipe although well-defined are singular in cylindrical coordinates, namely are not strictly defined on y = 0, but they should be considered only as the limit as y goes to zero. This implies that the conditions (2.31b) and (2.31c) must hold.

Equation (4.3) shows that $P = P(z)$ only, and thus $\partial P/\partial z \equiv P^{\prime}(z)$ is an unknown function that must be determined; this is a major feature of the classic lubrication theory at the limit of a vanishing small aspect ratio. It is also interesting that the corresponding equation for ${\tau _{\theta \theta }}$, (4.4c), is decoupled from the remaining equations. Thus, one can solve (4.1), (4.2) and (4.4a,b,d) to find the velocity components U and V, the pressure gradient $P^{\prime}(z)$ and the viscoelastic extra-stresses ${\tau _{zz}},\;{\tau _{yy}}$ and ${\tau _{yz}}$, and then (4.4c) to find ${\tau _{\theta \theta }}$ for completeness. This decoupling is a consequence of the lubrication approximation according to which ${\tau _{\theta \theta }}$ in the momentum balance is neglected, and the axisymmetry of the flow field.

We proceed by reporting the solution of the lubrication equations at the entrance region of the pipe, the solution of the lubrication equations in the varying region of the pipe for a Newtonian fluid at creeping flow conditions, and by providing information about the discontinuity of the viscoelastic extra-stress when the polymeric fluid exits the entrance region and enters the varying region. Most importantly, starting from the general expressions derived in § 3, we derive the simplified formulae for the total average pressure drop, the first normal stress difference and the Trouton ratio at the classic lubrication limit. We also use (4.1)–(4.3) to derive the individual contributions to the total mechanical energy of the system; this leads to an alternative decomposition of the average pressure drop in the pipe.

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 (2023)).

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 2022; Sialmas & Housiadas 2024). 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 (2024c), 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).

5. Trouton ratio based on the exact, velocity-based formula (3.18b)

Using the analytical solution for the extra-stress tensor at the entrance region, (4.5), i.e. the information that the polymer molecules are unstretched at the origin of the coordinate system ($\kern0.08em y = z = 0$), the general formula for the first normal stress difference at the lubrication limit is derived directly from (3.18b) because as mentioned in § 3, ${\varepsilon ^2}$ does not appear explicitly in (3.18b). Thus, ${N_1}$ can be found provided that the velocity profile on the axis of symmetry, $u = u(z)$, is known. Of course, this requires the solution of the full governing equations albeit approximation(s) of u can be utilized too. For instance, for a Newtonian fluid ($De = 0$) and creeping flow ($R{e_m} = 0$), and using (4.6) to find the velocity on the axis of symmetry, ${u_N}(z) \equiv {U_N}(y = 0,z) = 4/{H^2}(z)$ and ${u^{\prime}_N}(z) ={-} 8H^{\prime}(z)/{H^3}(z)$, reduces (3.18b) as follows:

(5.1)\begin{equation}{N_{1,N}} = 3{u^{\prime}_N}(z) ={-} 24H^{\prime}(z)/{H^3}(z).\end{equation}

The approximation of the velocity profile with the corresponding Newtonian one is a common approach in the literature for hyperbolic geometries and many fluids with different rheological behaviour such as viscoelastic Boger-type fluids (James & Roos 2021), inelastic shear thinning-fluids (Ober et al. 2013; Pérez-Salas et al. 2019) or viscoelastic fluids with power-law viscous behaviour (James & Tripathi 2023). Although the first normal stress difference shown in (5.1) appears to depend on the axial coordinate, due to the dependence of H on z, for the hyperbolic pipe described by (2.30), $H^{\prime}/{H^3} ={-} ({\varLambda ^2} - 1)/2$ yielding a simple formula for ${N_{1,N}}$ which depends on the geometric ratio $\varLambda \equiv h_0^\ast{/}h_f^\ast $ only, i.e. ${N_{1,N}} = 12({\varLambda ^2} - 1)$. Therefore, for $\varLambda \ne 1$, the Trouton ratio, defined as $\textrm{Tr} \equiv {N_1}/u^{\prime}(z) = {N_1}/4({\varLambda ^2} - 1)$, gives the well-established result for a Newtonian fluid under homogeneous uniaxial extension,

(5.2)\begin{equation}\textrm{T}{\textrm{r}_N} = 3,\quad \varLambda \ne 1.\end{equation}

For the viscoelastic case, one can use the Newtonian velocity profile in (3.18b), to find an approximation of the Trouton ratio; this merely corresponds to the uncoupled case ($\eta = 0$) at the lubrication limit (${\varepsilon ^2} \to 0$). Assuming the hyperbolic pipe, i.e. considering (2.30) for the shape function and substituting in (3.18b) yields $\textrm{Tr} = \textrm{Tr}(z,\eta ,\varLambda ,D{e_m})$,

(5.3a) \begin{align}\textrm{Tr} &= 3(1 - \eta ) \nonumber\\ &\quad + \frac{\eta }{{D{e_m}}}\!\left(\!{\underbrace{{\frac{{1 - 2D{e_m}{{(1 + ({\varLambda ^2} - 1)z)}^{2 - 1/D{e_m}}}}}{{1 - 2D{e_m}}}}}_{{{C_{zz}}(0,z)}} - \underbrace{{\frac{{1 + D{e_m}(1 + ({\varLambda ^2} - 1)z{)^{ - 1 - 1/D{e_m}}}}}{{1 + D{e_m}}}}}_{{{C_{yy}}(0,z)}}}\!\right)\!,\end{align}

or, after some rearrangement to show that the formula is not singular as $D{e_m} \to 0$, we find

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

where, assuming hereafter $\varLambda \ne 1$, the modified Deborah number $D{e_m} \equiv 4De({\varLambda ^2} - 1)$ has been used. Here $D{e_m}$ is a new dimensional group which combines both viscoelastic and geometric effects and plays an important role in the solution of the major quantities of interest. For a contracting pipe, $\varLambda$ is larger than one and thus $D{e_m}$ is positive. Note that the first term in parenthesis on the right-hand side of (5.3a) is merely ${C_{zz}}(0,z)$, while the second is ${C_{yy}}(0,z) = {C_{\theta \theta }}(0,z)$; for completeness, we also report that $\varphi (z) = {(1 + ({\varLambda ^2} - 1)z)^{ - 1/D{e_m}}}$. In contrast to the Trouton ratio for a Newtonian fluid given by (5.1), the formulae for a viscoelastic fluid given by (5.3a,b) depends on the axial coordinate, z. Evaluating (5.3a,b) at the exit of the pipe, we find $\textrm{T}{\textrm{r}_{ex}} \equiv \textrm{Tr}(1,\eta ,\varLambda ,D{e_m})$, respectively,

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

and

(5.4b) \begin{align}\textrm{T}{\textrm{r}_{ex}} = \underbrace{{3\left( {1 + \eta \frac{{D{e_m}(1 + 2D{e_m})}}{{(1 - 2D{e_m})(1 + D{e_m})}}} \right)}}_{{\textrm{T}{r_h}}} - \eta \underbrace{{\left( {\frac{{2{\varLambda ^{2(2 - 1/D{e_m})}}}}{{1 - 2D{e_m}}} + \frac{{{\varLambda ^{ - 2(1 + 1/D{e_m})}}}}{{1 + D{e_m}}}} \right)}}_{{\textrm{T}{r_{nh}}}}.\end{align}

Based on the quantities indicated in (5.4b) $\textrm{T}{\textrm{r}_{ex}} = \textrm{T}{\textrm{r}_h} - \eta \,\textrm{T}{\textrm{r}_{nh}}$ where $\textrm{T}{\textrm{r}_h},\textrm{T}{\textrm{r}_{nh}} > 0$. It is not difficult to recognize that $\textrm{T}{\textrm{r}_h}$ is the solution for homogeneous extension, while $- \eta \,\textrm{T}{\textrm{r}_{nh}}$ is a negative contribution that can be seen as a non-homogeneous term. The formula for $\textrm{T}{\textrm{r}_h}$ is derived from the Oldroyd-B model assuming that a pure steady uniaxial elongation is imposed (Bird et al. 1987; Tanner 2000; Housiadas 2017; Housiadas & Beris 2024b). Therefore, the Trouton ratio of the fluid at the exit of the hyperbolic pipe is always less than the corresponding Touton ratio for steady homogenous extension. It is very interesting that although $\textrm{T}{\textrm{r}_h}$ is valid and bounded only in the range $0 \le D{e_m} < 1/2$ (the well-known singularity of the Oldroyd-B model for steady homogeneous extension at a finite Deborah number), the formula for $\textrm{T}{\textrm{r}_{ex}} = \textrm{T}{\textrm{r}_h} - \eta \,\textrm{T}{\textrm{r}_{nh}}$ is well-defined and finite for any $D{e_m} > 0$. Similar formulae have also been derived and discussed for the steady viscoelastic flow with slip along the walls of a hyperbolic planar channel by Housiadas & Beris (2024b,c).

In figure 2(a), (5.3a) or (5.3b), is shown as function of z for De_m = 0.3 (black solid line), 0.5 (dashed red line) and 0.7 (dotted blue line) for $\varLambda $ = 3 and η = 4/10. A monotonic increase of the Trouton ratio with the distance from the inlet is depicted, and the same holds with the increase of $\varLambda $. Also, figure 2(b) shows the Trouton ratio at z = 1, i.e. $\textrm{T}{\textrm{r}_{ex}}$, in logarithmic scale as a function of $D{e_m}$, for η = 4/10 and $\varLambda $ = 2 (solid black line), 4 (dashed red line) and 8 (dotted blue line). A smooth enhancement of the Trouton ration is observed, although at high enough modified Deborah number a slight decrease is also seen.

Figure 2. The Trouton ratio, (5.3a) or (5.3b), for a fluid with η = 4/10 based on the Newtonian velocity profile at the classic lubrication limit (i.e. (4.6)) evaluated at y = 0: (a) versus z; (b) evaluated at z = 1 versus De_m.

We also comment on the peculiar behaviour of ${C_{zz}}(0,z)$ as $D{e_m} \to 1/2$. In this case, ${C_{zz}}(0,z)$ is defined only as a limit which is given by

(5.5)\begin{equation}\lim_{D{e_m} \to 1/2} {C_{zz}}(0,z) = 1 + 2\,\textrm{ln}(1 + ({\varLambda ^2} - 1)z).\end{equation}

Note, however, that ${C_{zz}}(0,z)$ is continuous and differentiable for any $D{e_m} > 0$ (including $D{e_m} = 1/2$). From (5.3a,b)–(5.5) we find the corresponding limits of $\textrm{Tr}$ and $\textrm{T}{\textrm{r}_{ex}}$, respectively,

(5.6)\begin{equation}\lim_{D{e_m} \to 1/2} \textrm{Tr} = 3 + \eta \left( {4\,\textrm{ln}(1 + ({\varLambda ^2} - 1)z) - \frac{2}{{3(1 + ({\varLambda ^2} - 1)z){)^3}}} - \frac{7}{3}} \right),\end{equation}

and

(5.7)\begin{equation}\lim_{D{e_m} \to 1/2} \textrm{T}{\textrm{r}_{ex}} = 3 + \eta \left( {8\,\textrm{ln}(\varLambda ) - \frac{2}{{3{\varLambda ^6}}} - \frac{7}{3}} \right).\end{equation}

We conclude this section by emphasizing that, generally, the Trouton ratio is defined only for pure extensional flow. In the hyperbolic tube studied here, the flow is purely extensional exclusively along the axis of symmetry of the pipe, y, making its calculation meaningful only at y = 0. Additional formulae for the Trouton ratio in series form follow in § 6, further comments and discussion are given in § 8, while advancements on the same issue for viscoelastic flows, with or without slip along the walls of a hyperbolic channel, can be found in Housiadas & Beris (2024a,b,c).

6. Low Deborah number analysis of the lubrication equations

Instead of solving the (4.1)–(4.4) for the primary field variables $\{ U,V,p,{\tau _{zz}},{\tau _{yz}},{\tau _{yy}},{\tau _{\theta \theta }}\}$, we introduce the stream function, $\varPsi$, which is defined with the aid of the two velocity components,

(6.1a,b)\begin{equation}U = \frac{1}{y}\frac{{\partial \varPsi }}{{\partial y}},\quad V ={-} \frac{1}{y}\frac{{\partial \varPsi }}{{\partial z}}.\end{equation}

Thus, the continuity equation, (4.1), is satisfied automatically and the pressure gradient can be eliminated by differentiating (4.2) with respect to y. Thus, (4.1)–(4.3) are replaced with the following equation:

(6.2)\begin{align} R{e_m}\dfrac{\partial }{{\partial y}}\left( {\dfrac{\textrm{D}}{{\textrm{D}t}}\left( {\dfrac{1}{y}\dfrac{{\partial \varPsi }}{{\partial y}}} \right)} \right) & = \dfrac{{(1 - \eta )}}{y}\left( {\dfrac{{{\partial^4}\varPsi }}{{\partial {y^4}}} - \dfrac{2}{y}\dfrac{{{\partial^3}\varPsi }}{{\partial {y^3}}} + \dfrac{3}{{{y^2}}}\dfrac{{{\partial^2}\varPsi }}{{\partial {y^2}}} - \dfrac{3}{{{y^3}}}\dfrac{{\partial \varPsi }}{{\partial y}}} \right)\nonumber\\ & \quad + \eta \left( {\dfrac{{{\partial^2}{\tau_{zz}}}}{{\partial z\partial y}} + \dfrac{{{\partial^2}{\tau_{yz}}}}{{\partial {y^2}}} + \dfrac{1}{y}\dfrac{{\partial {\tau_{yz}}}}{{\partial y}} - \dfrac{{{\tau_{yz}}}}{{{y^2}}}} \right), \end{align}

where

(6.3)\begin{align}\frac{\partial }{{\partial y}}& \left( {\frac{\textrm{D}}{{\textrm{D}t}}\left( {\frac{1}{y}\frac{{\partial \varPsi }}{{\partial y}}} \right)} \right)\nonumber\\ &\quad = \frac{1}{{{y^2}}}\left( {\left( {\frac{3}{{{y^2}}}\frac{{\partial \varPsi }}{{\partial z}} - \frac{1}{y}\frac{{\partial \varPsi }}{{\partial y}}\frac{{\partial \varPsi }}{{\partial z}} + \frac{{{\partial^3}\varPsi }}{{\partial y\partial {z^2}}}} \right)\frac{{\partial \varPsi }}{{\partial y}} + \left( {\frac{3}{y}\frac{{{\partial^2}\varPsi }}{{\partial {y^2}}} - \frac{{{\partial^3}\varPsi }}{{\partial {y^3}}}} \right)\frac{{\partial \varPsi }}{{\partial z}}} \right).\end{align}

Since both U and V (and their gradients) are expressed in terms of derivatives of the stream function only, a reference value for $\varPsi$ can be chosen arbitrary. Substituting U in the total mass balance, (2.31d), gives $\varPsi (H(z),z) - \varPsi (0,z) = 1$ which indicates that one of $\varPsi (0,z)$ and $\varPsi (H(z),z)$ can be set conveniently. Also, considering the boundary conditions at the wall and the restrictions needed so that the governing equations are well-defined on the axis of symmetry of the pipe, the final auxiliary conditions for the stream function, imposed for any $0 \le z \le 1$, are summarized as follows:

(6.4a–c) \begin{align}\left.\begin{aligned}\varPsi (0,z) &= \frac{{\partial \varPsi }}{{\partial y}}(0,z) = \frac{{{\partial ^3}\varPsi }}{{\partial {y^3}}}(0,z) = 0,\quad \varPsi (H(z),z) = 1,\\ \frac{{\partial \varPsi }}{{\partial z}}(H(z),z) &= \frac{{\partial \varPsi }}{{\partial y}}(H(z),z) = 0.\end{aligned}\right\}\end{align}

The conditions approaching the symmetry axis, can also be expressed as $\varPsi = O(\kern0.08em {y^2})$ as $y \to 0$ (Sisavath, Jing & Zimmerman 2001).

For a Newtonian fluid (De = 0) and creeping flow conditions ($R{e_m} = 0$) the analytical solution for the stream function in the varying region of the pipe, $0 \le z \le 1$, is

(6.5)\begin{equation}{\varPsi _N} = \frac{{{y^2}}}{{{H^2}}}\left( {2 - \frac{{{y^2}}}{{{H^2}}}} \right).\end{equation}

It can be trivially confirmed that (6.5) satisfies all conditions given in (6.4), as well as that $0 \le {\varPsi _N} \le 1$. At the entrance region, $z \le 0$, the shape function is constant, i.e. $H(z) = 1$, and (6.5) reduces to ${\varPsi _{en}} = {y^2}(2 - {y^2})$ where we have used the subscript ‘en’ in accordance with (4.5).

For the solution of (5.2) and (4.4a–d) we assume a regular perturbation scheme in terms of the Deborah number,

(6.6)\begin{equation}f \approx {f_N} + \sum_{k = 1}^\infty {{f_k}\,D{e^k}} \quad 0 < De \ll 1,\end{equation}

where $f = \varPsi ,\;{\tau _{zz}},\;{\tau _{yy}}$, ${\tau _{yz}}$ and ${\tau _{\theta \theta }}$. The series are substituted into the governing equations and accompanied boundary conditions yielding a sequence of equations which are solved analytically up to $O(D{e^8})$ with the aid of MATHEMATICA software (Wolfram 2023). This quite demanding task requires advanced techniques and substantial computer resources to be achieved; it is also more difficult than the planar case solved and presented recently by Housiadas & Beris (2023, 2024a). The $O(D{e^0}) = O(1)$ equations correspond to the equations for a simple Newtonian fluid, while the effect of viscoelasticity is built with the addition of the $O(D{e^k})$ terms, k = 1, 2, …, 8. Worth mentioning is also the fact that with this method of solution, inlet (at z = 0) and/or outlet (at z = 1) conditions for the stream function and the viscoelastic extra-stresses are not needed and cannot be imposed. Inlet (or initial) conditions, however, are required when the equations are solved numerically; more details are given in § 6. It will become clear below that the calculation of many terms in the asymptotic perturbation series is of great importance because it provides adequately information capable of demonstrating the convergence of the major quantities of interest.

The $O(De)$ solution is

(6.7)\begin{equation}\left. {\begin{array}{@{}c@{}} {{\varPsi _1} = 0,\quad {\tau_{zz,1}} = \dfrac{{128}}{{{H^6}}}{{\bar{y}}^2},\quad {\tau_{yz,1}} = \dfrac{{128}}{{{H^3}}}\left( {\dfrac{{H^{\prime}}}{{{H^3}}}} \right)(2{{\bar{y}}^2} - 1)\bar{y},}\\[12pt] {{\tau_{yy,1}} = 32(5 - 18{{\bar{y}}^2} + 17{{\bar{y}}^4}){{\left( {\dfrac{{H^{\prime}}}{{{H^3}}}} \right)}^2} - 32{{({{\bar{y}}^2} - 1)}^2}\dfrac{{H^{\prime\prime}}}{{{H^5}}},}\\[12pt] {{\tau_{\theta \theta ,1}} = 160{{({{\bar{y}}^2} - 1)}^2}{{\left( {\dfrac{{H^{\prime}}}{{{H^3}}}} \right)}^2} - 32{{({{\bar{y}}^2} - 1)}^2}\dfrac{{H^{\prime\prime}}}{{{H^5}}},} \end{array}} \right\}\end{equation}

where we have used the normalized radial coordinate $\bar{y} = y/H(z)$ for brevity. Since ${\varPsi _1} = 0$, the $O(De)$ velocity field is zero as previously predicted by Boyko & Stone (2022) and Housiadas & Beris (2023, 2024a) for the planar case, and according to the theorem of Tanner & Pipkin (1969). Note that (6.7) is given in a suitable form in terms of the quantities $H^{\prime}/{H^3}$ and $H^{\prime\prime}/{H^5}$ which for the hyperbolic function given by (2.30) are constants, i.e. $- ({\varLambda ^2} - 1)/2$ and $3{({\varLambda ^2} - 1)^2}/4$, respectively. This implies that ${H^6}{\tau _{zz,1}}$, ${H^3}{\tau _{yz,1}}$, ${\tau _{yy,1}}$ and ${\tau _{\theta \theta ,1}}$ depend on $\bar{y}$ only, i.e. there is not any axial evolution of the extra-stresses.

The $O(D{e^2})$ solution is

(6.8a)\begin{gather}{\varPsi _2} = \eta {\bar{y}^2}(\hat{\varPsi }_2^{(2)}(z) + {\bar{y}^2}\hat{\varPsi }_2^{(4)}(z) + {\bar{y}^4}\hat{\varPsi }_2^{(6)}(z) + {\bar{y}^6}\hat{\varPsi }_2^{(8)}(z)),\end{gather}

(6.8b)\begin{gather}{\tau _{zz,2}} = \frac{{3072}}{{{H^6}}}\left( {\frac{{H^{\prime}}}{{{H^3}}}} \right)(1 - {\bar{y}^2}){\bar{y}^2},\end{gather}

(6.8c)\begin{gather}{\tau _{yz,2}} = \frac{{256}}{{{H^3}}}\bar{y}(\hat{\tau }_{yz,2}^{(1)}(z) + {\bar{y}^3}\hat{\tau }_{yz,2}^{(3)}(z) + {\bar{y}^5}\hat{\tau }_{yz,2}^{(5)}(z)),\end{gather}

(6.8d)\begin{gather}{\tau _{yy,2}} = 128(\hat{\tau }_{yy,2}^{(0)}(z) + {\bar{y}^2}\hat{\tau }_{yy,2}^{(2)}(z) + {\bar{y}^4}\hat{\tau }_{yy,2}^{(4)}(\bar{z}) + {\bar{y}^6}\hat{\tau }_{yy,2}^{(6)}(z)),\end{gather}

(6.8e)\begin{gather}{\tau _{\theta \theta ,2}} = 128(\hat{\tau }_{\theta \theta ,2}^{(0)}(z) + {\bar{y}^2}\hat{\tau }_{\theta \theta ,2}^{(2)}(z) + {\bar{y}^4}\hat{\tau }_{\theta \theta ,2}^{(4)}(z) + {\bar{y}^6}\hat{\tau }_{\theta \theta ,2}^{(6)}(z)),\end{gather}

where $\hat{\varPsi }_2^{(j)}$, $\hat{\tau }_{yz,2}^{(j)}$, $\hat{\tau }_{zz,2}^{(j)}$ and $\hat{\tau }_{\theta \theta ,2}^{(j)}$ are functions of the shape function and its derivative(s) with respect to the axial coordinate. The precise form of these functions is provided in the Appendix A. The higher-order solutions are too long to be printed here.

Finally, we confirm that for a hyperbolic pipe, at any order of the perturbation scheme, i.e. for k = 0,1,2,…,8, ${\varPsi _k}$, ${H^6}{\tau _{zz,k}}$, ${H^3}{\tau _{yz,k}}$, ${\tau _{yy,k}}$ and ${\tau _{\theta \theta ,k}}$ depend on $\bar{y}$ only (see the Appendix A for k = 2), as well as that on the axis of symmetry of the pipe ${\tau _{zz,k}}(0,z) = {\tau _{yz,k}}(0,z) = 0$.

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 (2022) and Housiadas & Beris (2023), albeit these authors did not report the corresponding formula; the latter can be found in Housiadas & Beris (2024a) 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é 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 2022; Housiadas & Beris 2023, 2024a). 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, 2023, 2024a). 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 (2022) and further studied recently up to $O({\varepsilon ^{20}})$ by Sialmas & Housiadas (2024), 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 1996). Pakdel & McKinley (1996) 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 (1996) investigated in detail the lid-driven cavity flow problem of two different ideal elastic fluids, while McKinley et al. (1996) extended the analysis of Pakdel & McKinley (1996) 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. (1996) 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 (1996) 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$).

7. Numerical solution of the lubrication equations

Equations (6.2) and (4.4a–d) consist of an initial (in the axial direction, z) and boundary (in the radial direction, y) value problem accompanied with the boundary conditions given in (6.4) and appropriate initial conditions. For the numerical solution of this system, first we introduce new coordinates $(\xi ,\zeta )$ that map the varying boundary of the flow domain into a fixed one,

(7.1a,b)\begin{equation}\xi ={-} 1 + \frac{{2y}}{{H(z)}},\quad \zeta = z\; \Leftrightarrow \;y = \frac{{\xi + 1}}{2}H(\zeta ),\quad z = \zeta .\end{equation}

Thus, the shape function H appears into the governing equations and the domain of definition of the lubrication equations becomes $\bar{\varOmega } = \{ (\xi ,\zeta )|- 1 < \xi < 1,\;0 < \zeta \le 1\}$. Based on (7.1) the first-order partial differential operators are transformed as follows:

(7.2a,b)\begin{equation}\frac{\partial }{{\partial y}} = \frac{2}{{H(\zeta )}}\frac{\partial }{{\partial \xi }},\quad \frac{\partial }{{\partial z}} = \frac{\partial }{{\partial \zeta }} - \frac{{(1 + \xi )H^{\prime}(\zeta )}}{{H(\zeta )}}\frac{\partial }{{\partial \xi }},\end{equation}

where we have used $\partial H/\partial y = 0$ to show that $H = H(\zeta )$ only. Using (7.2) we find the higher-order derivatives with respect to y and/or z in terms of derivatives with respect to $\xi$ and $\zeta$. Substituting the transformed derivatives and the rescaled components of the conformation tensor into (6.2) and (4.4a–d) gives the final equations for $\varPsi$, ${c_{zz}}$, ${c_{yy}}$ and ${c_{yz}}$ defined on $\bar{\varOmega }$; notice that we prefer to solve for the conformation tensor components instead of the components of the viscoelastic extra-stress tensor. Based on the stream function and the new coordinates, the material derivative is transformed as follows:

(7.3)\begin{equation}\frac{\textrm{D}}{{\textrm{D}t}} \equiv \frac{4}{{(1 + \xi ){H^2}}}\left( {\frac{{\partial \varPsi }}{{\partial \xi }}\frac{\partial }{{\partial \zeta }} - \frac{{\partial \varPsi }}{{\partial \zeta }}\frac{\partial }{{\partial \xi }}} \right).\end{equation}

Similarly, the velocity components are given as

(7.4a,b)\begin{equation}U = \frac{4}{{(1 + \xi ){H^2}}}\frac{{\partial \varPsi }}{{\partial \xi }},\quad V = \frac{2}{H}\left( {\frac{{H^{\prime}}}{H}\frac{{\partial \varPsi }}{{\partial \xi }} - \frac{1}{{(1 + \xi )}}\frac{{\partial \varPsi }}{{\partial \zeta }}} \right).\end{equation}

Finally, the velocity gradients are

(7.5)\begin{gather}\frac{{\partial U}}{{\partial z}} = \frac{4}{{(1 + \xi ){H^2}}}\left( {\frac{{{\partial^2}\varPsi }}{{\partial \xi \,\partial \zeta }} - \frac{{H^{\prime}}}{H}\frac{{\partial \varPsi }}{{\partial \xi }} - \frac{{H^{\prime}}}{H}\frac{{{\partial^2}\varPsi }}{{\partial {\xi^2}}}(1 + \xi )} \right),\end{gather}

(7.6)\begin{gather}\frac{{\partial U}}{{\partial y}} = \frac{8}{{(1 + \xi ){H^3}}}\left( {\frac{{{\partial^2}\varPsi }}{{\partial {\xi^2}}} - \frac{1}{{(1 + \xi )}}\frac{{\partial \varPsi }}{{\partial \xi }}} \right),\end{gather}

(7.7)\begin{gather}\begin{aligned}[b]\frac{{\partial V}}{{\partial z}} &={-} \frac{2}{{(1 + \xi )H}}\left\{ \frac{{{\partial^2}\varPsi }}{{\partial {\zeta^2}}} - \left( {\left( {\frac{{H^{\prime\prime}}}{H} - \frac{{2{{H^{\prime}}^2}}}{{{H^2}}}} \right)\frac{{\partial \varPsi }}{{\partial \xi }} + \frac{{2H^{\prime}}}{H}\frac{{{\partial^2}\varPsi }}{{\partial \zeta \,\partial \xi }}} \right)(1 + \xi )\right.\\ &\quad \left. +\, \frac{{{{H^{\prime}}^2}}}{{{H^2}}}\frac{{{\partial^2}\varPsi }}{{\partial {\xi^2}}}{{(1 + \xi )}^2} \right\}.\end{aligned}\end{gather}

Note that $\partial V/\partial y$ is eliminated with the aid of the continuity equation, i.e. $\partial V/\partial y ={-} \partial U/\partial z - V/y$. As far as the initial conditions are concerned, the solution at the entrance region, i.e. for $z \le 0$, is imposed for the stream function,

(7.8a,b)\begin{equation}\varPsi (y,0) \equiv {\varPsi _{en}} = {y^2}(2 - {y^2}),\quad \frac{{\partial \varPsi }}{{\partial z}}(y,0) \equiv \frac{{\partial {\varPsi _{en}}}}{{\partial z}}(y,0) = 0,\end{equation}

and, similarly, for the rescaled components of the conformation tensor (for ${\tau _{zz}}$, ${\tau _{yz}}$ and ${\tau _{yy}}$ see § 4.1):

(7.9a–c)\begin{equation}{c_{zz}}(y,0) = 128D{e^2}{y^2},\quad {c_{yz}}(y,0) ={-} 8De\,y,\quad {c_{yy}}(y,0) = 1.\end{equation}

We emphasize that (7.8b) requires special treatment in the new coordinate system (see (7.2) with $H(0) = 1$)

(7.10)\begin{equation}\frac{{\partial \varPsi }}{{\partial \zeta }}(\xi ,0) = H^{\prime}({0^ + })(1 + \xi ){\varPsi ^{\prime}_{en}}(\xi ).\end{equation}

The integration along $\zeta$ is performed by discretizing the computational domain in N + 1 equidistant grid points, ${\zeta _n} = n\Delta \zeta = n/N$, $n = 0,1,2,3, \ldots ,N$. We use the A-stable and first-order accurate finite-difference formula for $\partial f/\partial \zeta$,

(7.11a)\begin{equation}\frac{{\partial f}}{{\partial \zeta }}(\xi ,{\zeta _n}) \equiv {\left. {\frac{{\partial f}}{{\partial \zeta }}} \right|_{(n)}} \approx \frac{{{f_{(n)}} - {f_{(n - 1)}}}}{{\mathrm{\Delta }\zeta }},\end{equation}

where $f = \varPsi$, ${c_{zz}}$, ${c_{yz}}$ and ${c_{yy}}$, while for ${\partial ^2}\varPsi /\partial {\zeta ^2}$ we use the first-order accurate formula,

(7.11b)\begin{equation}\frac{{{\partial ^2}\varPsi }}{{\partial {\zeta ^2}}}(\xi ,{\zeta _n}) \equiv {\left. {\frac{{{\partial^2}\varPsi }}{{\partial {\zeta^2}}}} \right|_{(n)}} \approx \frac{1}{{\mathrm{\Delta }\zeta }}\left( {{{\left. {\frac{{\partial \varPsi }}{{\partial \zeta }}} \right|}_{(n)}} - {{\left. {\frac{{\partial \varPsi }}{{\partial \zeta }}} \right|}_{(n - 1)}}} \right).\end{equation}

In the above formulae, $n \ge 1$ and a subscript in parenthesis denotes value at the corresponding grid $\zeta$-point, i.e. ${f_{(n)}} \equiv f(\xi ,{\zeta _n})$. This is achieved by developing a pseudospectral method. The field variables are given as series of Chebyshev orthogonal polynomials in terms of $\xi$, as follows:

(7.12)\begin{equation}{f_{(n)}} \equiv f(\xi ,{\zeta _n}) = \sum_{k = 0}^M {{{\hat{f}}_{2k}}({\zeta _n}){T_{2k}}(\xi )} ,\quad f = \varPsi ,{c_{zz}},{c_{yy}},{c_{yz}}.\end{equation}

The Chebyshev orthogonal polynomials are defined as ${T_k}(\xi ) = \textrm{cos}(k\,\textrm{co}{\textrm{s}^{ - 1}}(\xi ))$, $k = 0,1,2, \ldots ,M$ in the domain $\xi \in [ - 1,1]$. The spectral coefficients, namely the quantities with a hat in (6.12), are calculated pseudospectrally based on the Gauss–Lobatto points (Hesthaven, Gottlieb & Gottlieb 2007) $\{ {\xi _j}\} _{j = 0}^M = \{ - \textrm{cos}(j{\rm \pi} /M)\} _{j = 0}^M$, at which the unknown physical values ${f_{j,(n)}} \approx f({\xi _j},{\zeta _n})$ are assigned. The partial derivatives of f with respect to $\xi$ at the Gauss–Lobatto points are evaluated effectively based on the differentiation matrices (Hesthaven, Gottlieb & Gottlieb 2007). Then, ${f_{j,(n)}}$ are calculated according to the discretized version of the governing equations and accompanied boundary conditions at the Gauss–Lobatto points. The resulting strongly nonlinear system of algebraic equations is solved iteratively using a Newton scheme at each n-step yielding the solution for ${f_{j,(n)}}$; this scheme converges quadratically (as it should) with an absolute relative criterion ${10^{ - 12}}$ within two or three iterations.

The accuracy of the solution calculated by our pseudospectral code is checked by monitoring the magnitude of the spectral coefficients for the stream function. In figure 7(a), $|{\Delta {{\hat{\varPsi }}_k}} |$ is shown in logarithmic plot for $k=1,2,\ldots,M$, where M = 15, for ζ = 0.2, ζ = 0.5 (at the middle of the pipe) and ζ = 1 (at the outlet) for the case with $\eta = 0.4$ and $D{e_m} = 0.45$. First, however, we report that ${\varPsi _{en}}$ and ${\varPsi _N}$ have the same form when are expressed in terms of the new coordinates

(7.13)\begin{equation}{\varPsi _{en}}(\xi ) = {\varPsi _N}(\xi ) = {\textstyle{1 \over {16}}}{(1 + \xi )^2}(7 - 2\xi - {\xi ^2}),\end{equation}

as well as that there is no dependence on the ζ-coordinate. From (7.13) one can trivially confirm that ${\varPsi _j}( - 1) = {\varPsi ^{\prime}_j}( - 1) = {\varPsi ^{\prime\prime\prime}_j}( - 1) = 0$, ${\varPsi _j}(1) = 1,$ ${\varPsi ^{\prime}_j}(1) = 0$ where j = ‘en’ or ‘N’, and also to calculate the corresponding spectral coefficients:

(7.14)\begin{equation}\{ {\hat{\varPsi }_{en,k}}\} _{k = 0}^4 = \{ {\hat{\varPsi }_{N,k}}\} _{k = 0}^4 = \left\{ {\frac{{61}}{{128}},\frac{9}{{16}},\frac{1}{{32}}, - \frac{1}{{16}}, - \frac{1}{{128}}} \right\}.\end{equation}

Figure 9(a) shows that the magnitude of the spectral coefficients of $\mathrm{\Delta }\varPsi \equiv \varPsi (\xi ,\zeta ) - {\varPsi _N}(\xi )$ drops fast, almost down to machine accuracy, revealing that the stream function is fully resolved. The very small magnitude of $|{\Delta {{\hat{\varPsi }}_k}} |$, $k = 0,1,2, \ldots ,15$, reveals that the stream function, and consequently the velocity field, only slightly changes compared with the Newtonian velocity field (or compared with the solution for the straight pipe at the entrance region). The difference $\Delta \varPsi$ is also shown as function of the transformed radial coordinate near the inlet (ζ = 0.2), at the middle of the pipe (ζ = 1/2) and at the outlet (ζ = 1) in figure 9(b).

Figure 9. Simulations results for η = 4/10 and De_m = 0.45. (a) The magnitude of the spectral coefficients for ΔΨ; circles, ζ = 0.2, squares, ζ = 0.5; diamonds, ζ = 1. (b) Here ΔΨ = Ψ − Ψ_N as a function of the transformed radial coordinate ξ.

In table 2, we compare the results obtained with our pseudospectral method with the accelerated analytical solutions. The excellent agreement between the numerical and analytical results reveals the validity and accuracy for both the theoretical and numerical results up to at least $D{e_m} \approx 1$. The same accuracy is also achieved for the individual contributions to the normalized pressure-drop (shown in figure 5), and for the Trouton ratio (shown in figure 7a).

Table 2. Comparison of $\Delta \varPi /\Delta {\varPi _N}$ calculated numerically and analytically (rounded in four significant digits) using the Padé [M/M] diagonal approximants with M = 1, 2, 3 and 4. The percentage relative absolute error is given in parenthesis; a star indicates error less than 0.01 %.

Last, in figure 10, we present results for the rescaled components of the conformation tensor, ${H^6}{c_{zz}}$, ${H^3}{c_{yz}}$ and ${c_{yy}}$ as functions of the transformed radial coordinate, $\xi$. It is seen that these quantities remain practically the same for any ζ > 0; the curves for ζ = 0 correspond to the initial condition at the inlet. In other words, there is no evolution of the rescaled components of the conformation tensor with the distance from the inlet. Indeed, this is fully consistent with the perturbation solution found here up to $O(De_m^8)$; see (6.7) and (6.8b–e) along with the corresponding functions given in the Appendix A. The results reveal that ${c_{zz}} \gg {c_{zz,N}} = 0$ and that the magnitude of ${c_{zz}}$ is much larger that the magnitudes of ${c_{yz}}$ and ${c_{yy}}$. These clearly indicate that ${c_{zz}}$ (and the corresponding extra-stress ${\tau _{zz}}$) is the leading contribution to the pressure drop due to the viscoelasticity of the fluid.

Figure 10. Simulations results for $\varLambda $ = 3 and De_m = 0.45 as function of the transformed radial coordinate, ξ: (a) H ⁶c_zz; (b) H ³c_yz; (c) c_yy. All curves for ζ > 0 (shown with dashed lines) collapse.

8. Discussion

We showed above that the high-order asymptotic formulae for the Trouton ratio, (6.33a) and (6.33b), requires special attention. The main reason is that for the hyperbolic geometry the normal component of the conformation tensor, c_zz, (or τ_zz) shows a peculiar behaviour at De_m = 1/2; consequently, this peculiarity affects the Trouton ratio at De_m = 1/2 too. Thus, it did not come as a surprise that the perturbation solutions and the corresponding transformed formulae diverge near that point, which clearly poses an upper bound for their validity. These results are closely related to the known deficiency of the Oldroyd-B model which predicts an infinite Trouton ratio in homogeneous uniaxial extensional flow at a finite extension rate which corresponds to $D{e_m} = 1/2$ (Bird et al. 1987; Tanner 2000; Housiadas 2017; Housiadas & Beris 2024b).

The explanation for this peculiar behaviour of the Trouton ratio at $D{e_m} = 1/2$ is given with the aid of the exact analytical solution for ${C_{zz}}(0,z)$ (see (5.3a)) and its limit as $D{e_m} \to 1/2$ given by (5.5). These solutions directly affect the Trouton ratio at the exit of the pipe and its limit as $D{e_m} \to 1/2$ given by (5.4b) and (5.6), respectively. It is important to mention, however, that the exact solutions for ${C_{zz}}(0,z)$ and $\textrm{Tr}$ are well-behaved for any $D{e_m} > 0$. The formula around which the exact solution for $\textrm{Tr}$ was developed, (3.18b) reveals the presence of a power law term with exponent proportional to the inverse of the original Deborah number (see (3.16)). This term is not amenable to a representation in terms of a regular perturbation series with respect to De, causing the series to diverge at infinity at the same modified Deborah number as the exact solution for pure homogeneous uniaxial extensional flow (i.e. at $D{e_m} = 1/2$).

One can therefore naturally ask the following question. Suppose we do not know the well-behaved for any $D{e_m} > 0$ exact solution, is it possible, using a regular perturbation scheme, to provide reasonable and/or accurate predictions for the Trouton ratio at De_m values larger than 1/2? Although a strict and formal answer to this important issue cannot be given, previous work by Housiadas (2017) on pure homogeneous elongational flows showed that using only a series developed asymptotically as $D{e_m} \to {0^ + }$, this is not possible. Even worse, the use of more realistic (and highly nonlinear with respect to the polymer extra-stress tensor) models like the Giesekus, FENE-P (Bird et al. 1987) and PTT (Tanner 2000) models does not circumvent or resolve this issue.

However, recent work by Housiadas (2023) revealed that a possible answer to this question is to consider both limits (low and high Deborah numbers), in deriving asymptotic solutions. Then, use of convergence acceleration methods (like, for instance, two-point Padé approximants) followed by the development of uniformly valid approximations, can resolve this problem. The performance of this method, which depends on the available number of analytical terms at both limits, can be amazing. Another possible answer to this question may be to consider advanced perturbation schemes, such as the exponential asymptotics presented by Kataoka & Akylas (2022), that may capture more efficiently the anomaly/peculiar behaviour of the constitutive model at De_m = 1/2. This issue requires further in-depth investigation because for a long time is believed that the regular perturbation schemes in terms of the Deborah or the Weissenberg numbers are the only ones that should be used for developing asymptotic solutions for weakly viscoelastic flows (Bird et al. 1987).

9. Conclusions

We studied theoretically the steady viscoelastic flow of an Oldroyd-B fluid under creeping conditions in a long axisymmetric tube with variable cross-section, focusing on the hyperbolic contracting pipe. First, we developed the general theoretical framework for the evaluation of the average pressure drop, required to maintain a constant flow rate through the pipe. Second, we proved that the velocity field along the axis of symmetry of the pipe is a pure uniaxial extensional field, and we found the exact analytical solution for the components of the constitutive model on the axis of symmetry in terms of the fluid velocity only. The exact solution for the constitutive model along with the observation that the flow velocity changes little with viscoelasticity along the axis of symmetry (at least in the limit of a small aspect ratio, less than 0.1 approximately) allowed for a robust and accurate evaluation of the first normal stress difference up to high Deborah numbers, through which the dimensionless elongational viscosity of the fluid (or Trouton ratio, $\textrm{Tr}$) was extracted. For a hyperbolic contracting pipe and approximating the velocity with that for a Newtonian fluid under the same conditions and geometry, we showed that the Trouton ratio increases with increasing the Deborah number, polymer viscosity ratio and contraction ratio. The calculation of the Trouton ratio based on the exact analytical solution of the constitutive model along the centreline, (3.18a,b), is a major achievement of this work. It not only provides the correct result for the Trouton ratio for any $D{e_m} > 0$, but also reveals the reason for the failure of the perturbation series for the Trouton ratio (and their corresponding transformed formulae) for $D{e_m} \ge 1/2$.

Further analytical progress was achieved by invoking the classic lubrication approximation and solving the final equations for an arbitrary shape function, using a high-order asymptotic scheme in terms of the original Deborah number, De, valid for any value of the polymer viscosity ratio, i.e. for $0 \le \eta \le 1$. In this case, the entrance and exit regions of the tube do not play any role in the analysis, and consequently they do not affect the results. The lubrication equations were solved analytically up to eighth-order in the Deborah number, resulting in analytical formulae for all the field variables, the average pressure drop, $\Delta \varPi$, and the Trouton ratio, $\textrm{Tr}$; we reiterate that the latter is defined and is meaningful only when the flow is purely extensional, i.e. only along the axis of symmetry of the pipe. For the hyperbolic case, we revealed that the high-order perturbation solutions, when reduced by their corresponding Newtonian values, namely $\Delta \varPi /\Delta {\varPi _N}$ and $\textrm{Tr}$ given by (6.15) and (6.33a), respectively, can be recast in terms of $D{e_m}$ and $\eta$ only. The convergence of the perturbation series was checked and enhanced by deriving transformed solutions using Padé diagonal approximants. The latter showed convergence for $\Delta \varPi /\Delta {\varPi _N}$, for any $D{e_m} > 0$, revealing a monotonic decrease with increasing $D{e_m}$ and/or $\eta$. On the contrary, the transformed solutions for $\textrm{Tr}$ converged only in the window $0 < D{e_m} < 1/2$, showing an increase of the Trouton ratio with increasing $D{e_m}$ and/or $\eta$, and diverging near $D{e_m} \approx 1/2$.

This deficiency/divergence was traced to the very peculiar dependence of the exact analytical solution for the Trouton ratio on $D{e_m}$ and z which changes continuously from an algebraic type dependence for $0 < D{e_m} < 1/2$, to a logarithmic type dependence at $D{e_m} = 1/2$, and back to the algebraic type for $D{e_m} = 1/2$. Therefore, the transformed perturbation solution(s) for the Trouton ratio can be used only for $0 < D{e_m} < 1/2$ while the new developed exact formula ((3.18a) or (3.18b)) allows us to get well-behaved and accurate predictions even for $D{e_m} = 1/2$ given of course that a good approximation for the fluid velocity is available. We were also able to show that even the Newtonian velocity profile can give a sufficiently accurate result, which fully justifies this approach followed by many researchers in the literature. Thus, the corresponding formula for the Trouton ratio given by (5.3a,b), valid for a long hyperbolic axisymmetric pipe, represents the best theoretical predictions for Boger-type fluids (i.e. viscoelastic fluids with negligible shear thinning behaviour). Moreover, comparison of (5.3a), or (5.3b), with the high-order perturbation solution for $\textrm{Tr}$, postprocessed with the Padé diagonal transformation, revealed that the formulae agree very well in the range $0 < D{e_m} < 0.3$, approximately. Note that these observations and predictions are in full accordance with the well-known singularity of the Oldroyd-B model in homogeneous uniaxial extensional flow at a finite rate of extension which corresponds to $D{e_m} = 1/2$.

Two different decompositions for the average pressure drop were also performed with the aid of the total force balance and the total mechanical energy of the flow system in the hyperbolic section of the pipe. Both decompositions revealed the individual contributions to the pressure drop, confirming the correctness of the high-order analytical solution and highlighting important features of the flow. It was shown that the viscous stresses at the wall and the viscoelastic contribution in the bulk increase the pressure drop, while the net viscoelastic contribution between the inlet and outlet is responsible for the decrease of the pressure drop with increasing the fluid viscoelasticity. It was also revealed that the energy loss due to viscous dissipation remains almost constant with the increase of Deborah number, as well as the dissipative nature of viscoelasticity, and the fact that the work done by the normal viscoelastic forces gives energy to the fluid facilitating the transport of the fluid through the pipe. We also emphasize that all contributions to the total force balance and total mechanical energy of the system converge when at least the first five terms in the asymptotic solutions are taken into account in order to generate the Padé [2/2] approximant. Using more terms in the series and generating higher-order Padé approximants results in no visible differences in the illustrated outcomes.

Finally, we mention that our theoretical predictions, the analytical formulae, as well as their transformed (accelerated) solutions were checked for their accuracy and validity through comparison with numerical results obtained using pseudospectral methods to solve the lubrication equations. The comparison revealed excellent agreement between the solutions clearly demonstrating the accuracy, robustness and suitability of the theoretical methods and techniques used here. In conclusion, our high-order asymptotic results in terms of De_m and our pseudospectral numerical results are valid and very accurate in the parameters range $0 \le \varepsilon \le 0.1$, $0 \le D{e_m} \le 1$ and $0 \le \eta \le 1$.