2.1. Oldroyd-B

The Oldroyd-B model fluid is the simplest viscoelastic fluid, combining a simple viscous stress and a simple elastic stress:

\begin{equation*} {\boldsymbol{\sigma}} =-p\boldsymbol{I} + 2\mu_0\boldsymbol{e} + G\boldsymbol{A}, \end{equation*}

where $\mu _0$ is a viscosity, $\boldsymbol {e}$ is the strain-rate of the flow, $G$ is an elastic modulus, and $\boldsymbol {A}$ is the deformation of the microstructure. The viscosity $\mu _0$ will be called the ‘solvent viscosity’, as it occurs in the bead-and-spring model of dilute polymer solutions. The microstructure deforms due to velocity differences in the flow as represented by the velocity gradient ${\boldsymbol {\nabla }}\boldsymbol {u}$, and simultaneously relaxes to the isotropic state on a relaxation time $\tau$:

\begin{equation*} \frac{{\rm D}\boldsymbol{A}}{{\rm D} t} = \boldsymbol{A}\boldsymbol{\cdot}{\boldsymbol{\nabla}}\boldsymbol{u} + {\boldsymbol{\nabla}}\boldsymbol{u}^{\rm T}\boldsymbol{\cdot}\boldsymbol{A} - \frac{1}{\tau}\,(\boldsymbol{A} - \boldsymbol{I}). \end{equation*}

The left-hand side and first two terms of the right-hand side are the so-called upper-convective time derivative $\mathop{\boldsymbol{A}}\limits^{\nabla}$ of the Oldroyd-B equation, appropriate for fibrous microstructures. A modified right-hand side $-\boldsymbol {A}\boldsymbol {\cdot }{\boldsymbol {\nabla }}\boldsymbol {u}^{\rm T} - {\boldsymbol {\nabla }}\boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {A}$ corresponds to the lower-convective time derivative $\mathop{\boldsymbol{A}}\limits^{\triangle}$, which produces the Oldroyd-A equation, appropriate to disk-like microstructures.

The simple fluid model has just three parameters, $\mu _0$, $G$ and $\tau$, to be fitted to experimental data. It represents fairly well many aspects, but not all, of the flow behaviour of some polymer solutions, such as Boger fluids. For these two reasons, the Oldroyd-B equation is a sensible choice for numerical calculations of flows. Finally, the elastic-dumbbell model of polymer solutions is governed by the Oldroyd-B equation.

In a uniform steady simple shear flow $\boldsymbol {u}=(\gamma y,0,0)$, the Oldroyd-B fluid has viscous-like shear stresses

\begin{equation*} \sigma_{xy} = \sigma_{yx} = (\mu_0 + G\tau)\gamma.\end{equation*}

In the steady state, the microstructural deformation $A_{xy} = \gamma \tau$ should be thought of as the deformation-rate $\gamma$ multiplied by the time $\tau$ of how long the material can remember being deformed before it relaxes. With the steady deformation proportional to the shear-rate $\gamma$, both the elastic and viscous shear stresses are proportional to $\gamma$, with the combination described by an effective viscosity $\mu _* = \mu _0 +G\tau$. We shall call this effective viscosity the steady shear viscosity of the Oldroyd-B fluid.

In addition to shear stresses, there are so-called normal stresses

\begin{equation*} \sigma_{xx} =-p + 2G(\gamma\tau)^2,\quad \sigma_{yy} = \sigma_{zz} =- p. \end{equation*}

The difference $\sigma _{xx}-\sigma _{yy} = 2G(\gamma \tau )^2$ corresponds to an elastic tension in the streamlines. In order for this tension in the streamlines to play a role in the lubrication dynamics, it is necessary for the normal stresses to be much larger than the shear stresses. We see that this is possible if the Weissenberg number $We = \gamma \tau$ is very large. We shall need $We$ to be as large as the ratio of the length to the width of the contraction. There is a question of whether the common experimental Boger fluids keep a quadratic variation in shear-rate of the normal stresses at the high shear-rates that we consider theoretically.

To help us to understand the response of the Oldroyd-B fluid as it flows through a contraction, it is useful to consider the start up and shut down of a simple shear flow $\boldsymbol {u} = (\gamma (t)\,y,0,0)$, where

\begin{equation*} \gamma(t) = \begin{cases} 0 & \text{for}\ t<0,\\ \gamma_1 & \text{for}\ 0 \le t < t_1, \\ 0 & \text{for}\ t_1 \le t, \end{cases} \end{equation*}

with the duration of the shear much longer than the relaxation time, $t_1\gg \tau$. The shear stress response is

\begin{equation*} \sigma_{xy}(t) = \begin{cases} 0 & \text{for}\ t<0,\\ \mu_0\gamma_1 + G\gamma_1\tau(1 - \exp({-t/\tau})) & \text{for}\ 0 \le t < t_1, \\ G\gamma_1\tau \exp({-(t-t_1)/\tau}) & \text{for}\ t_1 \le t. \end{cases} \end{equation*}

We see that when the shearing is switched on, there is an immediate viscous stress $\mu _0\gamma _1$. As the microstructure begins to deform, an elastic component of the stress builds up. After a few relaxation times $t \gg \tau$, the deformation reaches the equilibrium value $\gamma _1\tau$, which is the deformation-rate $\gamma _1$ times the memory time $\tau$, and the shear stress takes a steady value $(\mu _0+G\tau )\gamma _1$. When the shearing is switched off at $t_1$, the viscous component $\mu _0\gamma _1$ immediately drops out, while the elastic component relaxes over several relaxation times.

The system of stresses in a simple shear flow can be understood by considering the deformation of a microstructure composed of fibres that deform like material line-elements. A fibre parallel to the flow, $\boldsymbol {\delta \ell } = (\delta \ell,0,0)$, will move with the flow without change in size or orientation. A fibre perpendicular to the flow, $\boldsymbol {\delta \ell } = (0,\delta \ell,0)$, will shear to $(\delta \ell \gamma t, \delta \ell, 0)$, i.e. the component perpendicular to the flow remains unchanged, while a component parallel to the flow grows linearly in time.

