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$.


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 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).
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), β = μ p /μ 0 1, where μ p is the polymer contribution to the total zero-shear-rate viscosity μ 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, μ 0 , is only slightly larger than the solvent viscosity, μ s (Remmelgas et al. 1999;Mokhtari et al. 2022).Furthermore, the limit β = μ p /μ 0 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.

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 , where h(z) , as illustrated in figure 1. Upstream of the contraction inlet (z = 0) there is an entry channel of height 2h 0 and length 0 , and downstream of the contraction outlet (z = ) there is an exit channel of height 2h and length .The fluid flow has velocity 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 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 , y), where the z and z = z − 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 (ξ, η) defined in § 2.3.
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 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 σ is 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 μ s , where E = (∇u + (∇u) T )/2 is the rate-of-strain tensor, and the last term, τ p , is the polymer contribution.We note that 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 τ p can be expressed in terms of the (symmetric) conformation tensor (or the deformation of the microstructure) A as (Bird et al. 1987;Larson 1988;Morozov & Spagnolie 2015) where G is the elastic modulus, λ is the relaxation time and μ p = Gλ is the polymer contribution to the shear viscosity at zero shear rate.It is also convenient to introduce the total zero-shear-rate viscosity μ 0 = μ s + μ p .The evolution equation for the deformation of the microstructure A of the Oldroyd-B model fluid is given at steady state as (Bird et al. 1987;Larson 1988;Morozov & Spagnolie 2015) (2.4)

Scaling analysis and non-dimensionalization We consider narrow configurations, in which h(z)
, 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) where we have introduced the aspect ratio of the configuration, which is assumed to be small the contraction ratio the viscosity ratios 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, λ, 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, λ, to the residence time in the contraction region, /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 = Wi, and for narrow geometries with 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 μ 0 = μ s + μ 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 Ãzz = O(De 2 ), Ãzy = O(De) and Ã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.

Dimensionless lubrication equations in Cartesian coordinates
Using the non-dimensionalization (2.5)-(2.9a,b), to the leading order in , the governing equations (2.1)-(2.4)take the form From (2.10c), it follows that P = P(Z), i.e. the pressure is independent of Y up to O( 2 ), consistent with the classical lubrication approximation.We note that the scaled Ãzz on the right-hand side of (2.10d) relaxes to 2 , which is neglected at the leading order in .
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 (ξ, η), as illustrated in As shown in Appendix A, the curvilinear coordinates (ξ, η) are orthogonal with a relative error of O( 4), i.e. ∇ξ • ∇η = O( 4 ).
Hereafter, we use u = ue ξ + ve η and A = A 11 e ξ e ξ + A 12 (e ξ e η + e η e ξ ) + A 22 e η e η to denote, respectively, the components of velocity and deformation of the microstructure in curvilinear coordinates (ξ, η).The corresponding non-dimensional velocity components in different coordinates are related through (see Appendix A) (2.12a,b) Similarly, the scaled conformation tensor components in different coordinates are related through (see Appendix A) (2.13a) (2.13b) (2.13c) Finally, we note that, since there is only a O( 2 ) difference between the ξ -and z-directions, for convenience, we continue to use Z rather than ξ in curvilinear coordinates.

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.14e) The corresponding boundary conditions on the velocity are 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 with H ≡ 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 ≡ H .

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) Substituting the expression for ∂ 2 U/∂η 2 from (2.14b) into (2.17),we obtain (2.18) which can be rearranged to yield the pressure gradient Integrating (2.19) with respect to Z from 0 to 1 provides the pressure drop P = P(0) − P(1) across the non-uniform region (2.20) Using integration by parts, (2.20) can be expressed as 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).

Low-β 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, β and De, in the lubrication limit ( 1).In this section, we derive analytical expressions for the velocity, conformation tensor and the q − p relation for the pressure-driven flow of a very dilute viscoelastic Oldroyd-B fluid, β = μ p /μ 0 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, β 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) 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( β) 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( β0 ) and the pressure drop up to O( β).
We note that the viscosity ratio β = μ p /μ 0 is related to the so-called concentration of the polymers c = μ p /μ s through β = c/(c + 1).Thus, at the leading order, the limits β 1 and c 1 are identical.

