Viscoelastic fluid flow in a slowly varying planar contraction: the role of finite extensibility on the pressure drop

Abstract

We analyse the steady viscoelastic fluid flow in slowly varying contracting channels of arbitrary shape and present a theory based on the lubrication approximation for calculating the flow rate–pressure drop relation at low and high Deborah ($De$) numbers. Unlike most prior theoretical studies leveraging the Oldroyd-B model, we describe the fluid viscoelasticity using a FENE-CR model and examine how the polymer chains’ finite extensibility impacts the pressure drop. We employ the low-Deborah-number lubrication analysis to provide analytical expressions for the pressure drop up to $O(De^4)$. We further consider the ultra-dilute limit and exploit a one-way coupling between the parabolic velocity and elastic stresses to calculate the pressure drop of the FENE-CR fluid for arbitrary values of the Deborah number. Such an approach allows us to elucidate elastic stress contributions governing the pressure drop variations and the effect of finite extensibility for all $De$. We validate our theoretical predictions with two-dimensional numerical simulations and find excellent agreement. We show that, at low Deborah numbers, the pressure drop of the FENE-CR fluid monotonically decreases with $De$, similar to the previous results for the Oldroyd-B and FENE-P fluids. However, at high Deborah numbers, in contrast to a linear decrease for the Oldroyd-B fluid, the pressure drop of the FENE-CR fluid exhibits a non-monotonic variation due to finite extensibility, first decreasing and then increasing with $De$. Nevertheless, even at sufficiently high Deborah numbers, the pressure drop of the FENE-CR fluid in the ultra-dilute and lubrication limits is lower than the corresponding Newtonian pressure drop.

1. Introduction

The ability to accurately predict the hydrodynamic features is at the core of understanding viscoelastic fluid flows. Such complex fluid flows may exhibit significantly different characteristics from Newtonian flows, even with a small concentration of polymer molecules present, giving rise to viscoelastic effects such as normal stress differences and extensional thickening (Bird et al. 1987; Steinberg 2021; Datta et al. 2022; Ewoldt & Saengow 2022).

One hydrodynamic feature that has received considerable attention in the fluid mechanics community is the relationship between the pressure drop $\triangle p$ and the flow rate $q$ in viscoelastic channel flows with spatially varying shapes. Over the years, the $q{-} \triangle p$ relation of viscoelastic fluid flows has been studied in different geometries, through numerical simulations (Szabo et al. 1997; Alves et al. 2003; Binding et al. 2006; Alves & Poole 2007; Zografos et al. 2020; Varchanis et al. 2022) and experimental measurements (Rothstein & McKinley 1999, 2001; Sousa et al. 2009; Ober et al. 2013; James & Roos 2021), and recently, via theoretical analysis (Pérez-Salas et al. 2019; Boyko & Stone 2022; Housiadas & Beris 2023, 2024c; Boyko et al. 2024; Hinch et al. 2024). For an overview of recent studies, the reader is referred to Boyko & Stone (2022) and Hinch et al. (2024).

The majority of previous numerical and experimental studies on the flow rate–pressure drop relation have focused on rapidly varying geometries with sharp corners, such as abrupt or hyperbolic contraction and contraction–expansion (constriction) channels (see e.g. Rothstein & McKinley 1999; Alves et al. 2003; Binding et al. 2006; Campo-Deaño et al. 2011; Keshavarz & McKinley 2016; Zografos et al. 2022). However, such rapidly varying geometries greatly complicate theoretical analysis. Therefore, to overcome this issue and enable asymptotic analysis, theoretical studies have considered instead a slowly varying geometry and exploited the narrowness of the geometry through the application of the lubrication theory (see e.g. Boyko & Stone 2022; Housiadas & Beris 2023,2024a, b, c, d). There have been numerous applications of lubrication theory to other viscoelastic fluid flows, such as thin films and tribology problems (Ro & Homsy 1995; Tichy 1996; Sawyer & Tichy 1998; Zhang et al. 2002; Saprykin et al. 2007; Ahmed & Biancofiore 2021; Gamaniel et al. 2021; Datt et al. 2022; Ahmed & Biancofiore 2023, 2024), as well as translation of a sphere near a rigid plane (Ardekani et al. 2007; Ruangkriengsin et al. 2024), and analysis of forces and torques acting on nearly touching spheres (Dandekar & Ardekani 2021). For example, Sari et al. (2024) recently explored the effect of fluid viscoelasticity in soft lubrication contacts using the Oldroyd-B model, considering elastohydrodynamic lubrication situations at low values of the Deborah number.

Using such a theoretical approach in conjunction with applying a perturbation expansion in powers of the Deborah number $De$ (see definition in $\S$ 2.1), Boyko & Stone (2022) studied the steady flow of an Oldroyd-B fluid in a slowly varying, arbitrarily shaped two-dimensional (2-D) channel and provided the expression for the $q{-} \triangle p$ relation up to $O(De^3)$ in the low-Deborah-number limit. Recently, Housiadas & Beris (2023) extended the analysis of Boyko & Stone (2022) to much higher asymptotic orders and provided analytical expressions for the pressure drop up to $O(De^8$ ) for different constitutive models, such as Oldroyd-B, Phan-Thien $-$ Tanner (PTT) (Phan-Thien & Tanner 1977; Phan-Thien 1978), Giesekus (Giesekus 1982) and a finitely extensible nonlinear elastic (FENE) model with the Peterlin approximation (FENE-P) (Bird et al. 1980, 1987). Their low-Deborah-number theoretical predictions for pressure drop using more complex constitutive models are very close to those of the Oldroyd-B model, showing a monotonic decrease in the scaled pressure drop with $De$ for the flow through a hyperbolic contraction (Housiadas & Beris 2023).

Recently, Hinch et al. (2024) and Boyko et al. (2024) analysed the flow of an Oldroyd-B fluid in a slowly varying 2-D channel in the high- $De$ limit using lubrication theory. Hinch et al. (2024) studied numerically the flow through a contraction, expansion and constriction for order-one Deborah numbers, and provided asymptotic solutions at high Deborah numbers. Boyko et al. (2024) studied the flow of the Oldroyd-B fluid in a slowly varying contraction considering the ultra-dilute limit, in which there is a one-way coupling between the Newtonian velocity and polymer stresses (Remmelgas et al. 1999; Moore & Shelley 2012; Li et al. 2019; Mokhtari et al. 2022). Such an approach allows for considerable theoretical progress beyond low $De$ , yielding semi-analytical expressions for the conformation tensor and pressure drop for arbitrary values of the Deborah number. For a contraction, Hinch et al. (2024) and Boyko et al. (2024) showed that the pressure drop of the Oldroyd-B fluid monotonically decreases with $De$ , scaling linearly with $De$ at high Deborah numbers, and identified two physical mechanisms responsible for the pressure drop reduction.

Although the Oldroyd-B model is the simplest viscoelastic model that combines viscous and elastic stresses and can be derived from kinetic theory, it has several shortcomings (Beris 2021; Hinch & Harlen 2021; Shaqfeh & Khomami 2021; Castillo-Sánchez et al. 2022; Stone et al. 2023). One well-known shortcoming of the Oldroyd-B model is that it allows the polymer chains, represented by elastic dumbbells, to be infinitely extensible (Bird et al. 1987). However, in reality, the polymer chains have a finite length. More importantly, theoretical and numerical predictions for the pressure drop reduction of an Oldroyd-B fluid in a contraction (Alves et al. 2003; Boyko & Stone 2022; Housiadas & Beris 2023; Boyko et al. 2024) are in contrast with the experiments showing a nonlinear increase in the pressure drop with $De$ for the flow of a Boger fluid through abrupt contraction–expansion and contraction geometries (Rothstein & McKinley 1999; Nigen & Walters 2002; Nigen & Walters 2002; Sousa et al. 2009). As pointed out by Hinch et al. (2024), this discrepancy might be resolved by using a more complex constitutive model that incorporates some form of extra dissipative stress instead of a simple Oldroyd-B model.

Different models, such as the FENE-CR model introduced by Chilcott & Rallison (1988) and the FENE-P model, incorporate the feature of finite extensibility through a nonlinear restoring force and include extra dissipation. Similar to the Oldroyd-B model, the FENE-CR model does not account for the shear-thinning effect and is suitable for describing constant shear-viscosity viscoelastic (Boger) fluids (James 2009). In contrast, the FENE-P model incorporates both the finite extensibility and the shear-thinning effect of viscoelastic fluids. Recently, Ahmed & Biancofiore (2023) used both the FENE-CR and FENE-P models to study the influence of viscoelastic effects in a thin lubricated contact with a boundary-driven flow in the low-Deborah-number limit.

There are several advantages of studying the FENE-CR model prior to the FENE-P model, particularly at high Deborah numbers. First, the FENE-CR model allows the study of elastic effects on the pressure drop without the influence of shear thinning in shear viscosity. Second, the FENE-CR model is more convenient for theoretical analysis. For example, in contrast to the conformation tensor components of the fully developed flow of a FENE-CR fluid in a straight channel, which have relatively simple expressions (see Appendix A), the corresponding expressions for the FENE-P fluid are more cumbersome (Cruz et al. 2005).

Nevertheless, it should be noted that at low $De$ , more complex constitutive models, such as PTT, Giesekus, FENE-P and FENE-CR, exhibit behaviour similar to Oldroyd-B due to the weak effect of additional microscopic features (Boyko & Stone 2024). Indeed, at low Deborah numbers, the PTT, Giesekus and FENE-P fluids showed only a slight difference in the pressure drop results compared with the Oldroyd-B fluid (Housiadas & Beris 2023). However, at high Deborah numbers, additional microscopic features, such as finite extensibility, become apparent and impact the elastic stresses (see e.g. Zografos et al. 2022). Therefore, one should anticipate significant differences between the predictions for the pressure drop of the Oldroyd-B and the more complex constitutive models, thus motivating further investigation.

In this work, we study the pressure-driven flow of the FENE-CR fluid in slowly varying, arbitrarily shaped, planar contracting channels using lubrication theory. Specifically, in the current work, we analyse the low-Deborah-number limit and the ultra-dilute limit, with the latter enabling us to explore arbitrary values of Deborah number. This is in contrast to Housiadas & Beris (2023), who considered the flow of a FENE-P fluid through a non-uniform channel at low $De$ . We first employ a perturbation expansion in powers of the Deborah number to calculate the non-dimensional pressure drop for the FENE-CR fluid up to $O(De^4)$ and then apply the Padé approximation (Housiadas 2017) to improve the convergence of the asymptotic series. We find that, at low Deborah numbers, the pressure drop of the FENE-CR fluid monotonically decreases with $De$ , similar to the Oldroyd-B and FENE-P fluid predictions.