The microstructure in an Oldroyd-B fluid is described by a second-order tensor $\boldsymbol {A}$. We think of this tensor as the tensor product of two vector fibres. Thus $A_{xx}=1$ corresponds to two fibres of unit length parallel to the flow, which are unchanged in the simple shear flow, so $A_{xx}=1$ remains constant. On the other hand, $A_{yy}=1$ corresponds to two fibres perpendicular to the flow, each of which will develop in the simple shear components parallel to the flow, $A_{xy}=A_{yx} = \gamma t$, while leaving the component perpendicular to the flow, $A_{yy}=1$, unchanged. Finally, these shear components $A_{xy}$ and $A_{yx}$ correspond to one fibre parallel and one fibre perpendicular. Through the shear, the perpendicular fibres develop components parallel to the flow, so producing $A_{xx} = 1 + 2(\gamma t)^2$.

The above discussion is for a uniform steady simple shear. While the flow through a contraction flow is a shearing flow with one layer of fluid sliding over another, it is not steady in the Lagrangian view moving with the fluid. As the flow accelerates, a fibre in the direction of the flow will stretch proportional to the increase in velocity. (In a steady flow, a material line-element evolves according to $\boldsymbol {u}\boldsymbol {\cdot }{\boldsymbol {\nabla }}\delta {\boldsymbol {\ell }} = \delta {\boldsymbol {\ell }}\boldsymbol {\cdot }{\boldsymbol {\nabla }}\boldsymbol {u}$. If the fibre $\delta {\boldsymbol {\ell }}$ starts parallel to $\boldsymbol {u}$, then the solution is $\delta {\boldsymbol {\ell }}\propto \boldsymbol {u}$.) Thus we can expect the component of $\boldsymbol {A}$ in the direction of the flow, corresponding to two fibres parallel to the flow, to increase in proportion to the square of the velocity. As fibres parallel to the flow stretch, fibres perpendicular to the flow must compress by the same proportion in order to conserve mass in a planar flow. (Axisymmetric flows would differ at this point.) Thus we can expect the component of $\boldsymbol {A}$ perpendicular to the flow to decrease proportional to the square of the velocity, while the shearing components of $\boldsymbol {A}$, corresponding to one fibre parallel and one fibre perpendicular, should be unchanged as the flow accelerates. Note that this discussion of the deformation of the microstructure has ignored its relaxation.

The previous paragraph gives the basis of the high-$De$ analysis in § 3.4.

The idea that the microstructure stretches like a material line-element in regions where the relaxation is negligible was suggested by Rallison & Hinch (Reference Rallison and Hinch1988) in § 5.2 when considering how dumbbells behave in rapid sink flows. This stretching was expressed by Hinch (Reference Hinch1993) at the beginning of § 3 as $\boldsymbol {A} = f(\psi )\,\boldsymbol {u}\boldsymbol {u}$, where $f$ is constant along a streamline, when considering flow around a sharp corner. Keiller (Reference Keiller1993a) added in his equation (18) the shear components, when improving a numerical method for entry flows. Finally, Renardy (Reference Renardy1994) in his equation (3) gave a decomposition of the deformation for steady planar flows into three components, parallel, shear and perpendicular to the flow. We shall return to this decomposition much later, in § 6.3.

2.2. Lubrication equations

We consider a flow in the $x$-direction through a planar contraction starting at $x=0$ and finishing at $x=\ell$. The width of the channel is $2h(x)$, with a straight centreline $y=0$ and rigid boundaries at $y = \pm h(x)$. There is a short inlet section in $x<0$ with a constant width $2h(0)$, and a long outlet section in $x>\ell$ with a constant width $2h(\ell )$. For the channel to be slowly varying, the width is taken to be much smaller than the length:

\begin{equation*} \epsilon = h(0)/\ell \ll 1. \end{equation*}

We shall work only to leading order in $\epsilon$, and not exhibit where any correction terms can occur. A second non-dimensional parameter of the geometry is the contraction ratio

\begin{equation*} H = h(0)/h(\ell). \end{equation*}

In our numerical studies, we shall take $H = 2^{1/2}$, $2$, $2^{3/2}$ and $4$. Between each case in this sequence, the shear-rate doubles and the normal stresses quadruple.

In the standard lubrication scalings, we non-dimensionalise the distance $x$ along the channel by the length $\ell$ of the contraction, and the distance $y$ across the channel as well as the width $h(x)$ by the width of the inlet $h_0= h(0)$. With a volume flux $2Q$, we non-dimensionalise the component $u$ of the flow along the channel by $Q/h_0$, and the much smaller component of the flow across the channel by $Q/\ell$. Time $t$ is non-dimensionalised by the residence time in the contraction $\ell h_0/Q$, and pressure $p$ by the pressure drop of a Newtonian viscous fluid with the solvent viscosity $\mu _0$ by $\mu _0 Q \ell /h_0^3$.

With these scalings, a non-dimensional parameter arises that measures the elasticity:

\begin{equation*} c = G\tau/\mu_0. \end{equation*}

In the bead-and-spring model of dilute polymer solutions, this $c$ would be the concentration of the polymers. For most of the numerical studies, we shall take $c=1$, which would not qualify as dilute.

We shall ignore inertia. This means taking the appropriate reduced Reynolds number to be negligibly small:

\begin{equation*} Re = \rho Q h_0/\mu_0 \ell \ll 1, \end{equation*}

where $\rho$ is the density of the fluid.

A final non-dimensional parameter is needed to measure the relaxation time $\tau$ in the constitutive equation. In most flows, there is no difference between a Weissenberg number and a Deborah number. An exception is in lubrication flows. The Weissenberg number $We = Q\tau /h_0^2$ measures the relaxation time compared with the local shear-rate, which we shall not use. As discussed in the previous section § 2.1, we need the Weissenberg number to be as large as $\epsilon ^{-1}$, in order to bring the tension in the streamlines into the leading-order dynamics. For this reason, we use the Deborah number

\begin{equation*} De = \epsilon\,We = Q\tau/\ell h_0. \end{equation*}

This Deborah number is the ratio of the relaxation time to the residence time in the contraction.