Velocity, conformation and pressure drop at the leading order in β
Substituting (3.1) into (2.14a)-(2.14b)and considering the leading order in β, the continuity and momentum equations take the form subject to the boundary conditions (3.3a-d)The solutions for the axial velocity U 0 and the pressure drop P 0 at the leading order are well known (see, e.g.Boyko & Stone 2022) Substituting (3.4a) into the continuity equation (3.2a) and using (3.3b), yields From (3.5), it follows that, in orthogonal curvilinear coordinates, the velocity in the η-direction is identically zero at O( β0 ), in contrast to the Cartesian coordinates where 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 β, the equations for the conformation tensor components, (2.14c)-(2.14e),simplify to subject to the boundary conditions Equations (3.6) represent a set of one-way coupled first-order semi-linear partial differential equations that can be solved first for Ã22,0 , followed by Ã12,0 and then for Ã11,0 .
Solving (3.6) together with (3.7), we obtain closed-form expressions for Ã22,0 , Ã12,0 and Ã11,0 for arbitrary values of De and the shape function H(Z) Ã22,0 where f (DeU 0 (Z, η)) is defined as It is worth noting that the right-hand sides of (3.8)-(3.10)depend on the product DeU 0 (Z, η) and are not functions of De and η separately.Furthermore, (3.8)-(3.10)clearly show that, while the distribution of Ã22,0 is set solely by the value at the beginning of the non-uniform region, the distribution of elastic shear and normal stresses, Ã12,0 and Ã11,0 , are coupled to the transverse normal stress Ã22,0 .In fact, the elastic normal stress Ã11,0 depends both on Ã12,0 and Ã22,0 .
From (3.8)-(3.10),one might think that the conformation tensor components diverge at the wall (η = ±1).However, using (3.6) and noting that U 0 = ∂U 0 /∂Z = 0 at η = ±1, it follows that, at the walls of the non-uniform channel, 18De 2 H(Z) 4 for all De. (3.12) In § § 3.1.1and 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.

Conformation tensor in the low-De limit
For De 1, we solve the equations iteratively for the conformation tensor components (3.6) to obtain (3.13c) 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 Ãzz , Ãzy and Ãyy up to O(De 2 ) in Cartesian coordinates.For example, using (2.13c) and (3.13), Ãyy can be expressed as

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 Ã22,0 and consider the core flow region.
For De 1, except close to the wall, (3.6a) reduces to Next, since Ã12,0 scales as O(De) while Ã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 so that elastic shear stresses preserve their value from the entry channel through the non-uniform region Ã12,0 (Z, η) = Ã12,0 (0, η) = −3Deη. (3.17) Finally, to determine Ã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 1, we expect the first and second terms to balance each other while the remaining terms are negligible, so that In fact, for De 1, there is a purely passive response of the microstructure, similar to a material line element, transported and deformed by the flow without relaxing.

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 η ∈ [0, 1], and introduce the rescaled inner-region coordinate 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ζ H(Z) Solving (3.22) together with (3.23), we obtain closed-form expressions for A 22 , A 12 and A 11 in the boundary-layer region where F(Z, ζ ) is defined as We note that solutions (3.24) satisfy the matching conditions between the inner and outer regions.Specifically,

Pressure drop at the first order in β
Equation (2.20) shows that the pressure drop depends on the elastic normal and shear stresses Ã11 and Ã12 , and thus, generally, requires the solution of the nonlinear viscoelastic problem.However, in the ultra-dilute limit, corresponding to β = μ p /μ 0 1, we can determine the pressure drop at O( β) 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( β) the pressure drop P 1 , or alternatively (3.27) Thus, for a given flow rate q, the dimensionless pressure drop P = p/(μ 0 q /2h 3 0 ), as a function of the shape function H(Z), the Deborah number De and the viscosity ratio β 1, up to O( β), is given by where the expressions for P 0 and 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 1 and β ∈ [0, 1] (Boyko & Stone 2022), the current result (3.28) applies to the limit of β 1, while allowing De = O(1).

