Flow of an Oldroyd-B fluid in a slowly varying contraction: theoretical results for arbitrary values of Deborah number in the ultra-dilute limit

Abstract

Pressure-driven flows of viscoelastic fluids in narrow non-uniform geometries are common in physiological flows and various industrial applications. For such flows, one of the main interests is understanding the relationship between the flow rate $q$ and the pressure drop $\Delta p$, which, to date, is studied primarily using numerical simulations. We analyse the flow of the Oldroyd-B fluid in slowly varying arbitrarily shaped, contracting channels and present a theoretical framework for calculating the $q-\Delta p$ relation. We apply lubrication theory and consider the ultra-dilute limit, in which the velocity profile remains parabolic and Newtonian, resulting in a one-way coupling between the velocity and polymer conformation tensor. This one-way coupling enables us to derive closed-form expressions for the conformation tensor and the flow rate–pressure drop relation for arbitrary values of the Deborah number ($De$). Furthermore, we provide analytical expressions for the conformation tensor and the $q-\Delta p$ relation in the high-Deborah-number limit, complementing our previous low-Deborah-number lubrication analysis. We reveal that the pressure drop in the contraction monotonically decreases with $De$, having linear scaling at high Deborah numbers, and identify the physical mechanisms governing the pressure drop reduction. We further elucidate the spatial relaxation of elastic stresses and pressure gradient in the exit channel following the contraction and show that the downstream distance required for such relaxation scales linearly with $De$.

1. Introduction

Viscoelastic fluid flows in non-uniform geometries consisting of contractions or expansions occur in physiological flows, e.g. arterial flows that may have such shape changes due to thrombus formation (Westein et al. 2013), and in various industrial applications (Pearson 1985). For such flows, one of the key interests is to understand the dependence of the pressure drop $\Delta p$ on the flow rate $q$. It is well known that adding even small amounts of polymer molecules in a Newtonian solvent may drastically change the hydrodynamic features of the flow of the solution due to polymer stretching, which generates elastic stresses in addition to viscous stresses (Bird, Armstrong & Hassager 1987; Alves, Oliveira & Pinho 2021; Steinberg 2021; Datta et al. 2022).

Pressure-driven flows of viscoelastic fluids and the corresponding flow rate–pressure drop relation have been studied extensively in various geometries, mainly through numerical simulations (Szabo, Rallison & Hinch 1997; Alves, Oliveira & Pinho 2003; Binding, Phillips & Phillips 2006; Alves & Poole 2007; Zografos et al. 2020; Varchanis et al. 2022) and experimental measurements (Rothstein & McKinley 1999, 2001; Sousa et al. 2009; Ober et al. 2013; James & Roos 2021). We refer the reader to overviews given recently by Boyko & Stone (2022) and Hinch, Boyko & Stone (2024).

In particular, the abrupt contraction and contraction–expansion channels have received much attention (Rothstein & McKinley 1999; Alves et al. 2003; Binding et al. 2006; Ferrás et al. 2020), and $4\,{:}\,1$ two-dimensional (2-D) and axisymmetric contraction flows have become benchmark flow problems in computational non-Newtonian fluid mechanics (Alves et al. 2021). Numerical simulations of viscoelastic fluid flow in these and other non-uniform geometries include a long downstream (exit) section to allow the stresses to reach their fully relaxed values (see, e.g. Debbaut, Marchal & Crochet 1988; Alves et al. 2003). This is because, once perturbed from their fully relaxed values, the elastic stresses require a long distance for spatial relaxation to enable stable and converged numerical solutions. For higher Deborah ($De$) or Weissenberg ($Wi$) numbers (see definitions in § 2.1), a longer downstream section is required (Keiller 1993).

Therefore, understanding the spatial relaxation of elastic stresses, velocity and pressure is of both fundamental and practical importance, as that determines the size of the computational domain (Alves et al. 2003). However, despite extensive study of viscoelastic channel flows, the spatial relaxation of stresses and pressure in these geometries is not well understood. As a result, the length of the exit channel is currently set somewhat arbitrarily, thus motivating the development of theory. Furthermore, in many applications, it is necessary to determine the total pressure drop over the configuration for a given flow rate, thus requiring us to account for the pressure drop in the entry and exit channels. However, most studies to date focused on the non-uniform region or close vicinity of the abrupt contraction and reported a suitably non-dimensionalized so-called Couette correction (or excess pressure drop), rather than the total non-dimensional pressure drop in the entire configuration (see, e.g. Rothstein & McKinley 1999; Alves et al. 2003; Binding et al. 2006), presumably due to the arbitrariness of the exit channel length in simulations.

One widely used approach to obtaining theoretical results in different viscoelastic fluid flow problems relies on considering the weakly viscoelastic limit by applying a perturbation expansion in powers of the Deborah or Weissenberg number, which are assumed to be small (see, e.g. Datt et al. 2017; Datt, Nasouri & Elfring 2018; Datt & Elfring 2019; Gkormpatsis et al. 2020; Dandekar & Ardekani 2021; Housiadas, Binagia & Shaqfeh 2021; Su et al. 2022). In particular, there have been many applications of such an expansion in conjunction with lubrication theory in studying thin films and tribology problems (Ro & Homsy 1995; Tichy 1996; Sawyer & Tichy 1998; Zhang, Matar & Craster 2002; Saprykin, Koopmans & Kalliadasis 2007; Ahmed & Biancofiore 2021; Gamaniel, Dini & Biancofiore 2021; Ahmed & Biancofiore 2023). Recently, we have applied lubrication theory and such an expansion in powers of $De$, developing a reduced-order model for the steady flow of an Oldroyd-B fluid in a slowly varying, arbitrarily shaped 2-D channel (Boyko & Stone 2022). We provided analytical expressions for the velocity and stress fields and the flow rate–pressure drop relation in the non-uniform region up to $O(De^2)$. We further exploited the reciprocal theorem (Boyko & Stone 2021, 2022) to obtain the flow rate–pressure drop relation at the next order, $O(De^3)$. Housiadas & Beris (2023) extended the low-Deborah-number lubrication analysis of Boyko & Stone (2022) to much higher asymptotic orders and provided analytical expressions for the pressure drop up to $O(De^8$).

However, the low-Deborah-number analysis cannot accurately capture the behaviour at high $De$ numbers where there are significant elastic stresses. Another approach to simplifying the governing equations while capturing the underlying physics at non-small Deborah numbers is to consider the ultra-dilute limit (Remmelgas, Singh & Leal 1999; Moore & Shelley 2012; Li, Thomases & Guy 2019; Mokhtari et al. 2022), $\tilde {\beta }=\mu _{p}/\mu _{0} \ll 1$, where $\mu _{p}$ is the polymer contribution to the total zero-shear-rate viscosity $\mu _{0}$ of the polymer solution. Physically, the ultra-dilute limit corresponds to a low concentration of polymer molecules in a Newtonian solvent, such that the viscosity of the polymer solution, $\mu _{0}$, is only slightly larger than the solvent viscosity, $\mu _{s}$ (Remmelgas et al. 1999; Mokhtari et al. 2022). Furthermore, the limit $\tilde {\beta }=\mu _{p}/\mu _{0} \ll 1$ is closely related to the diluteness criterion of a constant shear-viscosity viscoelastic Boger fluid (Moore & Shelley 2012). In the ultra-dilute limit, the flow field approximated as Newtonian creates elastic stresses that are not coupled back to change the flow. These elastic stresses can then be used to find the correction to the velocity and pressure fields due to fluid viscoelasticity, even at high Deborah numbers. Previous studies used this approach to determine the structure of the stress distribution in the flow around a cylinder (Renardy 2000), a sphere (Moore & Shelley 2012) and arrays of cylinders (Mokhtari et al. 2022), as well as in stagnation (Becherer, Van Saarloos & Morozov 2009; Van Gorder, Vajravelu & Akyildiz 2009) and cross-slot (Remmelgas et al. 1999) flows.

In this work, we continue our theoretical studies (Boyko & Stone 2022; Hinch et al. 2024) of the pressure-driven flow of the Oldroyd-B fluid in slowly varying, arbitrarily shaped, narrow channels. In contrast to Boyko & Stone (2022), who focused only on the flow through a non-uniform channel in the low-Deborah-number limit, and Hinch et al. (2024), who studied numerically the flow through a contraction, expansion and constriction for order-one Deborah numbers, and also provided an asymptotic description at high Deborah numbers, the current work examines the ultra-dilute limit and arbitrary values of the Deborah number. Specifically, we analyse the flow of the Oldroyd-B fluid in a contracting geometry and the relaxation of the elastic stresses and pressure in the exit channel. We apply the lubrication approximation and use a one-way coupling between the velocity and polymer stresses to derive semi-analytical expressions for the conformation tensor in the contraction and the exit channel for arbitrary values of the Deborah number in the ultra-dilute limit. These semi-analytical expressions allow us to calculate the pressure drop and elucidate the relaxation of the elastic stresses and pressure in the exit channel for all $De$. We provide analytical expressions for the conformation tensor and the pressure drop in the high-Deborah-number limit, which are consistent with recent results of Hinch et al. (2024), thus complementing our previous low-Deborah-number lubrication analysis (Boyko & Stone 2022). Furthermore, we analyse the viscoelastic boundary layer near the walls at high Deborah numbers and derive the boundary-layer asymptotic solutions. Given the well-known lack of accuracy and convergence difficulties associated with the high-Weissenberg-number problem in numerical simulations (Owens & Phillips 2002; Alves et al. 2021), our analytical and semi-analytical results for the ultra-dilute limit, valid at high Deborah numbers, are of fundamental importance as they may serve to validate simulation predictions or be compared with experimental measurements to understand more about the applicability of model constitutive equations.

2. Problem formulation and governing equations

We analyse the incompressible steady flow of a viscoelastic fluid in a slowly varying and symmetric 2-D contraction of height $2h(z)$ and length $\ell$, where $h(z)\ll \ell$, as illustrated in figure 1. Upstream of the contraction inlet ($z=0$) there is an entry channel of height $2h_{0}$ and length $\ell _{0}$, and downstream of the contraction outlet ($z=\ell$) there is an exit channel of height $2h_{\ell }$ and length $\ell _{\ell }$. The fluid flow has velocity $\boldsymbol {u}$ and pressure distribution $p$, which are induced by an imposed flow rate $q$ (per unit depth). Our primary interest is to determine the pressure drop $\Delta p$ over the contraction region and the spatial relaxation of pressure and elastic stresses in the exit channel. For our analysis, we shall employ two different systems of coordinates. The first is Cartesian coordinates $(z,y)$ and $(z_{\ell },y)$, where the $z$ and $z_{\ell }=z-\ell$ axes lie along the symmetry midplane of the channel (dashed-dotted line) and $y$ is in the direction of the shortest dimension. The second one is orthogonal curvilinear coordinates $(\xi,\eta )$ defined in § 2.3.

Figure 1. Schematic illustration of the 2-D configuration consisting of a slowly varying and symmetric contraction of height $2h(z)$ and length $\ell$ ($h\ll \ell$). The contraction is connected to two long straight channels of height $2h_{0}$ and $2h_{\ell }$, respectively, up- and downstream and contains a viscoelastic fluid steadily driven by the imposed flow rate $q$.

We consider low-Reynolds-number flows so that the fluid motion is governed by the continuity equation and Cauchy momentum equations in the absence of inertia

(2.1a,b)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u}=0,\quad\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\sigma}=\boldsymbol{0}.\end{equation}

To describe the viscoelastic behaviour of the fluid, we use the Oldroyd-B constitutive model (Oldroyd 1950), which represents the most simple combination of viscous and elastic stresses and is used widely to describe the flow of viscoelastic Boger fluids, characterized by a constant shear viscosity. The Oldroyd-B equation can be derived from microscopic principles by modelling polymer molecules as elastic dumbbells, which follow a linear Hooke's law for the restoring force as they are advected and stretched by the flow. The corresponding stress tensor $\boldsymbol {\sigma }$ is

(2.2)\begin{equation} \boldsymbol{\sigma}={-}p\boldsymbol{\mathsf{I}}+2\mu_{s}\boldsymbol{\mathsf{E}}+\boldsymbol{\tau}_{\!p},\end{equation}