In order to promote the non-Newtonian tension in the streamlines to a leading-order role, we scale differently the different components of the deformation of the microstructure: $A_{yy}$ by $\epsilon ^0$, $A_{xy}$ by $\epsilon ^{-1}$, and $A_{xx}$ by $\epsilon ^{-2}$.

With these non-dimensionalisations, the governing equations become the following. The conservation of mass is

\begin{equation*} {\frac{\partial u}{\partial x}} + {\frac{\partial v}{\partial y}} = 0. \end{equation*}

As in all thin-film flows, the $y$-component of the momentum equation gives at leading order that the pressure is constant across the flow, i.e. $p=p(x,t)$. The conservation of momentum in the $x$-direction is then

\begin{equation*} 0 =- \frac{{\rm d} p}{{\rm d}\kern0.06em x} + {\frac{\partial ^2u}{\partial y^2}} +\frac{c}{De}\left({\frac{\partial A_{xx}}{\partial x}} + {\frac{\partial A_{xy}}{\partial y}}\right)\!. \end{equation*}

Finally, the Oldroyd-B equations are

\begin{equation*} \begin{aligned} &\left({\frac{\partial }{\partial t}} + u\,{\frac{\partial }{\partial x}} + v\,{\frac{\partial }{\partial y}}\right) A_{xx} -2eA_{xx} - 2\gamma_1A_{xy} =-\frac{1}{De}\,A_{xx}, \notag\\ &\left({\frac{\partial }{\partial t}} + u\,{\frac{\partial }{\partial x}} + v\,{\frac{\partial }{\partial y}}\right) A_{xy} -\gamma_2A_{xx}\ \ - \gamma_1A_{yy} =-\frac{1}{De}\,A_{xy}, \notag\\ &\left({\frac{\partial }{\partial t}} + u\,{\frac{\partial }{\partial x}} + v\,{\frac{\partial }{\partial y}}\right) A_{yy} -2\gamma_2A_{xy} +2eA_{yy} =-\frac{1}{De}\,(A_{yy}-1), \end{aligned} \end{equation*}

where the shear-rates and extension-rate are

\begin{equation*} \gamma_1 = {\frac{\partial u}{\partial y}}, \quad \gamma_2 = {\frac{\partial v}{\partial x}}, \quad e = {\frac{\partial u}{\partial x}} =- {\frac{\partial v}{\partial y}}. \end{equation*}

These governing equations are equivalent to those written down before in the papers cited at the end of § 1. Compared with those earlier works, we have not included negligibly small terms in the small $\epsilon$, and have used a microstructural representation of the Oldroyd-B rheology. Note that the scaled $A_{xx}$ relaxes to $\epsilon ^2$, which is neglected in the first of the Oldroyd-B equations. Note that in lubrication approximations, the hoop-stress is neglected, so there can be no purely elastic instability due to curved streamlines.

2.3. Curvilinear coordinates

For numerical calculations, it is convenient to map the geometry of the contraction onto a unit square

\begin{equation*} 0 \le x \le 1,\quad 0 \le \eta=\frac{y}{h(x)} \le 1. \end{equation*}

The coordinate lines $x=\rm {const.}$ and $\eta =\rm {const.}$ are not orthogonal in the $xy$-plane, due to the slope $h'(x)$. Transforming partial differential equations is much easier with orthogonal curvilinear coordinates. Because the slope of the boundary is $O(\epsilon )$ small in the slowly varying geometry, only a small displacement in the $x$-direction is necessary to create an orthogonal system

\begin{equation*} x(\xi,\eta) = \xi + \epsilon^2 \tfrac14 (h^2(\xi))'(1-\eta^2) + O(\epsilon^4),\quad y(\xi,\eta) = h(\xi)\,\eta.\end{equation*}

In the previous subsection, we used $u$ for the downstream $x$-component of velocity, and $v$ for the cross-stream $y$-component. Henceforth, we change to use $u$ for the velocity in the downstream $\xi$-direction, and $v$ for the flow in the cross-stream $\eta$-direction. There is a negligible $O(\epsilon ^2)$ difference between the flow $u$ in the $\xi$- and $x$-directions. However, the flow in the $y$-direction is greater by $uh'$ than the flow $v$ in the $\eta$-direction, because the small slope $h'$ is multiplied by a large downstream velocity $u$. In a similar way, we shall use $A_{11}$, $A_{12}$ and $A_{22}$ for the components of the microstructure in the $\xi \xi$-, $\xi \eta$- and $\eta \eta$-directions. As with the velocity, there is a negligible $O(\epsilon ^2)$ difference between $A_{11}$ and $A_{xx}$, and $O(1)$ differences between the other curvilinear and Cartesian components.

In these curvilinear coordinates, the governing equations become for mass and momentum

(2.1)$$\begin{gather} \frac1{h}\,{\frac{\partial (hu)}{\partial x}} + \frac1{h}\,{\frac{\partial v}{\partial \eta}} = 0, \end{gather}$$

(2.2)$$\begin{gather}0 =- \frac{{\rm d} p}{{\rm d}\kern0.06em x} + \frac{1}{h^2}\,{\frac{\partial ^2u}{\partial \eta^2}} + \frac{c}{De}\left(\frac{1}{h}\,{\frac{\partial (hA_{11})}{\partial x}} + \frac{1}{h}\,{\frac{\partial A_{12}}{\partial \eta}}\right)\!, \end{gather}$$

and for Oldroyd-B

(2.3a)$$\begin{align}\left({\frac{\partial }{\partial t}} + u\,{\frac{\partial }{\partial x}} + \frac{v}{h}\,{\frac{\partial }{\partial \eta}}\right) A_{11} -2eA_{11} - 2\gamma_1A_{12} &=-\frac{1}{De}\,A_{11}, \end{align}$$

(2.3b)$$\begin{align}\left({\frac{\partial }{\partial t}} + u\,{\frac{\partial }{\partial x}} + \frac{v}{h}\,{\frac{\partial }{\partial \eta}}\right) A_{12} -\gamma_2A_{11} \ \ - \gamma_1A_{22} &=-\frac{1}{De}\,A_{12}, \end{align}$$