Pressure drop at O( β) in the low-De limit
To calculate the pressure drop 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( β) is where [ Ã11,0 ] Z=0 Z=1 = Ã11,0 (0, η) − Ã11,0 (1, η).The elastic shear stress (SS) contribution to the pressure drop at O( β) is 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 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.

Pressure drop at O( β) in the high-De limit
To calculate the pressure drop 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( β) are for De 1.
(3.33a,b) Substituting (3.33) into (3.27)yields the pressure drop at O( β) in the high-De limit so that the total pressure drop across the non-uniform channel in the high-De limit is for De 1. (3.35) 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 Ã12 = −3Deη/H 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 β 1 0 H(Z) −1 dZ, is smaller than the steady Poiseuille value of 3 β 1 0 H(Z) −3 dZ, further reducing the pressure drop.Finally, we note that the result (3.35) also holds for the expansion geometry H > 1, in which the two physical mechanisms mentioned above lead to an increase in the pressure drop.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 where L = / is the dimensionless length, H = H(Z = 1) = h /h 0 is the dimensionless height of the exit channel, Z = Z − 1, Ã11,0 and Ã12,0 are given in (B4) and (B5) and [ Ã11,0 ] Z =0 Z =L = Ã11,0 (Z = 0, η) − Ã11,0 (Z = L, η).It should be noted that we can express the first-order contribution P ,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) 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 ≡ H .Under this assumption, (4.2) clearly shows that the first-order contribution P ,1 is independent of L since the steady-state values of Ã11,0 , Ã12,0 and Ã22,0 depend solely on the η coordinate.Note, however, that the total pressure in the exit channel depends on L via P = 3L/H 3 + β P ,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 for De 1, (4.3) 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( βDe −1 ) correction to the pressure drop in the exit channel for De 1, as noted by Hinch et al. (2024).

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 where H = H(1)/H(0) = h /h 0 is the ratio of the exit to entry heights; for the contracting geometry we have H < 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 = 0.5 and β = 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 1 and β 1.To obtain the semi-analytical solutions for given values of De and H , we first used MATLAB's routine 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 Z = 10 −4 and η = 0.005.We then used MATLAB's routine trapz to calculate the pressure drop, (3.28) and (4.1), for a contraction and exit channel, respectively.

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 η = 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 = 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.
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 Z /[3De(1 − η 2 )]) (see (B3)-(B5)), the elastic stresses relax over slightly different length scales, with the shortest relaxation distance required for Ã22,0 and the longest for Ã11,0 .The latter behaviour is High-De asymptote η = 0.5 η = 0.5 associated with the nature of the coupling between the elastic stresses so that Ã11,0 depends both on Ã12,0 on Ã22,0 , while Ã12,0 depends only on Ã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 2 = 4, the transverse normal stress is squashed by a factor of H 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 η = 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 η = y/H 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, Ã12,0 (Z = 1, η)/(−3Deη/H 2 ) and Ã11,0 (Z = 1, η)/(18De 2 η 2 /H 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 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.
To provide insight into this viscoelastic boundary layer, we replot in figure 4  Furthermore, examining (3.8)-(3.10),we infer that their right-hand sides are not a function of De and η separately but depend on the product DeU 0 (Z, η).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) Ã12,0 (Z = 1, η)/(−3Deη/H 2 ) and (b) Ã11,0 (Z = 1, η)/(18De 2 η 2 /H 4 ) minus H 2 , divided by the factor 1 − H 2 , as a function of DeU 0 (Z = 1, η) for De = 0.5, 1 and H = 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.

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 ) 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 into (2.19)yields the pressure gradient in the exit channel Noting that in the exit channel U 0 = (3/2H )(1 − η 2 ) and dU 0 /dη = −(3/H )η, and using the expression for U 0 ∂ Ã11,0 /∂Z from (B2c), (5.2) can be written as where the right-hand side is independent of β. 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 (dP/dZ + 3/H 3 )/ β as a function of the downstream distance Z for De = 0.02, 0.2, 1 and 2. (b) Scaled pressure gradient (dP/dZ + 3/H 3 )/ β as a function of the rescaled downstream distance 2H Z /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 Z /3De).All calculations were performed using H = 0.5.
We present in figure 6(a) the relaxation of the scaled pressure gradient (dP/dZ + 3/H 3 )/ β as a function of the downstream distance Z 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).
Recalling that the elastic stresses relax exponentially over a distance proportional to (3De/2H ), we replot in figure 6(b) the scaled pressure gradient, (5.3), as a function of the rescaled downstream distance 2H Z /3De in a log−linear plot.As a result, all curves become parallel to the green dashed line 100 exp(−2H Z /3De), thus confirming that the pressure gradient relaxes over a length scale ∼(3De/2H ), 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 PG 1 % , is approximately where we obtain that the prefactor 5.3 ± 0.5 is weakly dependent on De throughout the investigated range of Deborah numbers.Equation (5.4) and the scaling 3De/2H indicate that, in the exit channel, the appropriate Deborah number is based on the exit height, i.e.De exit = λq/2h = De/H .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 PG 1 % ≈ 239 ± 23 for De exit = De/H = 30, which should be contrasted with 250 of Debbaut et al. (1988), who studied numerically the flow through the planar 4 : 1 contraction.

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 P = p/(μ 0 q /2h 3 0 ) in the contraction as a function of De = λq/(2 h 0 ) for H = 0.5 and β = 0.05.For further clarification, figure 7(b) shows the first-order contribution P 1 = p 1 /(μ 0 q /2h 3 0 ) as a function of De = λq/(2 h 0 ), which is independent of β.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 = 0.5, β = 0.05 and = 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).
It is evident that the semi-analytical solution for the pressure drop in the contraction approaches the high-De asymptotic solution for De 0.4 and linearly decreases with the Deborah number.First, such an agreement for De / 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 P in the exit channel as a function of De for H = 0.5, β = 0.05, and L = 50.For De = 2, a long exit channel of L 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 P ,1 as a function of De, which is independent of β.In contrast to the total pressure drop P , the first-order contribution P ,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 ≡ H .The inset in figure 8(a) shows a comparison of our semi-analytical predictions (black dots) and finite-element simulation results (grey triangles) for P − P ,0 = β P ,1 as a function of De for H = 0.5, β = 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 1. Similar to the contraction, the pressure drop in the exit channel linearly decreases with De for De 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 β P ,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.

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.

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 − 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 As expected, (B1) simply represents the solution for the velocity and pressure drop of a Newtonian fluid with a constant viscosity μ 0 flowing in a straight channel of (non-dimensional) height H and length L. Substituting (B1a) into (3.6),we obtain the governing equations for the conformation tensor components in the exit channel at the leading order in β Equations (B2), similar to (3.6), represent a set of one-way coupled first-order semi-linear partial differential equations that can be solved first for Ã22,0 , followed by Ã12,0 and then for Ã11,0 .The solution of these equations is ) where are the reference distributions of the conformation tensor components at the outlet (Z = 1) of the non-uniform channel that can be obtained from (3.8), (3.9) and (3.10).
We note that, under the assumption of a fully developed flow in the entire exit channel so that U(η) = (3/2H )(1 − η 2 ), the governing equations for the conformation tensor components (B2) and their solution (B3)-(B5) are valid not only at O( β0 ) but for arbitrary values of β.
Finally, we note that the components of the conformation tensor at the walls of the exit channel (η = ±1) are given in (3.12), with H(Z) ≡ H . Thus, the conformation tensor components at the walls of the exit channel attain their fully relaxed values without spatial development.