To elucidate the pressure drop behaviour at high $De$ , we consider the ultra-dilute limit of small polymer concentration and leverage a one-way coupling between the parabolic velocity and polymer stresses to calculate the pressure drop for arbitrary values of the Deborah number. Such an approach allows us to study the elastic stress contributions governing the pressure drop variations and the effect of finite extensibility for all $De$ . We show that, at high Deborah numbers, in contrast to a linear pressure drop reduction of the Oldroyd-B fluid, the pressure drop of the FENE-CR fluid exhibits a non-monotonic variation, first decreasing and then increasing with $De$ . Nevertheless, in the ultra-dilute limit, the pressure drop of the FENE-CR fluid is lower than the corresponding Newtonian pressure drop even at sufficiently high Deborah numbers. We validate our theoretical predictions with 2-D finite-volume numerical simulations and find excellent agreement. However, as expected, at sufficiently high $De$ , our 2-D finite-volume numerical simulations, implementing the log-conformation formulation, suffer from accuracy and convergence difficulties due to the high-Weissenberg-number problem (Owens & Phillips 2002; Alves et al. 2021). Therefore, we believe that our theoretical results for the FENE-CR fluid in the ultra-dilute limit, valid at high Deborah numbers, are of fundamental importance for validating simulation predictions and advancing our understanding of viscoelastic channel flows.

Table 1 presents a summary of previous theoretical and numerical work on the pressure drop of viscoelastic fluids in slowly varying contraction geometries, highlighting the novelty of this study. It is evident from table 1 that all prior studies incorporating the effect of finite extensibility considered only low values of $De$ and predicted a monotonic pressure drop reduction, similar to the Oldroyd-B behaviour. In contrast, our study shows that the finite extensibility leads to the non-monotonic pressure drop behaviour at high $De$ . Furthermore, we believe that this study is an important step in understanding and resolving the theoretical/experimental discrepancy for the pressure drop of a Boger fluid through contraction geometries, discussed in the previous paragraphs.

Table 1. Theoretical and numerical studies of the pressure-driven flows of viscoelastic fluids in slowly varying contraction geometries that employed lubrication theory and focused on understanding the pressure drop behaviour.

2. Problem formulation and governing equations

Figure 1. Schematic illustration of the planar configuration consisting of a slowly varying and symmetric contraction of height $2h(z)$ and length $\ell$ ( $h\ll \ell$ ). Upstream of the contraction inlet, there is an entry channel of height $2h_{0}$ and length $\ell _{0}$ , and downstream of the contraction outlet, there is an exit channel of height $2h_{\ell }$ and length $\ell _{\ell }$ . The imposed flow rate $q$ results in a viscoelastic fluid flow through the geometry, and we aim to determine the pressure drop $\triangle p$ across the contraction region. We have indicated the qualitatively expected extension of polymers as the fluid flows through the contraction since the extension affects the fluid response in the FENE-CR description.

We study the incompressible steady flow of a viscoelastic fluid in a slowly varying and symmetric planar channel of height $2h(z)$ and length $\ell$ , where $h \ll \ell$ , as shown in figure 1. Motivated by the geometries used in previous experimental and numerical studies (see e.g. Szabo et al. 1997; Rothstein & McKinley 1999; Alves et al. 2003; Alves & Poole 2007; Campo-Deaño et al. 2011; Ober et al. 2013; Zografos et al. 2020; Boyko & Stone 2022; Boyko et al. 2024; Hinch et al. 2024), we assume that the inlet ( $z=0$ ) and outlet ( $z=\ell$ ) of the contraction are connected to two long straight channels of height $2h_{0}$ and $2h_{\ell }$ , and length $\ell _{0}$ and $\ell _{\ell }$ , respectively. We consider the fluid motion with the pressure distribution $p$ and velocity $\boldsymbol {u} =(u_z, u_y)$ induced by an imposed flow rate $q$ (per unit depth). Our primary interest in this work is to examine the pressure drop $\triangle p$ of a viscoelastic fluid over the contraction region at low and high Deborah numbers while incorporating the finite extensibility of polymer chains.

We consider low-Reynolds-number flows and neglect the fluid inertia. In this creeping flow limit, the governing equations are the continuity and momentum equations

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

Here, the stress tensor $\boldsymbol {\sigma }$ can be expressed as

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

where $-p{\boldsymbol{I}}$ is the pressure contribution, $2\mu _s{\boldsymbol{E}}$ is the viscous stress contribution of a Newtonian solvent with a constant viscosity $\mu _s$ , where ${\boldsymbol{E}}=(\boldsymbol {\nabla }\boldsymbol {u}+(\boldsymbol {\nabla }\boldsymbol {u})^{\mathrm {T}})/2$ is the rate-of-strain tensor, and $\boldsymbol {\tau}_{\!p}$ is the polymer contribution to the stress tensor.

To describe the viscoelastic rheology of the fluid, we use the FENE-CR model introduced by Chilcott & Rallison (1988). In contrast to the Oldroyd-B constitutive equation (Oldroyd 1950), the FENE-CR constitutive model considers polymer molecules as dumbbells with a finite extensibility $L$ relative to their value at equilibrium. However, the FENE-CR model does not account for the shear-thinning effect, which can be captured using the FENE-P model (Bird et al. 1987). For the FENE-CR model, the polymer contribution to the stress tensor $\boldsymbol {\tau}_{\!p}$ can be expressed in terms of the symmetric conformation tensor (or the deformation of the microstructure) ${\boldsymbol{A}}$ as (Chilcott & Rallison 1988; Alves et al. 2021)

(2.3) \begin{equation} \boldsymbol \tau _{\!p}=\frac {\mu _p}{\lambda }F({\boldsymbol{A}})({\boldsymbol{A}}-{\boldsymbol{I}}\,), \end{equation}

where $\mu _p$ is the polymer contribution to the shear viscosity at zero shear rate and $\lambda$ is the relaxation time. We also introduce the total zero-shear-rate viscosity $\mu _0=\mu _s+\mu _p$ .

The function $F({\boldsymbol{A}})$ in (2.3) accounts for the finite extensibility of polymers represented by elastic dumbbells and is modelled using the Warner spring function (Warner 1972),

(2.4) \begin{equation} F({\boldsymbol{A}})=\frac {1}{1-(\mathrm {tr}{\boldsymbol{A}})/L^{2}}, \end{equation}

where $\mathrm {tr}{\boldsymbol{A}}$ denotes the trace of the conformation tensor ${\boldsymbol{A}}$ .

At a steady state, the conformation tensor of the FENE-CR model satisfies (Chilcott & Rallison 1988)

(2.5) \begin{equation} \boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {\nabla }{\boldsymbol{A}}-(\boldsymbol {\nabla }\boldsymbol {u})^{\mathrm {T}}\boldsymbol {\cdot }{\boldsymbol{A}}-{\boldsymbol{A}}\boldsymbol {\cdot }(\boldsymbol {\nabla }\boldsymbol {u})=-\frac {F({\boldsymbol{A}})}{\lambda }({\boldsymbol{A}}-{\boldsymbol{I}}\,). \end{equation}

For large values of $L$ , the function $F({\boldsymbol{A}})$ tends to 1, so that the FENE-CR model, given in (2.3) and (2.5), reduces to the steady form of the Oldroyd-B constitutive equation.

2.1. Non-dimensionalisation

We analyse the viscoelastic fluid flow through a narrow slowly varying channel, in which the channel height is much smaller than the channel length, $h(z)\ll \ell$ . Therefore, for the non-dimensionalisation of the viscoelastic flow problem, we introduce dimensionless variables based on the lubrication theory (Tichy 1996; Zhang et al. 2002; Saprykin et al. 2007; Ahmed & Biancofiore 2021; Boyko & Stone 2022; Boyko & Stone 2022; Boyko et al. 2024),

(2.6a) \begin{equation} Z=\frac {z}{\ell }, \qquad Y=\frac {y}{h_{\ell }}, \qquad U_z=\frac {u_z}{u_c}, \qquad U_y=\frac {u_y}{\epsilon u_c}, \end{equation}

(2.6b) \begin{equation} P=\frac {p}{\mu _0 u_c \ell /h^2_{\ell }}, \qquad \triangle P=\frac {\triangle p}{\mu _0 u_c \ell /h^2_{\ell }}, \qquad H(Z)=\frac {h(z)}{h_{\ell }}, \end{equation}

(2.6c) \begin{equation} {\tilde A}_{zz}=\epsilon ^2 A_{zz}, \qquad {\tilde A}_{yz}=\epsilon A_{yz}, \qquad {\tilde A}_{yy}=A_{yy}, \end{equation}

(2.6d) \begin{equation} \mathcal {T}_{p,zz}=\frac {\epsilon ^2\ell }{\mu _0 u_c}\tau _{p,zz}, \qquad \mathcal {T}_{p,yz}=\frac {\epsilon \ell }{\mu _0 u_c}\tau _{p,yz},\qquad \mathcal {T}_{p,yy}=\frac {\ell }{\mu _0 u_c}\tau _{p,yy}, \end{equation}

where $u_c=q/2h_{\ell }$ is the characteristic velocity scale, $q$ is the imposed flow rate per unit depth and $h_{\ell }$ is the half-height at $z=\ell$ . In addition, we introduce the aspect ratio of the configuration, which is assumed to be small,

(2.7) \begin{equation} \epsilon =\frac {h_{\ell }}{\ell }\ll 1, \end{equation}

the contraction ratio,

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

the viscosity ratios,

(2.9) \begin{equation} \tilde \beta =\frac {\mu _p}{\mu _s+\mu _p}=\frac {\mu _p}{\mu _0} \quad \text {and} \quad \beta =1-\tilde \beta =\frac {\mu _s}{\mu _0}, \end{equation}

and the Deborah and Weissenberg numbers,

(2.10) \begin{equation} De=\frac {\lambda u_c}{\ell } \quad \text {and} \quad Wi=\frac {\lambda u_c}{h_{\ell }}. \end{equation}

Finally, we note that the fluid inertia is negligible, provided the reduced Reynolds number is small,

(2.11) \begin{equation} \epsilon Re=\epsilon \frac {\rho u_c h_{\ell }}{\mu _0}=\frac {\rho q h_{\ell }}{2\mu _0\ell }\ll 1, \end{equation}

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

Note that we have defined both the Deborah and Weissenberg numbers. Although the Deborah and Weissenberg numbers are equivalent in many steady flows, in lubrication flows, they have different orders of magnitude due to the two distinct length scales. The Deborah number $De$ is the ratio of the relaxation time of the fluid, $\lambda$ , to the residence time in the non-uniform region, $\ell /u_c$ (Tichy 1996; Zhang et al. 2002; Saprykin et al. 2007; Ahmed & Biancofiore 2021; Boyko & Stone 2022; Ahmed & Biancofiore 2023; Housiadas & Beris 2023; Boyko et al. 2024; Hinch et al. 2024). The Weissenberg number $Wi$ is the product of the relaxation time of the fluid, $\lambda$ , and the characteristic shear rate of the flow, $u_c/h_{\ell }$ , and is related to the Deborah number through $De=\epsilon Wi$ . Therefore, for lubrication flows in narrow geometries with $\epsilon \ll 1$ , the Deborah number can be small, but $Wi=O(1)$ . In addition to the Deborah number $De=\lambda q/(2\ell h_{\ell })$ based on the exit height, we can introduce the Deborah number $De_{entry}=\lambda q/(2\ell h_{0})$ based on the entry height; the two Deborah numbers are related through $De_{entry}=De/H_0$ .

