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. Reference Bird, Armstrong and Hassager1987; Steinberg Reference Steinberg2021; Datta et al. Reference Datta2022; Ewoldt & Saengow Reference Ewoldt and Saengow2022).
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. Reference Szabo, Rallison and Hinch1997; Alves et al. Reference Alves, Oliveira and Pinho2003; Binding et al. Reference Binding, Phillips and Phillips2006; Alves & Poole Reference Alves and Poole2007; Zografos et al. Reference Zografos, Hartt, Hamersky, Oliveira, Alves and Poole2020; Varchanis et al. Reference Varchanis, Tsamopoulos, Shen and Haward2022) and experimental measurements (Rothstein & McKinley Reference Rothstein and McKinley1999, Reference Rothstein and McKinley2001; Sousa et al. Reference Sousa, Coelho, Oliveira and Alves2009; Ober et al. Reference Ober, Haward, Pipe, Soulages and McKinley2013; James & Roos Reference James and Roos2021), and recently, via theoretical analysis (Pérez-Salas et al. Reference Pérez-Salas, Sánchez, Ascanio and Aguayo2019; Boyko & Stone Reference Boyko and Stone2022; Housiadas & Beris Reference Housiadas and Beris2023, Reference Housiadas and Beris2024c; Boyko et al. Reference Boyko, Hinch and Stone2024; Hinch et al. Reference Hinch, Boyko and Stone2024). For an overview of recent studies, the reader is referred to Boyko & Stone (Reference Boyko and Stone2022) and Hinch et al. (Reference Hinch, Boyko and Stone2024).
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 Reference Rothstein and McKinley1999; Alves et al. Reference Alves, Oliveira and Pinho2003; Binding et al. Reference Binding, Phillips and Phillips2006; Campo-Deaño et al. Reference Campo-Deaño, Galindo-Rosales, Pinho, Alves and Oliveira2011; Keshavarz & McKinley Reference Keshavarz and McKinley2016; Zografos et al. Reference Zografos, Afonso and Poole2022). 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 Reference Boyko and Stone2022; Housiadas & Beris Reference Housiadas and Beris2023,Reference Housiadas and Beris2024a,Reference Housiadas and Beris b,Reference Housiadas and Beris c,Reference Housiadas and Beris d). There have been numerous applications of lubrication theory to other viscoelastic fluid flows, such as thin films and tribology problems (Ro & Homsy Reference Ro and Homsy1995; Tichy Reference Tichy1996; Sawyer & Tichy Reference Sawyer and Tichy1998; Zhang et al. Reference Zhang, Matar and Craster2002; Saprykin et al. Reference Saprykin, Koopmans and Kalliadasis2007; Ahmed & Biancofiore Reference Ahmed and Biancofiore2021; Gamaniel et al. Reference Gamaniel, Dini and Biancofiore2021; Datt et al. Reference Datt, Kansal and Snoeijer2022; Ahmed & Biancofiore Reference Ahmed and Biancofiore2023, Reference Ahmed and Biancofiore2024), as well as translation of a sphere near a rigid plane (Ardekani et al. Reference Ardekani, Rangel and Joseph2007; Ruangkriengsin et al. Reference Ruangkriengsin, Brandão, Mahapatra, Boyko and Stone2024), and analysis of forces and torques acting on nearly touching spheres (Dandekar & Ardekani Reference Dandekar and Ardekani2021). For example, Sari et al. (Reference Sari, Putignano, Carbone and Biancofiore2024) 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 (Reference Boyko and Stone2022) 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 (Reference Housiadas and Beris2023) extended the analysis of Boyko & Stone (Reference Boyko and Stone2022) 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 Reference Phan-Thien and Tanner1977; Phan-Thien Reference Phan-Thien1978), Giesekus (Giesekus Reference Giesekus1982) and a finitely extensible nonlinear elastic (FENE) model with the Peterlin approximation (FENE-P) (Bird et al. Reference Bird, Dotson Paul and Johnson1980, Reference Bird, Armstrong and Hassager1987). 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 Reference Housiadas and Beris2023).
Recently, Hinch et al. (Reference Hinch, Boyko and Stone2024) and Boyko et al. (Reference Boyko, Hinch and Stone2024) 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. (Reference Hinch, Boyko and Stone2024) 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. (Reference Boyko, Hinch and Stone2024) 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. Reference Remmelgas, Singh and Leal1999; Moore & Shelley Reference Moore and Shelley2012; Li et al. Reference Li, Thomases and Guy2019; Mokhtari et al. Reference Mokhtari, Latché, Quintard and Davit2022). 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. (Reference Hinch, Boyko and Stone2024) and Boyko et al. (Reference Boyko, Hinch and Stone2024) 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 Reference Beris2021; Hinch & Harlen Reference Hinch and Harlen2021; Shaqfeh & Khomami Reference Shaqfeh and Khomami2021; Castillo-Sánchez et al. Reference Castillo-Sánchez, Jovanović, Kumar, Morozov, Shankar, Subramanian and Wilson2022; Stone et al. Reference Stone, Shelley and Boyko2023). 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. Reference Bird, Armstrong and Hassager1987). 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. Reference Alves, Oliveira and Pinho2003; Boyko & Stone Reference Boyko and Stone2022; Housiadas & Beris Reference Housiadas and Beris2023; Boyko et al. Reference Boyko, Hinch and Stone2024) 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 Reference Rothstein and McKinley1999; Nigen & Walters Reference Nigen and Walters2002; Nigen & Walters Reference Nigen and Walters2002; Sousa et al. Reference Sousa, Coelho, Oliveira and Alves2009). As pointed out by Hinch et al. (Reference Hinch, Boyko and Stone2024), 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 (Reference Chilcott and Rallison1988) 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 Reference James2009). In contrast, the FENE-P model incorporates both the finite extensibility and the shear-thinning effect of viscoelastic fluids. Recently, Ahmed & Biancofiore (Reference Ahmed and Biancofiore2023) 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. Reference Cruz, Pinho and Oliveira2005).
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 Reference Boyko and Stone2024). 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 Reference Housiadas and Beris2023). 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. Reference Zografos, Afonso and Poole2022). 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 (Reference Housiadas and Beris2023), 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 Reference Housiadas2017) 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 Reference Owens and Phillips2002; Alves et al. Reference Alves, Oliveira and Pinho2021). 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.
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
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. Reference Szabo, Rallison and Hinch1997; Rothstein & McKinley Reference Rothstein and McKinley1999; Alves et al. Reference Alves, Oliveira and Pinho2003; Alves & Poole Reference Alves and Poole2007; Campo-Deaño et al. Reference Campo-Deaño, Galindo-Rosales, Pinho, Alves and Oliveira2011; Ober et al. Reference Ober, Haward, Pipe, Soulages and McKinley2013; Zografos et al. Reference Zografos, Hartt, Hamersky, Oliveira, Alves and Poole2020; Boyko & Stone Reference Boyko and Stone2022; Boyko et al. Reference Boyko, Hinch and Stone2024; Hinch et al. Reference Hinch, Boyko and Stone2024), 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
\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
\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 (Reference Chilcott and Rallison1988). In contrast to the Oldroyd-B constitutive equation (Oldroyd Reference Oldroyd1950), 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. Reference Bird, Armstrong and Hassager1987). 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 Reference Chilcott and Rallison1988; Alves et al. Reference Alves, Oliveira and Pinho2021)
\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 Reference Warner1972),
\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 Reference Chilcott and Rallison1988)
\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 Reference Tichy1996; Zhang et al. Reference Zhang, Matar and Craster2002; Saprykin et al. Reference Saprykin, Koopmans and Kalliadasis2007; Ahmed & Biancofiore Reference Ahmed and Biancofiore2021; Boyko & Stone Reference Boyko and Stone2022; Boyko & Stone Reference Boyko and Stone2022; Boyko et al. Reference Boyko, Hinch and Stone2024),
\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}
\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}
\begin{equation} {\tilde A}_{zz}=\epsilon ^2 A_{zz}, \qquad {\tilde A}_{yz}=\epsilon A_{yz}, \qquad {\tilde A}_{yy}=A_{yy}, \end{equation}
\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,
\begin{equation} \epsilon =\frac {h_{\ell }}{\ell }\ll 1, \end{equation}
the contraction ratio,
\begin{equation} H_{0}=\frac {h_{0}}{h_{\ell }}, \end{equation}
the viscosity ratios,
\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,
\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,
\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 Reference Tichy1996; Zhang et al. Reference Zhang, Matar and Craster2002; Saprykin et al. Reference Saprykin, Koopmans and Kalliadasis2007; Ahmed & Biancofiore Reference Ahmed and Biancofiore2021; Boyko & Stone Reference Boyko and Stone2022; Ahmed & Biancofiore Reference Ahmed and Biancofiore2023; Housiadas & Beris Reference Housiadas and Beris2023; Boyko et al. Reference Boyko, Hinch and Stone2024; Hinch et al. Reference Hinch, Boyko and Stone2024). 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
\begin{align} \frac {\partial U_z}{\partial Z}+\frac {\partial U_y}{\partial Y}=0,\qquad\qquad\qquad\qquad\qquad \end{align}
\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}
\begin{align} \frac {\partial P}{\partial Y}=0, \qquad\qquad\qquad\qquad\qquad\qquad\end{align}
\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}
\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}
\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
\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 Reference Ahmed and Biancofiore2023; Housiadas & Beris Reference Housiadas and Beris2023). Although we consider
$\epsilon \ll 1$
, since the realistic values of
$L^2$
are typically large (see e.g. Remmelgas et al. (Reference Remmelgas, Singh and Leal1999); Rothstein & McKinley (Reference Rothstein and McKinley1999)), we may have
$L^2 \epsilon ^2=O(1)$
.
The corresponding boundary conditions on the velocity are
\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)
\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}
\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
\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 Reference Boyko and Stone2022)
\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 Reference Boyko and Stone2021, Reference Boyko and Stone2022). 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$
,
\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 (Reference Boyko and Stone2022), 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 Reference Black and Denn1976; Boyko & Stone Reference Boyko and Stone2022; Housiadas & Beris Reference Housiadas and Beris2023). For example, considering the leading and first order in
$De$
, we find
\begin{equation} \tilde A_{zz,0}=0,\qquad \tilde A_{yz,0}=0,\qquad \tilde A_{yy,0}=1, \end{equation}
\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
\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}
\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 (Reference Housiadas and Beris2023).
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
\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 Reference Housiadas and Beris2023).
Finally, at the fourth order in
$De$
, the resulting expression for the pressure drop of the FENE-CR fluid is
\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.
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 (Reference Housiadas and Beris2023) 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 (Reference Housiadas and Beris2023). 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)$
,
\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 Reference Hinch1991; Housiadas Reference Housiadas2017; Housiadas & Beris Reference Housiadas and Beris2023),
It should be noted that Housiadas & Beris (Reference Housiadas and Beris2023) 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. (Reference Boyko, Hinch and Stone2024) and Hinch et al. (Reference Hinch, Boyko and Stone2024), 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. Reference Remmelgas, Singh and Leal1999; Moore & Shelley Reference Moore and Shelley2012; Li et al. Reference Li, Thomases and Guy2019; Mokhtari et al. Reference Mokhtari, Latché, Quintard and Davit2022; Boyko et al. Reference Boyko, Hinch and Stone2024; Hinch et al. Reference Hinch, Boyko and Stone2024), 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. Reference Boyko, Hinch and Stone2024; Hinch et al. Reference Hinch, Boyko and Stone2024)
\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
\begin{align} U_{z}=U-\epsilon ^{2}\eta H'(\xi ) V, \qquad U_{y}=\eta H'(\xi ) U+V, \,\,\end{align}
\begin{align} \tilde {A}_{zz}=\tilde {A}_{11}+O(\epsilon ^{2}),\qquad\qquad\qquad\quad \end{align}
\begin{align} \tilde {A}_{zy}=\tilde {A}_{12}+\eta H'(\xi )\tilde {A}_{11}+O(\epsilon ^{2}), \qquad\qquad\end{align}
\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. Reference Boyko, Hinch and Stone2024).
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. Reference Boyko, Hinch and Stone2024; Hinch et al. Reference Hinch, Boyko and Stone2024) take the form
\begin{align} \frac {\partial (HU)}{\partial Z}+\frac {\partial V}{\partial \eta }=0,\qquad\qquad\qquad\qquad\qquad \end{align}
\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}
\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}
\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}
\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
\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
\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
\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}
\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
\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
\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. (Reference Boyko, Hinch and Stone2024), 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 }$
,
\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}
\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}
\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. Reference Boyko, Hinch and Stone2024). Furthermore, Boyko et al. (Reference Boyko, Hinch and Stone2024) and Hinch et al. (Reference Hinch, Boyko and Stone2024) 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
\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
\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)$
,
\begin{equation} \triangle P_0=\frac {3}{4}\left(1+H^{-1}_{0}\right) \left(1+H^{-2}_{0}\right), \end{equation}
\begin{equation} \triangle P_1=-\frac {9}{2}\tilde \beta \left(1-H^{-4}_{0}\right), \end{equation}
\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}
\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}
\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.
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 (Reference Housiadas and Beris2023), 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 Reference Housiadas and Beris2023). 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 (Reference Housiadas and Beris2023) 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 (Reference Housiadas and Beris2023) 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.
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$
.
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. Reference Boyko, Hinch and Stone2024; Hinch et al. Reference Hinch, Boyko and Stone2024), 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. Reference Debbaut, Marchal and Crochet1988; Keiller Reference Keiller1993; Alves et al. Reference Alves, Oliveira and Pinho2003; Boyko et al. Reference Boyko, Hinch and Stone2024), 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 Reference Owens and Phillips2002; Alves et al. Reference Alves, Oliveira and Pinho2021). 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$
.
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. (Reference Hinch, Boyko and Stone2024) and Boyko et al. (Reference Boyko, Hinch and Stone2024). 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$
.
(
$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. Reference Zografos, Afonso and Poole2022). 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. (Reference Hinch, Boyko and Stone2024), under the lubrication approximation, the hoop stress is neglected, so purely elastic instability cannot arise due to curved streamlines.
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)$
.
(
$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
\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}
\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. (Reference Boyko, Hinch and Stone2024) and Hinch et al. (Reference Hinch, Boyko and Stone2024) for the Oldroyd-B model as well as to the low-
$De$
results of Housiadas & Beris (Reference Housiadas and Beris2023) 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. Reference Nyström, Tamaddon-Jahromi, Stading and Webster2012; Zografos et al. Reference Zografos, Afonso and Poole2022). 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 Reference Rothstein and McKinley1999; Nigen & Walters Reference Nigen and Walters2002; Nigen & Walters Reference Nigen and Walters2002; Sousa et al. Reference Sousa, Coelho, Oliveira and Alves2009). 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. Reference Phan-Thien, Manero and Leal1984; Boyko & Stone Reference Boyko and Stone2024), 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 Reference Keiller1993; Rothstein & McKinley Reference Rothstein and McKinley2001; Boyko et al. Reference Boyko, Hinch and Stone2024). 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.
Supplementary material.
Supplementary material is available at http://doi.org/10.1017/jfm.2025.142. Supplementary material includes the MATHEMATICA file containing the explicit expressions for the velocity, conformation tensor components and the pressure drop in the low-Deborah-number limit up to
$O(De^4)$
.
Funding.
B.M. and E.B. acknowledge the support by grant no. 2022688 from the US-Israel Binational Science Foundation (BSF) and grant no. 1942/23 from the Israel Science Foundation (ISF). H.A.S. acknowledges the support from grant no. CBET-2246791 from the United States National Science Foundation (NSF). E.B. acknowledges the support from the Israeli Council for Higher Education Yigal Alon Fellowship.
Declaration of interests.
The authors report no conflict of interest.
Appendix A. A fully developed flow of a FENE-CR fluid in a straight channel
Consider a steady and fully developed flow of a FENE-CR fluid in a straight and long channel of non-dimensional height
$2H_0$
. Under the assumption of a fully developed flow, we have
$ U_y\equiv 0$
, so that the governing equations (2.12) simplify to
\begin{align} \frac {\mathrm {d} P}{\mathrm {d} Z}=(1-\tilde \beta )\frac {\mathrm {d}^2 U_z}{\mathrm {d} Y^2}+\frac {\tilde \beta }{De} \frac {\mathrm {d}(\mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})\tilde A_{yz})}{\mathrm {d} Y}, \end{align}
\begin{align} 2\frac {\mathrm {d} U_z}{\mathrm {d} Y}\tilde A_{yz}=\frac {\mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})}{De}\tilde {A}_{zz},\qquad\quad \end{align}
\begin{align} \frac {\partial U_z}{\partial Y}\tilde A_{yy}=\frac {\mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})}{De}\tilde {A}_{yz}, \qquad\quad\end{align}
\begin{align} \tilde {A}_{yy}=1. \qquad\qquad\qquad\end{align}
Substituting (A1c) into (A1a) yields
\begin{equation} \frac {\mathrm {d} P}{\mathrm {d} Z}=\frac {\mathrm {d}^2 U_z}{\mathrm {d} Y^2}. \end{equation}
Solving for the velocity
$U_z$
subject to (2.14), we obtain a parabolic profile
\begin{equation} U_{z}(Y)=\frac {3}{2}\frac {H_{0}^2-Y^2}{H_{0}^3}. \end{equation}
Next, substituting (A1c) into (A1b), and using (A1d) and (2.13) leads to the nonlinear algebraic equation for
$\tilde A_{zz}$
\begin{equation} 2De^2\left (1-\frac {\tilde A_{zz}}{L^2\epsilon ^2}\right )^2\left (\frac {\partial U_z}{\partial Y}\right )^2=\tilde A_{zz}. \end{equation}
The corresponding solution of (A4) is
\begin{equation} \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}. \end{equation}
Combining (A1c) and (A5) provides the expression for
$\tilde A_{yz}$
:
\begin{equation} \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}. \end{equation}
Finally, we note that considering the limit
$L^2\epsilon ^2 \to \infty$
and using (A5), (A6), we obtain the corresponding expressions for the conformation tensor components of the Oldroyd-B fluid:
\begin{equation} \tilde {A}_{zz}=\frac {18De^{2}}{H_{0}^{6}}Y^{2},\qquad \tilde {A}_{yz}=-\frac {3De}{H_{0}^{3}}Y,\qquad \tilde{A}_{yy}=1. \end{equation}
Appendix B. Non-dimensional pressure drop across the contraction
The integral mass conservation along the channel (2.14
$d$
) sets the local value of the pressure gradient and allows one to calculate the pressure drop without solving for the velocity field. Integrating by parts the integral constraint (2.14
$d$
) and using (2.14
$a$
) and (2.14
$c$
), we obtain (Boyko et al. Reference Boyko, Hinch and Stone2024; Hinch et al. Reference Hinch, Boyko and Stone2024)
\begin{equation} 1=\int ^{H(Z)}_0 U_{z}\mathrm {d}Y=-\int _{0}^{H(Z)}Y\frac {\partial U_z}{\partial Y}\mathrm {d}Y=-\frac {1}{2}\int _{0}^{H(Z)}(H(Z)^2-Y^{2})\frac {\partial ^{2}U_z}{\partial Y^{2}}\mathrm {d}Y. \end{equation}
Substituting the expression for
$\partial ^{2}U_z/\partial Y^{2}$
from the momentum equation (2.12
$b$
) into (B1) and rearranging provides an expression for the pressure gradient,
\begin{align} \frac {\mathrm {d}P}{\mathrm {d}Z}=&-\frac {3(1-\tilde {\beta })}{H(Z)^{3}}\nonumber\\&+\frac {3 \tilde {\beta }}{2De H(Z)^{3}}\int _{0}^{H(Z)}(H(Z)^2-Y^{2})\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 ]\mathrm {d}Y. \end{align}
Next, integrating (B2) with respect to
$Z$
from 0 to 1 yields the pressure drop
$\triangle P=P(0)-P(1)$
across the non-uniform region
\begin{eqnarray} &&\triangle P=3(1-\tilde {\beta }) \int _{0}^{1} \frac {\mathrm {d} Z}{H(Z)^{3}} \nonumber \\ &&\,\,-\frac {3 \tilde {\beta }}{2De }\int _{0}^{1}\frac {1}{ H(Z)^{3}}\int _{0}^{H(Z)}(H(Z)^2-Y^{2})\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 ]\mathrm {d}Y \mathrm {d} Z.\quad\;\;\;\;\; \end{eqnarray}
Finally, using integration by parts, (B3) can be expressed as in the form given in (2.16).
Appendix C. Low-Deborah-number lubrication analysis: detailed derivation
We here provide details of the derivation of the analytical expressions for the pressure drop of the FENE-CR fluid in the low-
$De$
limit up to
$O(De^4)$
.
Before proceeding to the asymptotic solution of the pressure drop, we expand
$\mathcal {F}(\,\tilde {{\boldsymbol{\!A}}}) \tilde A_{zz}$
,
$\mathcal {F}(\,\tilde {{\boldsymbol{\!A}}}) \tilde A_{yz}$
and
$\mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})(\tilde A_{yy}-1)$
into perturbation series in
$De\ll 1$
. Specifically, using (2.13), (3.1), (3.2) and noting that
$\tilde A_{zz,1}=0$
, we obtain
\begin{equation} \mathcal {F}(\,\tilde {{\boldsymbol{\!A}}})\tilde A_{zz}=De^2\tilde A_{zz,2}+De^3\tilde A_{zz,3}+De^4 \left [\tilde A_{zz,4}+\frac {(\tilde A_{zz,2})^2}{L^2 \epsilon ^2}\right ]+O(De^5), \end{equation}
\begin{eqnarray} \mathcal {F}(\,\tilde {{\boldsymbol{\!A}}}) \tilde A_{yz}&=&De\tilde A_{yz,1}+De^2\tilde A_{yz,2}+De^3 \left [\tilde A_{yz,3}+\frac {\tilde A_{yz,1} \tilde A_{zz,2}}{L^2 \epsilon ^2}\right ] \nonumber \\ &&+De^4 \left [\tilde A_{yz,4}+\frac {\tilde A_{yz,2}\tilde A_{zz,2}+\tilde A_{yz,1}\tilde A_{zz,3}}{L^2 \epsilon ^2}\right ]+O(De^5), \end{eqnarray}
\begin{eqnarray} \mathcal {F}(\,\tilde {{\boldsymbol{\!A}}}) (\tilde A_{yy}-1)&=&De\tilde A_{yy,1}+De^2\tilde A_{yy,2}+De^3 \left [\tilde A_{yy,3}+\frac {\tilde A_{yy,1} \tilde A_{zz,2}}{L^2 \epsilon ^2}\right ] \nonumber \\ &&+De^4 \left [\tilde A_{yy,4}+\frac {\tilde A_{yy,2}\tilde A_{zz,2}+\tilde A_{yy,1}\tilde A_{zz,3}}{L^2 \epsilon ^2}\right ]+O(De^5). \end{eqnarray}
C.1. Leading-order solution for the pressure drop of a FENE-CR fluid
Substituting (3.1) into (2.12) and considering the leading order in
$De$
, and using (3.3) and (C1), we obtain
\begin{align} \frac {\partial U_{z,0}}{\partial Z}+\frac {\partial U_{y,0}}{\partial Y}=0,\qquad \frac {\mathrm {d} P_0}{\mathrm {d} Z}=(1-\tilde \beta )\frac {\partial ^2 U_z}{\partial Y^2}+\tilde \beta \frac {\partial \tilde A_{yz,1}}{\partial Y}=\frac {\partial ^2 U_{z,0}}{\partial Y^2}, \end{align}
subject to the boundary conditions
\begin{align} U_{z,0}(H(Z),Z)=0, \; U_{y,0}(H(Z),Z)=0, \; \frac {\partial U_{z,0}}{\partial Y}(0,Z)=0, \; \int ^{H(Z)}_0 U_{z,0}(Y,Z)\mathrm {d}Y=1. \end{align}
As expected, (C2
$b$
) is the classical momentum equation of the Newtonian fluid with a constant viscosity
$\mu _0$
. The leading-order solutions, previously derived by Boyko & Stone (Reference Boyko and Stone2022), are given as
\begin{align} U_{z,0}=\frac {3}{2}\frac {H(Z)^2-Y^2}{H(Z)^3}, \quad U_{y,0}=\frac {3}{2}\frac {H'(Z)Y(H(Z)^2-Y^2)}{H(Z)^4}, \quad \triangle P_{0}=3\int _{0}^{1}\frac {\mathrm {d}Z}{H(Z)^{3}},\nonumber\\ \end{align}
where primes indicate derivatives with respect to
$Z$
.
C.2. First-order solution for the pressure drop of a FENE-CR fluid
Substituting (3.1) and (C1) into (2.12) and considering the first order in
$De$
, we obtain
\begin{align} \frac {\partial U_{z,1}}{\partial Z}+\frac {\partial U_{y,1}}{\partial Y}=0,\qquad \frac {\mathrm {d} P_1}{\mathrm {d} Z}=(1-\tilde \beta )\frac {\partial ^2 U_{z,1}}{\partial Y^2}+\tilde \beta \left ( \frac {\partial \tilde A_{zz,2}}{\partial Z}+\frac {\partial \tilde A_{yz,2}}{\partial Y}\right ),\,\, \end{align}
\begin{align} 2\frac {\partial U_{z,0}}{\partial Y}\tilde A_{yz,1}=\tilde A_{zz,2}, \qquad\qquad\qquad\qquad\qquad\end{align}
\begin{align} U_{z,0}\frac {\partial \tilde A_{yz,1}}{\partial Z}+U_{y,0}\frac {\partial \tilde A_{yz,1}}{\partial Y}-\frac {\partial U_{z,0}}{\partial Y}\tilde A_{yy,1}-\frac {\partial U_{z,1}}{\partial Y}=-\tilde A_{yz,2}, \qquad\quad\end{align}
\begin{align} U_{z,0}\frac {\partial \tilde A_{yy,1}}{\partial Z}+U_{y,0}\frac {\partial \tilde A_{yy,1}}{\partial Y}-2\frac {\partial U_{y,0}}{\partial Z}\tilde A_{yz,1}-2\frac {\partial U_{ y,0}}{\partial Y}\tilde A_{yy,1}-2\frac {\partial U_{y,1}}{\partial Y}=-\tilde A_{yy,2}. \end{align}
These governing equations are supplemented by the boundary conditions
\begin{align} U_{z,1}(H(Z),Z)=0, \ U_{y,1}(H(Z),Z)=0, \ \frac {\partial U_{z,1}}{\partial Y}(0,Z)=0, \ \int ^{H(Z)}_0 U_{z,1}(Y,Z)\mathrm {d}Y=0. \end{align}
At the first order in
$De$
, the dimensionless governing equations for the FENE-CR fluid are equivalent to those of the Oldroyd-B fluid. Thus, from (C5), it follows that the expressions for the velocity and pressure drop at
$O(De)$
as well as
$\tilde A_{zz,2}$
,
$\tilde A_{yz,2}$
,
$\tilde A_{yy,2}$
are identical for the FENE-CR and Oldroyd-B fluids, and are given by (Boyko & Stone Reference Boyko and Stone2022)
\begin{align} \tilde {A}_{zz,2}=\frac {18 Y^2}{H(Z)^6}, \qquad \tilde {A}_{yz,2}=\frac {18 Y \left (2 Y^2-H(Z)^2\right ) H'(Z)}{H(Z)^7},\qquad \end{align}
\begin{align} \tilde {A}_{yy,2}=\frac {9}{2}\frac {4\left (-2Y^{2}+H(Z)^{2}\right )^{2}H'(Z)^{2}-H(Z)H''(Z)\left (Y^{2}-H(Z)^{2}\right )^{2}}{H(Z)^{8}}. \end{align}
We note that the FENE-CR and Oldroyd-B fluids exhibit a second-order fluid behaviour at
$O(De)$
, so that the velocity field remains Newtonian, i.e.
$U_{z,1}=U_{y,1}=0$
, following the theorem of Tanner and Pipkin (Tanner Reference Tanner1966; Tanner & Pipkin Reference Tanner and Pipkin1969).
C.3. Second-order solution for the pressure drop of a FENE-CR fluid
At the second order,
$O(De^2)$
, the governing equations (2.12) yield
\begin{equation} \frac {\partial U_{z,2}}{\partial Z}+\frac {\partial U_{y,2}}{\partial Y}=0, \end{equation}
\begin{equation} \frac {\mathrm {d} P_2}{\mathrm {d} Z}=(1-\tilde \beta )\frac {\partial ^2 U_{z,2}}{\partial Y^2}+\tilde \beta \left [ \frac {\partial \tilde A_{zz,3}}{\partial Z}+\frac {\partial }{\partial Y}\left (\tilde A_{yz,3}+\frac {\tilde A_{yz,1} \tilde A_{zz,2}}{L^2 \epsilon ^2}\right )\right ], \end{equation}
\begin{equation} U_{z,0}\frac {\partial \tilde A_{zz,2}}{\partial Z}+U_{y,0}\frac {\partial \tilde A_{zz,2}}{\partial Y}-2 \frac {\partial U_{z,0}}{\partial Z}\tilde A_{zz,2}-2 \frac {\partial U_{z,0}}{\partial Y}\tilde A_{yz,2}=-\tilde A_{zz,3}, \end{equation}
\begin{align} & U_{z,0}\frac {\partial \tilde A_{yz,2}}{\partial Z}+U_{y,0}\frac {\partial \tilde A_{yz,2}}{\partial Y}-\frac {\partial U_{y,0}}{\partial Z}\tilde A_{zz,2}-\frac {\partial U_{z,0}}{\partial Y}\tilde A_{yy,2}-\frac {\partial U_{z,2}}{\partial Y}=-\tilde A_{yz,3}-\frac {\tilde A_{yz,1}\tilde A_{zz,2}}{L^2 \epsilon ^2},\nonumber \\[2pt] & \end{align}
\begin{align} U_{z,0}\frac {\partial \tilde A_{yy,2}}{\partial Z}+U_{y,0}\frac {\partial \tilde A_{yy,2}}{\partial Y}&-2\frac {\partial U_{y,0}}{\partial Z}\tilde A_{yz,2}-2\frac {\partial U_{ y,0}}{\partial Y}\tilde A_{yy,2}-2\frac {\partial U_{y,2}}{\partial Y}\nonumber\\=&-\tilde A_{yy,3}-\frac {\tilde A_{yy,1}\tilde A_{zz,2}}{L^2 \epsilon ^2}, \end{align}
where we have used the expressions
$\tilde A_{zz,1}=0$
,
$U_{z,1}=0$
and
$U_{y,1}=0$
. The governing equations (C8) are subject to the boundary conditions
\begin{align} U_{z,2}(H(Z),Z)=0, \; U_{y,2}(H(Z),Z)=0, \; \frac {\partial U_{z,2}}{\partial Y}(0,Z)=0, \; \int ^{H(Z)}_0 U_{z,2}(Y,Z)\mathrm {d}Y=0. \end{align}
We note that the evolution equation for
$\tilde A_{zz,3}$
, given in (C8c), is the same for the FENE-CR and Oldroyd-B fluids. In contrast, the evolution equations for
$\tilde A_{yz,3}$
and
$\tilde A_{yy,3}$
, given in (C8d) and (C8e), are different for the two fluids due to additional terms for the FENE-CR fluid, which depend on
$L^2 \epsilon ^2$
. Nevertheless, similar to the first order, the expressions for the velocity and pressure drop at
$O(De^2)$
are the same for the FENE-CR and Oldroyd-B fluids. This can be seen by substituting (C8d) into the last term on the right-hand side of the momentum equation (C8b), thus clearly showing that the velocity and pressure are independent of
$L^2 \epsilon ^2$
at
$O(De^2)$
.
The resulting expressions for
$U_{z,2}$
,
$U_{y,2}$
and
$\tilde A_{zz,3}$
,
$\tilde A_{yz,3}$
,
$\tilde A_{yy,3}$
are readily found using MATHEMATICA software, but they are rather lengthy and, thus, not presented here. As the
$\tilde A_{yz,3}$
and
$\tilde A_{yy,3}$
for the FENE-CR fluid are coupled to
$L^2 \epsilon ^2$
, we expect the pressure drop to depend on the finite extensibility at the next order,
$O(De^3)$
. We show this dependence in the following subsection.
C.4. Third-order solution for the pressure drop of a FENE-CR fluid
Substituting (3.1) and (C1) into (2.12) and considering the third order in
$De$
, we obtain
\begin{equation} \frac {\partial U_{z,3}}{\partial Z}+\frac {\partial U_{y,3}}{\partial Y}=0, \end{equation}
\begin{align} \frac {\mathrm {d} P_3}{\mathrm {d} Z}&=(1-\tilde \beta )\frac {\partial ^2 U_{z,3}}{\partial Y^2}\nonumber \\ &\quad + \tilde \beta \left [ \frac {\partial }{\partial Z}\left (\tilde A_{zz,4}+\frac {(\tilde A_{zz,2})^2}{L^2 \epsilon ^2}\right )+\frac {\partial }{\partial Y}\left (\tilde A_{yz,4}+\frac {\tilde A_{yz,2}\tilde A_{zz,2}+\tilde A_{yz,1}\tilde A_{zz,3}}{L^2 \epsilon ^2}\right )\right ], \end{align}
\begin{eqnarray} U_{z,0}\frac {\partial \tilde A_{zz,3}}{\partial Z}&+&U_{y,0}\frac {\partial \tilde A_{zz,3}}{\partial Y}-2 \frac {\partial U_{z,0}}{\partial Z}\tilde A_{zz,3}-2 \frac {\partial U_{z,0}}{\partial Y}\tilde A_{yz,3} \nonumber \\ &-&2 \frac {\partial U_{z,2}}{\partial Y}\tilde A_{yz,1}=-\left (\tilde A_{zz,4}+\frac {(\tilde A_{zz,2})^2}{L^2 \epsilon ^2}\right ), \end{eqnarray}
\begin{eqnarray} U_{z,0}\frac {\partial \tilde A_{yz,3}}{\partial Z}&+&U_{z,2}\frac {\partial \tilde A_{yz,1}}{\partial Z}+U_{y,0}\frac {\partial \tilde A_{yz,3}}{\partial Y}+U_{y,2}\frac {\partial \tilde A_{yz,1}}{\partial Y}-\frac {\partial U_{y,0}}{\partial Z}\tilde A_{zz,3}-\frac {\partial U_{z,0}}{\partial Y}\tilde A_{yy,3} \nonumber \\ &-&\frac {\partial U_{z,2}}{\partial Y}\tilde A_{yy,1}-\frac {\partial U_{z,3}}{\partial Y}=-\left (\tilde A_{yz,4}+\frac {\tilde A_{yz,2}\tilde A_{zz,2}+\tilde A_{yz,1}\tilde A_{zz,3}}{L^2 \epsilon ^2}\right ), \end{eqnarray}
\begin{align} &U_{z,0}\frac {\partial \tilde A_{yy,3}}{\partial Z}+U_{z,2}\frac {\partial \tilde A_{yy,1}}{\partial Z}+U_{y,0}\frac {\partial \tilde A_{yy,3}}{\partial Y}+U_{y,2}\frac {\partial \tilde A_{yy,1}}{\partial Y}-2\frac {\partial U_{y,0}}{\partial Z}\tilde A_{yz,3}-2\frac {\partial U_{y,2}}{\partial Z}\tilde A_{yz,1}\nonumber \\ &-2\frac {\partial U_{y,0}}{\partial Y}\tilde A_{yy,3}-2\frac {\partial U_{y,2}}{\partial Y}\tilde A_{yy,1}-2\frac {\partial U_{y,3}}{\partial Y}=-\left (\tilde A_{yy,4}+\frac {\tilde A_{yy,2}\tilde A_{zz,2}+\tilde A_{yy,1}\tilde A_{zz,3}}{L^2 \epsilon ^2}\right ), \end{align}
where we have used the expressions
$\tilde A_{zz,1}=0$
,
$U_{z,1}=0$
and
$U_{y,1}=0$
. The governing equations (C10) are subject to the boundary conditions
\begin{align} U_{z,3}(H(Z),Z)=0, \; U_{y,3}(H(Z),Z)=0, \; \frac {\partial U_{z,3}}{\partial Y}(0,Z)=0, \; \int ^{H(Z)}_0 U_{z,3}(Y,Z)\mathrm {d}Y=0. \end{align}
First, we integrate (C10b) twice with respect to
$Y$
and apply the boundary conditions (C11
$a$
) and (C11
$c$
) to obtain the expression for
$U_{z,3}(Y,Z)$
that involves the pressure gradient
$\mathrm {d}P_3/\mathrm {d}Z$
. The resulting expression is lengthy and thus not shown here. To determine
$\mathrm {d}P_3/\mathrm {d}Z$
, we use the integral constraint (C11
$d$
), leading to
\begin{align} \frac {\mathrm {d}P_3}{\mathrm {d}Z} &=\frac {10692 \tilde \beta H'(Z)}{35 L^2 \epsilon ^2 H(Z)^9} \nonumber \\ &\quad+\frac {216 \tilde \beta }{35}\left [(\tilde \beta -8)\frac {H'''(Z)}{H(Z)^7}+ (110-13 \tilde \beta ) \frac {H'(Z) H''(Z)}{H(Z)^8}+24 (\tilde \beta -9) \frac {H'(Z)^3}{H(Z)^9}\right ]. \end{align}
Integrating (C12) with respect to
$Z$
from 0 to 1 provides an expression for the pressure drop of the FENE-CR fluid at
$O(De^3)$
given in (3.5).
C.5. Fourth-order solution for the pressure drop of a FENE-CR fluid
To calculate the pressure drop at the next order,
$O(De^4)$
, we use the expression (2.16), which resembles the result of an application of the reciprocal theorem (Boyko & Stone Reference Boyko and Stone2021, Reference Boyko and Stone2022), and requires only the knowledge of velocity and conformation tensor components from the previous orders. At
$O(De^4)$
, the expression for the pressure drop
$\triangle P_4$
takes the form
\begin{eqnarray} \triangle P_4&=&\tilde \beta \int ^{H(0)}_0[\mathcal {G}_{zz,5}\hat U_z]_{Z=0} \mathrm {d}Y-\tilde \beta \int ^{H(1)}_0[\mathcal {G}_{zz,5}\hat U_z]_{Z=1} \mathrm {d}Y \nonumber \\ &&+\tilde \beta \int ^1_0\int ^{H(Z)}_0\left (\mathcal {G}_{zz,5}\frac {\partial \hat U_z}{\partial Z}+\mathcal {G}_{yz,5}\frac {\partial \hat U_z}{\partial Y}\right ) \mathrm {d}Y \mathrm {d}Z, \end{eqnarray}
where
$\mathcal {G}_{zz,5}$
and
$\mathcal {G}_{yz,5}$
are given by
\begin{eqnarray} \mathcal {G}_{zz,5}&=&-U_{z,0}\frac {\partial \tilde A_{zz,4}}{\partial Z}-U_{z,2}\frac {\partial \tilde A_{zz,2}}{\partial Z}-U_{y,0}\frac {\partial \tilde A_{zz,4}}{\partial Y}-U_{y,2}\frac {\partial \tilde A_{zz,2}}{\partial Y}+2 \frac {\partial U_{z,0}}{\partial Z}\tilde A_{zz,4} \nonumber \\ &&+2 \frac {\partial U_{z,2}}{\partial Z}\tilde A_{zz,2}+2 \frac {\partial U_{z,0}}{\partial Y}\tilde A_{yz,4}+2 \frac {\partial U_{z,2}}{\partial Y}\tilde A_{yz,2}+2 \frac {\partial U_{z,3}}{\partial Y}\tilde A_{yz,1}, \end{eqnarray}
\begin{align} \mathcal {G}_{yz,5}= &-U_{z,0}\frac {\partial \tilde A_{yz,4}}{\partial Z}-U_{z,2}\frac {\partial \tilde A_{yz,2}}{\partial Z}-U_{z,3}\frac {\partial \tilde A_{yz,1}}{\partial Z}-U_{y,0}\frac {\partial \tilde A_{yz,4}}{\partial Y}-U_{y,2}\frac {\partial \tilde A_{yz,2}}{\partial Y} \nonumber \\ &-U_{y,3}\frac {\partial \tilde A_{yz,1}}{\partial Y} +\frac {\partial U_{y,0}}{\partial Z}\tilde A_{zz,4}+\frac {\partial U_{y,2}}{\partial Z}\tilde A_{zz,2}+\frac {\partial U_{z,0}}{\partial Y}\tilde A_{yy,4}+\frac {\partial U_{z,2}}{\partial Y}\tilde A_{yy,2} \nonumber \\ &+\frac {\partial U_{z,3}}{\partial Y}\tilde A_{yy,1}+\frac {\partial U_{z,4}}{\partial Y}. \end{align}
We note that, because of the integral constraint
$\int _{0}^{H(Z)}U_{z,4}\mathrm {d}Y=0$
, the last term appearing in (C14b),
$\partial U_{z,4}/\partial Y$
, satisfies
\begin{equation} \int ^{H(Z)}_0\frac {\partial U_{z,4}}{\partial Y}\frac {\partial \hat U_z}{\partial Y} \mathrm {d}Y =-\int ^{H(Z)}_0 U_{z,4}\frac {\partial ^2 \hat U}{\partial Y^2} \mathrm {d}Y =-\frac {\mathrm {d}\hat P}{\mathrm {d} Z}\int ^{H(Z)}_0U_{z,4} \mathrm {d}Y=0, \end{equation}
and thus, this term does not contribute to the pressure drop, since it is identically zero. Therefore, the expressions for
$\mathcal {G}_{zz,5}$
and
$\mathcal {G}_{yz,5}$
depend on the solution from the previous orders, and we can calculate the fourth-order pressure drop
$\triangle P_4$
using the results of the leading-, first-, second- and third-order viscoelastic problems. The resulting expression for
$\triangle P_4$
for the FENE-CR fluid is given in (3.6).
For completeness, in the supplementary material available at http://doi.org/10.1017/jfm.2025.142, we provide the MATHEMATICA file containing the explicit expressions for the velocity, conformation tensor components and the pressure drop in the low-Deborah-number limit up to
$O(De^4)$
.
Appendix D. Details of numerical simulations using OpenFOAM
In this appendix, we describe the numerical procedure used to solve the system of nonlinear governing equations (2.1)–(2.5) for the viscoelastic fluid flow. In addition to the FENE-CR fluid, we also consider the Oldroyd-B fluid for comparison and validation. We have performed two-dimensional finite-volume simulations using an open-source framework OpenFOAM (Jasak et al. Reference Jasak, Jemcov and Tukovic2007) integrated with viscoelastic flow solver RheoTool (Pimenta & Alves Reference Pimenta and Alves2017). We use the log-conformation method to calculate the polymer stress tensor by solving the equations for the logarithm of the conformation tensor
$\boldsymbol {\Theta }$
instead of
$\boldsymbol \tau _p$
(Pimenta & Alves Reference Pimenta and Alves2017; Habla et al. Reference Habla, Tan, Haßlberger and Hinrichsen2014; Kumar et al. Reference Kumar, Aramideh, Browne, Datta and Ardekani2021; Kumar & Ardekani Reference Kumar and Ardekani2021). The log-conformation approach delays the numerical instability caused by high Deborah/Weissenberg numbers by maintaining the positive definiteness of the conformation tensor (Fattal & Kupferman Reference Fattal and Kupferman2004, Reference Fattal and Kupferman2005). The details of the numerical implementation and the code validation are given in prior studies (see e.g. Pimenta & Alves Reference Pimenta and Alves2017; Favero et al. Reference Favero, Secchi, Cardozo and Jasak2010).
In our simulations, we impose the no-slip and no-penetration boundary conditions along the wall,
$y= \pm h(z)$
, and a fully developed unidirectional Poiseulille velocity profile at the entrance (
$z=-\ell _{0}$
) and exit (
$z=\ell _{\ell }$
). In addition, we specify a null value of polymeric stress tensor and zero-gradient of pressure at the channel entrance. At the channel wall, we impose a linear extrapolation for polymer stresses and zero-normal gradient for pressure (Pimenta & Alves Reference Pimenta and Alves2017). At the exit (
$z=\ell _{\ell }$
), we use a zero-gradient boundary condition for polymer stresses and prescribe a constant value for pressure,
$p=0$
. Finally, we calculate the pressure drop along the centreline between the inlet (
$z=0$
) and outlet (
$z=\ell$
) of the contraction, i.e.
$\triangle p=p(y=0,z=0)-p(y=0,z=\ell )$
, eliminating the entrance and exit effects.
We summarise in table 3 the values of physical and geometrical parameters used in the numerical simulations. We consider a geometry with an inlet-to-outlet ratio
$H_0=h_0/h_\ell =4$
and an aspect ratio
$\epsilon =h_\ell /\ell =0.02$
, and explore two different polymer-to-total viscosity ratios:
$\tilde {\beta }=\mu _p/\mu _0=0.4$
and
$\tilde {\beta }=\mu _p/\mu _0=0.05$
, where the latter corresponds to the ultra-dilute limit. In all simulations, we keep
$\ell _0=\ell$
and
$\ell _\ell =5\ell$
.
Values of physical and geometrical parameters used in the two-dimensional numerical simulations of the pressure-driven flow of the FENE-CR fluid in a hyperbolic contracting channel.
To study the effect of Deborah numbers in the case of the FENE-CR fluid, we mainly set the finite extensibility parameter to
$L^2=1250$
, corresponding to
$L^2\epsilon ^2=0.5$
, and change the value of relaxation time
$\lambda$
, while keeping the values of all other parameters. When investigating the effect of finite extensibility
$L^2\epsilon ^2$
on the pressure drop, we change the value of
$L^2$
and set different
$\lambda$
corresponding to different values of
$De$
, while keeping the values of all other parameters. Similarly, when analysing the pressure drop at different viscosity ratios
$\mu _p/\mu _0$
, we change the value of
$\mu _p$
and
$\mu _s$
, while setting
$\mu _0=1$
Pa s and keeping the values of all other parameters. We note that the effect of fluid inertia is negligible in our simulations because the reduced Reynolds number
$\epsilon Re=(h_{\ell }/\ell )\rho u_{c}h_{\ell }/\mu _{0}= 10^{-8}$
is very small. Eventually, we use the transient rheoFoam solver (Pimenta & Alves Reference Pimenta and Alves2017) for simulations, and once the residuals of the variables
$\boldsymbol {u}$
,
$p$
and
$\boldsymbol {\Theta }$
become less than
$10^{-6}$
, we terminate the simulation and calculate the pressure drop. We non-dimensionalise the time
$t$
using the residence time in the contraction
$t_c=\ell /u_c=1$
s. Typical non-dimensional values of the time step are
$\triangle T= 10^{-4}$
for the low-
$De$
simulations and a reduced time step
$\triangle T= 10^{-5}$
for the high-
$De$
simulations.
To assess the grid sensitivity, we have performed tests by considering three different mesh resolutions (total number of node points is 75 672, 114 882, and 139 482) at four different Deborah numbers (
$De$
= 1, 2, 3 and 4), and established grid independence with a maximum relative error of 0.3 % for the pressure drop. We have also carried out numerical simulations without the log-conformation approach and found an excellent agreement with the log-conformation results.
In addition, we cross-validate our OpenFOAM results for Oldroyd-B and FENE-CR fluids with those obtained from the finite-element software COMSOL Multiphysics (version 6.2, COMSOL AB, Stockholm, Sweden). The details of the numerical implementation in COMSOL are given by Boyko & Stone (Reference Boyko and Stone2022) for the Oldroyd-B fluid. To simulate the FENE-CR fluid in COMSOL, we impose the polymer stress distribution corresponding to the Poiseuille flow at the entrance.
Comparison of simulation results obtained from OpenFOAM and COMSOL for the pressure drop 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 })$
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. Black squares and red crosses represent the COMSOL simulation results for the Oldroyd-B and FENE-CR fluids. Cyan dotted and solid black lines represent the low-
$De$
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 9 the scaled pressure drop
$\triangle P/\triangle P_0$
for the Oldroyd-B and FENE-CR fluids as a function of
$De$
in a contracting channel for
$L^2\epsilon ^2=0.5$
(panel
$a$
) and
$L^2\epsilon ^2=0.25$
(panel
$b$
), with
$H_0=4$
and
$\tilde \beta =0.4$
. Grey triangles and purple circles represent the OpenFOAM simulation results for the Oldroyd-B and FENE-CR fluids. Black squares and red crosses represent the COMSOL simulation results 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. In COMSOL simulations of the Oldroyd-B fluid, we could not obtain converged results beyond
$De=0.45$
. In contrast, using OpenFOAM, we have performed simulations up to
$De=4$
with no difficulties, thus achieving the high-
$De$
limit. We encountered no convergence issues when running simulations with the FENE-CR model in OpenFOAM (up to
$De=4$
) and COMSOL (up to
$De=0.8$
). Clearly, for both Oldroyd-B and FENE-CR fluids, there is excellent agreement between the simulation results obtained from OpenFOAM and COMSOL. In particular, for the Oldroyd-B fluid, the maximum relative error is 1.3 % at
$De=0.45$
. Similarly, for the FENE-CR fluid, we find a maximum relative error of 0.4 % and 0.3 % at
$De=0.8$
, corresponding to
$L^2\epsilon ^2=0.5$
and
$L^2\epsilon ^2=0.25$
, respectively.
Open access