where the first term on the right-hand side of (2.2) is the pressure contribution, the second term is the viscous stress contribution of a Newtonian solvent with a constant viscosity $\mu _{s}$, where $\boldsymbol{\mathsf{E}}=(\boldsymbol {\nabla }\boldsymbol {u}+(\boldsymbol {\nabla }\boldsymbol {u})^{\mathrm {T}})/2$ is the rate-of-strain tensor, and the last term, $\boldsymbol {\tau }_{p}$, is the polymer contribution. We note that $\boldsymbol{\mathsf{I}}$ in (2.2) is the identity tensor and T is the transpose operator on a tensor.

For the Oldroyd-B model, the polymer contribution to the stress tensor $\boldsymbol {\tau }_{p}$ can be expressed in terms of the (symmetric) conformation tensor (or the deformation of the microstructure) $\boldsymbol{\mathsf{A}}$ as (Bird et al. 1987; Larson 1988; Morozov & Spagnolie 2015)

(2.3)\begin{equation} \boldsymbol{\tau}_{\!p}=G(\boldsymbol{\mathsf{A}}-\boldsymbol{\mathsf{I}})=\frac{\mu_{p}}{\lambda}(\boldsymbol{\mathsf{A}}-\boldsymbol{\mathsf{I}}),\end{equation}

where $G$ is the elastic modulus, $\lambda$ is the relaxation time and $\mu _{p}=G \lambda$ is the polymer contribution to the shear viscosity at zero shear rate. It is also convenient to introduce the total zero-shear-rate viscosity $\mu _0=\mu _s+\mu _p$.

The evolution equation for the deformation of the microstructure $\boldsymbol{\mathsf{A}}$ of the Oldroyd-B model fluid is given at steady state as (Bird et al. 1987; Larson 1988; Morozov & Spagnolie 2015)

(2.4)\begin{equation} \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{\mathsf{A}}-(\boldsymbol{\nabla}\boldsymbol{u})^{\mathrm{T}} \boldsymbol{\cdot}\boldsymbol{\mathsf{A}}-\boldsymbol{\mathsf{A}}\boldsymbol{\cdot}(\boldsymbol{\nabla}\boldsymbol{u})={-}\frac{1}{\lambda}(\boldsymbol{\mathsf{A}}-\boldsymbol{\mathsf{I}}).\end{equation}

2.1. Scaling analysis and non-dimensionalization

We consider narrow configurations, in which $h(z)\ll \ell$, $h_{0}$ is the half-height at $z=0$, and $u_{c}=q/2h_{0}$ is the characteristic velocity scale set by the cross-sectionally averaged velocity. We introduce non-dimensional variables based on lubrication theory (Tichy 1996; Zhang et al. 2002; Saprykin et al. 2007; Ahmed & Biancofiore 2021; Boyko & Stone 2022)

(2.5a)\begin{gather} Z=\frac{z}{\ell},\quad Y=\frac{y}{h_{0}}, \quad U_{z}=\frac{u_{z}}{u_{c}},\quad U_{y}=\frac{u_{y}}{\epsilon u_{c}}, \end{gather}

(2.5b)\begin{gather}P=\frac{p}{\mu_{0}u_{c}\ell/h_{0}^{2}},\quad\Delta P=\frac{\Delta p}{\mu_{0}u_{c}\ell/h_{0}^{2}},\quad H=\frac{h}{h_{0}}, \end{gather}

(2.5c)\begin{gather}\tilde{A}_{zz}=\epsilon^{2}A_{zz},\quad\tilde{A}_{zy}=\epsilon A_{zy},\quad\tilde{A}_{yy}=A_{yy}, \end{gather}

where we have introduced the aspect ratio of the configuration, which is assumed to be small

(2.6)\begin{equation} \epsilon=\frac{h_{0}}{\ell}\ll1,\end{equation}

the contraction ratio

(2.7)\begin{equation} H_{\ell}=\frac{h_{\ell}}{h_{0}},\end{equation}

the viscosity ratios

(2.8a,b)\begin{equation} \tilde{\beta}=\frac{\mu_{p}}{\mu_{s}+\mu_{p}}=\frac{\mu_{p}}{\mu_{0}}\quad\mbox{and}\quad\beta=1-\tilde{\beta}=\frac{\mu_{s}}{\mu_{0}},\end{equation}

and the Deborah and Weissenberg numbers

(2.9a,b)\begin{equation} De=\frac{\lambda u_{c}}{\ell}\quad\mbox{and}\quad Wi=\frac{\lambda u_{c}}{h_{0}}.\end{equation}

For lubrication flows through the narrow geometries that we consider, there is a difference between the Deborah and Weissenberg numbers because of the two distinct length scales. The Weissenberg number $Wi$ is the product of the relaxation time scale of the fluid, $\lambda$, and the characteristic shear rate of the flow, $u_{c}/h_{0}$. On the other hand, the Deborah number $De$ is the ratio of the relaxation time, $\lambda$, to the residence time in the contraction region, $\ell /u_{c}$, or alternatively, the product of the relaxation time and the characteristic extensional rate of the flow (Tichy 1996; Zhang et al. 2002; Saprykin et al. 2007; Ahmed & Biancofiore 2021). The Deborah and Weissenberg numbers are related through $De=\epsilon Wi$, and for narrow geometries with $\epsilon \ll 1$, $De$ can be small while keeping $Wi =O(1)$.

Similar to our previous study (Boyko & Stone 2022), we non-dimensionalize the pressure using the total zero-shear-rate viscosity $\mu _{0}=\mu _{s}+\mu _{p}$. However, for convenience, we non-dimensionalize the height based on the entry height rather than the exit height. In addition, unlike our previous study, we do not scale the deformation of the microstructure with $De^{-1}$. Our current scaling is consistent with a fully developed unidirectional flow of an Oldroyd-B fluid in a straight channel, which yields $\tilde {A}_{zz}=O(De^2)$, $\tilde {A}_{zy}=O(De)$ and $\tilde {A}_{yy}=O(1)$; see (2.10d)–(2.10f) and (2.16). This scaling is convenient when considering arbitrary and large values of the Deborah number.

Note that, in both Hinch et al. (2024) and here, the channel height is $2h$, but the total flow rate per unit depth in the former is $2q$, whereas in this work it is $q$ as in Boyko & Stone (2022). All results are compatible because the variables used for the non-dimensionalization are the same, i.e. the expressions for the characteristic velocity, characteristic pressure and the Deborah number are the same.

2.2. Dimensionless lubrication equations in Cartesian coordinates

Using the non-dimensionalization (2.5)–(2.9a,b), to the leading order in $\epsilon$, the governing equations (2.1)–(2.4) take the form

(2.10a)\begin{gather} \frac{\partial U_{z}}{\partial Z}+\frac{\partial U_{y}}{\partial Y}=0, \end{gather}

(2.10b)\begin{gather}\frac{\partial P}{\partial Z}=(1-\tilde{\beta})\frac{\partial^{2}U_{z}}{\partial Y^{2}}+\frac{\tilde{\beta}}{De}\left(\frac{\partial\tilde{A}_{zz}}{\partial Z}+\frac{\partial\tilde{A}_{zy}}{\partial Y}\right), \end{gather}

(2.10c)\begin{gather}\frac{\partial P}{\partial Y}=0, \end{gather}

(2.10d)\begin{gather}U_{z}\frac{\partial \tilde{A}_{zz}}{\partial Z}+U_{y}\frac{\partial \tilde{A}_{zz}}{\partial Y}-2\frac{\partial U_{z}}{\partial Z}\tilde{A}_{zz}-2\frac{\partial U_{z}}{\partial Y}\tilde{A}_{zy}={-}\frac{1}{De}\tilde{A}_{zz}, \end{gather}

(2.10e)\begin{gather}U_{z}\frac{\partial \tilde{A}_{zy}}{\partial Z}+U_{y}\frac{\partial \tilde{A}_{zy}}{\partial Y}-\frac{\partial U_{y}}{\partial Z}\tilde{A}_{zz}-\frac{\partial U_{z}}{\partial Y}\tilde{A}_{yy}={-}\frac{1}{De}\tilde{A}_{zy}, \end{gather}

(2.10f)\begin{gather}U_{z}\frac{\partial \tilde{A}_{yy}}{\partial Z}+U_{y}\frac{\partial \tilde{A}_{yy}}{\partial Y}-2\frac{\partial U_{y}}{\partial Z}\tilde{A}_{zy}-2\frac{\partial U_{y}}{\partial Y}\tilde{A}_{yy}={-}\frac{1}{De}(\tilde{A}_{yy}-1). \end{gather}

From (2.10c), it follows that $P=P(Z)$, i.e. the pressure is independent of $Y$ up to $O(\epsilon ^2)$, consistent with the classical lubrication approximation. We note that the scaled $\tilde {A}_{zz}$ on the right-hand side of (2.10d) relaxes to $\epsilon ^2$, which is neglected at the leading order in $\epsilon$.

2.3. Orthogonal curvilinear coordinates for a slowly varying geometry

For our theoretical analysis, it is convenient to transform the geometry of the contraction from the Cartesian coordinates ($Z,Y$) to curvilinear coordinates ($\xi,\eta$), as illustrated in figure 2, with the mapping (Hinch et al. 2024)

(2.11a,b)\begin{equation} \xi=Z-\frac{1}{2}\epsilon^{2}\frac{H'(Z)}{H(Z)}(H(Z)^{2}-Y^{2})+O(\epsilon^{4}), \quad\eta=\frac{Y}{H(Z)}.\end{equation}

As shown in Appendix A, the curvilinear coordinates ($\xi,\eta$) are orthogonal with a relative error of $O(\epsilon ^{4})$, i.e. $\boldsymbol {\nabla }\xi \boldsymbol {\cdot }\boldsymbol {\nabla }\eta =O(\epsilon ^{4})$.

Figure 2. Schematic illustration of the orthogonal curvilinear coordinates ($\xi,\eta$) for a slowly varying geometry. The coordinate $\xi$ is constant along vertical grid lines, and $\eta$, defined in (2.11a,b), is constant along the curves going from left to right.

Hereafter, we use $\boldsymbol {u}=u\boldsymbol {e}_{\xi }+\upsilon \boldsymbol {e}_{\eta }$ and $\boldsymbol{\mathsf{A}}=A_{11}\boldsymbol {e}_{\xi }\boldsymbol {e}_{\xi }+A_{12}(\boldsymbol {e}_{\xi } \boldsymbol {e}_{\eta }+\boldsymbol {e}_{\eta }\boldsymbol {e}_{\xi })+A_{22} \boldsymbol {e}_{\eta }\boldsymbol {e}_{\eta }$ to denote, respectively, the components of velocity and deformation of the microstructure in curvilinear coordinates ($\xi,\eta$). The corresponding non-dimensional velocity components in different coordinates are related through (see Appendix A)

(2.12a,b)\begin{equation} U_{z}=U-\epsilon^{2}\eta H'(\xi) V, \quad U_{y}=\eta H'(\xi) U+V.\end{equation}

Similarly, the scaled conformation tensor components in different coordinates are related through (see Appendix A)

(2.13a)\begin{gather} \tilde{A}_{zz}=\tilde{A}_{11}+O(\epsilon^{2}), \end{gather}

(2.13b)\begin{gather}\tilde{A}_{zy}=\tilde{A}_{12}+\eta H'(\xi)\tilde{A}_{11}+O(\epsilon^{2}), \end{gather}

(2.13c)\begin{gather}\tilde{A}_{yy}=\tilde{A}_{22}+2\eta H'(\xi)\tilde{A}_{12}+\eta^2 (H'(\xi))^2\tilde{A}_{11}+O(\epsilon^{2}). \end{gather}

Finally, we note that, since there is only a $O(\epsilon ^{2})$ difference between the $\xi$- and $z$-directions, for convenience, we continue to use $Z$ rather than $\xi$ in curvilinear coordinates.

2.4. Dimensionless lubrication equations in orthogonal curvilinear coordinates

Using the mapping (2.11a,b), the governing equations (2.10) take the form in curvilinear coordinates (Hinch et al. 2024)

(2.14a)\begin{gather} \frac{\partial(HU)}{\partial Z}+\frac{\partial V}{\partial\eta}=0, \end{gather}