B.1.1. Conformation tensor in the exit channel at low De numbers
At low Deborah numbers, we use (3.13) to obtain the reference distributions of the conformation tensor components at the beginning of the exit channel where, for a smooth geometry, we have assumed that H (1) = H (1) = 0.
Substituting (B6) into (B3), we obtain explicit expressions for the spatial relaxation of the conformation tensor components in the exit channel for De 1 18De 2 H (1)  Z =L = Ã11,0 (Z = 0, η) − Ã11,0 (Z = L, η).The three terms on the right-hand side of (B11) represent, respectively, the Newtonian solvent stress contribution, the elastic normal stress contribution and the elastic shear stress contribution to the pressure drop.
It is possible to express the first-order contribution P ,1 in terms of the difference between the conformation tensor components at the beginning and end of the exit channel.First, integrating (B2a) and (B2b) with respect to Z from L to 0, we obtain Substituting (B15) into (B11) provides the alternative expression for P ,1 Under the assumption that L is such that the elastic stresses reach their fully relaxed values by the end of the exit channel, (B16) shows that the first-order contribution P ,1 is independent of L since the steady-state values of Ã11,0 , Ã12,0 and Ã22,0 depend solely on the η coordinate.
Figure 1.Schematic illustration of the 2-D configuration consisting of a slowly varying and symmetric contraction of height 2h(z) and length (h).The contraction is connected to two long straight channels of height 2h 0 and 2h , respectively, up-and downstream and contains a viscoelastic fluid steadily driven by the imposed flow rate q.
.8a,b) and the Deborah and Weissenberg numbers De = λu c and Wi = λu c h 0 .(2.9a,b) Figure 2. Schematic illustration of the orthogonal curvilinear coordinates (ξ, η) for a slowly varying geometry.The coordinate ξ is constant along vertical grid lines, and η, defined in (2.11a,b), is constant along the curves going from left to right.
(b,d) the elastic shear (b) and axial normal stresses (d) at the end of the contraction as a function of the rescaled coordinate ζ = De(1 − η) 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)