2.2. Non-dimensional governing equations in Cartesian coordinates and pressure drop

Substituting the non-dimensional variables (2.6)–(2.10) into the governing equations (2.1)–(2.5) and considering the leading order in $\epsilon$ , we obtain

(2.12a) \begin{align} \frac {\partial U_z}{\partial Z}+\frac {\partial U_y}{\partial Y}=0,\qquad\qquad\qquad\qquad\qquad \end{align}

(2.12b) \begin{align} \frac {\partial P}{\partial Z}=(1-\tilde \beta )\frac {\partial ^2 U_z}{\partial Y^2}+\frac {\tilde \beta }{De}\left (\frac {\partial (\mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})\tilde A_{zz})}{\partial Z}+ \frac {\partial (\mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})\tilde A_{yz})}{\partial Y}\right ),\quad \end{align}

(2.12c) \begin{align} \frac {\partial P}{\partial Y}=0, \qquad\qquad\qquad\qquad\qquad\qquad\end{align}

(2.12d) \begin{align} U_z\frac {\partial \tilde A_{zz}}{\partial Z}+U_y\frac {\partial \tilde A_{zz}}{\partial Y}-2\frac {\partial U_z}{\partial Z}\tilde A_{zz}-2\frac {\partial U_z}{\partial Y}\tilde A_{yz}=-\frac {\mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})}{De}\tilde {A}_{zz}, \quad\end{align}

(2.12e) \begin{align} U_z\frac {\partial \tilde A_{yz}}{\partial Z}+U_y\frac {\partial \tilde A_{yz}}{\partial Y}-\frac {\partial U_y}{\partial Z}\tilde A_{zz}-\frac {\partial U_z}{\partial Y}\tilde A_{yy}=-\frac {\mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})}{De}\tilde {A}_{yz},\qquad \end{align}

(2.12f) \begin{align} U_z\frac {\partial \tilde A_{yy}}{\partial Z}+U_y\frac {\partial \tilde A_{yy}}{\partial Y}-2\frac {\partial U_y}{\partial Z}\tilde A_{yz}-2\frac {\partial U_y}{\partial Y}\tilde A_{yy}=-\frac {\mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})}{De}(\tilde {A}_{yy}-1), \end{align}

where

(2.13) \begin{equation} \mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})=\frac {1}{1-\dfrac {1}{L^2 \epsilon ^2}(\tilde {A}_{zz}+\epsilon ^2 \tilde {A}_{yy})}\approx \frac {1}{1-\tilde {A}_{zz}/( L^2 \epsilon ^2)}. \end{equation}

From the $y$ -momentum equation, (2.12 $c$ ), it follows that $P=P(Z)+O(\epsilon ^2)$ , i.e. the pressure is constant across a cross-section but varies along the $z$ -direction. Under the non-dimensionalisation (2.6c), the right-hand side of (2.12 $d$ ) becomes $-(\mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})/De)$ $(\tilde {A}_{zz}-\epsilon ^2)$ . Thus, at the leading order in $\epsilon$ , we have $-(\mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})/De)\tilde {A}_{zz}$ .

For lubrication flows through the slowly varying geometries that we consider, (2.13) clearly indicates that the finite extensibility is governed by the dimensionless parameter $L^2 \epsilon ^2$ rather than $L^2$ (Ahmed & Biancofiore 2023; Housiadas & Beris 2023). Although we consider $\epsilon \ll 1$ , since the realistic values of $L^2$ are typically large (see e.g. Remmelgas et al. (1999); Rothstein & McKinley (1999)), we may have $L^2 \epsilon ^2=O(1)$ .

The corresponding boundary conditions on the velocity are

(2.14a–d) \begin{align} U_{z}(H(Z),Z)=0, \quad U_{y}(H(Z),Z)=0, \quad \frac {\partial U_{z}}{\partial Y}(0,Z)=0, \quad \int ^{H(Z)}_0 U_{z}(Y,Z)\mathrm {d}Y=1. \end{align}

These boundary conditions respectively represent the no-slip and no-penetration 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 flow of a FENE-CR fluid in the straight entry channel, given by the Poiseuille velocity profile, and the corresponding conformation tensor (see the derivation in Appendix A)

(2.15a) \begin{align} \tilde A_{zz}=L^2\epsilon ^2+L^3\epsilon ^3\frac {L\epsilon -\sqrt {L^2\epsilon ^2+72De^2 Y^2/H^6_0}}{36 De^2 Y^2/H^6_0},\quad \end{align}

(2.15b) \begin{align} \tilde {A}_{yz}=L\epsilon \frac {L\epsilon -\sqrt {L^2\epsilon ^2+72De^2 Y^2/H^6_0}}{12De Y/H^3_0} \quad \text {and} \quad \tilde {A}_{yy}=1. \end{align}

Following the steps outlined in Appendix B, the non-dimensional pressure drop $\triangle P=P(0)-P(1)$ across the non-uniform region can be expressed as

(2.16) \begin{align} \triangle P&= (1-\tilde {\beta })\triangle \hat {P}+\frac {\tilde {\beta }}{De }\int _{0}^{H(0)}[\mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})\tilde A_{zz}\hat {U}_{z}]_{Z=0}\mathrm {d}Y-\frac {\tilde {\beta }}{De}\int _{0}^{H(1)}[\mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})\tilde A_{zz}\hat {U}_{z}]_{Z=1}\mathrm {d}Y \nonumber \\ &\quad +\frac {\tilde {\beta }}{De }\int _{0}^{1}\int _{0}^{H(Z)}\mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})\tilde A_{zz}\frac {\partial \hat {U}_{z}}{\partial Z}\mathrm {d}Y\mathrm {d}Z+\frac {\tilde {\beta }}{De }\int _{0}^{1}\int _{0}^{H(Z)}\mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})\tilde A_{yz}\frac {\partial \hat {U}_{z}}{\partial Y}\mathrm {d}Y\mathrm {d}Z. \end{align}

Here, the function $\mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})$ is given in (2.13), and $\triangle \hat {P}$ and $\hat {U}_{z}$ are the corresponding pressure drop and axial velocity of a Newtonian fluid given by (Boyko & Stone 2022)

(2.17a,b) \begin{align} \triangle \hat {P}=3\int _{0}^{1}\frac {\mathrm {d}Z}{H(Z)^{3}}, \qquad \hat {U}_{z}=\frac {3}{2}\frac {H(Z)^2-Y^2}{H(Z)^3}. \end{align}

Equation (2.16) represents the expression for the non-dimensional pressure drop previously obtained from an application of the reciprocal theorem in a slowly varying channel (Boyko & Stone 2021, 2022). The first term on the right-hand side of (2.16) represents the contribution of the Newtonian solvent to the pressure drop. The second and third terms represent the contribution of the elastic normal stresses at the inlet and outlet of the non-uniform channel. Finally, the fourth and fifth terms represent the contribution of the elastic normal stresses and elastic shear stresses within the non-uniform channel.

3. Low-Deborah-number lubrication analysis

In this section, we employ the low-Deborah-number lubrication analysis to derive asymptotic expressions for the velocity, conformation tensor and pressure drop of a weakly viscoelastic FENE-CR fluid up to $O(De^4)$ . To this end, we expand the velocity, pressure drop and conformation tensor components into perturbation series in the Deborah number $De\ll 1$ ,

(3.1) \begin{equation} \begin{pmatrix} U_z\\ U_y\\ P\\ \tilde A_{zz}\\ \tilde A_{yy}\\ \tilde A_{yz} \end{pmatrix} = \begin{pmatrix} U_{z,0}\\ U_{y,0}\\ P_0\\ \tilde A_{zz,0}\\ \tilde A_{yy,0}\\ \tilde A_{yz,0} \end{pmatrix}+De \begin{pmatrix} U_{z,1}\\ U_{y,1}\\ P_1\\ \tilde A_{zz,1}\\ \tilde A_{yy,1}\\ \tilde A_{yz,1} \end{pmatrix}+De^2\begin{pmatrix} U_{z,2}\\ U_{y,2}\\ P_2\\ \tilde A_{zz,2}\\ \tilde A_{yy,2}\\ \tilde A_{yz,2} \end{pmatrix}+\ldots \:. \end{equation}

As noted by Boyko & Stone (2022), in the weakly viscoelastic and lubrication limits, $De\ll 1$ and $\epsilon \ll 1$ , it is sufficient to apply the boundary conditions on the velocity (2.14) to find the flow field, conformation tensor components and pressure drop at each order in $De$ . Indeed, the iterative structure of the solution eliminates the need to use the boundary condition (2.15) on the conformation tensor (Black & Denn 1976; Boyko & Stone 2022; Housiadas & Beris 2023). For example, considering the leading and first order in $De$ , we find

(3.2) \begin{equation} \tilde A_{zz,0}=0,\qquad \tilde A_{yz,0}=0,\qquad \tilde A_{yy,0}=1, \end{equation}

(3.3) \begin{equation} \tilde A_{zz,1}=0,\qquad \tilde A_{yz,1}=\frac {\partial U_{z,0}}{\partial Y}, \qquad \tilde A_{yy,1}=2\frac {\partial U_{y,0}}{\partial Y}. \end{equation}

In Appendix C, we provide a detailed derivation of the expressions for the pressure drop of the FENE-CR fluid in the low- $De$ limit up to $O(De^4)$ . We obtain that the expressions for the pressure drop at the leading, first and second order in $De$ are the same for the FENE-CR and Oldroyd-B fluids, and are given by

(3.4a,b) \begin{align} \triangle P_{0}=3\int _{0}^{1}\frac {\mathrm {d}Z}{H(Z)^{3}}, \qquad \triangle P_1 = \frac {9}{2}\tilde \beta \left (\frac {1}{H(0)^4}-\frac {1}{H(1)^4} \right ), \end{align}