(2.14b)\begin{gather}\frac{\mathrm{d}P}{\mathrm{d}Z}=(1-\tilde{\beta})\frac{1}{H^{2}}\frac{\partial^{2}U}{\partial\eta^{2}}+\frac{\tilde{\beta}}{De}\left(\frac{1}{H}\frac{\partial(H\tilde{A}_{11})}{\partial Z}+\frac{1}{H}\frac{\partial\tilde{A}_{12}}{\partial\eta}\right), \end{gather}

(2.14c)\begin{gather}U\frac{\partial\tilde{A}_{11}}{\partial Z}+\frac{V}{H}\frac{\partial\tilde{A}_{11}}{\partial\eta}-2\frac{\partial U}{\partial Z}\tilde{A}_{11}-\frac{2}{H}\frac{\partial U}{\partial\eta}\tilde{A}_{12}={-}\frac{1}{De}\tilde{A}_{11}, \end{gather}

(2.14d)\begin{gather}U\frac{\partial\tilde{A}_{12}}{\partial Z}+\frac{V}{H}\frac{\partial\tilde{A}_{12}}{\partial\eta}-H\frac{\partial}{\partial Z}\left(\frac{V}{H}\right)\tilde{A}_{11}-\frac{1}{H}\frac{\partial U}{\partial\eta}\tilde{A}_{22}={-}\frac{1}{De}\tilde{A}_{12}, \end{gather}

(2.14e)\begin{gather}U\frac{\partial\tilde{A}_{22}}{\partial Z}+\frac{V}{H}\frac{\partial\tilde{A}_{22}}{\partial\eta}-2H\frac{\partial}{\partial Z}\left(\frac{V}{H}\right)\tilde{A}_{12}+2\frac{\partial U}{\partial Z}\tilde{A}_{22}={-}\frac{1}{De}(\tilde{A}_{22}-1). \end{gather}

The corresponding boundary conditions on the velocity are

(2.15a–d)\begin{equation} U(Z,1)=0,\quad V(Z,1)=0,\quad \frac{\partial U} {\partial\eta}(Z,0)=0,\quad H(Z)\int_{0}^{1}U(Z,\eta)\,\mathrm{d}\eta=1,\end{equation}

which represent, respectively, the no-slip and no-penetration boundary conditions along the channel walls, the symmetry boundary condition at the centreline and the integral mass conservation along the channel. In addition, we assume a fully developed unidirectional Poiseuille flow in the straight entry channel and the corresponding deformation of the microstructure

(2.16a–c)\begin{equation} \tilde{A}_{11}=\frac{18De^{2}}{H^{4}}\eta^{2},\quad \tilde{A}_{12}={-}\frac{3De}{H^{2}}\eta,\quad \tilde{A}_{22}=1,\end{equation}

with $H \equiv 1$ at the entrance. We also assume that, far downstream in the exit channel, the deformation of the microstructure attains a fully relaxed value, given by (2.16) with $H \equiv H_{\ell }$.

2.5. Pressure drop across the non-uniform region in the lubrication limit

In this subsection, we show that one can calculate the pressure drop without solving directly for the velocity field. To this end, we first integrate by parts the integral constraint (2.15d), repeatedly, using (2.15a) and (2.15c), e.g. (Hinch et al. 2024)

(2.17) \begin{align} \frac{1}{H(Z)}&=\int_{0}^{1}U\,\mathrm{d}\eta=\underset{0}{\underbrace{\left.\eta U\right|_{0}^{1}}}-\int_{0}^{1}\eta\frac{\partial U}{\partial\eta}\mathrm{d}\eta=\underset{0}{\underbrace{\left.\frac{1}{2}(1-\eta^{2})\frac{\partial U}{\partial\eta}\right|_{0}^{1}}}-\frac{1}{2}\int_{0}^{1}(1-\eta^{2})\frac{\partial^{2}U}{\partial\eta^{2}}\mathrm{d}\eta.\end{align}

Substituting the expression for $\partial ^{2}U/\partial \eta ^{2}$ from (2.14b) into (2.17), we obtain

(2.18)\begin{equation} -\frac{1-\tilde{\beta}}{H(Z)^3}=\frac{1}{2}\int_{0}^{1}(1-\eta^{2})\left[\frac{\mathrm{d}P}{\mathrm{d}Z}-\frac{\tilde{\beta}}{De}\left(\frac{1}{H}\frac{\partial(H\tilde{A}_{11})}{\partial Z}+\frac{1}{H}\frac{\partial\tilde{A}_{12}}{\partial\eta}\right)\right]\mathrm{d}\eta,\end{equation}

which can be rearranged to yield the pressure gradient

(2.19)\begin{equation} \frac{\mathrm{d}P}{\mathrm{d}Z}={-}\frac{3(1-\tilde{\beta})}{H(Z)^{3}}+\frac{3 \tilde{\beta}}{2De}\int_{0}^{1}(1-\eta^{2})\left[\frac{1}{H(Z)}\frac{\partial(H(Z)\tilde{A}_{11})}{\partial Z}+\frac{1}{H(Z)}\frac{\partial\tilde{A}_{12}}{\partial\eta}\right]\mathrm{d}\eta.\end{equation}

Integrating (2.19) with respect to $Z$ from 0 to 1 provides the pressure drop $\Delta P=P(0)-P(1)$ across the non-uniform region

(2.20)\begin{align} \Delta P&=3(1-\tilde{\beta})\int_{0}^{1}\frac{\mathrm{d}Z}{H(Z)^{3}} \nonumber\\ &\quad -\frac{3\tilde{\beta}}{2De}\int_{0}^{1}\int_{0}^{1}(1-\eta^{2})\left[\frac{1}{H(Z)} \frac{\partial(H(Z)\tilde{A}_{11})}{\partial Z}+\frac{1}{H(Z)}\frac{\partial\tilde{A}_{12}}{\partial\eta}\right]\mathrm{d}\eta\,\mathrm{d}Z. \end{align}

Using integration by parts, (2.20) can be expressed as