Figure 7 .
Figure 7. Non-dimensional pressure drop for the Oldroyd-B fluid in a contracting channel in the ultra-dilute limit.(a) Dimensionless pressure drop P = p/(μ 0 q /2h 3 0 ) as a function of De = λq/(2 h 0 ) for β = 0.05.(b) First-order contribution P 1 = p 1 /(μ 0 q /2h 3 0 ) to the dimensionless pressure drop as a function of De = λq/(2 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 = 0.5.

)B. 2 .
Pressure drop in the exit channel at the first order in βUsing (2.21) and (3.27), the expressions for the pressure drop at O( β), P ,1 and the total pressure drop in the exit channel up to O( β), where Ã11,0 and Ã12,0 are given in (B4) and (B5) and [ Ã11,0 ] Z =0 (B14), the last term on the right-hand side of (B11) can be expressed

B. 2 . 1 .L
Pressure drop in the exit channel at O( β) in the low-De limitTo calculate the pressure drop P in the exit channel at low Deborah numbers, we use (B7b)-(B7c) and (B10).The elastic normal stress contribution to P ,1 is neglected terms multiplying exp(−2H L/[3De(1 − η 2 )]) stress and shear stress contributions, (B17) and (B20), provides the expression for the pressure drop at O( β) in the low-De limit P ,total pressure drop in the exit channel in the low-De limit is ) shows that for a smooth contraction with H (1) = H (1) = 0, the first non-vanishing viscoelastic contribution to the pressure drop in the exit channel at low Deborah numbers is only at O(De 3 ) as the O(De) and O(De 2 ) contributions are identically zero.B.2.2.Pressure drop in the exit channel at O( β) in the high-De limitTo calculate the pressure drop P in the exit channel at high Deborah numbers, we use (B9b)-(B9c) and (B10).The elastic normal stress contribution to P ,1 is Ã12,0 dZ , after neglecting terms multiplying exp(−2H L/ [3De(1 − η 2 )]) ≈ 0, is given as 0 L Ã12,0 dZ ≈ 3DeL H 2 η + 9De 2 (H 2 − 1) H 3 η(1 − η 2 ) for De 1. (B25)Combining the normal stress and shear stress contributions, (B23) and (B24), provides the expression for the pressure drop at O( β) in the high-De limit total pressure drop in the exit channel in the high-De limit is