(3.4c) \begin{align} \triangle P_2=\frac {324}{35}\tilde \beta \int ^1_0\left (\frac {14 H'(Z)^2}{H(Z)^7} - \frac {3 H''(Z)}{H(Z)^6} \right ) \mathrm {d}Z.\qquad \end{align}

Interestingly, unlike the FENE-CR fluid, the pressure drop of the FENE-P fluid is different from the Oldroyd-B case at $O(De^2)$ and depends on finite extensibility through $L^2\epsilon ^2$ , as recently shown by Housiadas & Beris (2023).

At the third order in $De$ , the pressure drop of the FENE-CR fluid is different from the Oldroyd-B fluid due to the finite extensibility and is given as

(3.5) \begin{align} \triangle P_{3}&=- \frac {2673\tilde {\beta }}{70L^2\epsilon ^2}\left (\frac {1}{H(0)^{8}}-\frac {1}{H(1)^{8}}\right ) \nonumber \\ &\quad +\frac {648\tilde {\beta }(9-\tilde {\beta })}{35}\left (\frac {H'(0)^{2}}{H(0)^{8}}-\frac {H'(1)^{2}}{H(1)^{8}}\right )-\frac {216\tilde {\beta }(8-\tilde {\beta })}{35}\left (\frac {H''(0)}{H(0)^{7}}-\frac {H''(1)}{H(1)^{7}}\right )\!.\end{align}

From (3.5), it follows that $\triangle P_{3}$ may increase, decrease or not change the total pressure drop of the FENE-CR fluid, depending on the geometry. For a contraction ( $H(0)\gt H(1)$ ), the first term, which depends on finite extensibility through $L^2\epsilon ^2$ and distinguishes the FENE-CR fluid from the Oldroyd-B fluid, leads to an increase in the pressure drop. However, for an expansion ( $H(0)\lt H(1)$ ), the first term leads to a decrease in the pressure drop, and for a constriction ( $H(1)=H(0)$ ), it does not contribute to the pressure drop. We also note that our expression for the pressure drop $\triangle P_3$ of the FENE-CR fluid is similar to the expression for the pressure drop of the FENE-P fluid at $O(De^3)$ , albeit a different number in the coefficient of the first term in (3.5) (Housiadas & Beris 2023).

Finally, at the fourth order in $De$ , the resulting expression for the pressure drop of the FENE-CR fluid is

(3.6) \begin{eqnarray} \triangle P_4&=&\frac {3888\tilde \beta (8 \tilde \beta +25)}{175L^2 \epsilon ^2}\left [\frac {H'(1)}{H(1)^{10}}-\frac {H'(0)}{H(0)^{10}} \right ]+\frac {648 \tilde \beta }{175L^2 \epsilon ^2}\int _0^1\left [a_{1}\frac {H''}{H^{10}} +a_{2}\frac {H'^2}{H^{11}} \right ]\mathrm {d}Z \nonumber \\ &&+\int _0^1 \left [a_3 \frac {H''^2}{H^9} +a_4\frac {H''' H'}{H^9}+a_5\frac {H^3 H''''}{H^{11}}+ a_6\frac {H'^4}{H^{11}} + a_7 \frac {H'^2 H''}{H^{10}} \right ]\mathrm {d}Z\nonumber \\ &&+\,a_8\left [\frac {H'''(0)}{H(0)^8}-\frac {H'''(1)}{H(1)^8} \right ]+a_9\left [\frac {H'(1) H''(1)}{H(1)^9}-\frac {H'(0) H''(0)}{H(0)^9}\right ]\nonumber \\ &&+\,a_{10}\left [\frac {H'(0)^3}{H(0)^{10}}-\frac {H'(1)^3} {H(1)^{10}}\right ] , \end{eqnarray}

where the coefficients $a_1, . . ., a_{10}$ are summarised in table 2.

Table 2. Coefficients appearing in the expression (3.6) for the fourth-order pressure drop $\triangle P_4$ of the FENE-CR fluid in a planar contracting channel.

The first two terms on the right-hand side of (3.6) depend on $L^2\epsilon ^2$ , and thus clearly distinguish the analytical prediction for $\triangle P_4$ of the FENE-CR fluid from the Oldroyd-B fluid. For the Oldroyd-B fluid, our analytical result for $\triangle P_4$ fully agrees with the solution of Housiadas & Beris (2023) when accounting for the differences in the non-dimensionalisation. However, as expected based on the previous orders, our expression for the fourth-order pressure drop of the FENE-CR fluid differs from the expression for the FENE-P fluid given by Housiadas & Beris (2023). Specifically, the first two terms in (3.6) that include $L^2\epsilon ^2$ appear in the fourth-order expressions for both FENE-CR and FENE-P fluids with different coefficients. Furthermore, the expression for the FENE-P fluid has an additional term of the form of $\int ^1_0 H(Z)^{-11} \mathrm {d}Z$ that depends on $1/L^4\epsilon ^4$ .

For a given flow rate, we have determined the dimensionless pressure drop $\triangle P=\triangle p/(\mu _{0}q\ell /2h_{\ell }^{3})$ of a FENE-CR fluid as a function of the shape function $H(Z)$ , the viscosity ratio $\tilde {\beta }$ , the parameter $L^2\epsilon ^2$ and the Deborah number $De$ up to $O(De^4)$ ,

(3.7) \begin{equation} \triangle P=\triangle P_{0}+De\triangle P_{1}+De^{2}\triangle P_{2}+De^{3}\triangle P_{3}+De^{4}\triangle P_{4}+O(\epsilon ^{2},De^{5}), \end{equation}

where the expressions for $\triangle P_{0}$ , $\triangle P_{1}$ , $\triangle P_{2}$ , $\triangle P_{3}$ and $\triangle P_{4}$ are given in (3.4 $a$ ), (3.4 $b$ ), (3.4 $c$ ), (3.5) and (3.6), respectively. Physically, the non-dimensional quantity $\triangle P=\triangle p/(\mu _{0}q\ell /2h_{\ell }^{3})$ represents the dimensionless flow resistance ( $\triangle p/q$ ) for a given geometry.

Having the low- $De$ asymptotic expressions for $\triangle P_{0}$ , $\triangle P_{1}$ , $\triangle P_{2}$ , $\triangle P_{3}$ and $\triangle P_{4}$ , we can improve the convergence of the asymptotic series (3.7) by using the diagonal Padé [2/2] approximation (Hinch 1991; Housiadas 2017; Housiadas & Beris 2023),

(3.8)

It should be noted that Housiadas & Beris (2023) extended the low-Deborah-number lubrication analysis to much higher asymptotic orders and provided analytical expressions for the pressure drop of the Oldroyd-B and FENE-P fluids up to $O(De^8)$ . Nevertheless, as shown for the Oldroyd-B fluid, the low- $De$ perturbation solutions obtained from the Padé approximations remain indistinguishable when adding more terms in the asymptotic series beyond $O(De^4)$ .

4. Low- $\tilde {\beta}$ lubrication analysis

In the previous section, we derived analytical expressions for the non-dimensional pressure drop of a FENE-CR fluid in a non-uniform channel of arbitrary shape $H(Z)$ in the low-Deborah-number limit, $De\ll 1$ . However, as pointed out by Boyko et al. (2024) and Hinch et al. (2024), the low-Deborah-number asymptotic analysis cannot accurately predict the pressure drop at high $De$ where there are significant elastic stresses.

In this section, we employ orthogonal curvilinear coordinates and consider the ultra-dilute limit, $\tilde {\beta }=\mu _{p}/\mu _{0}\ll 1$ (Remmelgas et al. 1999; Moore & Shelley 2012; Li et al. 2019; Mokhtari et al. 2022; Boyko et al. 2024; Hinch et al. 2024), which allows us to analyse the pressure drop and conformation tensor at high Deborah numbers.

4.1. Orthogonal curvilinear coordinates for a slowly varying geometry

For our low- $\tilde {\beta }$ lubrication analysis, we first transform the geometry of the contraction from the Cartesian coordinates ( $Z,Y$ ) to orthogonal curvilinear coordinates ( $\xi ,\eta$ ) with the mapping (Boyko et al. 2024; Hinch et al. 2024)

(4.1) \begin{equation} \xi =Z-\frac {1}{2}\epsilon ^{2}\frac {H'(Z)}{H(Z)}(H(Z)^{2}-Y^{2})+O(\epsilon ^{4}), \qquad \eta =\frac {Y}{H(Z)}, \end{equation}

and use $U\boldsymbol{e}_{\xi}+V\boldsymbol{e}_{\eta}$ and $\tilde A_{11}\boldsymbol{e}_{\xi}\boldsymbol{e}_{\xi}+\tilde A_{12}(\boldsymbol{e}_{\xi}\boldsymbol{e}_{\eta}+\boldsymbol{e}_{\eta}\boldsymbol{e}_{\xi})+\tilde A_{22}\boldsymbol{e}_{\eta}\boldsymbol{e}_{\eta}$ to denote the components of velocity and conformation tensor in curvilinear coordinates ( $\xi ,\eta$ ).

The corresponding components of the non-dimensional velocity field and conformation tensor in different coordinates are related through

(4.2a) \begin{align} U_{z}=U-\epsilon ^{2}\eta H'(\xi ) V, \qquad U_{y}=\eta H'(\xi ) U+V, \,\,\end{align}

(4.2b) \begin{align} \tilde {A}_{zz}=\tilde {A}_{11}+O(\epsilon ^{2}),\qquad\qquad\qquad\quad \end{align}

(4.2c) \begin{align} \tilde {A}_{zy}=\tilde {A}_{12}+\eta H'(\xi )\tilde {A}_{11}+O(\epsilon ^{2}), \qquad\qquad\end{align}

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

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

4.2. Non-dimensional governing equations in orthogonal curvilinear coordinates

Using the mapping (4.1), the governing equations (2.12), (2.13) and the corresponding boundary conditions (2.14), (2.15) in curvilinear coordinates (Boyko et al. 2024; Hinch et al. 2024) take the form

(4.3a) \begin{align} \frac {\partial (HU)}{\partial Z}+\frac {\partial V}{\partial \eta }=0,\qquad\qquad\qquad\qquad\qquad \end{align}

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

(4.3c) \begin{align} U\frac {\partial \tilde {A}_{11}}{\partial Z}+\frac {V}{H}\frac {\partial \tilde {A}_{11}}{\partial \eta }-2\frac {\partial U}{\partial Z}\tilde {A}_{11}-\frac {2}{H}\frac {\partial U}{\partial \eta }\tilde {A}_{12}=-\frac {\mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})}{De}\tilde {A}_{11}, \qquad\end{align}

(4.3d) \begin{align} U\frac {\partial \tilde {A}_{12}}{\partial Z}+\frac {V}{H}\frac {\partial \tilde {A}_{12}}{\partial \eta }-H\frac {\partial }{\partial Z}\left (\frac {V}{H}\right )\tilde {A}_{11}-\frac {1}{H}\frac {\partial U}{\partial \eta }\tilde {A}_{22}=-\frac {\mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})}{De}\tilde {A}_{12}, \end{align}

(4.3e) \begin{align} U\frac {\partial \tilde {A}_{22}}{\partial Z}+\frac {V}{H}\frac {\partial \tilde {A}_{22}}{\partial \eta }-2H\frac {\partial }{\partial Z}\left (\frac {V}{H}\right )\tilde {A}_{12}+2\frac {\partial U}{\partial Z}\tilde {A}_{22}=-\frac {\mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})}{De}(\tilde {A}_{22}-1), \end{align}

where

(4.4) \begin{equation} \mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})=\frac {1}{1-\dfrac {1}{L^2 \epsilon ^2}( \tilde {A}_{11}+\epsilon ^2 \tilde {A}_{22})}\approx \frac {1}{1-\tilde {A}_{11}/(L^2 \epsilon ^2)}, \end{equation}

subject to the boundary conditions

(4.5a–d) \begin{align} U(Z,1)=0,\quad V(Z,1)=0,\quad \frac {\partial U} {\partial \eta }(Z,0)=0,\quad H(Z)\int _{0}^{1}U(Z,\eta )\mathrm {d}\eta =1, \end{align}

and

(4.6a) \begin{align} \tilde A_{11}(0,\eta )=L^2\epsilon ^2+L^3\epsilon ^3\frac {L\epsilon -\sqrt {L^2\epsilon ^2+72De^2 \eta ^2/H^4_0}}{36 De^2 \eta ^2/H^4_0}, \qquad\end{align}

(4.6b) \begin{align} \tilde {A}_{12}(0,\eta )=L\epsilon \frac {L\epsilon -\sqrt {L^2\epsilon ^2+72De^2 \eta ^2/H^4_0}}{12De \eta /H^2_0} \quad \text {and} \quad \tilde {A}_{22}(0,\eta )=1. \end{align}

Following similar steps as in Appendix B and using the integral constraint (4.5 $d$ ), the non-dimensional pressure drop can be expressed in curvilinear coordinates as

(4.7) \begin{align} \triangle P= \underset {\mathrm {Solvent\,stress\,contribution}}{\underbrace {3(1-\tilde {\beta })\int _{0}^{1}\frac {\mathrm {d}Z}{H(Z)^{3}}}} +\underset {\mathrm {Elastic\,shear\,stress\,contribution}}{\underbrace {-\frac {3\tilde {\beta }}{De}\int _{0}^{1}\left [\frac {1}{H(Z)}\int _{0}^{1}\eta \mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})\tilde {A}_{12}\mathrm {d}\eta \right ]\mathrm {d}Z}} \qquad\qquad\nonumber \\ +\underset {\mathrm {Elastic\,normal\,stress\,contribution}}{\underbrace {\frac {3\tilde {\beta }}{2De}\left (\int _{0}^{1}(1-\eta ^{2})\left [\mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})\tilde {A}_{11}\right ]^{0}_{1}\mathrm {d}\eta -\int _{0}^{1}\left [\frac {H'(Z)}{H(Z)}\left (\int _{0}^{1}(1-\eta ^{2})\mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})\tilde {A}_{11}\mathrm {d}\eta \right )\right ]\mathrm {d}Z\right )}}, \end{align}

where $[\mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})\tilde {A}_{11}]^{0}_{1}=\mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})\tilde {A}_{11}|_{Z=0}-\mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})\tilde {A}_{11}|_{Z=1}$ .

Equation (4.7) represents the pressure drop in curvilinear coordinates and is an analogue of (2.16), written in Cartesian coordinates. The first term on the right-hand side of (4.7) represents the viscous contribution of the Newtonian solvent to the pressure drop. The second term represents the contribution of the elastic shear stresses and the last term represents the contribution of the elastic normal stresses to the pressure drop.

4.3. Velocity, conformation and pressure drop in the ultra-dilute limit

Next, we consider the ultra-dilute limit, $\tilde {\beta } \ll 1$ , representing a one-way coupling between the velocity and pressure fields and the conformation tensor. At the leading order in $\tilde {\beta }$ , the velocity field of the FENE-CR fluid is parabolic, similar to Newtonian and Oldroyd-B fluids, and is given as

(4.8a,b) \begin{align} U=\frac {3}{2}\frac {1}{H(Z)}(1-\eta ^{2})\qquad \hbox {and}\qquad V\equiv 0. \end{align}

We note that in orthogonal curvilinear coordinates, the velocity in the $\eta$ -direction is identically zero at $O(\tilde {\beta }^{0})$ , in contrast to the Cartesian coordinates where $U_{y}=(3/2)H'(Z)Y(H(Z)^{2}-Y^{2})/H(Z)^{4}$ . As pointed out by Boyko et al. (2024), the latter fact significantly simplifies the theoretical analysis.

Substituting (4.8) into (4.3c) $-$ (4.3e) and using (4.4), we obtain the simplified equations for the conformation tensor components of the FENE-CR fluid at leading order in $\tilde {\beta }$ ,

(4.9a) \begin{align} U\frac {\partial \tilde {A}_{22}}{\partial Z}+2\frac {\partial U}{\partial Z}\tilde {A}_{22}=-\frac {1}{De}\frac {1}{1-\tilde {A}_{11}/(L^2 \epsilon ^2)}(\tilde {A}_{22}-1), \quad\end{align}

(4.9b) \begin{align} U\frac {\partial \tilde {A}_{12}}{\partial Z}-\frac {1}{H}\frac {\partial U}{\partial \eta }\tilde {A}_{22}=-\frac {1}{De}\frac {1}{1-\tilde {A}_{11}/(L^2 \epsilon ^2)}\tilde {A}_{12}, \qquad\end{align}

(4.9c) \begin{align} U\frac {\partial \tilde {A}_{11}}{\partial Z}-2\frac {\partial U}{\partial Z}\tilde {A}_{11}-\frac {2}{H}\frac {\partial U}{\partial \eta }\tilde {A}_{12}=-\frac {1}{De}\frac {1}{1-\tilde {A}_{11}/(L^2 \epsilon ^2)}\tilde {A}_{11}, \end{align}

where $U$ is given in (4.8 $a$ ).

Equation (4.9) represents a set of coupled first-order semi-linear partial differential equations that should be solved to obtain $\tilde {A}_{22}$ , $\tilde {A}_{12}$ and $\tilde {A}_{11}$ for the FENE-CR fluid. When $L^2\epsilon ^2 \to \infty$ , corresponding to the Oldroyd-B fluid, (4.9) reduces to a set of one-way coupled equations, allowing us to derive semi-analytical expressions for the conformation tensor for arbitrary values of the Deborah number in the ultra-dilute limit (Boyko et al. 2024). Furthermore, Boyko et al. (2024) and Hinch et al. (2024) provided analytical expressions for the conformation tensor and the pressure drop of the Oldroyd-B fluid in the high-Deborah-number limit. In particular, the pressure drop of the Oldroyd-B fluid across the non-uniform channel in the high- $De$ limit is

(4.10) \begin{equation} \triangle P=\underset {\mathrm {Solvent\,stress}}{\underbrace {3(1-\tilde {\beta })\int _{0}^{1}\frac {\mathrm {d}Z}{H(Z)^{3}}}}+\underset {\mathrm {Elastic\,shear\,stress}}{\underbrace {3\tilde {\beta }H_0^{-2}\int _{0}^{1}\frac {\mathrm {d}Z}{H(Z)}}}+\underset {\mathrm {Elastic\,normal\,stress}}{\underbrace {\frac {9}{5}\tilde {\beta }De\left(H_0^{-4}-H_{0}^{-2}\right)}}\;\mathrm {for}\; De\gg 1. \end{equation}

The coupling between equations in (4.9) greatly complicates the analytical progress, particularly in the high- $De$ asymptotic limit for the FENE-CR fluid. Nevertheless, examining the expressions in (4.9), we observe that for a given value of $\eta \in [0,1]$ , (4.9) represent a set of first-order ordinary differential equations for $\tilde {A}_{22}$ , $\tilde {A}_{12}$ and $\tilde {A}_{11}$ of the FENE-CR fluid. Therefore, we solve numerically the coupled equations (4.9) subject to the boundary conditions (4.6) using MATLAB’s routine ode45 and obtain the distribution of $\tilde {A}_{22}$ , $\tilde {A}_{12}$ and $\tilde {A}_{11}$ in a contraction for different values of $De$ and $H_0$ in the limit of $\tilde {\beta } \ll 1$ . Typical values of the grid size are $\triangle Z = 10^{-3}$ and $\triangle \eta = 0.005$ . Once $\tilde {A}_{11}$ and $\tilde {A}_{12}$ are determined, we use MATLAB’s routine trapz to calculate the pressure drop (4.7) in a contraction.

5. Results

In this section, we present our theoretical results for the pressure drop and elastic stresses of the FENE-CR fluid as developed in the previous sections. We also validate the predictions of our theoretical model against the two-dimensional numerical simulations with the finite-volume software OpenFOAM. The details of the numerical procedure are provided in Appendix D. For comparison and validation, in addition to the FENE-CR fluid, we show the results for the Oldroyd-B fluid.

As an illustrative example, we consider a hyperbolic contracting channel of the form

(5.1) \begin{equation} H(Z)=\frac {H_{0}}{(H_{0}-1)Z+1}, \end{equation}

where $H_{0}=h_0/h_\ell$ is the ratio of the heights at the inlet and outlet; for the contracting geometry, we have $H_0\gt 1$ . The present study focuses on the contraction ratio $H_{0}=h_0/h_\ell =4$ .

5.1. Pressure drop at low Deborah numbers

In this subsection, we elucidate the pressure drop behaviour of the FENE-CR fluid at low Deborah numbers using our analytical predictions and OpenFOAM simulation results. In addition, we present the pressure drop of the Oldroyd-B and FENE-P fluids, thus highlighting how the finite extensibility (without the influence of shear thinning) incorporated by the FENE-CR model impacts pressure drop.

For the planar hyperbolic contracting channel (5.1), using (3.4 $a$ ), (3.4 $b$ ), (3.4 $c$ ), (3.5) and (3.6), we obtain analytical expressions for the pressure drop contributions of the FENE-CR fluid up to $O(De^4)$ ,

(5.2a) \begin{equation} \triangle P_0=\frac {3}{4}\left(1+H^{-1}_{0}\right) \left(1+H^{-2}_{0}\right), \end{equation}

(5.2b) \begin{equation} \triangle P_1=-\frac {9}{2}\tilde \beta \left(1-H^{-4}_{0}\right), \end{equation}

(5.2c) \begin{equation} \triangle P_2=\frac {648}{35}\tilde \beta \left(1-H^{-1}_{0}\right)^2\left(1+H^{-1}_{0}\right)\left(1+H^{-2}_{0}\right),\end{equation}

(5.2d) \begin{equation} \triangle P_3=\frac {2673}{70L^2 \epsilon ^2}\tilde \beta \left(1-H^{-8}_{0}\right)-\frac {216}{35}\tilde \beta (11-\tilde \beta )\left(1-H^{-1}_{0}\right)^3\left(1+H^{-1}_{0}\right)\left(1+H^{-2}_{0}\right), \end{equation}

(5.2e) \begin{align} \triangle P_4=&-\frac {162}{35L^2 \epsilon ^2}\tilde \beta (32 \tilde \beta +139)\left(1-H^{-1}_{0}-H^{-8}_{0} +H^{-9}_{0}\right) \nonumber \\ &+\frac {324}{13475}\tilde \beta (840 \tilde \beta ^2-3351 \tilde \beta +9800)\left(1-H^{-1}_{0}\right)^4\left(1+H^{-1}_{0}\right)\left(1+H^{-2}_{0}\right). \end{align}

Using (5.2) in conjunction with (3.8), we obtain the Padé approximation for the pressure drop. Note that for $L^2\epsilon ^2 \to \infty$ , we recover the Oldroyd-B limit. In this case, the first terms in (5.2d) and (5.2e), which are dependent on $L^2\epsilon ^2$ , vanish.

Figure 2. Non-dimensional pressure drop at low Deborah numbers for the Oldroyd-B and FENE-CR fluids in a contracting channel described by (5.1). ( $a$ – $d$ ) Scaled pressure drop $\triangle P/\triangle P_0$ as a function of $De=\lambda q/(2\ell h_{\ell })$ for ( $a$ ) $L^2\epsilon ^2=10$ , ( $b$ ) $L^2\epsilon ^2=5$ , ( $c$ ) $L^2\epsilon ^2=0.5$ and ( $d$ ) $L^2\epsilon ^2=0.1$ . Grey triangles and purple circles respectively represent the OpenFOAM simulation results for the Oldroyd-B and FENE-CR fluids. Grey solid and green dashed lines represent the fourth-order asymptotic solutions for the Oldroyd-B and FENE-CR fluids, given by (5.2a) $-$ (5.2e). Cyan dotted and solid black lines represent the Padé approximation (3.8) applied to the fourth-order asymptotic solutions for the Oldroyd-B and FENE-CR fluids. All calculations were performed using $H_0=4$ and $\tilde \beta =0.4$ .

We present in figure 2 the scaled pressure drop $\triangle P/\triangle P_0$ as a function of $De=\lambda q/(2\ell h_{\ell })$ for the Oldroyd-B and FENE-CR fluids in a contracting channel for different values of $L^2\epsilon ^2$ . Grey triangles and purple circles represent the OpenFOAM simulation results for the Oldroyd-B and FENE-CR fluids obtained from calculating the pressure drop along the centreline. Grey solid and green dashed lines represent the fourth-order asymptotic solutions for the Oldroyd-B and FENE-CR fluids. Cyan dotted and black solid lines respectively represent the Padé approximation (3.8) applied to the fourth-order asymptotic solutions for the Oldroyd-B and FENE-CR fluids.

First, we observe that the fourth-order asymptotic solutions (grey solid and green dashed lines) cannot accurately capture the pressure drop except for very low values of $De$ , consistent with results of Housiadas & Beris (2023), indicating that the asymptotic series has a very small radius of convergence. Nevertheless, when using the Padé approximation to accelerate the convergence of the asymptotic series, we find that our analytical predictions for the pressure drop are in excellent agreement with numerical simulations for both Oldroyd-B and FENE-CR fluids. For example, even for $L^2\epsilon ^2=0.1$ , where the Padé approximation slightly overpredicts the pressure drop of the FENE-CR fluid, the relative error is approximately 5 % for up to $De = 0.5$ .

Second, it is evident that, at low Deborah numbers, the dimensionless pressure drop of both Oldroyd-B and FENE-CR fluids monotonically decreases with $De$ , similar to Giesekus and FENE-P fluids (Housiadas & Beris 2023). Furthermore, as expected, for $L^2\epsilon ^2=10$ and $L^2\epsilon ^2=5$ , the pressure drop behaviour of both fluids is almost indistinguishable. However, when the finite extensibility becomes more apparent, i.e. as $L^2\epsilon ^2$ decreases, the FENE-CR model predicts a higher dimensionless pressure drop than the Oldroyd-B model, as shown in figure 2( $d$ ).

It is of particular interest to compare and contrast our predictions for the pressure drop of the FENE-CR fluid with the recent low- $De$ results of Housiadas & Beris (2023) for the FENE-P fluid. Such a comparison of the non-dimensional pressure drop is shown in figure 3 for Oldroyd-B, FENE-CR and FENE-P fluids in a contracting channel for $L^2\epsilon ^2=0.5$ and $0.25$ . Blue dashed lines represent the Padé approximation (3.8) applied to the fourth-order asymptotic solutions obtained from Housiadas & Beris (2023) for the FENE-P fluid, when accounting for the differences in characteristic scales. Similar to the Oldroyd-B and FENE-CR fluids, the dimensionless pressure drop of the FENE-P fluid monotonically decreases with $De$ at low Deborah numbers. Furthermore, as expected, the FENE-P fluid shows a higher pressure drop than the Oldroyd-B fluid due to the effects of finite extensibility. However, due to the shear-thinning effects, the resulting pressure drop of the FENE-P fluid is lower than that of the FENE-CR fluid.

Figure 3. Comparison of non-dimensional pressure drop at low Deborah numbers for the Oldroyd-B, FENE-CR and FENE-P fluids in a contracting channel. ( $a{,}b$ ) Scaled pressure drop $\triangle P/\triangle P_0$ as a function of $De=\lambda q/(2\ell h_{\ell })$ for ( $a$ ) $L^2\epsilon ^2=0.5$ and ( $b$ ) $L^2\epsilon ^2=0.25$ . Grey triangles and purple circles represent the OpenFOAM simulation results for the Oldroyd-B and FENE-CR fluids, respectively. Cyan dotted, solid black and dashed blue lines represent the Padé approximation (3.8) applied to the fourth-order asymptotic solutions for the Oldroyd-B, FENE-CR and FENE-P fluids. All calculations were performed using $H_0=4$ and $\tilde \beta =0.4$ .

Although our low- $De$ analysis using the Padé approximation predicts well the pressure drop at low Deborah numbers, it cannot accurately capture the pressure drop behaviour at high Deborah numbers. To this end, in the next subsections, we employ numerical simulations and the low- $\tilde \beta$ lubrication analysis.

5.2. Pressure drop and elastic stresses at high Deborah numbers

In this subsection, we study and contrast the elastic stresses and pressure drop of the Oldroyd-B and FENE-CR fluids across the contraction at high Deborah numbers. Specifically, we first consider high Deborah numbers up to $De= 4$ using the OpenFOAM simulations and validate the predictions of our low- $\tilde \beta$ lubrication analysis. Then, we employ the low- $\tilde \beta$ lubrication analysis to study the behaviour of the elastic stresses and pressure drop at sufficiently high Deborah numbers up to $De=20$ .

Figure 4. Non-dimensional pressure drop at high Deborah numbers for the Oldroyd-B and FENE-CR fluids in a contracting channel. ( $a,b$ ) Scaled pressure drop $\triangle P/\triangle P_0$ as a function of $De=\lambda q/(2\ell h_{\ell })$ (or $De_{entry}=\lambda q/(2\ell h_{0})$ ) for ( $a$ ) $\tilde \beta =0.4$ and ( $b$ ) $\tilde \beta =0.05$ . Grey triangles and purple circles represent the OpenFOAM simulation results for the Oldroyd-B and FENE-CR fluids. Black dots and grey crosses in ( $b$ ) represent the results of the low- $\tilde \beta$ lubrication analysis for the Oldroyd-B and FENE-CR fluids. Cyan dotted and solid black lines represent the low- $De$ Padé approximation (3.8) for the Oldroyd-B and FENE-CR fluids. Red dashed lines represent the high- $De$ asymptotic solution (4.10) for the Oldroyd-B fluid. All calculations were performed using $H_0=4$ and $L^2\epsilon ^2=0.5$ .

First, in figure 4( $a,b$ ), we present the scaled pressure drop $\triangle P/\triangle P_0$ of the Oldroyd-B and FENE-CR fluids in the contraction as a function of $De=\lambda q/(2\ell h_{\ell })$ for $\tilde \beta =0.4$ (panel $a$ ) and $\tilde \beta =0.05$ (panel $b$ ) , with $L^2\epsilon ^2 = 0.5$ . Grey triangles and purple circles respectively represent the OpenFOAM simulation results for Oldroyd-B and FENE-CR fluids. Black dots and grey crosses respectively represent the results of the low- $\tilde \beta$ lubrication analysis for the Oldroyd-B and FENE-CR fluids. Cyan dotted and solid black lines represent the low- $De$ Padé approximation (3.8) for the Oldroyd-B and FENE-CR fluids. Red dashed lines represent the high- $De$ asymptotic solution (4.10) for the Oldroyd-B fluid in the ultra-dilute limit. As both the Deborah number $De=\lambda q/(2\ell h_{\ell })$ based on the exit height and the Deborah number $De_{entry}=\lambda q/(2\ell h_{0})$ based on the entry height are used in the literature, we present our results both as a function of $De$ and $De_{entry}$ .

Consistent with the previous studies (Boyko et al. 2024; Hinch et al. 2024), the pressure drop of the Oldroyd-B fluid monotonically decreases with $De$ and scales linearly with $De$ at high Deborah numbers for $\tilde \beta =0.05$ , corresponding to the ultra-dilute limit, as represented by the red dashed line in figure 4( $b$ ). Furthermore, there is excellent agreement between the predictions of the low- $\tilde \beta$ lubrication analysis with $\tilde \beta =0.05$ and the OpenFOAM simulations. In particular, for the Oldroyd-B fluid, the relative error at $De =2$ is 0.2 %. Nevertheless, as expected, for $\tilde \beta =0.4$ (figure 4 $a$ ), the high- $De$ asymptotic solution (4.10) for the Oldroyd-B fluid in the ultra-dilute limit does not accurately capture the slope of the OpenFOAM simulations due to the deviations in the flow velocity from the parabolic profile when .

In contrast to a monotonic pressure drop reduction with $De$ observed for the Oldroyd-B fluid, the pressure drop of the FENE-CR fluid levels off to a plateau at high Deborah numbers for both $\tilde \beta =0.4$ and 0.05, with a slight increase for $De\gtrsim 3$ . Understanding this non-monotonic pressure drop variation for the FENE-CR fluid necessitates analysing higher Deborah numbers. We note that the presented OpenFOAM simulations for the FENE-CR fluid are in the range of $0\leqslant De\leqslant 4$ . Performing simulations at higher Deborah numbers requires a longer downstream (exit) section to allow the elastic stresses to reach their fully relaxed values (see e.g. Debbaut et al. 1988; Keiller 1993; Alves et al. 2003; Boyko et al. 2024), thus significantly increasing the computational time. Furthermore, above a certain high $De$ , we expect our OpenFOAM simulations to suffer from accuracy and convergence difficulties associated with the high-Weissenberg-number problem (Owens & Phillips 2002; Alves et al. 2021). Indeed, for the Oldroyd-B fluid with $\tilde \beta =0.05$ , we cannot obtain reliable results above $De \approx 2$ .

Therefore, instead of carrying out computationally expensive simulations, we use the low- $\tilde \beta$ lubrication analysis considering the ultra-dilute limit, which is considerably faster and allows us to access the behaviour of the elastic stresses and pressure drop at arbitrary values of $De$ . Such an approach is strongly supported by the excellent agreement between the pressure drop predictions of the low- $\tilde \beta$ lubrication analysis and the OpenFOAM simulation results, as shown in figure 4( $b$ ). Specifically, for the FENE-CR fluid, we find a relative error of approximately 0.3 % for up to $De = 4$ .

Figure 5. Streamwise variation of elastic stresses of the FENE-CR fluid on $\eta =0.5$ in a contracting channel in the ultra-dilute limit. ( $a$ – $i$ ) Elastic normal and shear stresses $\mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})\tilde {A}_{11}$ and $\mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})\tilde {A}_{12}$ , scaled by their entry values, as a function of $Z$ for different values of $De$ and $L^2\epsilon ^2$ . Solid lines represent the results of the low- $\tilde \beta$ lubrication analysis. Cyan dotted lines in panel ( $a$ – $c$ ) represent the low- $De$ asymptotic solutions for the FENE-CR fluid. Red dashed lines in panel ( $g$ ) represent the high- $De$ asymptotic solutions for the Oldroyd-B fluid. All calculations were performed using $H_{0}=4$ .

Before investigating the pressure drop behaviour at higher Deborah numbers, it is of particular interest to elucidate the spatial variation of the elastic stresses. We present in figure 5 the streamwise variation of the elastic normal and shear stresses of the FENE-CR fluid, scaled by their entry values, on $\eta =0.5$ in a contracting channel in the ultra-dilute limit for different values of $De$ and $L^2\epsilon ^2$ . As expected, for $L^2\epsilon ^2=50$ , we recover the Oldroyd-B behaviour previously studied by Hinch et al. (2024) and Boyko et al. (2024). Specifically, we find that, at low Deborah numbers ( $De=0.05$ , figure 5 $a$ ), the elastic shear and normal stresses increase by a factor of $H_0^2=16$ and $H_0^4=256$ , respectively, by the end of contraction. In contrast, at high Deborah numbers ( $De=5$ , figure 5 $g$ ), the elastic shear stress $\mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})\tilde {A}_{12}$ preserves its entry value and the elastic normal stress $\mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})\tilde {A}_{11}$ increases by a factor of $H_0^2=16$ .

It is evident from figure 5( $a$ – $c$ ) that at $De=0.05$ , the elastic shear stress weakly depends on the finite extensibility parameter $L^2\epsilon ^2$ , where the magnitude of the elastic normal stress decreases as $L^2\epsilon ^2$ is reduced from 50 to 0.005. At higher Deborah numbers, $De=0.5$ and $De=5$ , we observe a trade-off between the axial component of the conformation tensor $\tilde {A}_{11}$ and the finite extensibility, incorporated by the nonlinear spring function $\mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})=(1-\tilde {A}_{11}/(L^2 \epsilon ^2))^{-1}$ . On the one hand, when $L^2\epsilon ^2$ is large (the Oldroyd-B limit), the dumbbell extension, as measured by $\mathrm {tr}\,\tilde {{\boldsymbol{\!A}}}\approx \tilde {A}_{11}$ , is large and $\mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})\approx 1$ . On the other hand, when $L^2\epsilon ^2$ is small, the dumbbell extension $\mathrm {tr}\,\tilde {{\boldsymbol{\!A}}}\approx \tilde {A}_{11}$ is also small but $\mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})$ can be large. Therefore, as shown in figures 5( $d$ – $f$ ) and 5( $g$ – $i$ ), for a sufficient large $De$ , the maximum value of elastic normal stress $\mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})\tilde {A}_{11}$ , achieved at the end of contraction, may exhibit a non-monotonic variation with $L^2\epsilon ^2$ (see also figure 8 $b$ ). For example, when $De=5$ , the maximum value of $\mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})\tilde {A}_{11}$ for $L^2\epsilon ^2=0.5$ is greater than the corresponding values for $L^2\epsilon ^2=50$ and $L^2\epsilon ^2=0.005$ . Furthermore, for $De=5$ , in contrast to the Oldroyd-B fluid where $\mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})\tilde {A}_{12}$ maintains its entry value, when the finite extensibility is significant, i.e. $L^2\epsilon ^2=0.5$ and $L^2\epsilon ^2=0.005$ , we observe a non-monotonic increase of elastic shear stress with axial position $Z$ .

Figure 6. ( $a$ ) Scaled pressure drop $\triangle P/\triangle P_0$ as a function of $De=\lambda q/(2\ell h_{\ell })$ for $\tilde \beta =0.05$ . Black dots and grey crosses represent the results of the low- $\tilde \beta$ lubrication analysis for the Oldroyd-B and FENE-CR fluids. Cyan dotted and solid black lines represent the low- $De$ Padé approximation (3.8) for the Oldroyd-B and FENE-CR fluids. The red dashed line represents the high- $De$ asymptotic solution (4.10) for the Oldroyd-B fluid. ( $b$ ) Elastic contributions to the non-dimensional pressure drop, scaled by $\tilde \beta$ , as a function of $De=\lambda q/(2\ell h_{\ell })$ in the ultra-dilute limit. Black circles and grey dots represent ultra-dilute predictions of the Oldroyd-B fluid for elastic shear and normal stress contributions. Black crosses and purple squares represent ultra-dilute predictions of the FENE-CR fluid for elastic shear and normal stress contributions. Red and black dashed lines represent the high- $De$ asymptotic solution of the Oldroyd-B fluid for elastic shear and normal stress contributions. All calculations were performed using $H_0=4$ and $L^2\epsilon ^2=0.5$ .

Next, we analyse the pressure drop variation at significantly higher Deborah numbers using our low- $\tilde \beta$ lubrication analysis. We present in figure 6( $a$ ) the scaled pressure drop $\triangle P/\triangle P_0$ of the Oldroyd-B and FENE-CR fluids in the contraction as a function of $De=\lambda q/(2\ell h_{\ell })$ for $L^2\epsilon ^2 = 0.5$ and $\tilde \beta =0.05$ , corresponding to the ultra-dilute limit. Black dots represent the results of the low- $\tilde \beta$ lubrication analysis, the cyan dotted line represents the low- $De$ Padé approximation (3.8) and the red dashed line represents the high- $De$ asymptotic solution (4.10) for the Oldroyd-B fluid. We observe excellent agreement between our low- and high- $De$ asymptotic solutions and the low- $\tilde \beta$ lubrication results. Moreover, somewhat surprisingly, from figure 6( $a$ ) it follows that the low- $De$ Padé approximation (3.8) captures fairly well the pressure drop reduction with $De$ for up to $De=2$ ( $De_ {entry}=0.5$ ) for both Oldroyd-B and FENE-CR fluids. More importantly, unlike a linear pressure drop reduction of the Oldroyd-B fluid at high Deborah numbers, the pressure drop of the FENE-CR fluid (grey crosses) exhibits a non-monotonic variation, first decreasing with $De$ , attaining a local minimum at $De \approx 2.8$ and then increasing with $De$ . Nevertheless, the non-dimensional pressure drop for the FENE-CR fluid in the ultra-dilute limit is lower than the corresponding Newtonian pressure drop, i.e. $\triangle P/\triangle P_0\lt 1$ , even for very high Deborah numbers. Such a non-monotonic variation in the pressure drop of the FENE-CR fluid is consistent with the previous numerical studies on the flow of the FENE-P fluid in 2-D abruptly contracting geometries (Zografos et al. 2022). It is well known that the FENE-P fluid incorporates both features of finite extensibility and shear thinning. However, we believe that the main source of the non-monotonic pressure drop behaviour is associated with the finite extensibility since the shear-thinning effect leads to the pressure drop reduction, as shown in figure 3, and thus is unlikely to reverse the trend, causing the pressure drop increase.

To probe deeper into the source of the non-monotonic variation of the pressure drop for the FENE-CR fluid, we present in figure 6( $b$ ) the elastic contributions to the non-dimensional pressure drop, scaled by $\tilde \beta$ , as a function of $De=\lambda q/(2\ell h_{\ell })$ in the ultra-dilute limit. Black circles and grey dots represent the elastic shear and normal stress contributions obtained from the low- $\tilde \beta$ lubrication analysis for the Oldroyd-B fluid. Black crosses and purple squares represent the elastic shear and normal stress contributions obtained from the low- $\tilde \beta$ lubrication analysis for the FENE-CR fluid. As expected, for the Oldroyd-B fluid, there is excellent agreement between our low- $\tilde \beta$ lubrication results and the high- $De$ asymptotic solution (4.10), represented by red and black dashed lines.

In contrast to the Oldroyd-B fluid, where the elastic normal stress contribution decreases with $De$ and scales linearly with $De$ at high Deborah numbers, for the FENE-CR fluid, we observe a non-monotonic variation. In particular, the elastic normal stress contribution of the FENE-CR fluid first decreases, attains a minimum at $De \approx 3.2$ and then increases with $De$ . Such an increase is associated with the dissipative effect of the finite extensibility. Despite this increase, figure 6( $b$ ) clearly shows that the elastic normal stress contribution of the FENE-CR fluid is negative at $De=20$ , leading to a reduction in the pressure drop, similar to the Oldroyd-B fluid. However, we find that, at $De \approx 118$ , the elastic normal stress contribution of the FENE-CR fluid becomes positive and then increases with $De$ . For $\tilde \beta =0.05$ , we have confirmed that up to $De=1000$ , this positive elastic normal stress contribution is too weak since it scales with $\tilde \beta$ , and thus cannot lead to the pressure drop enhancement above the Newtonian value $\triangle P_0$ . Note that we have assumed steady flows, so further investigation is necessary to determine if there might be flow instabilities at these high Deborah numbers. Nevertheless, as pointed out by Hinch et al. (2024), under the lubrication approximation, the hoop stress is neglected, so purely elastic instability cannot arise due to curved streamlines.

Figure 7. Influence of the finite extensibility on the non-dimensional pressure drop of the FENE-CR fluid in a contracting channel. ( $a{,}b$ ) Scaled pressure drop $\triangle P/\triangle P_0$ as a function of the finite extensibility parameter $L^2\epsilon ^2$ for ( $a$ ) low and ( $b$ ) high Deborah numbers. Triangles in ( $a$ ) represent the OpenFOAM simulation results. Dots represent the results obtained from the low- $\tilde {\beta }$ lubrication analysis. Dash-dotted lines represent the low- $De$ Padé approximation (3.8) applied up to the fourth-order asymptotic solution. Cyan dotted lines represent the low- $L^2\epsilon ^2$ asymptotic solution, corresponding to the Newtonian limit. Red dashed lines represent the high- $L^2\epsilon ^2$ asymptotic solution, corresponding to the Oldroyd-B limit. All calculations were performed using $H_0=4$ and $\tilde \beta =0.05$ .

The elastic shear stress contribution of the FENE-CR fluid also exhibits a non-monotonic variation with the Deborah number. It first decreases, attains a minimum at $De \approx 1.2$ and then approaches a plateau at high Deborah numbers. Such a non-monotonic variation of the elastic normal and shear stress contributions rationalises the non-monotonic pressure drop behaviour, shown in figure 6( $a$ ). Similar to the Oldroyd-B fluid, the elastic shear stress contribution of the FENE-CR fluid is independent of $De$ at high Deborah numbers, but with a constant value higher than for the Oldroyd-B fluid, due to the dissipative effect of the finite extensibility. This higher value of elastic shear stress contribution leads to an even greater increase in the pressure drop of the FENE-CR fluid compared with the Oldroyd-B fluid.

5.3. Assessing the effect of the finite extensibility on the pressure drop

In the previous subsections, we analysed the pressure drop variation with the Deborah number $De$ and the viscosity ratio $\tilde {\beta }$ , mainly considering the finite extensibility parameter $L^2\epsilon ^2=0.5$ . In this subsection, we study how the finite extensibility parameter $L^2\epsilon ^2$ impacts the pressure drop.

First, in figure 7( $a,b$ ), we present the variation of the scaled pressure drop $\triangle P/\triangle P_0$ as a function of $L^2\epsilon ^2$ for the FENE-CR fluid in a contracting channel for low and high Deborah numbers, with $\tilde \beta =0.05$ . Triangles and dots respectively represent the results of the OpenFOAM simulations and low- $\tilde {\beta }$ lubrication analysis. Dash-dotted lines represent the low- $De$ Padé approximation (3.8) applied up to the fourth-order asymptotic solution. Cyan dotted and red dashed lines represent the low- and high- $L^2\epsilon ^2$ asymptotic solutions, corresponding to the Newtonian and Oldroyd-B limits, respectively.

At low Deborah numbers, it is evident from figure 7( $a$ ) that the pressure drop monotonically decreases with increasing $L^2\epsilon ^2$ . Clearly, there is excellent agreement between our low- $De$ asymptotic solutions based on the Padé approximation, the OpenFOAM simulation results and the predictions of the low- $\tilde {\beta }$ lubrication analysis. Consistent with the low- $De$ Padé approximation (3.8), for small values of $L^2\epsilon ^2$ , the pressure drop becomes independent of $De$ , approaching the Newtonian limit for all values of $De$ , represented by cyan dotted lines. As expected, for large values of $L^2\epsilon ^2$ , the pressure drop approaches the Oldroyd-B limit, represented by red dashed lines.

Next, we consider the variation in pressure drop with $L^2\epsilon ^2$ at high Deborah numbers, as shown in figure 7( $b$ ). At high Deborah numbers, the pressure drop shows Newtonian and Oldroyd-B asymptotic behaviour for $L\epsilon \ll 1$ and $L\epsilon \gg 1$ , similar to the low- $De$ limit. However, in contrast to low Deborah numbers, at high Deborah numbers $De=2$ and 3, pressure drop exhibits a strong non-monotonic behaviour with $L^2\epsilon ^2$ . Specifically, we observe that the pressure drop first decreases and then increases with $L^2\epsilon ^2$ approaching the Oldroyd-B limit, with the transition occurring at $L^2\epsilon ^2 =O(1)$ .

Figure 8. ( $a$ ) Elastic contributions to the non-dimensional pressure drop of the FENE-CR fluid, scaled by $\tilde \beta$ , as a function of the finite extensibility parameter $L^2\epsilon ^2$ for $De=3$ in the ultra-dilute limit. Black crosses and purple squares represent the elastic shear and normal stress contributions obtained from the low- $\tilde {\beta }$ lubrication analysis. Cyan and grey dotted lines represent the low- $L^2\epsilon ^2$ asymptotic solution for the elastic shear and normal stress contributions, corresponding to the Newtonian limit. Red and black dashed lines represent the high- $L^2\epsilon ^2$ asymptotic solution (4.10) for the elastic shear and normal stress contributions, corresponding to the Oldroyd-B limit at high $De$ . ( $b$ ) Elastic normal stress $\mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})\tilde {A}_{11}(Z,\eta =0.7)$ as a function of $Z$ for $De=3$ and $L^2\epsilon ^2=0.45$ (dotted line), $L^2\epsilon ^2=2.8$ (dashed line) and $L^2\epsilon ^2=100$ (solid line). All calculations were performed using $H_0=4$ .

To provide further insight into the pressure drop dependence on the finite extensibility $L^2\epsilon ^2$ for a given $De$ , we study the relative importance of elastic contributions to the pressure drop. The elastic contributions to the non-dimensional pressure drop across the contraction, scaled by $\tilde {\beta }$ , as a function of $L^2\epsilon ^2$ are shown in figure 8( $a$ ) for $De=3$ . Black crosses and purple squares represent the elastic shear and normal stress contributions obtained from the low- $\tilde {\beta }$ lubrication analysis. Red and black dashed lines represent the high- $L^2\epsilon ^2$ asymptotic solution (4.10) for the elastic shear and normal stress contributions, corresponding to the Oldroyd-B limit at high $De$ .

For small values of $L^2\epsilon ^2$ , the elastic normal stress contribution to the pressure drop approaches zero, while the elastic shear stress contribution approaches an order-one Newtonian value. We rationalise this behaviour by noting from (4.4) and (4.6a), (4.6b) that at the beginning of the contraction, in the low- $L^2\epsilon ^2$ limit, we have

(5.3) \begin{equation} \tilde {A}_{11}=L^2\epsilon ^2 -\frac {H_0^2 L^3\epsilon ^3}{3 \sqrt {2} De \eta }+O(L^4\epsilon ^4), \!\quad\! \tilde {A}_{12}=-\frac {L\epsilon }{\sqrt {2}}+\frac {H_0^2 L^2\epsilon ^2 }{12 De\eta }+O(L^3\epsilon ^3) \!\!\quad \text {for} \; L\epsilon \ll 1, \end{equation}

(5.4) \begin{equation} \mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})\tilde {A}_{11}=\frac {3 \sqrt {2} L\epsilon De }{H_0^2}\eta +O(L^2\epsilon ^2), \quad \mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})\tilde {A}_{12}=-\frac {3 De }{H_0^2}\eta \quad \text {for} \; L\epsilon \ll 1. \end{equation}

This result is valid for all $De$ . Therefore, for $L\epsilon \ll 1$ , the elastic normal stress $\mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})\tilde {A}_{11}$ scales as $O(L\epsilon De )$ and the elastic shear stress $\mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})\tilde {A}_{12}$ scales as $O(De)$ . Using (4.7), the latter scaling arguments imply that the elastic normal stress contribution to the pressure drop scales as $O(L\epsilon )$ and thus is negligible when $De=O(1)$ . However, the elastic shear stress has a Newtonian contribution, which is independent of $De$ , as shown in figure 8( $a$ ).

Furthermore, we observe that, while the elastic shear stress contribution monotonically decreases with increasing $L^2\epsilon ^2$ , the elastic normal stress contribution exhibits a non-monotonic variation with $L^2\epsilon ^2$ . Thus, the non-monotonic behaviour of the pressure drop, shown in figure 7( $b$ ) for $De=2$ and 3 at $L^2\epsilon ^2 =O(1)$ , arises due to the elastic normal stress contribution. Such a non-monotonic variation with $L^2\epsilon ^2$ for a given $De$ can be attributed to the trade-off between the axial component of the conformation tensor $\tilde {A}_{11}$ and the finite extensibility $L^2 \epsilon ^2$ through $\mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})=(1-\tilde {A}_{11}/(L^2 \epsilon ^2))^{-1}$ , as discussed in $\S$ 5.2. For example, as shown in figure 8( $b$ ), for a given $De$ , the elastic normal stress $\mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})\tilde {A}_{11}$ can exhibit similar spatial variations for small (dotted line) and large (solid line) values of $L^2\epsilon ^2$ , rationalising the non-monotonic behaviour of elastic normal stress contribution to the pressure drop.

6. Concluding remarks

In this work, we studied the flow of a FENE-CR fluid in slowly varying contracting channels at low and high Deborah numbers. Employing the low-Deborah-number lubrication analysis, we provided analytical expressions for the non-dimensional pressure drop for the FENE-CR fluid up to $O(De^4)$ and applied the Padé approximation to improve the convergence of the asymptotic series. To understand the pressure drop behaviour of the FENE-CR fluid at high Deborah numbers, we considered the ultra-dilute limit of small polymer concentration and exploited the one-way coupling between the parabolic velocity and elastic stresses to calculate the pressure drop for arbitrary values of $De$ . We further compared and contrasted the predictions of the FENE-CR model to the recent results of Boyko et al. (2024) and Hinch et al. (2024) for the Oldroyd-B model as well as to the low- $De$ results of Housiadas & Beris (2023) for the FENE-P model. We validated our theoretical results for the dimensionless pressure drop in a contracting channel with 2-D finite-volume numerical simulations for both Oldroyd-B and FENE-CR fluid and found excellent agreement.

At low Deborah numbers, the pressure drop of the FENE-CR fluid monotonically decreases with $De$ , as shown in figure 2, similar to the predictions of the Oldroyd-B and FENE-P fluids. However, at high Deborah numbers, unlike a linear pressure drop reduction of the Oldroyd-B fluid, the pressure drop of the FENE-CR fluid exhibits a non-monotonic variation, first decreasing and then increasing with $De$ . Note that the pressure drop for the FENE-CR fluid remains lower than the corresponding Newtonian pressure drop even for very high Deborah numbers, as shown in figure 6( $a$ ). We identified two causes for such pressure drop variation of the FENE-CR fluid (see figure 6 $b$ ). The first cause is the elastic normal stress contribution to the pressure drop, which becomes less negative as $De$ increases at high Deborah numbers due to the dissipative effect of the finite extensibility. The second cause is the contribution of elastic shear stresses, which is higher compared with the Oldroyd-B fluid, again owing to the dissipative effect of the finite extensibility.

In general, the pressure drop of the FENE-CR fluid increases compared with the Oldroyd-B fluid as the finite extensibility becomes more apparent (when $L^2\epsilon ^2$ decreases). Nevertheless, for very small values of $L^2\epsilon ^2$ , the pressure drop of FENE-CR fluid becomes independent of $De$ and approaches the Newtonian value. Specifically, when $L^2\epsilon ^2 \ll 1$ , the elastic normal stress contribution vanishes while the elastic shear stress contribution shows a Newtonian behaviour for all $De$ (see figure 7 and figure 8 $a$ ).

Our theoretical framework, based on lubrication theory and the ultra-dilute limit, allows us to study the behaviour of the elastic stresses and pressure drop of a FENE-CR fluid at sufficiently high Deborah numbers. We emphasise that we are currently unable to achieve these high values of the Deborah number using finite-volume or finite-element simulations. We, therefore, believe that our theoretical results for the FENE-CR fluid in the ultra-dilute limit, valid at all $De$ , are of fundamental interest and can be helpful for simulation validation and enhancing our understanding of viscoelastic channel flows.

The theoretical predictions of the non-monotonic pressure drop behaviour of the FENE-CR fluid in a contraction are consistent with the previous numerical studies on contraction geometries (see e.g. Nyström et al. 2012; Zografos et al. 2022). However, these predictions are in contrast with the experimental results showing a nonlinear increase in the pressure drop with $De$ above the Newtonian pressure drop value for the flow of a Boger fluid through abrupt axisymmetric contraction and contraction–expansion geometries (Rothstein & McKinley 1999; Nigen & Walters 2002; Nigen & Walters 2002; Sousa et al. 2009). Our results with the FENE-CR model that incorporates the feature of finite extensibility cannot resolve this contradiction. Thus, as a future research direction, it is interesting to study more complex elastic dumbbell models that account for additional microscopic features of realistic polymer chains, such as the conformation-dependent friction coefficient and the conformation-dependent non-affine deformation (Phan-Thien et al. 1984; Boyko & Stone 2024), and to elucidate their effect on the pressure drop.

Finally, we note that, in this work, we have focused on studying the pressure drop across the contraction region. However, numerical simulations and experimental set-ups include a long downstream (exit) section to allow the stresses to reach their fully relaxed values (Keiller 1993; Rothstein & McKinley 2001; Boyko et al. 2024). Therefore, one interesting extension of the present work is to study the spatial relaxation of elastic stresses, velocity and pressure of viscoelastic fluids in the exit channel using the FENE-CR model and more complex constitutive equations.