(2.3c)$$\begin{align}\left({\frac{\partial }{\partial t}} + u\,{\frac{\partial }{\partial x}} + \frac{v}{h}\,{\frac{\partial }{\partial \eta}}\right) A_{22} -2\gamma_2A_{12} +2eA_{22} &=-\frac{1}{De}\,(A_{22}-1), \end{align}$$

with shear-rates and extension-rate

(2.4)\begin{equation} \gamma_1 = \frac{1}{h}\,{\frac{\partial u}{\partial \eta}},\quad \gamma_2 = h\,{\frac{\partial }{\partial x}}\left(\frac{v}{h} \vphantom{\frac{\partial }{\partial t}}\right)\!,\quad e = {\frac{\partial u}{\partial x}} =- \frac{1}{h}\,{\frac{\partial v}{\partial \eta}}- \frac{uh'}{h}. \end{equation}

There are boundary and inlet conditions. On the velocity, we impose no-slip on the upper rigid boundary,

(2.5)\begin{equation} u = v = 0\quad \text{on}\ \eta = 1, \end{equation}

and symmetry across the centreline,

(2.6)\begin{equation} {\frac{\partial u}{\partial \eta}} = v = 0\quad \text{on}\ \eta = 0. \end{equation}

The constitutive equation is hyperbolic and requires knowledge of the microstructural state at the entry. We assume a Poiseuille flow in the straight entry channel, and the microstructure is in the steady state of the simple shear there:

(2.7)\begin{equation} u = \frac{3}{2h}\,(1-\eta^2),\quad A_{22} = 1,\quad A_{12} =-3\,De\,\eta/h^2,\quad A_{11} = 18\,De^2\,\eta^2/h^4, \end{equation}

with $h=1$ at the entry. The factor $3/2$ in the velocity ensures that the volume flux is $1$ in the half-width of the channel $0\le \eta \le 1$.

The vanishing of the cross-stream velocity $v$ on the centreline and upper boundary means that the downstream volume flux is constant along the contraction, and we have chosen to normalise this flux to be unity:

\begin{equation*} h \int_0^1 u\,{\rm d}\eta = 1,\quad \text{independent of $x$.}\end{equation*}

Maintaining the flux to be constant sets the local value of the pressure gradient. The last expression can be twice integrated by parts. Using the boundary conditions (2.5) and (2.6), we have

\begin{equation*} -h\int_0^1 \frac12\,(1-\eta^2)\,{\frac{\partial ^2u}{\partial \eta^2}}\,{\rm d}\eta = 1.\end{equation*}

Substituting the expression for $\partial ^2u/\partial \eta ^2$ from the momentum equation (2.2) and performing a couple of integrals, we find an expression for the pressure gradient:

(2.8)\begin{equation} \frac{{\rm d} p}{{\rm d}\kern0.06em x} =-\frac{3}{h^3} + \frac{c}{h\,De} \int_0^1\left(\frac32\,(1-\eta^2)\,{\frac{\partial (hA_{11}) }{\partial x}} + 3\eta A_{12}\right){\rm d}\eta. \end{equation}

An important advantage of the lubrication approximation is that the pressure is determined locally with no need to solve a numerically expensive Poisson problem.

2.4. Numerical method

The contraction runs from $x=0$ to $x=1$ in the non-dimensionalisation. A short entry channel and a long exit channel, both of constant width, are added to the contraction. There is no upstream influence in the governing physics (2.1) to (2.3), so theoretically no entry channel is required. However, it is convenient for coding the numerical method and for controlling the effects of a small numerical diffusion to have a short entry. On the other hand, the exit channel must be long to allow the elastic stresses to relax to their steady values. We shall discuss in § 4.1 how long the exit channel must be, but typical values range from $x_\infty =5$ to $x_\infty =20$. A single simple shape is used for the contraction. To allow second-order numerical accuracy, a smooth shape is used with vanishing slope at the start and end:

(2.9)\begin{equation} h(x) = \begin{cases} 1 & \text{for }{-0.5} < x < 0, \\ 1- (1-H^{-1})x^2(2 - x)^2 & \text{for }0< x < 1,\\ H^{-1} & \text{for }1 < x < x_\infty. \end{cases} \end{equation}

Here, $H$ is the contraction ratio, the ratio of the entry to exit widths. This shape is shown in figure 3 for $H=2$.

The numerical approach is at each instant to solve for the flow $(u,v)$ given the instantaneous values of the elastic stress $A_{11}$, $A_{12}$ and $A_{22}$, and then to use this flow to time step the elastic stresses until they reach a steady state.

Given the elastic stresses and the pressure gradient from (2.8), the momentum equation (2.2) is integrated with respect to $\eta$, starting at the centreline $\eta =0$ where $\partial u/\partial \eta = 0$, to find the shear-rate $\partial u/\partial \eta$. This shear-rate is then integrated with respect to $\eta$, starting at the upper boundary $\eta =1$ where $u=0$, to find the velocity $u$. Note that these integrations at each downstream section are independent of one another. Finally, with $u(x,\eta )$ known, the cross-stream velocity $v$ is found by integrating the mass conservation (2.1) with respect to $\eta$, starting at the centreline where $v=0$.

Given the flow $u(x,\eta )$ and $v(x,\eta )$, the Oldroyd equations (2.3) are time-stepped to find the microstructure $A_{11}$, $A_{12}$ and $A_{22}$ at the next time step. For initial conditions of the microstructure, we take the expressions in (2.7) corresponding to a steady uniform Poiseuille flow in a channel of width $h(x)$.

Equispaced finite differences are used. Second-order spatial central differencing is used, except for the $u$ advection, which is first-order upwinded because downstream advection is important at high Deborah numbers. Upwinding is not used for the small $v$ advection. When calculating some fluxes, a second-order correction is applied to produce a fourth-order result for the flux. The time stepping is made by first-order forward differencing. Typical values of the grid size and time step are $\delta x=0.0208$, $\delta \eta = 0.0333$ and $\delta t = 10^{-3}$. Random spot checks are made to confirm an accuracy of around three figures. Run times are around 15 seconds on a laptop, considerably faster than solving the full non-lubrication equations by finite elements. At high Deborah numbers, a numerical instability occurred, which could be delayed by decreasing $\delta \eta$ while keeping $\delta x$ fixed, as suggested by Keiller (Reference Keiller1992b). Also, for stability, a smaller $\delta t$ was needed at very small Deborah numbers $De$. The time to come to a steady state was typically $t=5$, shorter for small Deborah numbers, and longer for high Deborah numbers with long exit channels.

4.1. Stress relaxation

Figure 4 in § 3.3 shows that most of the relaxation of the stress occurs in the exit channel, and takes a long distance. We study stress relaxation in this subsection. Debbaut et al. (Reference Debbaut, Marchal and Crochet1988) noted in some finite-element calculations that a long downstream section of some 250 radii in length was required to calculate the Couette correction reliably at $De=30$. Keiller (Reference Keiller1993b) examined the spatial decay of steady perturbations of plane Poiseuille flow for the Oldroyd-B equation, finding continuous spectra of eigenvalues. Far downstream, the stress relaxed exponentially on a length scale of the maximum velocity multiplied by the relaxation time.

In figure 9 we plot the decay of the difference between the pressure gradient $\textrm {d} p/{\textrm {d}\kern 0.06em x}$ and its eventual steady value $-3(1+c)H^3$ as a function of distance downstream for various contraction ratios and Deborah numbers. The log-linear plot shows an exponential approach to the steady value. Scaling the distance with the local Deborah number $H\,De$ brings the curves parallel, corresponding to a decay on a length scale $\frac 32 H\,De$. This is the result above of Keiller (Reference Keiller1993b). To obtain better than two-figure accuracy in the Couette correction, one must therefore extend the exit channel to $7H\,De$. This estimate of the necessary length of the exit channel of $210$ for $H\,De=30$ is consistent with the 250 of Debbaut et al. (Reference Debbaut, Marchal and Crochet1988).

Figure 9. Relaxation of the pressure gradient in the exit channel. The difference between the pressure gradient $\textrm {d} p/{\textrm {d}\kern 0.06em x}$ and its final steady value $- 3(1+c)H^3$ as a function of the distance downstream $x$ divided by the local Deborah number, $H\,De$, for $De=0.3$ with $H= 2^{1/2}, 2, 2^{3/2}$, and for $De=0.5$ with $H=2$. The black line is $1000\exp (-2x/(3H\,De))$.

Knowing that the pressure gradient approaches its steady value exponentially with a known exponential rate, one can shorten the numerical length of the exit channel and add a simple extrapolation to the Couette correction. This reduction of the computation was not employed because the calculations took less than a minute on a laptop.

4.2. Pressure drop

To calculate the pressure drop in the exit channel, we return to the expression (2.8) for the pressure gradient. This expression can be integrated from the start of the exit channel at $x=1$ to some distance down the channel, to $x=L$. Because the width of the channel, $h(x)= 1/H$, is constant in the exit section, it can be taken outside the partial derivative, outside the integral, and then cancelled, to yield

\begin{equation*} \Delta p = \frac{3(L-1)}{h^3} - \frac{c}{De} \int_0^1 \frac32\,(1-\eta^2)\,[A_{11}]_1^L\,{\rm d}\eta - \frac{c}{De\,h}\int_0^1 3\eta \left(\int_1^L A_{12}\,{{\rm d}\kern0.06em x}\right){\rm d}\eta.\end{equation*}

Thus we need only the starting and finishing values of $A_{11}$, and the integral down the channel of $A_{12}$. Note that we do not set $h=1/H$ in the first half of this subsection, because we shall use the same analysis for the exit of the constriction where $h=1$.

While the shear stress term $A_{12}$ is $O(De)$ smaller than the normal stress $A_{11}$, its relaxation takes place over a long $O(De)$ distance, so its contribution to the pressure drop is of a similar magnitude. We need the Oldroyd-B equations for how $A_{12}$ relaxes down the channel.

At this stage we make an approximation that the velocity has the same constant parabolic profile all the way down the exit channel; see figure 3. This approximation is asymptotically correct in the low concentration limit $c\ll 1$.

With $u$ independent of $x$, the cross-flow component of velocity vanishes, $v=0$. The Oldroyd-B equations needed, (2.3c) and (2.3b), are then

\begin{align*} u\,{\frac{\partial A_{22}}{\partial x}} &= -\ \kern0.12em\frac{1}{De}\,(A_{22}-1),\\ u\,{\frac{\partial A_{12}}{\partial x}} &= \frac{1}{h}\,{\frac{\partial u}{\partial \eta}}\,A_{22} \,-\, \frac{1}{De}\,A_{12}. \end{align*}

While these are simple to solve, we do not need the evolution of the relaxation, we just need an integral. Integrating down the flow, using $u$, and so $\partial u/\partial \eta$, are constant along the flow,

\begin{align*} u\,[A_{22}]_1^L &= -\ \kern0.12em\frac{1}{De}\int_1^L(A_{22}-1)\,{{\rm d}\kern0.06em x}, \\ u\,[A_{12}]_1^L &=\frac{1}{h}\,{\frac{\partial u}{\partial \eta}}\left(L-1 + \int_1^L(A_{22}-1)\,{{\rm d}\kern0.06em x}\right) \,-\,\frac{1}{De}\int_1^LA_{12}\,{{\rm d}\kern0.06em x}. \end{align*}

Hence the required integral of $A_{12}$ is

\begin{equation*} \int_1^L A_{12}\,{{\rm d}\kern0.06em x} =-De\,u\,[A_{12}]_1^L + \frac{1}{h}\,{\frac{\partial u}{\partial \eta}}\left(De\,(L-1) -De^2\,u\,[A_{22}]_1^L\right)\!. \end{equation*}

Substituting the expressions for the flow $u$ and the shear-rate $\partial u/\partial \eta$, we have an expression for the pressure drop in the exit channel:

(4.1) \begin{align} \Delta p &= \frac{3(1+c)(L-1)}{h^3} - \frac{c}{De}\int_0^1 \frac32\,(1-\eta^2)\,[A_{11}]_1^L\,{\rm d}\eta \nonumber\\ &\quad + \frac{c}{h^2}\int_0^1\frac92\,\eta(1-\eta^2)\,[A_{12}]_1^L\,{\rm d}\eta - \frac{c\,De}{h^4}\int_0^1\frac{27}{2}\,\eta^2(1-\eta^2)\,[A_{22}]_1^L\,{\rm d}\eta. \end{align}

This result is valid for all $De$. It assumes that the velocity profile takes a constant parabolic form.

At the far end of the exit channel, we have the steady Poiseuille values for all $De$:

\begin{equation*} A_{11} = \frac{18\,De^2}{h^4(L)}\,\eta^2, \quad A_{12} =-\frac{3\,De}{h^2(L)}\,\eta, \quad A_{22} = 1. \end{equation*}

In the low-$De$ limit, § 3.2 gives

\begin{equation*} [A_{22}]_1^L = O(De^2),\quad [A_{12}]_1^L = O(De^3)\quad\text{and}\quad [A_{11}]_1^L = O(De^4), \end{equation*}

using our smooth geometry $h(x)$ with zero slope at the end of the contraction, $h'(1)=0$. Hence the pressure drop is little changed from the steady Poiseuille value:

\begin{equation*} \Delta p = 3(1+c)(L-1)H^3 + O(De^3). \end{equation*}

In the high-$De$ limit, we have (3.4a,b) at the start of the exit channel,

\begin{equation*} A_{11} = \frac{18\,De^2}{h^2(1)}\,\eta^2, \quad A_{12} = -3\,De\,\eta, \quad A_{22} = h^2(1), \end{equation*}

so

(4.2)\begin{equation} \Delta p = {3(1+c)(L-1)}{H^3} - \tfrac{36}{5}c\,De\,(H^4-H^2). \end{equation}

There are equal contributions of $18/5$ from the normal and elastic shear stresses. The normal stresses are stronger far downstream compared with the entry to the exit channel. These normal stresses pull the fluid along, so requiring less pressure to drive the flow. The shear stress is below the steady Poiseuille value, so exerts less resistance until relaxed, again requiring less pressure drop. The shear stress contribution of $18/5$ is made up of two equal parts, arising from $A_{12}$ and $A_{22}$. In the contraction at high $De$, $A_{12}$ remains constant while $A_{22}$ is squashed by the same factor of $H^2$ that $A_{11}$ is stretched.

At high $De$, the net rate of working against the normal stresses at the beginning and end of the exit channel is $\frac {18}{5}c\,De\,(H^4-H^2)$, and the net out-flux of elastic energy is half this. Stress relaxation is, however, more complicated in the exit channel, and contributes more to the pressure drop.

The expressions (3.4a,b) used above are for the elastic stresses at high $De$ in the fast central part of the flow. Near to the boundary, there is a layer of thickness $O(1/De)$ with larger normal and shear stresses. Examining the integrals in (4.1), we see that the boundary layer will contribute only small corrections $O(cH^4/De)$.

Combining the pressure drop across the contraction (3.5) with the extra pressure drop from the stresses relaxing in the exit channel (4.2) at high $De$, we predict the total extra pressure drop (non-dimensionalised by $\mu _0Q\ell /h^3_0$)

(4.3)\begin{equation} \Delta p = \Delta p_1 + c\,\Delta p_2 - \frac{9 c\,De}{5}\,(H^2- 1)(4H^2+1). \end{equation}

To test this prediction, we plot in figure 10 the numerically calculated pressure drop across the contraction and exit channel minus the first two terms in (4.3), and then rescaled by $(H^2-1)(4H^2+1)$.

Figure 10. The extra pressure drop through a contraction including the extra from the exit channel minus the first two terms in (4.3), i.e. $\Delta p - \Delta p_1 - c\,\Delta p_2$, scaled by $(H^2 - 1)(4H^2+1)$, for $c=1$ and $H = 2^{1/2}$, $2$, $2^{3/2}$ and $4$. The black line is the high-$De$ prediction $-\frac {9}{5}c\,De$.

The curves in figure 10 for each of the four contraction ratios seem to be straight lines with the predicted slope $-\frac {9}{5}c$ and with offsets varying between 0 and 0.2. Other than a possible small reduction in the slope in $De<0.1$, the high-$De$ prediction (4.3) seems to be applicable for all $De$. We note that at low $De$, the exit channel contributes only an $O(De^3)$ correction, while the pressure drop through the contraction (3.3) varies as $-4.5 c\,De\,(H^4-1)$. This low-$De$ variation differs little from the $-7.2c\,De\,(H^4-H^2)$ from the exit channel at high $De$. At high $De$, the exit channel contributes a greater reduction in the pressure drop than the contraction itself, by a factor $4H^2$.

With the viscoelastic change in the pressure drop (4.3) decreasing linearly with the Deborah number, there must be a worry of the total pressure drop becoming negative above a critical Deborah number. That would yield a mechanical instability. Fortunately, the reduction $7.2cH^4\,De$ attains 90 % of this value only at a distance $3.5H\,De$ down the exit channel, over which distance the Newtonian pressure drop is $10.5(1+c)H^4\,De$. Hence the total pressure drop cannot become negative.

6.1. Low $De$

In (3.3), we gave an expression for the $O(De)$ reduction in the pressure drop at low $De$. It involves the difference between the width of the channel at the beginning and end of the flow. This difference vanishes for a constriction, so there should be no change in the pressure drop at $O(De)$. Boyko & Stone (Reference Boyko and Stone2022) gave two further terms in an expansion for small $De$ in their (4.11) and (4.14). The structure of these further terms suggests that the reduction in the pressure drop across a constriction should be an even function of $De$, should be proportional to $(H^4-1)$, where $H$ is the constriction ratio, and should depend on the effective Deborah number for the narrowest width, $De\,H$. In figure 13, we plot the reduction in the pressure drop with this scaling as a function of $De\,H$. We see that the results for different constriction ratios follow the same quadratic variation in $De\,H<0.2$. The pressure drops in figure 13 are for just the constriction, not including the contribution from the exit channel, and they continue to decrease at the higher $De$ plotted. This should be contrasted with the levelling off of the pressure drops when the contribution from the exit channel is included in figure 12. This difference is clearest in the case $H=2^{1/2}$.

Figure 13. The low Deborah number behaviour through a constriction. The scaled reduction in the pressure drop through just the constriction, $(\Delta p - 2(1+c)\Delta p_1)/(H^4-1)$, as a function of the effective Deborah number at the narrowest point $De\,H$, for $c=1$ and $H = 2^{1/2}$, $2$, $2^{3/2}$ and $4$.

6.2. Problems at high $De$

Our theory for high $De$ worked well in the contraction and the expansion, giving good predictions for $De>0.3$. The constriction is not so kind, possibly requiring $De>20$ for good predictions.

The high-$De$ theory above was based on the observation in figure 6 that there was little relaxation in the contraction. Without relaxation, the microstructure $A_{11}$ varies along the accelerating streamlines in proportion to the square of the velocity, while $A_{22}$ varies with the inverse of the square of the velocity, and $A_{12}$ remains constant along the streamline. Applying those variations to a constriction, one would predict that

1. (i) the normal stresses are symmetric about the minimum gap and so contribute nothing to the pressure drop;

2. (ii) the elastic shear stresses contribute as in the contraction and in the expansion;

3. (iii) from (i) and (ii), the total pressure drop should be $2(\Delta p_1 + c\,\Delta p_2)$;

4. (iv) at the end of the constriction, all the stresses return to their values at the start, so the exit channel contributes nothing.

However, we have already seen from comparing figures 12 and 13 that the exit channel does contribute significantly.

In figure 14, we revisit the streamwise variations of the microstructure in figure 6 now plotted through the constriction from $x=0$ to $x=2$. Only the streamlines in the central core $0\le \eta \le 0.6$ are included. In the first contraction part from $x=0$ to $x=1$, the normal stresses increase by a factor of just over $H^2=4$, corresponding approximately to the square of the increase in the velocity, the elastic shear stresses start approximately constant but are beginning to increase by $x=1$, and $A_{22}$ drops to only 0.5 at $x=1$ rather than the inverse of the square of the velocity $0.25$. In the second expansion part from $x=1$ to $x=2$, there are clear signs of relaxation. The normal stresses do not return to their entry values but end at least 50 % higher, while the elastic shear stresses and $A_{22}$ end at twice their entry values.

Figure 14. Streamwise variations on constant $\eta = 0.1\ {\rm in\ steps\ of}\ 0.1\ {\rm to}\ 0.6$ of the elastic normal stress $NS$, elastic contribution to the shear stress $SS$, and $A_{22}$, all scaled by their entry values, for $H= 2$, $De=0.5$, $c=1.0$.

It should be noted that the Deborah number was defined to be the ratio of the residence time from $x=0$ to $x=1$ to the relaxation time. Now that the constriction extends to $x=2$, the effective Deborah number has been halved, from $De=0.5$ to $De=0.25$ in figure 14, so that it should not be unexpected that relaxation can no longer be neglected.

The higher normal stresses at the end of the constriction compared with the entry reduce the pressure drop through the constriction. However, those same higher normal stresses at the end of the constriction compared with far downstream in the exit channel lead to a compensating increase in the pressure drop through the exit channel. On the other hand, the higher (i.e. higher in magnitude, as they are negative) elastic shear stresses at the end of the constriction increase the pressure drop in both the constriction and the exit channel. Finally, (4.1) shows that $A_{22}$ starting at the exit channel with values exceeding those far downstream also leads to an increase in the pressure drop. For $De=0.5$ and independent of the constriction ratio $H= 2, 2^{3/2}, 4$, approximately 60 % of the extra pressure drop in the exit channel comes from the normal stresses, 25 % from the elastic shear stresses, and 15 % from $A_{22}$.

There are two problems to constructing a theory to predict the behaviour at high $De$ in a constriction. First, the leading $O(De)$ effect in the contraction and expansion came from the normal stresses, and this effect vanishes in the constriction. Therefore, we are considering corrections to the leading-order theory, and these are more subtle. Within the constriction, our high-$De$ theory gives a first correction from the elastic shear stress not having time to build up to the higher values, predicting a pressure drop $2(\Delta p_1 + c\,\Delta p_2)$ at high $De$ compared with $2(1+c)\,\Delta p_1$ at low $De$. For $De=0.5$ and $H=2$, this predicts a 30 % reduction in the pressure drop, whereas our numerical calculations give only a 19 % reduction. Clearly, second and further corrections must be contributing at the not-so-large $De=0.5$. While the leading-order $O(De)$ term from normal stresses dominated, these correction terms lead to small offsets, as observed in figures 10 and 11.

The second problem comes from $A_{22}$. So far, it has made only an incidental appearance in the calculation of the relaxation of stresses in the exit channel, and we have not yet discussed its role, which was minor in the contraction and expansion, but becomes major in the constriction. The microstructure $A_{22}$ corresponds to two fibres perpendicular to the flow. In an accelerating flow without relaxation, each fibre is compressed as the streamlines come together, so $A_{22}$ varies with the inverse of the square of the velocity. But in a constriction, relaxation cannot be ignored. A shearing of the microstructure $A_{22}$ produces the shear component $A_{12}$. Thus the larger than expected values of $A_{22}$ at the end of the constriction produce larger shear stresses, and hence a larger pressure drop in the exit channel.

6.3. An asymptotic theory for high Deborah numbers

The high-$De$ theory presented so far was based on physical insights. Such an approach is sufficient for the leading-order terms. For the correction terms, we need a more organised asymptotic approach.

Progress can be made in the low-$c$ limit, where the velocity profile remains parabolic with a local magnitude proportional to $1/h(x)$. In our curvilinear coordinates, the cross-stream component of the velocity then vanishes. This removes several terms in the Oldroyd-B equations. The flow, shear-rate and extension-rate become

\begin{equation*}u(x,\eta) = \frac{F(\eta)}{h(x)},\quad \gamma(x,\eta) = \frac{f(\eta)}{h^2(x)}\quad \text{and}\quad e(x,\eta) =-\frac{h'\,F(\eta)}{h^2(x)}, \end{equation*}

with

\begin{equation*}F(\eta) = \tfrac32(1-\eta^2) \quad\text{and}\quad f(\eta) =-3\eta. \end{equation*}

Using the known leading-order behaviour, we make a multiplicative substitution

\begin{equation*}A_{11} = 2\,De^2\,f^2h^{-2}b_{11},\quad A_{12} = De\,f\kern0.06em b_{12}, \quad A_{22} = h^2 b_{22}. \end{equation*}

The Oldroyd-B equations (2.3) then become

(6.1) \begin{align} F\,{\frac{\partial b_{11}}{\partial x}} &=-\frac{h}{De}\,(b_{11}-b_{12}),\quad\ \,\text{with}\ b_{11}(0) = 1, \nonumber\\ F\,{\frac{\partial b_{12}}{\partial x}} &=-\frac{h}{De}\,(b_{12}-b_{22}),\quad\ \,\text{with}\ b_{12}(0) = 1, \nonumber\\ F\,{\frac{\partial b_{22}}{\partial x}} &=-\frac{h}{De}\left(b_{22} -\frac{1}{h^2}\right)\!,\quad\text{with}\ b_{22}(0) = 1, \end{align}

similar to the equations given by Renardy (Reference Renardy1994). These are valid for arbitrary $De$. For $De\gg 1$ and in the central core of the flow where $F(\eta )\gg 1/De$, we solve iteratively for

\begin{align*}b_{22} &= 1 + \frac{1}{De\,F}\,I_1 - \frac{1}{(De\,F)^2}\,I_2 + \frac{1}{(De\,F)^{3}}\,I_3 +O((De\,F)^{-4}),\\ b_{12} &= 1 + \frac{1}{(De\,F)^2}\,I_2 - \frac{1}{(De\,F)^{3}}\,I_3 + O((De\,F)^{-4}),\\ b_{11} &=1 + \frac{1}{(Dev\,F)^{3}}\,I_3 + O((De\,F)^{-4}), \end{align*}

where

\begin{gather} I_1 = \int_0^x\frac{1-h^2(x')}{h(x')}\,{{\rm d}\kern0.06em x}',\quad I_2 = \int_0^x h(x')\int_0^{x'} \frac{1-h^2(x'')}{h(x'')}\,{{\rm d}\kern0.06em x}''\,{{\rm d}\kern0.06em x}' \notag\\ \text{and}\quad I_3 = \int_0^x h(x')\int_0^{x'} h(x'')\int_0^{x''} \frac{1-h^2(x''')}{h(x''')}\,{{\rm d}\kern0.06em x}'''\,{{\rm d}\kern0.06em x}''\,{{\rm d}\kern0.06em x}'. \notag \end{gather}

Using these expressions in (2.8) for the pressure gradient in the constriction from $x=0$ to $x=2$, we find that the pressure drop through the constriction itself is

\begin{equation*} \Delta p = 2(\Delta p_1 + c\,\Delta p_2) + O(De^{-1}). \end{equation*}

The correction comes from the boundary layer next to the wall, $\eta = 1 - O(De^{-1})$. Using the expressions evaluated at $x=2$ in (4.1) for the extra pressure drop from relaxation in the exit channel, we find

\begin{equation*} \Delta p = c\,\Delta p_3 + O(De^{-1}), \end{equation*}

with again the correction coming the boundary layer next to the wall. Here

(6.2)\begin{equation} \Delta p_3 = 3\int_0^2 \frac{1-h^2(x)}{h(x)}\,{{\rm d}\kern0.06em x}, \end{equation}

whose values are given in table 1. Note that while $A_{22}$ differs from its equilibrium value far downstream by only $O(1/De)$, the long $O(De)$ relaxation distance produces an $O(1)$ effect in the pressure drop. Hence the pressure drop across the constriction, plus the extra contribution from relaxation in the exit channel, is predicted to be

(6.3)\begin{equation} 2\,\Delta p_1 + 2c\,\Delta p_2 + c\,\Delta p_3 . \end{equation}

In figure 15, we plot, as a function of the Deborah number, the pressure drop including the extra contribution from the exit channel minus the above predicted result, all scaled by $H^{3/2}$. The scaling is arbitrary, chosen to bring the curves together. If our predictions were correct, the curves should plateau to zero at high $De$. They seem to plateau, but not to zero. For a constriction ratio $H=2$, the pressure drop reduces by 13 % from $De=0$ to $De=0.5$, see figure 12, whereas we predict a 20 % reduction; and for the larger constriction ratio $H=4$, the pressure drop reduces by 28 %, whereas we predict 39 %. We cannot explain this discrepancy. It may be that $De=0.5$ is not sufficiently high. It may be that $c=1$ is not small, and the velocity profile changes significantly.

Figure 15. The high Deborah number behaviour through a constriction. The scaled pressure drop, including the extra drop from the exit channel, $(\Delta p - 2(\Delta p_1 + c\,\Delta p_2) -c\,\Delta p_3)/H^{3/2}$ as a function of $De$, for $c=1$ and $H = 2^{1/2}$, $2$, $2^{3/2}$ and $4$.

As an initial investigation into the discrepancy, we have examined the contribution of the relaxation of $A_{22}$ to the pressure drop in the exit channel. We calculate $A_{22}(x,\eta )$ at $x=2$ by integrating numerically the low-$c$ Oldroyd-B equation (6.1) along the flow from $x=0$ to $x=2$ at fixed $\eta$. Then integrating across the flow, we evaluate the last integral in (4.1) for the pressure drop in the exit channel. The results are plotted in figure 16 as a function of $1/De$ for four different constriction ratios. We see that as $1/De\to 0$, all curves tend to the predicted values of $\Delta p_3$ given in table 1. However, a large value $De=20$ is required to come within 10 % of the high-$De$ limit. At $De=0.5$, the values are typically only 15 % of the limit values. Our numerical method for solving the lubrication equations struggles beyond $De=1$, so we are currently unable to access the high $De$ needed to study the constrictions further.

Figure 16. The contribution of $A_{22}$ to the extra pressure drop in the exit channel following a constriction as a function of $1/De$, for constriction ratios $H=2^{1/2}$, $2$, $2^{3/2}$ and $4$. The horizontal lines are the limiting values $\Delta p_3$ given in table 1.