(2.21)\begin{align} \Delta P&= 3(1-\tilde{\beta})\int_{0}^{1}\frac{\mathrm{d}Z}{H(Z)^{3}}+\frac{3 \tilde{\beta}}{2De}\int_{0}^{1}(1-\eta^{2})\left[\tilde{A}_{11}(0,\eta)-\tilde{A}_{11}(1,\eta)\right] \mathrm{d}\eta \nonumber\\ &\quad -\frac{3\tilde{\beta}}{2De}\int_{0}^{1}\left[\frac{H'(Z)}{H(Z)}\left(\int_{0}^{1}(1-\eta^{2})\tilde{A}_{11}\, \mathrm{d}\eta\right)\right]\mathrm{d}Z \nonumber\\ &\quad -\frac{3\tilde{\beta}}{De}\int_{0}^{1}\left[\frac{1}{H(Z)}\int_{0}^{1}\eta\tilde{A}_{12}\,\mathrm{d}\eta \right]\mathrm{d}Z, \end{align}

where prime indicates a derivative with respect to $Z$.

Equation (2.21) resembles the result of an application of the reciprocal theorem previously derived for the pressure drop of the flow of an Oldroyd-B fluid in a slowly varying channel (Boyko & Stone 2021, 2022). The first term on the right-hand side of (2.21) represents the viscous contribution of the Newtonian solvent to the pressure drop. The second term represents the contribution of the elastic normal stress difference at the inlet and outlet of the non-uniform channel. The third term represents the contribution of the elastic normal stresses that arise due to the spatial variations in the channel shape, which is a contribution that is absent in a straight channel. Finally, the last term represents the elastic contribution due to shear stresses within the fluid domain of the non-uniform channel. It should be noted that we do not assume a priori the particular shape of the channel $H(Z)$ but rather consider a flow in a slowly varying channel of arbitrary shape $H(Z)$.

3. Low-$\tilde {\beta }$ lubrication analysis in a slowly varying region

In the previous section, we obtained the dimensionless equations (2.14), which are governed by the two non-dimensional parameters, $\tilde {\beta }$ and $De$, in the lubrication limit ($\epsilon \ll 1$). In this section, we derive analytical expressions for the velocity, conformation tensor and the $q-\Delta p$ relation for the pressure-driven flow of a very dilute viscoelastic Oldroyd-B fluid, $\tilde {\beta }=\mu _{p}/\mu _{0}\ll 1$ in a slowly varying channel of arbitrary shape $H(Z)$.

In contrast to our previous study that employed a low-Deborah-number lubrication analysis (Boyko & Stone 2022), in this work, we assume $De=O(1)$ and consider the ultra-dilute limit, $\tilde {\beta }\ll 1$ (see Remmelgas et al. 1999; Moore & Shelley 2012; Li et al. 2019; Mokhtari et al. 2022). To this end, we seek solutions of the form

(3.1)\begin{equation} \left(\begin{array}{c} U\\ V\\ P\\ \tilde{A}_{11}\\ \tilde{A}_{12}\\ \tilde{A}_{22} \end{array}\right)=\left(\begin{array}{c} U_{0}\\ V_{0}\\ P_{0}\\ \tilde{A}_{11,0}\\ \tilde{A}_{12,0}\\ \tilde{A}_{22,0} \end{array}\right)+\tilde{\beta}\left(\begin{array}{c} U_{1}\\ V_{1}\\ P_{1}\\ \tilde{A}_{11,1}\\ \tilde{A}_{12,1}\\ \tilde{A}_{22,1} \end{array}\right)+O(\epsilon^{2},\tilde{\beta}^{2}).\end{equation}

The ultra-dilute limit represents a one-way coupling between the velocity and pressure fields and the deformation of the microstructure (polymer stresses or conformation tensor). At leading order, the velocity and pressure are Newtonian, and the deformation of the microstructure (i.e. polymer stresses) arises from this Newtonian flow. Accordingly, the velocity and pressure at $O(\tilde {\beta })$ arise due to leading-order polymer stresses. In the next subsections, we provide closed-form asymptotic expressions for the velocity field and conformation tensor components at $O(\tilde {\beta }^{0})$ and the pressure drop up to $O(\tilde {\beta })$.

We note that the viscosity ratio $\tilde {\beta }=\mu _{p}/\mu _{0}$ is related to the so-called concentration of the polymers $c=\mu _{p}/\mu _{s}$ through $\tilde {\beta }=c/(c+1)$. Thus, at the leading order, the limits $\tilde {\beta }\ll 1$ and $c\ll 1$ are identical.

3.1. Velocity, conformation and pressure drop at the leading order in $\tilde {\beta }$

Substituting (3.1) into (2.14a)–(2.14b) and considering the leading order in $\tilde {\beta }$, the continuity and momentum equations take the form

(3.2a,b)\begin{equation} \frac{\partial(HU_{0})}{\partial Z}+\frac{\partial V_{0}}{\partial\eta}=0 \quad\hbox{and}\quad \frac{\mathrm{d}P_0}{\mathrm{d}Z}=\frac{1}{H^{2}}\frac{\partial^{2}U_{0}}{\partial\eta^{2}}, \end{equation}

subject to the boundary conditions

(3.3a–d)\begin{equation} U_{0}(Z,1)=0,\quad V_{0}(Z,1)=0,\quad \frac{\partial U_{0}}{\partial\eta}(Z,0)=0,\quad H(Z)\int_{0}^{1}U_{0}(Z,\eta)\,\mathrm{d}\eta=1.\end{equation}

The solutions for the axial velocity $U_0$ and the pressure drop $\Delta P_{0}$ at the leading order are well known (see, e.g. Boyko & Stone 2022)

(3.4a,b)\begin{equation} U_{0}=\frac{3}{2}\frac{1}{H(Z)}(1-\eta^{2})\quad\hbox{and}\quad \Delta P_{0}=3\int_{0}^{1}\frac{\mathrm{d}Z}{H(Z)^{3}}.\end{equation}

Substituting (3.4a) into the continuity equation (3.2a) and using (3.3b), yields

(3.5)\begin{equation} V_{0}\equiv0.\end{equation}

From (3.5), it follows that, in orthogonal curvilinear coordinates, the velocity in the $\eta$-direction is identically zero at $O(\tilde {\beta }^{0})$, in contrast to the Cartesian coordinates where $U_{y,0}=(3/2)H'(Z)Y(H(Z)^{2}-Y^{2})/H(Z)^{4}$. As we shall see, this fact significantly simplifies the theoretical analysis and allows us to derive closed-form expressions for the components of the conformation tensor.

Using (3.5), at leading order in $\tilde {\beta }$, the equations for the conformation tensor components, (2.14c)–(2.14e), simplify to

(3.6a)\begin{gather} U_{0}\frac{\partial\tilde{A}_{22,0}}{\partial Z}+2\frac{\partial U_{0}}{\partial Z}\tilde{A}_{22,0}={-}\frac{1}{De}(\tilde{A}_{22,0}-1), \end{gather}

(3.6b)\begin{gather}U_{0}\frac{\partial\tilde{A}_{12,0}}{\partial Z}-\frac{1}{H}\frac{\partial U_{0}}{\partial\eta}\tilde{A}_{22,0}={-}\frac{1}{De}\tilde{A}_{12,0}, \end{gather}

(3.6c)\begin{gather}U_{0}\frac{\partial\tilde{A}_{11,0}}{\partial Z}-2\frac{\partial U_{0}}{\partial Z}\tilde{A}_{11,0}-\frac{2}{H}\frac{\partial U_{0}}{\partial\eta}\tilde{A}_{12,0}={-}\frac{1}{De}\tilde{A}_{11,0}, \end{gather}

subject to the boundary conditions

(3.7a–c)\begin{equation} \tilde{A}_{11,0}(0,\eta)=18De^{2}\eta^{2},\quad\tilde{A}_{12,0}(0,\eta)={-}3De\eta, \quad\tilde{A}_{22,0}(0,\eta)=1.\end{equation}

Equations (3.6) represent a set of one-way coupled first-order semi-linear partial differential equations that can be solved first for $\tilde {A}_{22,0}$, followed by $\tilde {A}_{12,0}$ and then for $\tilde {A}_{11,0}$.

Solving (3.6) together with (3.7), we obtain closed-form expressions for $\tilde {A}_{22,0}$, $\tilde {A}_{12,0}$ and $\tilde {A}_{11,0}$ for arbitrary values of $De$ and the shape function $H(Z)$

(3.8)\begin{gather} \frac{\tilde{A}_{22,0}}{H(Z)^{2}}=\exp(\,{f(DeU_{0}(Z,\eta))}) \left[1+\int_{0}^{Z}\exp({-f(DeU_{0}(\tilde{Z},\eta))}) \frac{1}{DeU_{0}(\tilde{Z},\eta)H(\tilde{Z})^2}\mathrm{d}\tilde{Z}\right], \end{gather}

(3.9)\begin{gather}\frac{\tilde{A}_{12,0}}{({-}3De\eta)}=\exp(\,{f(DeU_{0}(Z,\eta))}) \left[1\!+\!\int_{0}^{Z}\exp({-f(DeU_{0}(\tilde{Z},\eta))}) \frac{\tilde{A}_{22,0}(\tilde{Z},\eta)}{DeU_{0}(\tilde{Z},\eta)H(\tilde{Z})^2}\mathrm{d}\tilde{Z}\right], \end{gather}

(3.10)\begin{align} &\frac{\tilde{A}_{11,0}}{18De^{2}\eta^{2}/H(Z)^2}\nonumber\\ &\quad =\exp(\,{f(DeU_{0}(Z,\eta))})\left[1+\int_{0}^{Z}\exp({-f(DeU_{0}(\tilde{Z},\eta))}) \frac{\tilde{A}_{12,0}(\tilde{Z},\eta)}{({-}3\eta De)DeU_{0}(\tilde{Z},\eta)}\mathrm{d}\tilde{Z}\right], \end{align}

where $f(DeU_{0}(Z,\eta ))$ is defined as

(3.11)\begin{equation} f(DeU_{0}(Z,\eta))={-}\int_{0}^{Z}\frac{1}{De U_{0}(\tilde{Z},\eta)}\mathrm{d}\tilde{Z}={-}\int_{0}^{Z}\frac{2H(\tilde{Z})}{3De(1-\eta^{2})} \mathrm{d}\tilde{Z}.\end{equation}

It is worth noting that the right-hand sides of (3.8)–(3.10) depend on the product $DeU_{0}(Z,\eta )$ and are not functions of $De$ and $\eta$ separately. Furthermore, (3.8)–(3.10) clearly show that, while the distribution of $\tilde {A}_{22,0}$ is set solely by the value at the beginning of the non-uniform region, the distribution of elastic shear and normal stresses, $\tilde {A}_{12,0}$ and $\tilde {A}_{11,0}$, are coupled to the transverse normal stress $\tilde {A}_{22,0}$. In fact, the elastic normal stress $\tilde {A}_{11,0}$ depends both on $\tilde {A}_{12,0}$ and $\tilde {A}_{22,0}$.

From (3.8)–(3.10), one might think that the conformation tensor components diverge at the wall ($\eta =\pm 1$). However, using (3.6) and noting that $U_{0}=\partial U_{0}/\partial Z=0$ at $\eta =\pm 1$, it follows that, at the walls of the non-uniform channel,

(3.12)\begin{equation} \tilde{A}_{22,0}^{wall}=1,\quad\tilde{A}_{12,0}^{wall}={\mp}\frac{3De}{H(Z)^{2}},\quad\tilde{A}_{11,0}^{wall}= \frac{18De^{2}}{H(Z)^{4}}\quad\mathrm{for\ all} \ De.\end{equation}

In §§ 3.1.1 and 3.1.2, we provide explicit expressions for the conformation tensor components in the low- and high-$De$ limits. We also note that the results shown in our figure 4(a,c) and the work of Hinch et al. (2024) suggest the existence of a viscoelastic boundary layer near the walls in the high-$De$ limit, which we analyse in § 3.1.3.

3.1.1. Conformation tensor in the low-$De$ limit

For $De\ll 1$, we solve the equations iteratively for the conformation tensor components (3.6) to obtain

(3.13a)\begin{align} \tilde{A}_{22,0}&=1+\frac{3DeH'}{H^{2}}(1-\eta^{2})+\frac{9De^{2}[4H'^{2}-HH'']}{2H^{4}}(1-\eta^{2})^{2}\nonumber\\ &\quad+ \frac{27De^{3}[24H'^{3}-13HH'H''+H^{2}H''']}{4H^{6}}(1-\eta^{2})^{3}, \end{align}

(3.13b)\begin{gather} \tilde{A}_{12,0}={-}\frac{3De}{H^{2}}\eta-\frac{18De^{2}H'}{H^{4}}\eta(1-\eta^{2}) - \frac{81De^{3}[4H'^{2}-HH'']}{2H^{6}}\eta(1-\eta^{2})^{2},\end{gather}

(3.13c)\begin{gather} \tilde{A}_{11,0}=\frac{18De^{2}}{H^{4}}\eta^{2}+\frac{162De^{3}H'}{H^{6}}\eta^{2}(1-\eta^{2}) +\frac{486De^{4}[4H'{}^{2}-HH'']}{H^{8}}\eta^{2}(1-\eta^{2})^{2}.\end{gather}

We note that the low-$De$ results (3.13) are consistent with our previous work (Boyko & Stone 2022), in which we provided explicit expressions for $\tilde {A}_{zz}$, $\tilde {A}_{zy}$ and $\tilde {A}_{yy}$ up to $O(De^2)$ in Cartesian coordinates. For example, using (2.13c) and (3.13), $\tilde {A}_{yy}$ can be expressed as $\tilde {A}_{yy}=1+3DeH'(Z)(H(Z)^2-3Y^2)/H(Z)^4+O(De^2)$, in agreement with (3.9a) in Boyko & Stone (2022).

3.1.2. Conformation tensor in the high-$De$ limit

We here provide the closed-form expressions for the conformation tensor components in the high-$De$ limit. We begin with the expression for $\tilde {A}_{22,0}$ and consider the core flow region.

For $De\gg 1$, except close to the wall, (3.6a) reduces to

(3.14)\begin{equation} U_{0}\frac{\partial\tilde{A}_{22,0}}{\partial Z}+2\frac{\partial U_{0}}{\partial Z}\tilde{A}_{22,0}=0,\end{equation}

whose solution subject to (3.7c) is

(3.15)\begin{equation} \tilde{A}_{22,0}(Z,\eta)=\tilde{A}_{22,0}(0,\eta)\frac{U_{0}(0,\eta)^{2}}{U_{0}(Z,\eta)^{2}}=H(Z)^{2}.\end{equation}

Next, since $\tilde {A}_{12,0}$ scales as $O(De)$ while $\tilde {A}_{22,0}$ is $O(1)$, within the core flow region in the high-$De$ limit we obtain that the first term in (3.6b) dominates over all the remaining terms

(3.16)\begin{equation} U_{0}\frac{\partial\tilde{A}_{12,0}}{\partial Z}=0,\end{equation}

so that elastic shear stresses preserve their value from the entry channel through the non-uniform region

(3.17)\begin{equation} \tilde{A}_{12,0}(Z,\eta)=\tilde{A}_{12,0}(0,\eta)={-}3De\eta.\end{equation}

Finally, to determine $\tilde {A}_{11,0}$, we note that the third and fourth terms in (3.6c) scale as $O(De)$, while the first and second terms are $O(De^{2})$. Thus, for $De\gg 1$, we expect the first and second terms to balance each other while the remaining terms are negligible, so that

(3.18)\begin{equation} U_{0}\frac{\partial\tilde{A}_{11,0}}{\partial Z}-2\frac{\partial U_{0}}{\partial Z}\tilde{A}_{11,0}=0.\end{equation}

Solving (3.18) subject to (3.7a) yields

(3.19)\begin{equation} \tilde{A}_{11,0}(Z,\eta)=\tilde{A}_{11,0}(0,\eta)\frac{U_{0}(Z,\eta)^{2}}{U_{0}(0,\eta)^{2}}=\frac{18De^{2}\eta^{2}}{H(Z)^{2}}.\end{equation}

In fact, for $De\gg 1$, there is a purely passive response of the microstructure, similar to a material line element, transported and deformed by the flow without relaxing.

The high-$De$ results (3.15), (3.17) and (3.19) can be also directly obtained from the closed-form solutions (3.8)–(3.10) by noting that, for $De\gg 1$, $\exp ({\pm f(DeU_{0}(Z,\eta ))})\approx 1$, and neglecting the $O(De^{-1})$ terms.

3.1.3. Boundary-layer analysis in the high-$De$ limit

In the previous section, we obtained analytical expressions for the components of the conformation tensor in the high-$De$ limit within the core flow region. However, these expressions do not hold near the walls, where a viscoelastic boundary layer of $O(De^{-1})$ thickness exists (Hinch et al. 2024). In this section, we analyse this boundary-layer region and provide boundary-layer equations and their closed-form solutions. To this end, we focus on the region $\eta \in [0,1]$, and introduce the rescaled inner-region coordinate

(3.20)\begin{equation} \zeta=De(1-\eta)=De\tilde{\eta}\quad\mathrm{for}\ \tilde{\eta}\ll1,\end{equation}

so that $De(1-\eta ^{2})=\zeta (2-\tilde {\eta })\approx 2\zeta$. Noting that, in the boundary layer, $\tilde {A}_{22,0}=O(1)$, $\tilde {A}_{12,0}=O(De)$ and $\tilde {A}_{11,0}=O(De^{2})$ (see (3.12)), to eliminate the dependence on $De$ in the governing equations and boundary conditions (3.7), we rescale $\tilde {A}_{22,0}$, $\tilde {A}_{12,0}$ and $\tilde {A}_{11,0}$, which are functions of $Z$ and $\zeta$, as

(3.21a–c)\begin{equation} \mathcal{A}_{22}=\frac{\tilde{A}_{22,0}}{H(Z)^{2}},\quad\mathcal{A}_{12}=\frac{\tilde{A}_{12,0}}{({-}3\eta De)},\quad\mathcal{A}_{11}=\frac{\tilde{A}_{11,0}}{18\eta^{2}De^{2}/H(Z)^{2}}.\end{equation}

Substituting (3.20) and (3.21a–c) into (3.6) and using (3.4a), we obtain the boundary-layer equations in the high-$De$ limit

(3.22a)\begin{gather} \frac{3\zeta}{H(Z)}\frac{\partial\mathcal{A}_{22}}{\partial Z}={-}\left(\mathcal{A}_{22}-\frac{1}{H(Z)^{2}}\right), \end{gather}

(3.22b)\begin{gather}\frac{3\zeta}{H(Z)}\frac{\partial\mathcal{A}_{12}}{\partial Z}={-}(\mathcal{A}_{12}-\mathcal{A}_{22}), \end{gather}

(3.22c)\begin{gather}\frac{3\zeta}{H(Z)}\frac{\partial\mathcal{A}_{11}}{\partial Z}={-}(\mathcal{A}_{11}-\mathcal{A}_{12}), \end{gather}

subject to the inlet conditions

(3.23a–c)\begin{equation} \mathcal{A}_{11}(0,\zeta)=1,\quad\mathcal{A}_{12}(0,\zeta)=1,\quad\mathcal{A}_{22}(0,\zeta)=1. \end{equation}

Solving (3.22) together with (3.23), we obtain closed-form expressions for $\mathcal {A}_{22}$, $\mathcal {A}_{12}$ and $\mathcal {A}_{11}$ in the boundary-layer region

(3.24a)\begin{gather} \mathcal{A}_{22}=\mathrm{e}^{\mathcal{F}(Z,\zeta)}\left[1+\int_{0}^{Z}\mathrm{e}^{-\mathcal{F}(\tilde{Z},\zeta)}\frac{1}{3\zeta H(\tilde{Z})}\mathrm{d}\tilde{Z}\right], \end{gather}

(3.24b)\begin{gather}\mathcal{A}_{12}=\mathrm{e}^{\mathcal{F}(Z,\zeta)}\left[1+\int_{0}^{Z}\mathrm{e}^{-\mathcal{F}(\tilde{Z},\zeta)}\frac{\mathcal{A}_{22}(\tilde{Z},\zeta)H(\tilde{Z})}{3\zeta}\mathrm{d}\tilde{Z}\right], \end{gather}

(3.24c)\begin{gather}\mathcal{A}_{11}=\mathrm{e}^{\mathcal{F}(Z,\zeta)}\left[1+\int_{0}^{Z}\mathrm{e}^{-\mathcal{F}(\tilde{Z},\zeta)}\frac{\mathcal{A}_{12}(\tilde{Z},\zeta)H(\tilde{Z})}{3\zeta}\mathrm{d}\tilde{Z}\right], \end{gather}

where $\mathcal {F}(Z,\zeta )$ is defined as

(3.25)\begin{equation} \mathcal{F}(Z,\zeta)={-}\frac{1}{3\zeta}\int_{0}^{Z}H(\tilde{Z})\,\mathrm{d}\tilde{Z}.\end{equation}

We note that solutions (3.24) satisfy the matching conditions between the inner and outer regions. Specifically, $\mathcal {A}_{22}|_{\zeta \rightarrow \infty }=[\tilde {A}_{22,0}^{core}/H(Z)^{2}]_{\eta =1}=1$, $\mathcal {A}_{12}|_{\zeta \rightarrow \infty }=[\tilde {A}_{12,0}^{core}/(-3\eta De)]_{\eta =1}=1$ and $\mathcal {A}_{11}|_{\zeta \rightarrow \infty }=[\tilde {A}_{11,0}^{core}/(18\eta ^{2}De^{2}/H(Z)^{2})]_{\eta =1}=1$.

3.2. Pressure drop at the first order in $\tilde {\beta }$

Equation (2.20) shows that the pressure drop depends on the elastic normal and shear stresses $\tilde {A}_{11}$ and $\tilde {A}_{12}$, and thus, generally, requires the solution of the nonlinear viscoelastic problem. However, in the ultra-dilute limit, corresponding to $\tilde {\beta }=\mu _{p}/\mu _{0}\ll 1$, we can determine the pressure drop at $O(\tilde {\beta })$ for arbitrary values of $De$ only with the knowledge of the velocity field and conformation tensor components at $O(1)$. Specifically, substituting (3.1) into (2.20) yields at $O(\tilde {\beta })$ the pressure drop $\Delta P_1$,

(3.26)\begin{align} \Delta P_1&={-}3\int_{0}^{1}\frac{\mathrm{d}Z}{H(Z)^{3}} \nonumber\\ &\quad -\frac{3}{2De}\int_{0}^{1}\int_{0}^{1}(1-\eta^{2})\left[\frac{1}{H(Z)}\frac{\partial(H(Z)\tilde{A}_{11,0})}{\partial Z}+\frac{1}{H(Z)}\frac{\partial\tilde{A}_{12,0}}{\partial\eta}\right]\mathrm{d}\eta\,\mathrm{d}Z, \end{align}

or alternatively

(3.27)\begin{align} \Delta P_{1}&={-}3\int_{0}^{1}\frac{ \mathrm{d}Z}{H(Z)^{3}}+\frac{3}{2De}\int_{0}^{1}(1-\eta^{2}) \left[\tilde{A}_{11,0}(0,\eta)-\tilde{A}_{11,0}(1,\eta)\right]\mathrm{d}\eta \nonumber\\ &\quad -\frac{3}{2De}\int_{0}^{1}\left[\frac{H'(Z)}{H(Z)}\left(\int_{0}^{1}(1-\eta^{2})\tilde{A}_{11,0}\,\mathrm{d}\eta\right) \right]\mathrm{d}Z\nonumber\\ &\quad -\frac{3}{De}\int_{0}^{1}\left[\frac{1}{H(Z)}\int_{0}^{1}\eta\tilde{A}_{12,0}\,\mathrm{d}\eta\right]\mathrm{d}Z. \end{align}

Thus, for a given flow rate $q$, the dimensionless pressure drop $\Delta P=\Delta p/(\mu _{0}q\ell /2h_{0}^{3})$, as a function of the shape function $H(Z)$, the Deborah number $De$ and the viscosity ratio $\tilde {\beta }\ll 1$, up to $O(\tilde {\beta })$, is given by

(3.28)\begin{equation} \Delta P=\Delta P_{0}(H(Z))+\tilde{\beta}\Delta P_{1}(De,H(Z))+O(\epsilon^{2},\tilde{\beta}^{2}),\end{equation}

where the expressions for $\Delta P_{0}$ and $\Delta P_{1}$ are given in (3.4b) and (3.27), respectively.

Notably, in contrast to our previous results for the pressure drop obtained in the weakly viscoelastic and lubrication limits with $De\ll 1$ and $\tilde {\beta }\in [0,1]$ (Boyko & Stone 2022), the current result (3.28) applies to the limit of $\tilde {\beta }\ll 1$, while allowing $De=O(1)$.

3.2.1. Pressure drop at $O(\tilde {\beta })$ in the low-$De$ limit

To calculate the pressure drop $\Delta P_{1}$ at low Deborah numbers in the non-uniform shape region, we use (3.13b)–(3.13c) and (3.27). The elastic normal stress (NS) contribution to the pressure drop at $O(\tilde {\beta })$ is

(3.29)\begin{align} \Delta P_{1}^{NS}&=\frac{3}{2De}\int_{0}^{1}(1\!-\!\eta^{2})\left[\tilde{A}_{11,0}\right]_{Z=1}^{Z=0}\, \mathrm{d}\eta -\frac{3}{2De}\int_{0}^{1}\left[\frac{H'(Z)}{H(Z)}\left(\int_{0}^{1} (1\!-\!\eta^{2})\tilde{A}_{11,0}\,\mathrm{d}\eta\right)\right]\mathrm{d}Z \nonumber\\ &= \frac{27}{10}De(1-H_{\ell}^{{-}4})\quad\mathrm{for}\ De\ll1, \end{align}

where $[\tilde {A}_{11,0}]_{Z=1}^{Z=0}=\tilde {A}_{11,0}(0,\eta )-\tilde {A}_{11,0}(1,\eta )$.

The elastic shear stress (SS) contribution to the pressure drop at $O(\tilde {\beta })$ is

(3.30)\begin{align} \Delta P_{1}^{SS}&={-}\frac{3}{De}\int_{0}^{1}\left[\frac{1}{H(Z)}\int_{0}^{1}\eta\tilde{A}_{12,0}\, \mathrm{d}\eta\right]\mathrm{d}Z \nonumber\\ &= 3\int_{0}^{1}\frac{\mathrm{d}Z}{H(Z)^{3}}+\frac{18}{10}De(1-H_{\ell}^{{-}4})\quad\mathrm{for}\ De\ll1. \end{align}

Substituting (3.29) and (3.30) into (3.27) provides the pressure drop at $O(\tilde {\beta })$ in the low-$De$ limit up to $O(De)$

(3.31)\begin{equation} \Delta P_{1}=\frac{9}{2}De(1-H_{\ell}^{{-}4})+O(De^2)\quad\mathrm{for}\ De \ll1,\end{equation}

so that the total pressure drop across the non-uniform channel in the low-$De$ limit, accounting for the leading-order effect of viscoelasticity, is

(3.32) \begin{align} \Delta P&=\underset{\textit{Solvent stress}}{\underbrace{3(1-\tilde{\beta})\int_{0}^{1} \frac{\mathrm{d}Z}{H(Z)^{3}}}}+\underset{\textit{Elastic shear stress}}{\underbrace{ 3\tilde{\beta}\int_{0}^{1}\frac{\mathrm{d}Z}{H(Z)^{3}}+\frac{18}{10}\tilde{\beta} De(1-H_{\ell}^{{-}4})}}+\underset{\textit{Elastic normal stress}}{\underbrace{ \frac{27}{10}\tilde{\beta}De(1-H_{\ell}^{{-}4})}} \nonumber\\ &=3\int_{0}^{1}\frac{\mathrm{d}Z}{H(Z)^{3}}+\frac{9}{2}\tilde{\beta}De(1-H_{\ell}^{{-}4})+ O(De^2)\quad\mathrm{for}\ De\ll1, \end{align}

in agreement with the results of our previous work (Boyko & Stone 2022). The three terms on the right-hand side of (3.32) represent, respectively, the Newtonian solvent stress contribution, the elastic shear stress contribution and the elastic normal stress contribution to the pressure drop.

3.2.2. Pressure drop at $O(\tilde {\beta })$ in the high-$De$ limit

To calculate the pressure drop $\Delta P_{1}$ at high Deborah numbers in the non-uniform region, we use (3.17), (3.19) and (3.27). The elastic normal and shear stress contributions to the pressure drop at $O(\tilde {\beta })$ are

(3.33a,b)\begin{equation} \Delta P_{1}^{NS}=\frac{9}{5}De(1-H_{\ell}^{{-}2})\quad\mathrm{and}\quad \Delta P_{1}^{SS}=3\int_{0}^{1}\frac{\mathrm{d}Z}{H(Z)}\quad\mathrm{for}\ De\gg1.\end{equation}

Substituting (3.33) into (3.27) yields the pressure drop at $O(\tilde {\beta })$ in the high-$De$ limit

(3.34)\begin{equation} \Delta P_{1}={-}3\int_{0}^{1}\frac{\mathrm{d}Z}{H(Z)^{3}}+3\int_{0}^{1} \frac{\mathrm{d}Z}{H(Z)}+\frac{9}{5}De(1-H_{\ell}^{{-}2})\quad\mathrm{for}\ De\gg1,\end{equation}

so that the total pressure drop across the non-uniform channel in the high-$De$ limit is

(3.35) \begin{equation} \Delta P=\underset{\textit{Solvent stress}}{\underbrace{3(1-\tilde{\beta}) \int_{0}^{1}\frac{\mathrm{d}Z}{H(Z)^{3}}}}+\underset{\textit{Elastic shear stress}}{ \underbrace{3\tilde{\beta}\int_{0}^{1}\frac{\mathrm{d}Z}{H(Z)}}}+\underset{\textit{Elastic normal stress}}{ \underbrace{\frac{9}{5}\tilde{\beta}De(1-H_{\ell}^{{-}2})}}\quad\mathrm{for}\ De\gg1.\end{equation}

Similar to the low-$De$ limit, for the contraction geometry, the last term, corresponding to the elastic normal stress contribution, leads to a decrease in the pressure drop, which is linear in the Deborah number. As noted by Hinch et al. (2024), the tension in the streamlines at the end of the contraction pulls the flow through the contraction, thus requiring less pressure to push. Furthermore, at high Deborah numbers, the elastic shear stresses are lower than the fully relaxed value $\tilde {A}_{12}=-3De\eta /H_{\ell }^{2}$ due to insufficient time (distance) to approach their fully relaxed value in the contraction. Thus, the elastic shear stress contribution to the pressure drop, $3\tilde {\beta }\int _{0}^{1}H(Z)^{-1}\,\mathrm {d}Z$, is smaller than the steady Poiseuille value of $3\tilde {\beta }\int _{0}^{1}H(Z)^{-3}\,\mathrm {d}Z$, further reducing the pressure drop. Finally, we note that the result (3.35) also holds for the expansion geometry $H_{\ell }>1$, in which the two physical mechanisms mentioned above lead to an increase in the pressure drop.

4. Low-$\tilde {\beta }$ lubrication analysis in the exit channel

In this section, we analyse the spatial relaxation of the elastic stresses and the pressure drop in the uniform exit channel. From examining the expressions (3.8)–(3.10) for the conformation tensor, when there are no longer shape changes, we expect the elastic stresses and the pressure in the exit channel to relax exponentially, with a strong dependence on $De^{-1}$. Thus, for higher Deborah numbers, a longer downstream section is required (Keiller 1993) for polymer relaxation, consistent with previous numerical simulations using the Oldroyd-B model (Debbaut et al. 1988; Alves et al. 2003).

Following similar steps as in the previous section, in Appendix B, we derive closed-form expressions for the conformation tensor and the pressure drop in the uniform exit channel for arbitrary values of the Deborah number. Furthermore, we provide analytical expressions for the conformation tensor and the pressure drop in the low- and high-$De$ limits. We summarize in table 1 the semi-analytical solutions and low- and high-$De$ asymptotic expressions for the deformation of the microstructure and the pressure drop of the Oldroyd-B fluid in a contraction and exit channel in the ultra-dilute limit derived in this work.

Table 1. A summary of the semi-analytical solutions and low- and high-$De$ asymptotic expressions for the deformation of the microstructure and the pressure drop of the Oldroyd-B fluid in a contraction and exit channel in the ultra-dilute limit.

In particular, we show that the total pressure drop in the exit channel can be expressed as

(4.1) \begin{align} \Delta P_{\ell}=\underset{\textit{Solvent stress}}{\underbrace{(1\!-\!\tilde{\beta})\frac{3L}{H_{\ell}^{3}}}}+\underset{\textit{Elastic normal stress}}{\underbrace{\frac{3\tilde{\beta}}{2De}\int_{0}^{1}(1\!-\!\eta^{2})\left[\tilde{A}_{11,0}\right]_{Z_{\ell}=L}^{Z_{\ell}=0}\,\mathrm{d}\eta}} +\underset{\textit{Elastic shear stress}}{\underbrace{\frac{3\tilde{\beta}}{DeH_{\ell}}\int_{0}^{1} \eta\left[\int_{L}^{0}\tilde{A}_{12,0}\,\mathrm{d}Z_{\ell}\right]\mathrm{d}\eta}},\end{align}

where $L=\ell _{\ell }/\ell$ is the dimensionless length, $H_{\ell }=H(Z=1)=h_{\ell }/h_{0}$ is the dimensionless height of the exit channel, $Z_{\ell }=Z-1$, $\tilde {A}_{11,0}$ and $\tilde {A}_{12,0}$ are given in (B4) and (B5) and $[\tilde {A}_{11,0}]_{Z_{\ell }=L}^{Z_{\ell }=0}=\tilde {A}_{11,0}(Z_{\ell }=0,\eta )-\tilde {A}_{11,0}(Z_{\ell }=L,\eta )$.

It should be noted that we can express the first-order contribution $\Delta P_{\ell, 1}$ in terms of the difference between the conformation tensor components at the beginning and end of the exit channel (see Appendix B and Hinch et al. 2024)

(4.2)\begin{align} \Delta P_{\ell,1}&=\frac{3}{2De}\int_{0}^{1}(1-\eta^{2})\left[\tilde{A}_{11,0}\right]_{Z_{\ell}=L}^{Z_{\ell}=0}\, \mathrm{d}\eta -\frac{9}{2H_{\ell}^2}\int_{0}^{1}\eta(1-\eta^2)\left[\tilde{A}_{12,0}\right]_{Z_{\ell}=L}^{Z_{\ell}=0}\, \mathrm{d}\eta \nonumber\\ &\quad +\frac{27De}{2H_{\ell}^{4}}\int_{0}^{1}\eta^2(1-\eta^2)\left[\tilde{A}_{22,0}\right]_{Z_{\ell}=L}^{Z_{\ell}=0}\,\mathrm{d}\eta. \end{align}

Hereafter, we assume that the length of the exit channel, $L$, is such that the elastic stresses reach their fully relaxed values by the end of the exit channel, given by (2.16) with $H \equiv H_{\ell }$. Under this assumption, (4.2) clearly shows that the first-order contribution $\Delta P_{\ell, 1}$ is independent of $L$ since the steady-state values of $\tilde {A}_{11,0}$, $\tilde {A}_{12,0}$ and $\tilde {A}_{22,0}$ depend solely on the $\eta$ coordinate. Note, however, that the total pressure in the exit channel depends on $L$ via $\Delta P_{\ell }=3L/H_{\ell }^{3}+\tilde {\beta }\Delta P_{\ell,1}$.

In addition, we show in Appendix B that the total pressure drop in the exit channel in the low- and high-$De$ limits is

(4.3)\begin{gather} \Delta P_{\ell}= \frac{3L}{H_{\ell}^{3}}-\frac{1728\tilde{\beta}De^{3}H''(1)}{35H_{\ell}^{7}}\quad\mathrm{for}\ De\ll1, \end{gather}

(4.4)\begin{gather}\Delta P_{\ell}= \frac{3L}{H_{\ell}^{3}}+\frac{36}{5}\tilde{\beta}De(H_{\ell}^{{-}2}-H_{\ell}^{{-}4})\quad\mathrm{for}\ De\gg1. \end{gather}

From (4.3) and (4.4), it follows that, similar to the contraction, the pressure drop in the exit channel decreases with $De$. Furthermore, the physical mechanisms responsible for the pressure drop reduction are the same in both the contraction and the exit channels.

The asymptotic result (4.4) is obtained using expressions (B9a)–(B9c), which hold in the high-$De$ limit within the core flow region. As discussed above, near the walls, there exists a viscoelastic boundary layer of thickness $O(De^{-1})$. Nevertheless, this boundary layer will contribute only a small $O(\tilde {\beta } De^{-1})$ correction to the pressure drop in the exit channel for $De\gg 1$, as noted by Hinch et al. (2024).

5. Results

In this section, we present the theoretical results for the pressure drop and conformation tensor distribution of the Oldroyd-B fluid in the ultra-dilute limit developed in §§ 3 and 4. As an illustrative example, we specifically consider the case of a smooth contraction of the form

(5.1)\begin{equation} H(Z)=1-(1-H_{\ell})Z^{2}(2-Z)^2\quad 0\leq Z\leq 1,\end{equation}

where $H_{\ell }=H(1)/H(0)=h_{\ell }/h_{0}$ is the ratio of the exit to entry heights; for the contracting geometry we have $H_{\ell }<1$. This contraction shape function is illustrated in figure 2 and satisfies $H'(0)=0$ and $H'(1)=H'''(1)=0$.

In this work, we present the results for $H_{\ell }=0.5$ and $\tilde {\beta }=0.05$. While the current study focuses only on one contraction ratio, in our previous work, we considered four contraction ratios, in which the elastic normal stresses vary by almost two decades (Hinch et al. 2024). In addition, figure 8 of our previous paper shows a 0.1 % difference between $c=0.1$ and $c=0.05$ for the pressure drop in the contraction at $De=0.8$. Nevertheless, our current analysis allows us to analyse slowly varying arbitrarily shaped channels provided $\epsilon \ll 1$ and $\tilde {\beta }\ll 1$. To obtain the semi-analytical solutions for given values of $De$ and $H_{\ell }$, we first used MATLAB's routine $\texttt{cumtrapz}$ to find the conformation tensor components, given in (3.8)–(3.10) and (B3)–(B5), for a contraction and exit channel. Typical values of the grid size were $\Delta Z = 10^{-4}$ and $\Delta \eta = 0.005$. We then used MATLAB's routine $\texttt{trapz}$ to calculate the pressure drop, (3.28) and (4.1), for a contraction and exit channel, respectively.

5.1. Streamwise variation of elastic stresses in the contraction and exit channel

We present in figure 3 the streamwise variation of the leading-order elastic stresses, scaled by their entry values, on $\eta =0.5$ in contraction and exit channels for $De=0.01$ (a,d), $De=0.1$ (b,e) and $De=1$ (c,f). As expected, for a small Deborah number of $De=0.01$, the elastic stresses achieve their downstream fully relaxed values by the end of contraction (figure 3a), and thus we observe very little variation in the relaxation along the exit channel (figure 3d). Consistent with the low-$De$ asymptotic solutions (3.13), represented by cyan dotted lines, for $H_{\ell }=0.5$, the elastic shear and axial normal stresses increase by a factor of 4 and 16, respectively, while the transverse normal stress preserves its entry value.

Figure 3. The streamwise variation of leading-order elastic stresses on $\eta =0.5$ in a smooth contraction and exit channel in the ultra-dilute limit. (a–c) Scaled elastic stresses $\tilde {A}_{11,0}/(18De^2\eta ^2)$, $\tilde {A}_{12,0}/(-3De\eta )$ and $\tilde {A}_{22,0}$ in the contraction as a function of $Z$ for ($a$) $De=0.01$, ($b$) $De=0.1$ and ($c$) $De=1$. (d–e) Scaled elastic stresses in the exit channel $\tilde {A}_{11,0}/(18De^2\eta ^2)$, $\tilde {A}_{12,0}/(-3De\eta )$ and $\tilde {A}_{22,0}$ as a function of $Z_{\ell }$ for ($d$) $De=0.01$, ($e$) $De=0.1$ and ($\,f$) $De=1$. Solid lines represent the semi-analytical solutions (3.8)–(3.10) (contraction) and (B3)–(B5) (exit channel). Cyan dotted lines represent the low-$De$ asymptotic solutions (3.13) (contraction) and (B7) (exit channel). Red dashed lines represent the high-$De$ asymptotic solutions (3.15), (3.17) and (3.19) (contraction) and (B9) (exit channel). All calculations were performed using $H_{\ell }=0.5$.

For the case of $De=0.1$, shown in figure 3(b,e), the elastic stresses do not have enough residence time to attain their downstream steady-state values in the contraction. Therefore, there is a significant spatial relaxation in the exit channel. Interestingly, although the relaxation in the exit channel is governed mainly by $\exp ({-2H_{\ell }Z_{\ell }/[3De(1-\eta ^{2})]})$ (see (B3)–(B5)), the elastic stresses relax over slightly different length scales, with the shortest relaxation distance required for $\tilde {A}_{22,0}$ and the longest for $\tilde {A}_{11,0}$. The latter behaviour is associated with the nature of the coupling between the elastic stresses so that $\tilde {A}_{11,0}$ depends both on $\tilde {A}_{12,0}$ on $\tilde {A}_{22,0}$, while $\tilde {A}_{12,0}$ depends only on $\tilde {A}_{22,0}$ (see (B3)–(B5)).

When $De=1$, it is evident from figure 3(c) that, at the end of the contraction, the axial normal stress increases by a factor of $1/H_{\ell }^2=4$, the transverse normal stress is squashed by a factor of $H_{\ell }^2=1/4$, and the elastic shear stress preserves its entry value. Figure 3(f) presents the spatial relaxation of the elastic stresses in the exit channel for $De=1$, clearly showing that a very long exit channel is required to attain the downstream fully relaxed values of all stresses ($L>16$ for $\eta =0.5$). Furthermore, we observe excellent agreement between the semi-analytical results (solid lines) and the high-$De$ asymptotic solutions (3.15), (3.17), (3.19) and (B9) (dashed red lines). Such an agreement for $De=1$ is consistent with recent results of Hinch et al. (2024), who found that the high-$De$ analysis works well for $De>0.4$.

The closed-form solutions for the conformation tensor components, (B3)–(B5), clearly show that the spatial relaxation of the elastic stresses in the exit channel strongly depends on the stresses at the end of the contraction ($Z = 1$). Therefore, it is of particular interest to elucidate the behaviour of the elastic stresses at the end of the contraction and the extent to which they are perturbed relative to their downstream fully relaxed values.

The solid lines in figure 4(a,c) present the elastic shear ($a$) and axial normal stresses ($c$) at the end of the contraction as a function of $\eta =y/H_{\ell }$ for $De=0.01,0.1,1$ and 10, scaled by their downstream fully relaxed values. For a small Deborah number of $De=0.01$, $\tilde {A}_{12,0}(Z=1,\eta )/(-3De\eta /H_{\ell }^2)$ and $\tilde {A}_{11,0}(Z=1,\eta )/(18De^2\eta ^2/H_{\ell }^4)$ only slightly differ from their downstream values, and this behaviour is well captured by the low-$De$ asymptotic solutions (3.13b)–(3.13c), represented by cyan dotted lines. As $De$ increases, the elastic stresses become considerably suppressed within the core flow region relative to their eventual relaxed values far downstream, and for $De=1$ and $De=10$, the elastic shear and axial normal stresses approach the high-$De$ asymptote of $H_{\ell }^2=1/4$, represented by red dashed lines. Furthermore, in the high-$De$ limit, we observe the presence of a viscoelastic boundary layer close to the walls, where the elastic stresses reach their downstream fully relaxed values.

Figure 4. The cross-stream variation of leading-order elastic shear and normal stresses at the end of the contraction in the ultra-dilute limit. (a,c) Scaled elastic shear and normal stresses at the end of the contraction, (a) $\tilde {A}_{12,0}(Z=1,\eta )/(-3De\eta /H_{\ell }^2)$ and (c) $\tilde {A}_{11,0}(Z=1,\eta )/(18De^2\eta ^2/H_{\ell }^4)$, as a function of $\eta$ for $De=0.01,0.1,1$ and 10, respectively; (b) $\tilde {A}_{12,0}(Z=1,\eta )/(-3De\eta /H_{\ell }^2)$ and (d) $\tilde {A}_{11,0}(Z=1,\eta )/(18De^2\eta ^2/H_{\ell }^4)$ as a function of the rescaled coordinate $\zeta =De(1-\eta )$ for $De=0.1,1$ and 10. Solid lines represent the semi-analytical solutions (3.9)–(3.10). Cyan dotted lines represent the low-$De$ asymptotic solutions (3.13b)–(3.13c). Red dashed lines represent the high-$De$ asymptotic solutions (3.17) and (3.19). Green dashed lines represent the boundary-layer solutions (3.24b)–(3.24c). All calculations were performed using $H_{\ell }=0.5$.

To provide insight into this viscoelastic boundary layer, we replot in figure 4(b,d) the elastic shear ($b$) and axial normal stresses ($d$) at the end of the contraction as a function of the rescaled coordinate $\zeta =De(1-\eta )$ for $De=0.1,1$ and 10 (see § 3.1.3). It is evident from figures 4(b) and 4(d) that this rescaling collapses the results for the different Deborah numbers onto the same curves, which are the boundary-layer asymptotic solutions (3.24b) and (3.24c) (green dashed lines). Clearly, for $De = 1$ and $De = 10$, which are graphically almost indistinguishable, there is excellent agreement between the semi-analytical results and the boundary-layer asymptotic solutions, thus confirming the thickness of a boundary layer as $O(De^{-1})$.

Furthermore, examining (3.8)–(3.10), we infer that their right-hand sides are not a function of $De$ and $\eta$ separately but depend on the product $DeU_{0}(Z,\eta )$. To test this prediction, we show in figure 5(a,b) the scaled elastic shear ($a$) and axial normal stresses ($b$) at the end of the contraction, ($a$) $\tilde {A}_{12,0}(Z=1,\eta )/(-3De\eta /H_{\ell }^2)$ and ($b$) $\tilde {A}_{11,0}(Z=1,\eta )/(18De^2\eta ^2/H_{\ell }^4)$ minus $H_{\ell }^2$, divided by the factor $1-H_{\ell }^2$, as a function of $DeU_{0}(Z=1,\eta )$ for $De=0.5, 1$ and $H_{\ell }=0.125,0.25,0.5$. We observe that the results for two different values of $De$ approximately collapse onto the same curve across three contraction ratios.

Figure 5. (a,b) Scaled elastic shear and normal stresses at the end of the contraction, ($a$) $\tilde {A}_{12,0}(Z=1,\eta )/(-3De\eta /H_{\ell }^2)$ and ($b$) $\tilde {A}_{11,0}(Z=1,\eta )/(18De^2\eta ^2/H_{\ell }^4)$ minus $H_{\ell }^2$, divided by the factor $1-H_{\ell }^2$, as a function of $DeU_{0}(Z=1,\eta )$ for $De=0.5, 1$ and $H_{\ell }=0.125,0.25$ and 0.5. This rescaling leads to an approximate collapse of the results on the single uniform curve for different Deborah numbers and contraction ratios.

5.2. Pressure gradient relaxation in the exit channel

It follows from figure 3(d–f) in the previous subsection that, as $De$ increases, there is a significant relaxation of the elastic stresses in the exit channel, which occurs over a long distance. Specifically, the elastic stresses relax exponentially over a distance which is proportional to the centreline velocity $(3/2H_{\ell })$ multiplied by the Deborah number $De$ (see (B3)–(B5)). For this reason, a longer downstream section is required at higher $De$.

In this subsection, we study the relaxation of the pressure gradient in the downstream section. Substituting $H(Z)=H_{\ell }$ into (2.19) yields the pressure gradient in the exit channel

(5.2)\begin{equation} \frac{\mathrm{d}P}{\mathrm{d}Z}={-}\frac{3(1-\tilde{\beta})}{H_{\ell}^{3}}+\frac{3 \tilde{\beta}}{2De}\int_{0}^{1}(1-\eta^{2})\frac{\partial\tilde{A}_{11,0}}{\partial Z}\mathrm{d}\eta+\frac{3 \tilde{\beta}}{H_{\ell}De}\int_{0}^{1}\eta\tilde{A}_{12,0}\,\mathrm{d}\eta+O(\tilde{\beta}^2).\end{equation}

Noting that in the exit channel $U_{0}=(3/2H_{\ell })(1-\eta ^{2})$ and $\mathrm {d} U_{0}/\mathrm {d}\eta =-(3/H_{\ell })\eta$, and using the expression for $U_{0}\partial \tilde {A}_{11,0}/\partial Z$ from (B2c), (5.2) can be written as

(5.3)\begin{equation} \left(\frac{\mathrm{d}P}{\mathrm{d}Z}+\frac{3}{H_{\ell}^{3}}\right)\frac{1}{\tilde{\beta}}= \frac{3}{H_{\ell}^{3}}-\frac{H_{\ell}}{De^2}\int_{0}^{1}\tilde{A}_{11,0}\,\mathrm{d}\eta-\frac{3 }{H_{\ell}De}\int_{0}^{1}\eta\tilde{A}_{12,0}\,\mathrm{d}\eta,\end{equation}

where the right-hand side is independent of $\tilde {\beta }$.

We present in figure 6(a) the relaxation of the scaled pressure gradient $(\mathrm {d}P/\mathrm {d}Z+3/H_{\ell }^3)/\tilde {\beta }$ as a function of the downstream distance $Z_{\ell }$ for $De=0.02,0.2,1$ and 2. Similar to elastic stresses, the scaled pressure gradient relaxes exponentially over the downstream distance, which significantly increases with $De$. Furthermore, we observe a good agreement between the low- and high-$De$ asymptotic solutions (cyan dotted and red dashed lines) and the semi-analytical results (solid lines).

Figure 6. The spatial relaxation of the pressure gradient for the Oldroyd-B fluid in the uniform exit channel of a contraction in the ultra-dilute limit. ($a$) Scaled pressure gradient $(\mathrm {d}P/\mathrm {d}Z+3/H_{\ell }^3)/\tilde {\beta }$ as a function of the downstream distance $Z_{\ell }$ for $De=0.02,0.2,1$ and 2. ($b$) Scaled pressure gradient $(\mathrm {d}P/\mathrm {d}Z+3/H_{\ell }^3)/\tilde {\beta }$ as a function of the rescaled downstream distance $2H_{\ell }Z_{\ell }/3De$ in a log$-$linear plot. Solid lines represent the semi-analytical solutions obtained from (5.3) using (B3)–(B5). Cyan dotted lines represent the low-$De$ asymptotic solutions obtained from (5.3) using (B7). Red dashed lines represent the high-$De$ asymptotic solutions obtained from (5.3) using (B9). The green dashed line is $100\exp ({-2H_{\ell }Z_{\ell }/3De})$. All calculations were performed using $H_{\ell }=0.5$.

Recalling that the elastic stresses relax exponentially over a distance proportional to $(3De/2H_{\ell })$, we replot in figure 6(b) the scaled pressure gradient, (5.3), as a function of the rescaled downstream distance $2H_{\ell }Z_{\ell }/3De$ in a log$-$linear plot. As a result, all curves become parallel to the green dashed line $100\exp ({-2H_{\ell }Z_{\ell }/3De})$, thus confirming that the pressure gradient relaxes over a length scale ${\sim }(3De/2H_{\ell })$, similar to the elastic stresses. More specifically, it follows from figure 6(b) that the downstream distance over which the scaled pressure gradient (PG) decays to 1 % of its maximum value, $L_{1\,\%}^{PG}$, is approximately

(5.4)\begin{equation} L_{1\,\%}^{PG}\approx (5.3 \pm 0.5)\times \frac{3De}{2H_{\ell}},\end{equation}

where we obtain that the prefactor $5.3 \pm 0.5$ is weakly dependent on $De$ throughout the investigated range of Deborah numbers. Equation (5.4) and the scaling $3De/2H_{\ell }$ indicate that, in the exit channel, the appropriate Deborah number is based on the exit height, i.e. $De_{exit}=\lambda q/2h_{\ell }\ell =De/H_{\ell }$.

We note that our estimate of the length of the downstream section, (5.4), is consistent with previous numerical studies on the viscoelastic flows in 2-D abrupt contractions (Debbaut et al. 1988; Alves et al. 2003). Specifically, (5.4) predicts $L_{1\,\%}^{PG}\approx 239 \pm 23$ for $De_{exit}=De/H_{\ell }=30$, which should be contrasted with 250 of Debbaut et al. (1988), who studied numerically the flow through the planar $4\,{:}\,1$ contraction.

5.3. Pressure drop in the contraction and exit channel

In this subsection, we study the pressure drop across the contraction and the exit channel. First, in figure 7(a) we present the non-dimensional pressure drop $\Delta P=\Delta p/(\mu _{0}q\ell /2h_{0}^{3})$ in the contraction as a function of $De=\lambda q/(2\ell h_{0})$ for $H_{\ell }=0.5$ and $\tilde {\beta }=0.05$. For further clarification, figure 7(b) shows the first-order contribution $\Delta P_{1}=\Delta p_{1}/(\mu _{0}q\ell /2h_{0}^{3})$ as a function of $De=\lambda q/(2\ell h_{0})$, which is independent of $\tilde {\beta }$. Black dots represent the semi-analytical solution (3.28), cyan dotted lines represent the low-$De$ asymptotic solution (3.32) and red dashed lines represent the high-$De$ asymptotic solution (3.35). Clearly, there is excellent agreement between our low- and high-$De$ asymptotic solutions and the semi-analytical results. We also validate the predictions of our semi-analytical and asymptotic results against the 2-D finite-element simulations with $H_{\ell }=0.5$, $\tilde {\beta }=0.05$ and $\epsilon =0.02$ (grey triangles), showing very good agreement. The details of the numerical implementation in the finite-element software COMSOL Multiphysics are provided in Boyko & Stone (2022).

Figure 7. Non-dimensional pressure drop for the Oldroyd-B fluid in a contracting channel in the ultra-dilute limit. ($a$) Dimensionless pressure drop $\Delta P=\Delta p/(\mu _{0}q\ell /2h_{0}^{3})$ as a function of $De=\lambda q/(2\ell h_{0})$ for $\tilde {\beta }=0.05$. ($b$) First-order contribution $\Delta P_{1}=\Delta p_{1}/(\mu _{0}q\ell /2h_{0}^{3})$ to the dimensionless pressure drop as a function of $De=\lambda q/(2\ell h_{0})$. Grey triangles in ($a$) represent the results of the finite-element simulation. Black dots represent the semi-analytical solution (3.28). Cyan dotted lines represent the low-$De$ asymptotic solution (3.32). Red dashed lines represent the high-$De$ asymptotic solution (3.35). All calculations were performed using $H_{\ell }=0.5$.

It is evident that the semi-analytical solution for the pressure drop in the contraction approaches the high-$De$ asymptotic solution for $De\gtrsim 0.4$ and linearly decreases with the Deborah number. First, such an agreement for $De\gg \!\!\!\!\!\!/~1$ is consistent with our results for the elastic stresses, shown in figure 3, and recent results of Hinch et al. (2024). Second, and more importantly, from the excellent agreement between the semi-analytical results and the high-$De$ asymptotic solution, based on the components of the conformation tensor within the core flow region, we conclude that the viscoelastic boundary layer near the walls makes a negligible contribution to the pressure drop in the contracting channel.

Next, in figure 8(a) we present the non-dimensional pressure drop $\Delta P_{\ell }$ in the exit channel as a function of $De$ for $H_{\ell }=0.5$, $\tilde {\beta }=0.05$, and $L=50$. For $De=2$, a long exit channel of $L\gtrsim 30$ is required to reach the full relaxation of the elastic stresses and pressure gradient, consistent with (5.4). Figure 8(b) shows the first-order contribution $\Delta P_{\ell,1}$ as a function of $De$, which is independent of $\tilde {\beta }$. In contrast to the total pressure drop $\Delta P_{\ell }$, the first-order contribution $\Delta P_{\ell,1}$ does not depend on $L$, as shown in (4.2), provided that $L$ is sufficiently long so that by the end of the exit channel the elastic stresses have achieved their fully relaxed values (2.16) with $H \equiv H_{\ell }$.

Figure 8. Non-dimensional pressure drop for the Oldroyd-B fluid in the exit channel of a contraction in the ultra-dilute limit. (a) Dimensionless pressure drop $\Delta P_{\ell }=\Delta p_{\ell }/(\mu _{0}q\ell /2h_{0}^{3})$ as a function of $De=\lambda q/(2\ell h_{0})$ for $\tilde {\beta }=0.05$ and $L=50$. (b) First-order contribution $\Delta P_{\ell,1}=\Delta p_{\ell,1}/(\mu _{0}q\ell /2h_{0}^{3})$ to the dimensionless pressure drop as a function of $De=\lambda q/(2\ell h_{0})$. Black dots represent the semi-analytical solutions (4.1) ($\Delta P_{\ell }$ in (a)) and (4.2) ($\Delta P_{\ell,1}$ in (b)). The cyan dotted curve represents the low-$De$ asymptotic solution (4.3). Red dashed lines represent the high-$De$ asymptotic solution (4.4). The inset in (a) shows a comparison of semi-analytical predictions (black dots) and finite-element simulation results (grey triangles) for $\Delta P_{\ell }-\Delta P_{\ell,0}=\tilde {\beta } \Delta P_{\ell,1}$ as a function of $De$ for $\tilde {\beta }=0.05$ and $L=5$. The inset in (b) shows $\Delta P_{\ell }-\Delta P_{\ell,0}=\tilde {\beta } \Delta P_{\ell,1}$ as a function of $De$ for $\tilde {\beta }=0.05$ in range of $1\leq De\leq 10$. All calculations were performed using $H_{\ell }=0.5$.

The inset in figure 8(a) shows a comparison of our semi-analytical predictions (black dots) and finite-element simulation results (grey triangles) for $\Delta P_{\ell }-\Delta P_{\ell,0}=\tilde {\beta } \Delta P_{\ell,1}$ as a function of $De$ for $H_{\ell }=0.5$, $\tilde {\beta }=0.05$ and $L=5$. We observe excellent agreement between the semi-analytical and numerical results. In addition, the low-$De$ asymptotic solution (cyan dotted curve) accurately captures the numerical results for $De<0.05$ and indicates that the pressure drop in the exit channel scales as $De^3$ for $De\ll 1$.

Similar to the contraction, the pressure drop in the exit channel linearly decreases with $De$ for $De\gtrsim 0.3$, as shown in figure 8. While our semi-analytical solution linearly diminishes with the slope of $-36/5$, as predicted by the high-$De$ asymptotic solution (red dashed lines), there is an offset between the two results for $\tilde {\beta } \Delta P_{\ell,1}$. In particular, for $De=0.4$, we have a non-negligible relative error of approximately 30 %. However, the inset in figure 8(b) shows that as $De$ increases, the agreement between our semi-analytical solution and the high-$De$ asymptotic prediction significantly improves, resulting in relative errors of only approximately 5 % and 1 % for $De=2$ and $De=10$, respectively.

We note that our theoretical approach, based on the ultra-dilute limit, allows us to study the behaviour of the elastic stresses and pressure drop at arbitrary values of $De$. In particular, we can predict the behaviour in the high-Deborah-number regime, for example, $De=2$ and even $De=10$, which we are currently unable to access via finite-element simulations. Note, however, that we have assumed steady flows, so further investigation would be required to assess whether there might be flow instabilities at higher $De$.

5.4. Different contributions to the pressure drop in the contraction and exit channel

In the previous subsection, we observed a monotonic reduction in the dimensionless pressure drop with increasing $De$ for an Oldroyd-B fluid flowing through the contraction and exit channel (figures 7 and 8). To understand the source of such pressure drop reduction, we elucidate the relative importance of elastic contributions to the pressure drop.

The elastic contributions to the non-dimensional pressure drop across the contraction and exit channel, scaled by $\tilde {\beta }$, as a function of $De$ are shown in figures 9(a) and 9(b), respectively. Black circles and grey dots represent the elastic shear and normal stress contributions obtained from the semi-analytical solutions (3.28) and (4.1). Cyan dotted and purple curves represent the elastic shear and normal stress contributions obtained from the low-$De$ asymptotic solutions (3.32) and (4.3). Red and black dashed lines represent the elastic shear and normal stress contributions obtained from the high-$De$ asymptotic solutions (3.35) and (4.4). As expected based on our previous results, we observe excellent agreement between our low- and high-$De$ asymptotic solutions and the semi-analytical predictions.

Figure 9. Elastic contributions to the non-dimensional pressure drop of the Oldroyd-B fluid, scaled by $\tilde {\beta }$, in (a) the contraction and (b) the exit channel in the ultra-dilute limit. Black circles and grey dots represent the semi-analytical solutions (3.28) (contraction) and (4.1) (exit channel) for elastic shear and normal stress contributions. Cyan dotted and purple curves represent the low-$De$ asymptotic solutions (3.32) (contraction) and (4.3) (exit channel) for elastic shear and normal stress contributions. Red and black dashed lines represent the high-$De$ asymptotic solutions (3.35) (contraction) and (4.4) (exit channel) for elastic shear and normal stress contributions. All calculations were performed using $H_{\ell }=0.5$ and $L=50$.

The first main source for the pressure drop reduction is the elastic normal stress contribution, which linearly decreases with $De$ in the contraction and exit channel at low and high Deborah numbers. As noted by Hinch et al. (2024), this is because the elastic normal stresses, which correspond to the tension in the streamlines, are higher at the end of the contraction (exit channel) compared with the beginning of the contraction (exit channel). These higher elastic normal stresses pull the fluid along and thus require less pressure to push.

The second main source for the pressure drop reduction is the decrease of elastic shear stress contribution with $De$ due to the long time (or long distance) required for the elastic shear stresses to approach their eventual relaxed values far downstream. As a result, the elastic shear stresses are lower than the fully relaxed value $\tilde {A}_{12}=-3De\eta /H_{\ell }^{2}$ (see figure 3), and their contribution to the pressure drop is smaller than the steady Poiseuille value of $3\tilde {\beta }\int _{0}^{1}H(Z)^{-3}\,\mathrm {d}Z$ (contraction) and $3\tilde {\beta }L/H_{\ell }^3$ (exit channel), thus reducing the pressure drop. At low Deborah numbers, such a decrease scales as $De$ and $De^3$ for a smooth contraction and exit channel, respectively. However, at high Deborah numbers, it approaches a constant asymptotic value of $3\tilde {\beta }\int _{0}^{1}H(Z)^{-1}\,\mathrm {d}Z$ for the contraction. For the exit channel, $\Delta P_{\ell,1}^{SS}$ linearly depends on the Deborah number since the relaxation of the elastic shear stresses occurs over the distance $L$, which scales linearly with $De$, as shown in (5.4).

6. Concluding remarks

In this work, we applied the lubrication approximation and considered the ultra-dilute limit to study the flow of an Oldroyd-B fluid in arbitrarily shaped contracting channels. Specifically, we exploited the one-way coupling between the parabolic velocity and polymer conformation tensor in the ultra-dilute limit to derive closed-form expressions for the microstructure deformation and the flow rate–pressure drop relation for arbitrary values of the Deborah number. We provided analytical expressions for the conformation tensor and the $q-\Delta p$ relation in the low- and high-Deborah-number limits for the contraction and exit channels, complementing the asymptotic results of Boyko & Stone (2022) and the analysis of Hinch et al. (2024) at any concentration. We further analysed the viscoelastic boundary layer of thickness $O(De^{-1})$, existing near the walls at high Deborah numbers, and derived the boundary-layer asymptotic solutions. We validated our semi-analytical and asymptotic results for the pressure drop in the smooth contraction and exit channels with 2-D finite-element numerical simulations and found excellent agreement.

For both contraction and exit channels, the pressure drop of an Oldroyd-B fluid monotonically decreases with increasing $De$ and scales linearly with $De$ at high Deborah numbers, as shown in figures 7 and 8. We identified two mechanisms for such pressure drop reduction (see figure 9). The first is higher elastic normal stresses at the end of the contraction and exit channels, relative to the corresponding entry values, that pull the fluid along and thus require less pressure to push. The second source for the pressure drop reduction is because, once perturbed from their upstream values, the elastic shear stresses require a long distance to approach their new downstream fully relaxed values, as shown in figure 3, so again reducing the pressure drop.

Our theoretical approach, which relies on lubrication theory and the ultra-dilute limit, allows us to study the behaviour of the elastic stresses and pressure drop of an Oldroyd-B fluid at arbitrary values of $De$. Our theory is not restricted to the case of 2-D contracting channels and can be utilized to study different slowly varying geometries, such as expansions and constrictions. The approach can also be extended to axisymmetric geometries. Furthermore, the theoretical framework we presented enables us to access sufficiently high Deborah numbers, which are difficult and sometimes impossible to study via numerical simulations due to the high-Weissenberg-number problem (Owens & Phillips 2002; Alves et al. 2021). We, therefore, believe that our analytical and semi-analytical results for the ultra-dilute limit are of fundamental importance as they may serve for simulation validation.

Finally, we note that our theoretical predictions for the pressure drop reduction of an Oldroyd-B fluid in a contraction are consistent with the previous numerical reports on 2-D abruptly contracting geometries (Aboubacar, Matallah & Webster 2002; Alves et al. 2003; Binding et al. 2006; Aguayo, Tamaddon-Jahromi & Webster 2008). However, these predictions are opposite to the experiments showing a nonlinear increase in the pressure drop with $De$ for the flow of a Boger fluid through abrupt axisymmetric contraction–expansion and contraction geometries (Rothstein & McKinley 1999, 2001; Nigen & Walters 2002; Sousa et al. 2009). As noted by Alves et al. (2003) and Hinch et al. (2024), this discrepancy might be attributed to the lack of dissipative effects in the Oldroyd-B model. Thus, as a future research direction, it is interesting to study more complex constitutive equations, such as a finitely extensible nonlinear elastic (FENE) model introduced by Chilcott & Rallison (1988) (FENE-CR) and a finitely extensible nonlinear elastic model with the Peterlin approximation (FENE-P), that incorporate dissipation and additional microscopic features of polymer solutions and understand how these features affect the pressure drop. We anticipate that even for a more complex constitutive model, the theoretical framework presented here will enable the development of a simplified, reduced-order theory, allowing us to study the behaviour at non-small Deborah numbers.