1. Introduction
Fluid mixtures comprising elastic material and subjected to fast flow processes exhibit strongly non-Newtonian traits. Examples of such fluid flow are commonly found in both biological and mechanical systems, such as capillary blood flow (Pandey et al. Reference Pandey, Karpitschka, Venner and Snoeijer2016), human synovial joints consisting of a lubricating fluid and rapid motion of the joint (Tichy & Bou-Saïd Reference Tichy and Bou-Saïd2008; Yousfi, Bou-Saïd & Tichy Reference Yousfi, Bou-Saïd and Tichy2013), tear films subjected to rapid movement of the eyelids (Dunn et al. Reference Dunn, Tichy, Urueña and Sawyer2013) and mechanical systems. The lubrication of these components has been around for several centuries and serves an important function in maintaining the performance of the machine elements. In thin-film hydrodynamic lubrication, a thin layer of fluid (liquid or gas) persists between two moving surfaces (bearings, gears, artificial joints, etc.).
The lubricant, when engineered, is often a mixture comprising a base Newtonian solvent which is then enhanced via the addition of polymeric additives (Mortier, Orszulik & Fox Reference Mortier, Orszulik and Fox2010). These additives are often necessary in counteracting the unavoidable effects due to high shear stresses and high pressure, such as (i) shear thinning (Ahmed & Biancofiore Reference Ahmed and Biancofiore2022) in which the shear stress causes a decline in viscosity, (ii) viscoplasticity, an irreversible decrease in lubricant viscosity versus the applied shear stress, (iii) piezoviscosity (Jeng, Hamrock & Brewe Reference Jeng, Hamrock and Brewe1986), an increase in the lubricant viscosity versus the film pressure, (iv) thermal degradation due to uncoiling of polymeric chains as the film temperature increases and also (v) cavitation of the film, the generation of a liquid–vapour mixture region due to a drop in film pressure below cavitation pressure (Biancofiore, Giacopini & Dini Reference Biancofiore, Giacopini and Dini2019; Gamaniel, Dini & Biancofiore Reference Gamaniel, Dini and Biancofiore2021). These effects are studied extensively in tribology to optimise the performance of machine elements (Ahmed & Biancofiore Reference Ahmed and Biancofiore2022). The emergence of viscoelasticity, arising purely due to a finitely elastic polymer additive, is somewhat overlooked owing to the nonlinear interaction of these phenomena which makes it difficult to isolate any single effect, at least experimentally.
In classical thin-film hydrodynamic lubrication (where the bounding surfaces are separated and assumed rigid) the film thickness,
$H_0$
, is some orders of magnitude less than the channel length,
$\ell _x$
(characterised by the ratio
$\epsilon = H_0 / \ell _x \ll 1$
). There are several operating modes depending on the motion of the surfaces, e.g. sliding, rolling, squeezing and a combination of thereof (Szeri Reference Szeri2010). The lubrication mode and the channel surface configuration unsurprisingly have a strong influence on the film pressure distribution and magnitude. For channel configurations where the bounding surfaces are in relative motion, i.e. a sliding configuration in which at least one surface has a fixed speed,
$U_0$
(such as bearings, meshing of gear teeth), the film pressure is generated only if some surface gradient is present. The strengths of these pressure gradients are measured via the bearing number, defined using the ambient pressure
$p_a$
,
$\varLambda = \eta _0 U_0 \ell _x / p_a H_0^2$
, i.e. the ratio between the viscous and pressure forces (Kundu, Cohen & Dowling Reference Kundu, Cohen and Dowling2015), which tends to be large leading to high pressure gradients. A combination of the two (large
$\varLambda$
and surface sliding) exposes the elastic polymeric additives to a high applied shear rate that persists along the entire channel preventing their otherwise relaxed state (when the polymeric stress decays to zero or reaches its Newtonian value). We focus primarily on the impact of viscoelasticity on the film pressure (and its ability to bear the applied load, the so-called load-carrying capacity) in the sliding case which has received comparatively little attention when compared with other nonlinear phenomena mentioned previously.
For the thin sliding channel, we have (i) a high shear strain rate,
$\dot {\gamma } \approx U_0 / H_0$
, and (ii) a residence time of the order of
$\ell _x / U_0$
, where
$U_0$
is some relative surface speed. When elastic additives are exposed to these two distinct time scales they tend to exhibit characteristics strongly differing from the classical Newtonian behaviour such as the onset of shear thinning. However, when the Deborah number
$\textit{De} = \lambda U_0 / \ell _x$
, i.e. the ratio of the polymer relaxation time
$\lambda \approx \eta _p / G$
(
$\eta _p$
is the polymer’s viscosity and
$G$
is the shear modulus) to the residence time scale
$\ell _x / U_0$
, is not negligible then viscoelastic effects also emerge. In addition, the Weissenberg number, i.e. the ratio of the polymer relaxation time to the shear time scale
$\textit{Wi} = \lambda \dot {\gamma }$
, tends to acquire large values since the shear strain rate can be large in practice. These two parameters when plotted against one another span the Pipkin space which characterises the problem (Tanner Reference Tanner2000), e.g. linearly viscoelastic (
$\textit{De} \ll 1$
,
$\textit{Wi} \ll 1$
) or viscometric (
$\textit{De} \ll 1$
,
$\textit{Wi} \gg 1$
). For thin-film lubrication they are related via the thin-film ratio
$\textit{De} = \epsilon \textit{Wi}$
. In these types of lubricated systems, the thin-film ratio prevents the Deborah number from being O(
$\textit{Wi}$
), but considering them as independent parameters, one can also place the flow in the
$\textit{De} \gg 1$
and
$\textit{Wi} \gg 1$
region in which the fluid acts as an elastic solid. In such a region, the validity of the fluid constitutive relations like Oldroyd-B is questionable (Housiadas & Beris Reference Housiadas and Beris2024b
).
Experimental evidence suggests that increasing the fluid elasticity has an observable impact, beyond the typical viscous effects. For example, in the case of sliding hydrodynamic lubrication (such as that in bearings and gears), it increases the lubricating film’s load-carrying capacity and thereby enhances the performance of the machine components. Tests conducted via a cone-and-plate rheometer have demonstrated a finite increase in the load-carrying capacity when viscoelasticity is present (Schuh et al. Reference Schuh, Lee, Allison and Ewoldt2017). In fact, when elastic compounds in the form of amino acid ionic liquids are added to water-based lubricants (which remain Newtonian even under extreme flow conditions), an increase in the first normal stress difference is found (Feng & Jabbarzadeh Reference Feng and Jabbarzadeh2024), corresponding to a positive net force. Similarly, a thin film squeezed between two parallel plates predicts an increased first normal stress difference and load-carrying capacity (Tichy & Winer Reference Tichy and Winer1978). However, owing to the small dimensions of the contact region, point-wise variation of the pressure is difficult to obtain via experimental approaches.
Numerical studies of the flow of polymeric material through a thin film enable visualisation of the pressure, velocity and stress distributions within the channel necessary for calculating the load-carrying capacity, flow rate and friction, respectively. A full-scale mathematical model comprises the mass and momentum conservation laws coupled with a valid constitutive relation for the polymer extra stress tensor,
$\boldsymbol{\tau }^*$
, for example, the Oldroyd-B model approximating an elastic material as an ideal Hookean spring element, being able to stretch several orders of magnitude in comparison with its equilibrium length. When the equations are solved via direct numerical simulations (DNS), it was discovered that a critical numerical instability prevents convergence for large values of
$\textit{De}$
or
$\textit{Wi}$
, termed the high-Weissenberg-number problem (HWNP) (Keunings Reference Keunings1986; Owens & Phillips Reference Owens and Phillips2002). Its manifestation was traced to a loss of symmetric positive-definite property of the conformation tensor,
$\unicode{x1D63E}$
, related to the polymeric stress via the relation
$\boldsymbol{\tau }^* = \eta _p / \lambda (\unicode{x1D63E} - \unicode{x1D644}\, )$
, as
$\textit{De}$
increases beyond a certain value. This problem improves, only temporarily, when the maximum possible extension of the polymer chains is restricted as in FENE-CR, and further by shear thinning, e.g. the FENE-P model (Alves, Oliveira & Pinho Reference Alves, Oliveira and Pinho2021; Zhang et al. Reference Zhang, Zhang, Wang, Li, Li and Li2023).
Extensive efforts in devising novel discretisation methods did not yield significant improvements as the degeneracy of the eigenvalues of
$\unicode{x1D63E}$
was traced to the use of polynomial-based approximations for the stress or conformation tensor. Specifically, beyond a certain critical
$\textit{De}$
(or
$\textit{Wi}$
), negative eigenvalues tend to emerge which are unphysical (Chakraborty et al. Reference Chakraborty, Bajaj, Yeo, Friend, Pasquali and Prakash2010). This was resolved to a great extent in the seminal work of Fattal & Kupferman (Reference Fattal and Kupferman2004) which employed a logarithmic transformation of
$\unicode{x1D63E}$
(the so-called logarithmic conformation representation) and its success was immediately demonstrated on some classical problems (Fattal & Kupferman Reference Fattal and Kupferman2005), later extended to more complicated and numerically challenging cases, such as the 4 : 1 contraction problem (Alves et al. Reference Alves, Oliveira and Pinho2021). Still, for some critical value of
$\textit{De}$
, simulation accuracy becomes questionable and eventually fails. The HWNP only having been delayed prompted further research and recent efforts have shown that it can be triggered also due to (i) a change in type of the numerical system (Renardy & Thomases Reference Renardy and Thomases2021) or (ii) improper interpolation of the eigenvalues of
$\unicode{x1D63E}$
(Zhang et al. Reference Zhang, Zhang, Wang, Li, Li and Li2023), but dampened via artificial diffusion methods (Fernandes et al. Reference Fernandes, Araujo, Ferrás and Nóbrega2017). Interestingly, for thin-film lubrication problems, the HWNP is not alleviated as the grid is refined by including more nodes/elements, contrary to earlier findings that suggest the opposite. This was detected very early on for simulations of journal bearings (Beris, Armstrong & Brown Reference Beris, Armstrong and Brown1986) in which a saturation of the stability versus mesh refinements was observed.
For thin-film boundary-driven lubrication, the use of DNS for predicting the pressure profile along the channel via a computationally efficient implementation of the logarithmic conformation representation in RheoTool (Alves et al. Reference Alves, Oliveira and Pinho2021) was limited to small values of
$\textit{De}$
(Ahmed & Biancofiore Reference Ahmed and Biancofiore2021, Reference Ahmed and Biancofiore2023). When combined with other nonlinear effects, such as film cavitation, the critical
$\textit{De}$
value was further reduced (Gamaniel et al. Reference Gamaniel, Dini and Biancofiore2021). These difficulties are alleviated via simplifications resulting from the application of lubrication theory (Szeri Reference Szeri2010) on the governing set of equations. This effort yields a set of reduced equations that can be solved to give the classical Reynolds equation for the Newtonian case. However, it does not simplify the constitutive relations. Efforts have therefore been directed at methods which attempt to further reduce the nonlinear coupling between pressure, velocity and polymer extra stress components. In this regard, the first major attempt was the generalised Reynolds equation (Johnson & Tevaarwerk Reference Johnson and Tevaarwerk1977), celebrated for its capacity in capturing cross-film variation of lubricant properties (shear thinning, thermal effects, liquid compressibility, etc.), modified to include the effects of elasticity. However, the problem was approached in a Newtonian-like manner only considering the effects of the non-Newtonian shear stress and ignoring the normal stress (since, for the Newtonian case, the latter is of the order of
$\epsilon ^2$
) (Wolff & Kubo Reference Wolff and Kubo1996). Furthermore, only a portion of the substantive derivative appearing in the upper-convected rate operator was kept for simplicity leading to an over-simplified prediction of the polymeric stress. Owing to these limitations, the elasticity of the lubricant film was not considered significant. Nonetheless, these early efforts, driven by a practical need to resolve the role of viscoelasticity in mechanical components, highlight the complexity of the problem.
The problem in obtaining a simplified analytical model lies in the nonlinear coupling of the velocity field (and its gradients) with the polymer extra stress tensor. Noting that the Deborah number (i) appears naturally in the constitutive relations and (ii) can be argued to be small for some cases, a perturbation approach in
$\textit{De}$
was applied which yielded an equation for the
$\textit{De}$
-order pressure depending only on the channel surface function,
$h$
, for an upper-convected Maxwell fluid (Tichy Reference Tichy1996) and also for a second-order fluid (Sawyer & Tichy Reference Sawyer and Tichy1998). An approach based on the reciprocal theorem yielded a third-order expansion in
$\textit{De}$
for purely pressure-driven flows for the Oldroyd-B model (Boyko & Stone Reference Boyko and Stone2022), comparing favourably with DNS over a wide range. Similarly, we find other studies for pressure-driven flows that leverage the perturbation analysis and demonstrate its effectiveness in predicting the relevant quantities, such as flow rate, pressure drop across the channel, etc. (Housiadas & Beris Reference Housiadas and Beris2024a
,
Reference Housiadas and Berisb
). These models generally express the desired quantity (usually pressure or flow rate) as a function of the surface height but a recent notable work by Boyko, Hinch & Stone (Reference Boyko, Hinch and Stone2024) achieves near-analytical expressions for the pressure and stress by (i) leveraging a curvilinear coordinate transformation and (ii) noting the velocity field is Newtonian when the polymer concentration is small. These assumptions allow for the elimination of the cross-film velocity component (in curvilinear space). In a related work, a numerical treatment of the thin-film reduced system reveals the stress relaxation mechanism and how it contributes to a pressure drop for contracting channels (Hinch et al. Reference Hinch, Boyko and Stone2024). More recently, an analytical solution for moderate values of the Deborah number was proposed for axisymmetric configurations (Sialmas & Housiadas Reference Sialmas and Housiadas2025). These simplified numerical models demonstrate the strongly contrasting behaviour of a viscoelastic fluid mainly owing to the generation of streamwise-normal stresses (depending quadratically on the applied shear strain rate for small
$\textit{De}$
) which are O(
$\epsilon ^2$
) for a Newtonian fluid, and are therefore negligible in a thin film. Furthermore, they reveal how the surface profile and its gradients (converging or diverging gaps) determine whether the extra stress will enhance or inhibit the pressure for problems neglecting spanwise effects.
The same measure of sensitivity in increasing viscoelasticity is not observed for the velocity field which is only weakly perturbed (Zheng et al. Reference Zheng, Xie, Chen, Liu, Yang, Xu and Huang2023). In certain special cases, the viscoelastic influence on the kinematics vanishes entirely and the flow field remains Newtonian (Phan-Thien & Tanner Reference Phan-Thien and Tanner1983; Phan-Thien & Tanner Reference Phan-Thien and Tanner1984; Phan-Thien et al. Reference Phan-Thien, Dudek, Boger and Tirtaatmadja1985) (see the Tanner and Pipkin theorem for a two-dimensional Cartesian configuration). In fact, for two-dimensional channels (representing the cross-section of an infinitely wide real channel) the elastic component of the velocity field vanishes when the channel inlet and outlet heights are equal and, also, when the polymeric viscosity,
$\eta _p$
, is small in comparison with that of the solvent,
$\eta _s$
(the ultra-dilute limit). Exploiting the weak dependence of
$\boldsymbol {u}$
on the Deborah number, it is possible to decouple the momentum equation from the constitutive relation, and obtain a Reynolds-type equation for the pressure, the so-called viscoelastic Reynolds (VR) approach (Ahmed & Biancofiore Reference Ahmed and Biancofiore2021). For sliding lubricated contacts, a straightforward numerical treatment shows a better prediction to the classical first-order perturbation versus the Deborah number and compares favourably with DNS. It predicts the emergence of a nonlinear trend in the load-carrying capacity versus
$\textit{De}$
, demonstrating not only the positive contribution of the polymer normal stress but a mechanism that differs significantly from the Newtonian case in which the primary contribution is only due to viscous effects (Tichy & Bou-Saïd Reference Tichy and Bou-Saïd2008).
The studies mentioned thus far focus on a two-dimensional approximation of real lubricated channels either arising from symmetry arguments or being infinite across the third dimension. However, practical problems involving the thin-film lubrication of two sliding surfaces deal with channels that have a finite spanwise width with respect to the length. These channels are arguably three-dimensional when the length scales are of the same order of magnitude, leading to a larger set of scalar-differential equations to be solved for the viscoelastic case. It is evident from the outset that the lack of analytical solutions even for the Newtonian case significantly increases the overall complexity of the problem. This necessitates the use of a numerical procedure and simplified reduced-order models.
Our aim is to understand the influence of a finite spanwise width on the polymeric response. To do so, we consider different length-to-width ratios to mimic different types of sliding channels, depicted in figure 1. (i) the contact region in cylindrical roller bearing elements has a constant cross-section and is very wide in comparison with its length (figure 1
a), (ii) ball bearings have equal dimensions at the contact point (figure 1
b), (iii) while the contact region in non-textured journal bearings that have no surface gradients along the width and the contact region is long in comparison with its width (figure 1
c), in contrast to surface textures like pockets or surface roughness. There are, therefore, two geometric considerations when dealing with three-dimensional problems, namely the surface gradients and the channel aspect ratio (length-to-width ratio). This paper extends the VR approach as a potential reduced-order modelling technique for the flow of a viscoelastic lubricant through a thin finite-width channel examining mainly the load-carrying capacity per unit width versus
$\textit{De}$
and the channel aspect ratio. In addition, we compare the results against a model linearised in
$\textit{De}$
to shed light on its usefulness for the three-dimensional case. Furthermore, we attempt to explain the relation between the observed trends in the load and the polymer stress distribution.
Figure 1. A schematic diagram indicating the characteristic streamwise (spanwise) length in blue (red) of (a) a cylindrical roller-element bearing (
$\ell _z \gt \ell _x$
for the rollers where the cross-section is denoted by the blue circle), (b) a ball bearing (
$\ell _z \approx \ell _x$
) and (c) a journal bearing (
$\ell _z \lt \ell _x$
). The green arrow is an indication of the direction of rotation. For the journal bearing (c), the contact area is the blue circle. The images were generated using Adobe Firefly.
Figure 2. A schematic diagram of a sliding lubricating channel with a stationary upper surface (red) and a flat lower surface (blue) moving at constant speed.
2. Problem formulation
In this section, we present the mathematical problem for the boundary-driven flow of a viscoelastic lubricant through a three-dimensional channel. In § 2.1, the governing equations along with the relevant dimensionless parameters characterising the flow are presented. In § 2.2, these equations are then simplified in two ways: (i) using classical perturbation technique and (ii) extending the VR method (Ahmed & Biancofiore Reference Ahmed and Biancofiore2021).
2.1. Governing equations
The flow of a viscoelastic polymer solution comprising a Newtonian base solvent and a dilute concentration of elastic polymeric additives through a channel with impermeable walls, depicted in figure 2, is mathematically described via the system of equations
\begin{align} &\qquad\qquad\qquad \boldsymbol{\nabla }^{*} \boldsymbol{\cdot }\boldsymbol{u}^{*} = 0, \end{align}
\begin{align}& \rho _0\frac {{\rm D} \boldsymbol{u}^{*}}{{\rm D} t^{*}} = -\boldsymbol{\nabla }^{*} p^{*} + \eta _s \boldsymbol{\nabla }^{*2}\boldsymbol{u}^{*} + \boldsymbol{\nabla }^{*} \boldsymbol{\cdot }\boldsymbol{\tau }^{*}, \end{align}
\begin{align}&\qquad\qquad\,\,\, \boldsymbol{\tau }^* + \lambda {\mathop {\boldsymbol{\tau}} \limits^{\triangledown}}{}^{*} = \eta _p \unicode{x1D63F}^{\kern1pt *}, \end{align}
where
$\rho ^{*}$
is the fluid density,
$\eta _s$
is the solvent viscosity,
$\eta _p$
is the polymer viscosity,
$\lambda$
is the polymer relaxation time,
$\boldsymbol {u}^{*}$
is the velocity field,
$t^{*}$
is the time,
$p^{*}$
is the fluid pressure,
$\boldsymbol{\tau }^{*}$
is the polymer extra stress tensor,
$\unicode{x1D63F}^{\kern1pt *} = \boldsymbol{\nabla }^{*}\boldsymbol{u}^{*} + \boldsymbol{\nabla }^{*}\boldsymbol{u}^{*T}$
is the deformation rate tensor,
$\unicode{x1D647} = \boldsymbol{\nabla }^{*}\boldsymbol{u}^{*}$
is the velocity gradient tensor and
${\mathop {\boldsymbol{\tau}} \limits^{\triangledown}}{}^{*} = ({{\rm D} \boldsymbol{\tau }^*}/{{\rm D}t^*}) -\unicode{x1D647}^*\boldsymbol{\tau }^* - \boldsymbol{\tau }^* \unicode{x1D647}^{*T}$
is the upper convected rate. Equation (2.1c
) is the Oldroyd-B constitutive relation (Bird, Armstrong & Hassager Reference Bird, Armstrong and Hassager1987).
Using the length and velocity scaling
\begin{align} x^* & = x \ell _x, \quad y^* = y H_0,\quad z^* = z \ell _z, \end{align}
\begin{align} u^* & = u U_0,\quad v^* = v \epsilon U_0,\quad w^* = w U_0, \end{align}
where
$\epsilon = {H_0}/{\ell _x}$
is the thin-film ratio, and a stress scaling (Tichy Reference Tichy1996; Li Reference Li2014; Ahmed & Biancofiore Reference Ahmed and Biancofiore2021; Boyko & Stone Reference Boyko and Stone2022; Housiadas & Beris Reference Housiadas and Beris2023),
\begin{align} & \tau _{xx} = \tau _{xx}^*\frac { h_0^2}{\eta _0 U_0 \ell _x},\quad \tau _{zx} = \tau _{zx}^*\frac {h_0^2}{\eta _0 U_0 \ell _x},\quad \tau _{xy} = \tau _{xy}^*\frac {h_0}{\eta _0 U_0}, \end{align}
\begin{align} & \tau _{yy} = \tau _{yy}^*\frac { \ell _x}{\eta _0 U_0},\quad \tau _{yz} = \tau _{yz}^*\frac { h_0}{\eta _0 U_0},\quad \tau _{zz} = \tau _{zz}^*\frac { h_0^2}{\eta _0 U_0 \ell _x}, \end{align}
in which
$\eta _0 = \eta _s + \eta _p$
is the reference viscosity and all the normal stress components are of the order of the film pressure, we arrive at a reduced system of equations:
\begin{align} &\qquad\qquad\qquad\qquad\quad\,\, \frac {\partial u }{\partial x} + \frac {\partial v}{\partial y} + a\frac {\partial w}{\partial z}= 0, \end{align}
\begin{align} &\qquad\qquad\qquad\,\, \frac {\partial p}{\partial x} = \beta \frac {\partial ^2 u}{\partial y^2} + \frac {\partial \tau _{xx}}{\partial x} + \frac {\partial \tau _{yx}}{\partial y} + a \frac {\partial \tau _{zx}}{\partial z}, \end{align}
\begin{align} &\qquad\qquad\qquad\qquad\qquad\qquad\,\, \frac {\partial p}{\partial y} = 0, \end{align}
\begin{align} &\qquad\qquad\qquad a\frac {\partial p}{\partial z} = \beta \frac {\partial ^2 w}{\partial y^2} + \frac {\partial \tau _{xz}}{\partial x} + \frac {\partial \tau _{yz}}{\partial y} + a\frac {\partial \tau _{zz}}{\partial z}, \end{align}
\begin{align} &\qquad \tau _{xx} + De \left (\frac {D\tau _{xx}}{Dt} -2\tau _{xy}\frac {\partial u}{\partial y} - 2\tau _{xx}\frac {\partial u}{\partial x} - 2a\frac {\partial u}{\partial z}\tau _{xz} \right ) = 0, \end{align}
\begin{align} \nonumber & \tau _{xy} + De \left (\frac {D\tau _{xy}}{Dt} -\tau _{yy}\frac {\partial u}{\partial y} - \tau _{xx}\frac {\partial v}{\partial x} - \frac {\partial u}{\partial x}\tau _{xy} - \frac {\partial v}{\partial y}\tau _{xy} - a\frac {\partial u}{\partial z}\tau _{yz}\right . \nonumber\\ &\quad\left . -\, a\frac {\partial v}{\partial z}\tau _{xz} \right ) = (1 - \beta )\frac {\partial u}{\partial y}, \end{align}
\begin{align} \nonumber &\qquad \tau _{xz} + De \left (\frac {D\tau _{xz}}{Dt} -\tau _{xz}\frac {\partial u}{\partial x} - \tau _{xx}\frac {\partial w}{\partial x} -\tau _{yz}\frac {\partial u}{\partial y} -\tau _{xy}\frac {\partial w}{\partial y}\right . \nonumber\\ &\qquad\quad\left . -a\tau _{zz}\frac {\partial u}{\partial z} -a\tau _{xz}\frac {\partial w}{\partial z} \right ) = 0, \end{align}
\begin{align} & \tau _{yy} + De \left (\frac {D\tau _{yy}}{Dt} -2\tau _{xy}\frac {\partial v}{\partial x} - 2\tau _{yy}\frac {\partial v}{\partial y} - 2a\frac {\partial v}{\partial z}\tau _{yz} \right ) = 2(1 - \beta )\frac {\partial v}{\partial y}, \end{align}
\begin{align} \quad \tau _{yz} & + De \left (\frac {D\tau _{yz}}{Dt} -\tau _{xz}\frac {\partial v}{\partial x} - \tau _{xy}\frac {\partial w}{\partial x} - \tau _{yz}\frac {\partial v}{\partial y} - \tau _{yy}\frac {\partial w}{\partial y} - a\tau _{zz}\frac {\partial v}{\partial z}\right . \nonumber \\ & \quad \left . -\, a\tau _{yz}\frac {\partial w}{\partial z}\right ) = (1 - \beta )\frac {\partial w}{\partial y}, \end{align}
\begin{align} & \tau _{zz} + De \left (\frac {D\tau _{yy}}{Dt} - 2\tau _{xz}\frac {\partial w}{\partial x} - 2\tau _{xy}\frac {\partial w}{\partial y} - 2a\frac {\partial w}{\partial z}\tau _{zz} \right ) = 0, \end{align}
where the Deborah number
$\textit{De}$
is defined as
$De=({\lambda U_0}/{\ell _x})$
,
$\beta$
is the solvent concentration in terms of viscosity (
$\beta = {\eta _s}/{\eta _0}$
),
$a = {\ell _x}/{\ell _z}$
is the channel aspect ratio (length-to-width ratio) and
$( {{\rm D}}/{{\rm D}t}) = u({\partial }/{\partial x}) + v( {\partial }/{\partial y}) + aw({\partial }/{\partial z})$
is the material derivative operator. Inherent in the derivation is the assumption of a vanishingly small Reynolds number (
$ \textit{Re} = ({\rho _0 U_0 H_0}/{\eta _0}) \ll 0$
). The system of (2.4) is similar to that obtained in Boyko & Christov (Reference Boyko and Christov2023) but with a different definition of the channel aspect ratio,
$a$
, that allows retrieving the two-dimensional case when applying the limit
$a\rightarrow 0$
in (2.4). In addition, values of
$a$
of order
$\epsilon ^{-1}$
cannot be used, since terms of the order of
$\epsilon ^2a^2$
have also been ignored when applying lubrication theory. This is true for both the Newtonian and the viscoelastic cases. A summary of the dimensionless parameters is given in table 1, including two additional parameters, the spanwise Deborah number
$ \textit{De}_s$
and the leakage Deborah number
$De_a$
, that are defined along the length of the curved region in the spanwise direction and the total spanwise length, respectively. These additional quantities help in explaining the effects due to viscoelasticity along the spanwise direction. A detailed discussion is given in § 3.3.2.
Table 1. The dimensionless parameters arising from the rescaling of the scalar governing equations. For the definition of
$W_s$
and
$W_a$
, see § 3.3.2.
The nonlinear coupled system yields analytical solutions in certain special cases, such as a constant shear flow (Oliveira Reference Oliveira2002). However, for most practical purposes it must be solved numerically which is challenging owing to the nonlinear coupling between the flow field and the polymer stress (intensifying with increasing
$\textit{De}$
), and the onset of the HWNP which prevents solutions at large
$\textit{De}$
(or
$\textit{Wi}$
).
For this work, the system of (2.4) and any reduced model obtained thereof are subject to the boundary conditions
\begin{align} &\qquad\qquad\quad u(x, 0, z) = 1,\quad u(x, h, z) = 0, \end{align}
\begin{align} &\qquad\qquad\qquad v(x, 0, z) = v(x, h, z) = 0, \end{align}
\begin{align} &\qquad\qquad\quad\,\,\,\, w(x, 0, z) = w(x, h, z) = 0, \end{align}
\begin{align} &\,\, p(0, z) = p(1,z) = p(x,1)= 0,\mbox{ and } \frac {\partial p}{\partial z}(x,0) = 0, \end{align}
\begin{align} & \frac {\partial \boldsymbol{\tau }}{\partial x}(0,y,z) = 0,\quad \frac {\partial \boldsymbol{\tau }}{\partial z}(x,y,0) = 0,\quad \frac {\partial \boldsymbol{\tau }}{\partial y}(x,0,z) = 0, \end{align}
which describe a configuration in which the lower channel wall (
$y=0$
) is sliding at a fixed speed and the upper wall (at
$y=h(x,z)$
) is stationary. We assume symmetry along the spanwise axis (
$z=0$
), and therefore employ a zero-gradient condition. These conditions imply gauge pressure at the open boundaries. See Tichy (Reference Tichy1996), Sawyer & Tichy (Reference Sawyer and Tichy1998) and Ahmed & Biancofiore (Reference Ahmed and Biancofiore2021) for a discussion, and a brief analysis covering the deviation from this condition is presented in § 3.4. Particularly, in three-dimensional contacts, assuming zero pressure at the lateral boundary (i.e.
$z=1$
) is common to model standard mechanical components, such as bearings, pads, etc. (Rastogi & Gupta Reference Rastogi and Gupta1991; Bertocchi et al. Reference Bertocchi, Dini, Giacopini, Fowell and Baldini2013; Çam et al. Reference Çam, Giacopini, Dini and Biancofiore2023). Note that we assume that the environment is pressurised to avoid the presence of cavitation, which is often found in mechanical components at ambient pressure (Dowson & Taylor Reference Dowson and Taylor1979; Gamaniel et al. Reference Gamaniel, Dini and Biancofiore2021; Çam et al. Reference Çam, Giacopini, Dini and Biancofiore2023).
2.2. Reduced-order models
The thin-film reduced system of equations is a nonlinear and fully coupled mathematical problem that does not readily yield any analytical solutions. In this section, we introduce two potential reduced-order models that decouple the momentum equations from the scalar constitutive equations, namely the linearised model (LIN) following a first-order perturbation in
$\textit{De}$
and the VR numerical approach which assumes a weakly perturbed velocity field (extended here to three-dimensional channels (Ahmed & Biancofiore Reference Ahmed and Biancofiore2021)).
2.2.1. Linearised viscoelastic Reynolds equation
The LIN is obtained by substituting the expansion
\begin{equation} \{p, \boldsymbol{u}, \boldsymbol{\tau }\} = \{p, \boldsymbol{u}, \boldsymbol{\tau }\}^{(0)} + De \{p, \boldsymbol{u}, \boldsymbol{\tau }\}^{(1)} + O(De^2) \end{equation}
into (2.8) and neglecting all terms of the order of
$De^2$
, where
$\{\boldsymbol{\cdot }\}^{(0)}$
and
$\{\boldsymbol{\cdot }\}^{(1)}$
denote the leading order and
$\textit{De}$
order, respectively. Collecting the relevant terms and simplifying, we find the leading-order or Newtonian system of equations:
\begin{align} &\qquad\quad \frac {\partial u^{(0)} }{\partial x} + \frac {\partial v^{(0)} }{\partial y} + \frac {\partial w^{(0)} }{\partial z}= 0, \end{align}
\begin{align}& \frac {\partial p^{(0)} }{\partial x} = \beta \frac {\partial ^2 u^{(0)}}{\partial y^2} + \frac {\partial \tau _{xx}^{(0)} }{\partial x} + \frac {\partial \tau _{yx}^{(0)} }{\partial y} + a \frac {\partial \tau _{zx}^{(0)} }{\partial z}, \end{align}
\begin{align} &\qquad\qquad\qquad\quad \frac {\partial p^{(0)}}{\partial y} = 0, \end{align}
\begin{align} & a\frac {\partial p^{(0)}}{\partial z} = \beta \frac {\partial ^2 w^{(0)} }{\partial y^2} + \frac {\partial \tau _{xz}^{(0)} }{\partial x} + \frac {\partial \tau _{yz}^{(0)} }{\partial y} + a\frac {\partial \tau _{zz}^{(0)} }{\partial z}, \end{align}
\begin{align} &\qquad\qquad\qquad\quad\,\, \tau _{xx}^{(0)} = 0, \end{align}
\begin{align} &\qquad\qquad\quad \tau _{xy}^{(0)} = (1 - \beta )\frac {\partial u^{(0)} }{\partial y}, \end{align}
\begin{align} &\qquad\qquad\qquad\quad\,\, \tau _{xz}^{(0)} = 0, \end{align}
\begin{align} &\qquad\qquad\quad \tau _{yy}^{(0)} = 2(1 - \beta )\frac {\partial v^{(0)} }{\partial y}, \end{align}
\begin{align} &\qquad\qquad\quad\, \tau _{yz}^{(0)} = (1 - \beta )\frac {\partial w^{(0)} }{\partial y}, \end{align}
\begin{align} &\qquad\qquad\qquad\qquad \tau _{zz}^{(0)} = 0, \\[6pt] \nonumber \end{align}
and the
$\textit{De}$
-order system of equations:
\begin{align} &\qquad\qquad\qquad\qquad\qquad\quad\,\, \frac {\partial u^{(1)} }{\partial x} + \frac {\partial v^{(1)} }{\partial y} + a\frac {\partial w^{(1)}}{\partial z}= 0, \end{align}
\begin{align} &\qquad\qquad\qquad\qquad \frac {\partial p^{(1)}}{\partial x} = \beta \frac {\partial ^2 u^{(1)} }{\partial y^2} + \frac {\partial \tau _{xx}^{(1)} }{\partial x} + \frac {\partial \tau _{yx}^{(1)} }{\partial y} + a \frac {\partial \tau _{zx}^{(1)} }{\partial z}, \end{align}
\begin{align} &\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad \frac {\partial p^{(1)} }{\partial y} = 0, \end{align}
\begin{align} &\qquad\qquad\qquad\quad a\frac {\partial p^{(1)} }{\partial z} = \beta \frac {\partial ^2 w^{(1)} }{\partial y^2} + \frac {\partial \tau _{xz}^{(1)} }{\partial x} + \frac {\partial \tau _{yz}^{(1)} }{\partial y} + a\frac {\partial \tau _{zz}^{(1)} }{\partial z }, \end{align}
\begin{align} &\qquad\quad \tau _{xx}^{(1)} = -De \left (\frac {{\rm D}\tau _{xx}^{(0)}}{{\rm D}t} -2\tau _{xy}^{(0)}\frac {\partial u^{(0)} }{\partial y} - 2\tau _{xx}^{(0)}\frac {\partial u^{(0)} }{\partial x} - 2a\frac {\partial u^{(0)} }{\partial z}\tau _{xz}^{(0)} \right )\!, \end{align}
\begin{align} \tau _{xy}^{(1)} & = -De \left (\frac {{\rm D}\tau _{xy}^{(0)}}{{\rm D}t} -\tau _{yy}^{(0)}\frac {\partial u^{(0)} }{\partial y} - \tau _{xx}^{(0)}\frac {\partial v^{(0)}}{\partial x} - \frac {\partial u^{(0)}}{\partial x}\tau _{xy}^{(0)} - \frac {\partial v^{(0)} }{\partial y}\tau _{xy}^{(0)} - a\frac {\partial u^{(0)}}{\partial z}\tau _{yz}^{(0)}\right . \nonumber \\ &\qquad\,\left . -\, a\frac {\partial v^{(0)}}{\partial z}\tau _{xz}^{(0)} \right ) + (1 - \beta )\frac {\partial u^{(1)} }{\partial y}, \end{align}
\begin{align} \nonumber &\qquad \tau _{xz}^{(1)} = -De \left (\frac {{\rm D}\tau _{xz}^{(0)} }{{\rm D}t} -\tau _{xz}^{(0)}\frac {\partial u^{(0)}}{\partial x} -\tau _{xx}^{(0)}\frac {\partial w^{(0)}}{\partial x} -\tau _{yz}^{(0)}\frac {\partial u^{(0)}}{\partial y} -\tau _{xy}^{(0)}\frac {\partial w^{(0)}}{\partial y} \right . \\ &\qquad\qquad\,\left . -\, a\tau _{zz}^{(0)}\frac {\partial u^{(0)}}{\partial z} -a\tau _{xz}^{(0)}\frac {\partial w^{(0)}}{\partial z} \right )\!, \end{align}
\begin{align} & \tau _{yy}^{(1)} = -De \left (\frac {{\rm D}\tau _{yy}}{{\rm D}t} -2\tau _{xy}^{(0)}\frac {\partial v^{(0)}}{\partial x} - 2\tau _{yy}^{(0)}\frac {\partial v^{(0)}}{\partial y} - 2a\frac {\partial v^{(0)}}{\partial z}\tau _{yz}^{(0)} \right ) + 2(1 - \beta )\frac {\partial v^{(1)}}{\partial y}, \end{align}
\begin{align} \tau _{yz}^{(1)} & = -De \left (\frac {{\rm D}\tau _{yz}^{(0)}}{{\rm D}t} -\tau _{xz}^{(0)}\frac {\partial v^{(0)}}{\partial x} - \tau _{xy}^{(0)}\frac {\partial w^{(0)}}{\partial x} - \tau _{yz}^{(0)}\frac {\partial v^{(0)}}{\partial y} - \tau _{yy}^{(0)}\frac {\partial w^{(0)}}{\partial y} - a\tau _{zz}^{(0)}\frac {\partial v^{(0)}}{\partial z} \right . \nonumber \\ &\quad\,\left . -\, a\tau _{yz}^{(0)}\frac {\partial w^{(0)}}{\partial z}\right ) + (1 - \beta )\frac {\partial w^{(1)}}{\partial y}, \end{align}
\begin{align} &\qquad\quad \tau _{zz}^{(1)} = -De \left (\frac {{\rm D}\tau _{zz}^{(0)}}{{\rm D}t} - 2\tau _{xz}^{(0)}\frac {\partial w}{\partial x} - 2\tau _{xy}^{(0)}\frac {\partial w}{\partial y} - 2a\frac {\partial w^{(0)}}{\partial z}\tau _{zz}^{(0)} \right)\!. \end{align}
The solution to the leading-order system (subject to the conditions (2.5)) gives the Newtonian pressure and the velocity field (Szeri Reference Szeri2010):
\begin{align} &\qquad\qquad\qquad \frac {\partial }{ \partial x}\left ( \frac {h^3}{12}\frac {\partial p^{(0)}}{\partial x} \right ) +a^2 \frac {\partial }{ \partial z}\left ( \frac {h^3}{12}\frac {\partial p^{(0)}}{\partial z} \right ) = \frac {1}{2}\frac {\partial h}{\partial x}, \end{align}
\begin{align} &\qquad\qquad\qquad\quad u^{(0)} = \frac {1}{2}\frac {\partial p^{(0)}}{\partial x}\left(y^2 - yh\right) + \left(1 - \frac {y}{h}\right)\!, \end{align}
\begin{align} &v^{(0)} = -\frac {y^3}{6}\left ( \frac {\partial ^2 p^{(0)}}{\partial x^2} + a^2 \frac {\partial ^2 p^{(0)}}{\partial z^2} \right ) + \frac {y^2}{2}\left ( \frac {\partial ^2 p^{(0)}}{\partial x^2}\frac {h}{2} + a^2 \frac {\partial ^2 p^{(0)}}{\partial z^2}\frac {h}{2} + \frac {1}{2}\frac {\partial p^{(0)}}{\partial x}\frac {\partial h}{\partial x}\right . \nonumber\\&\qquad\,\,\left . +\, \frac {a^2}{2}\frac {\partial p^{(0)}}{\partial z}\frac {\partial h}{\partial z}\right )\!, \end{align}
\begin{align} &\qquad\qquad\qquad\qquad\qquad w^{(0)} = \frac {a}{2}\frac {\partial p^{(0)}}{\partial z}(y^2 - yh). \end{align}
It is important to note that (2.9) do not yield an analytical solution for even the simplest surface profiles. The matter is made worse when factoring in discontinuities along the surface, such as a textured pocket (Schuh et al. Reference Schuh, Lee, Allison and Ewoldt2017; Çam et al. Reference Çam, Giacopini, Dini and Biancofiore2023). There are two limiting cases for (2.9a
): (i)
$a \ll 1$
, describing an infinitely wide channel, and (ii)
$a \gg 1$
, which models a slender channel.
The solution to the
$\textit{De}$
-order system requires additional effort owing to the increase in the number of equations and terms thereof. The final system is a lengthy set of equations for the pressure and velocity components. For the sake of brevity, we present only the
$\textit{De}$
-order pressure equation:
\begin{align} &\qquad\qquad\qquad \frac {\partial }{\partial x}\left ( \frac {h^3}{12}\frac {\partial p^{(1)}}{\partial x}\right ) + a^2 \frac {\partial }{\partial z}\left ( \frac {h^3}{12}\frac {\partial p^{(1)}}{\partial z}\right ) = \mathcal{F} + \mathcal{G}_1 + \mathcal{G}_2, \end{align}
\begin{align} \mathcal{F} & = (1-\beta ) \left ( \frac {h^5}{48} \left ( \frac {\partial p^{(0)}}{\partial x} \frac {\partial ^3 p^{(0)}}{\partial x^3} + \frac {\partial ^2 p^{(0)}}{\partial x^2} \right ) + \frac {7 h^4}{48} \frac {\partial h}{\partial x} \frac {\partial p^{(0)}}{\partial x} \frac {\partial ^2 p^{(0)}}{\partial x^2}+\frac {h^4}{48} \frac {\partial ^2 h}{\partial x^2} \frac {\partial p^{(0)}}{\partial x} \right .\nonumber\\&\quad\left . \times \left ( \frac {\partial p^{(0)}}{\partial x}\right )^2 + \frac {h^3}{12} \left ( \frac {\partial h}{\partial x} \right )^2 \left ( \frac {\partial p^{(0)}}{\partial x} \right )^2 - \frac {h^3}{24} \frac {\partial ^3 p^{(0)}}{\partial x^3} - \frac {h^2}{8} \frac {\partial h}{\partial x} \frac {\partial ^2 p^{(0)}}{\partial x^2} - \frac {1}{12} \frac {\partial ^2 h}{\partial x^2}\right )\!, \end{align}
\begin{align} & \!\! \mathcal{G}_1 = a^{2} (1-\beta ) \left ( \frac {h^5}{48} \frac {\partial p^{(0)}}{\partial x} \frac {\partial ^3 p^{(0)}}{\partial z^2 \partial x} + \frac {h^5}{40} \frac {\partial ^2 p^{(0)}}{\partial x^2} \frac {\partial ^2 p^{(0)}}{\partial z^2} + \frac {h^5}{48} \frac {\partial p^{(0)}}{\partial z} \frac {\partial ^3 p^{(0)}}{\partial x^2 \partial z}\right .\nonumber \\ & \!\!\!\! +\frac {h^5}{60} \left ( \frac {\partial ^2 p^{(0)}}{\partial x \partial z}\right )^2 + \frac {h^4}{12} \frac {\partial h}{\partial x} \frac {\partial p^{(0)}}{\partial x} \frac {\partial ^2 p^{(0)}}{\partial z^2} + \frac {h^4}{16} \frac {\partial h}{\partial x} \frac {\partial p^{(0)}}{\partial z} \frac {\partial ^2 p^{(0)}}{\partial x\partial z } - \frac {h^3}{24} \frac {\partial h}{\partial x} \frac {\partial ^3 p^{(0)}}{\partial z^2 \partial x} \nonumber \\ & \!\!\!\! \left . +\, \frac {h^4}{16} \frac {\partial h}{\partial z} \frac {\partial p^{(0)}}{\partial x} \frac {\partial ^2 p^{(0)}}{\partial x \partial z} + \frac {h^4}{12} \frac {\partial h}{\partial z} \frac {\partial ^2 p^{(0)}}{\partial x^2} \frac {\partial p^{(0)}}{\partial z} + \frac {h^3}{6} \frac {\partial h}{\partial x} \frac {\partial h}{\partial z} \frac {\partial p^{(0)}}{\partial x} \frac {\partial p^{(0)}}{\partial z} - \frac {h^2}{8}\frac {\partial p^{(0)}}{\partial z} \frac {\partial ^2 p^{(0)}}{\partial x \partial z} \! \right )\!, \nonumber \\ & \mathcal{G}_2 = a^4 (1-\beta ) \left ( \frac {1}{48} h^{5} \frac {\partial p^{(0)}}{\partial z} \frac {\partial ^{3} {p^{(0)}}}{\partial z^{3}} + \frac {1}{48} h^{5} \left(\frac {\partial ^{2} {p^{(0)}} }{\partial z^{2}} \right)^{2} + \frac {7}{48} h^{4} \frac {\partial h}{\partial z} \frac {\partial {p^{(0)}}}{\partial z} \frac {\partial ^{2} {p^{(0)}}}{\partial z^{2}}\right . \nonumber \\&\qquad\left . +\,\frac {1}{48} h^{4} \frac {\partial ^{2} h}{\partial z^{2}} \left(\frac {\partial p^{(0)}}{\partial z} \right)^{2} + \frac {1}{12} h^{3} \left(\frac {\partial h}{\partial z} \right)^{2} \left(\frac {\partial {p^{(0)}}}{\partial z} \right)^{2} \right )\!, \end{align}
where the function
$\mathcal{F}$
represents the streamwise forcing and
$\mathcal{G}_1$
(
$\mathcal{G}_2$
) represents the order-
$a^2$
(order-
$a^4$
) spanwise forcing functions. We define also
$\mathcal{G}$
as the sum of
$\mathcal{G}_1$
and
$\mathcal{G}_2$
,
$\mathcal{G}=\mathcal{G}_1+\mathcal{G}_2$
. The expressions in both the systems, i.e. the leading order and the
$\textit{De}$
order, involve higher powers of
$a$
. Furthermore, third-order and mixed partial derivatives of the leading-order pressure emerge. If
$a \rightarrow 0$
, then the linearised system of equations for a two-dimensional problem are easily retrieved.
2.2.2. Viscoelastic Reynolds equation
An alternative numerical approach, the VR method, significantly alleviates these issues. In this numerical treatment, we enforce the perturbation only in the velocity field, which implies a weak dependence of the flow kinematics due the presence of the polymers (Tanner Reference Tanner1969). This has been reported for several cases involving thin-film lubrication. Hence, we assume
$\boldsymbol{u} = \boldsymbol{u}^{(0)} + De \boldsymbol{u}^{(1)} + O(De^2)$
. This assumption implies that the velocity field responds linearly to increasing viscoelasticity (as
$\textit{De}$
increases), whereas the stress components, in particular along the streamwise- and spanwise-normal directions, can vary strongly and naturally, following their inherent nonlinearity. More importantly, with VR method we decouple the momentum equation from the stress constitutive relation, since the velocity field is known a priori from LIN. The Reynolds-type equation including the polymeric stress components is obtained (see Appendix A):
\begin{align} &\qquad\qquad\qquad\frac {\partial }{\partial x} \left ( \frac {H^3}{12} \frac {\partial p}{\partial x} \right ) + a^2 \frac {\partial }{\partial z} \left ( \frac {H^3}{12} \frac {\partial p}{\partial z} \right ) = \frac {\partial \mathcal{F}}{\partial x} + \frac {\partial \mathcal{G}}{\partial z} , \\[-12pt] \nonumber \end{align}
\begin{align} & \mathcal{F} = \frac {\beta }{2}H + \frac {H^3}{2}\overline {\left[(Y-Y^2)\frac {\partial \tau _{xx}}{\partial x}\right]} + \frac {H^2}{2}\overline {\left[(2Y-1)\tau _{xy}\right]} + a\frac {H^3}{2}\overline {\left[(Y-Y^2)\frac {\partial \tau _{xz}}{\partial z}\right]}, \\[-12pt] \nonumber \end{align}
\begin{align} &\quad \mathcal{G} = a^2\frac {H^3}{2}\overline {\left[(Y-Y^2)\frac {\partial \tau _{zz}}{\partial z}\right]} + a\frac {H^2}{2}\overline {\left[(2Y-1)\tau _{zy}\right]} + a\frac {H^3}{2}\overline {\left[(Y-Y^2)\frac {\partial \tau _{xz}}{\partial x}\right]}, \\[6pt] \nonumber \end{align}
where
$\overline {[\ldots ]} = \int _{0}^{1} [\ldots ] {\rm d}Y$
and
$Y = {y}/{h(x,z)}$
. For an extremely wide channel,
$a \ll 1$
, (2.11) reduce to the two-dimensional case. In addition, if
$h(0, z) = h(1, z)$
, then
$\boldsymbol{u}^{(1)} = 0$
and, based on the expressions obtained from the application of the reciprocal theorem (Boyko & Stone Reference Boyko and Stone2021), the higher-order expansions for
$\textit{De}$
(for a two-dimensional pressure-driven flow problem) also vanish, greatly simplifying the ensuing numerical procedure.
We resort to a numerical treatment via the finite-difference method (FDM) of the leading-order system (2.8), the
$\textit{De}$
-order system (2.10) and the VR equations (2.11) subject to the boundary conditions given by (2.5). Details regarding the numerical procedure are presented in Ahmed & Biancofiore (Reference Ahmed and Biancofiore2023). Briefly, we employ a curvilinear transformation for the velocity given in Appendix B (see Boyko et al. (Reference Boyko, Hinch and Stone2024) for a detailed treatment of all physical quantities), allowing a flux-conserving implicit treatment of the stress advection, presented in Appendix C, and finally the relevant differential operators in the three-dimensional case are summarised in Appendix D. The solution procedure begins by solving the leading-order system (2.8) first, followed by the constitutive relation for the polymeric stress (2.7
e
–
j
) and, finally, retrieving the pressure via the VR equations (2.11). For the stress components, we terminate the iterative process when a tolerance of
$||\tau _{\textit{ij}}^{(k+1)} - \tau _{\textit{ij}}^{(k)}|| \lt 10^{-12}$
is reached, where
$k$
is the iteration level. We report that the eigenvalues of the conformation tensor were checked to ensure that it remains symmetric-positive-definite. A violation of this condition for very large values of
$\textit{De}$
led to convergence failure and is not pursued in this work.
Similarly to the two-dimensional cases (Ahmed & Biancofiore Reference Ahmed and Biancofiore2021, Reference Ahmed and Biancofiore2023), we validate the VR model versus the DNS for the three-dimensional case. An excellent agreement between the two was reached for
$De=0.04$
. Details regarding the numerical validation are presented in Appendix E.
3. Results and discussion
In this section, we examine the influence of viscoelasticity on the load-carrying capacity per unit width
$F_\ell ^* = 2 ({1}/{\ell _z}) \int _{0}^{\ell _x}\int _{0}^{\ell _z} \sigma _{yy}^* \,{\rm d}x^* \,{\rm d}z^*$
which, upon re-scaling using (2.2) and (2.3), reduces to
$F_\ell = \int _{0}^{1}\int _{0}^{1} p \,{\rm d}x \,{\rm d}z$
(since
$\tau _{yy}$
is O(
$\epsilon ^2$
), while
$p$
is O(
$1$
)). The dependency of
$F_\ell$
on the polymer elasticity measured via
$\textit{De}$
and the additional influence of the channel’s spanwise variation measured via the aspect ratio
$a$
are examined using the VR approach, and the
$\textit{De}$
-order perturbed model (LIN). The choice of the channel surface geometry is given in § 3.1, and some considerations of the different boundary conditions for the stress are provided in § 3.4. The solvent concentration in terms of viscosity is set to
$\beta =0.8$
throughout the paper.
3.1. Channel geometry
The success of the numerical solution of the reduced constitutive relations (2.8) and (2.11) is strongly dependent on the channel geometry. For sharply varying surfaces, such as pocketed textures and contractions, convergence of the numerical procedure becomes a challenge (Schuh et al. Reference Schuh, Lee, Allison and Ewoldt2017). However, the main mechanisms involved in a viscoelastic flow can also be studied via a smoothly varying channel surface profile, avoiding several numerical hurdles (Hinch et al. Reference Hinch, Boyko and Stone2024).
We consider two surface configurations: a spanwise-varying surface modelled via a symmetric Gaussian distribution that allows controlling the spread of the curved region and its depth,
\begin{equation} h = 1 - d \exp {\bigg ( -\frac {( x - X_c)^2}{s^2}-\frac {(z - Z_c)^2}{s^2 a^2} \bigg ) }, \end{equation}
where
$d$
is the depth,
$s$
is the spread along the spanwise and streamwise directions and
$X_c$
and
$Z_c$
represent the coordinates of the centre. Note that in this work
$s$
,
$X_c$
and
$Z_c$
are always set to
$s = 0.1$
,
$X_c = 0.5$
and
$Z_c = 0$
. Equation (3.1) represents a surface protrusion that is sensitive to the channel width and has gradients along both the streamwise (sliding) and spanwise directions, and is depicted in figure 3(a–c).
The extruded case approximates channels that are simply an extrusion of the cross-section along the streamwise direction, i.e. the
$(x,y)$
plane, making the surface height
$h$
independent of
$z$
. These types of channels are found in cylindrical roller bearings, slider bearings, line contacts, etc. The surface height for this case is obtained by modifying (3.1):
\begin{equation} h = 1 - d \exp {\bigg ({-}\frac {(x - X_c)^2}{s^2} \bigg ) }, \end{equation}
and the resulting surface profiles are depicted in figure 3(d–f).
Figure 3. The channel surface height variation for three different values of the aspect ratio (
$a = \ell _x / \ell _z$
),
$a = 1/2$
,
$a = 1$
and
$a=2$
, for (a–c) the spanwise-varying (3.1) and (d–f) the extruded (3.2). The channel depth is
$d = 0.2$
, while the spread is
$s = 0.1$
. Note that there is a mass exchange along the spanwise boundary indicated via the red line.
Figure 4. The variation in the load-carrying capacity per unit width, predicted by the
$\textit{De}$
-order LIN model and the VR approach, for (a) the spanwise-varying channel and (b) the extruded surface versus the Deborah number for three different aspect ratios, using
$ d=0.2$
.
3.2. Comparison between LIN and VR approaches
Typically, the pressure will increase when the viscosity of the lubricant is increased, leading to larger shear stresses and, consequently, higher pressures to overcome the increased resistance to flow. As such, a purely Newtonian fluid (
$\textit{De} = 0$
) gives an antisymmetric pressure distribution along
$x$
such that
$\int _{0}^{1}p\,{\rm d}x = 0$
with
$p \gt 0$
(
$p \lt 0$
) for converging (diverging) sections of the channel. Contrary to this, the viscoelastic mixture (for constant viscosity) delivers a net positive load. In figure 4, we compare the load-carrying capacity versus the Deborah number for different channel aspect ratios using
$h=h(x)$
, using
$d=0.2$
and
$\beta = 0.8$
. The simplified models (LIN and VR) predict an expected increase in the net force as we enhance the polymer elastic contribution (Ahmed & Biancofiore Reference Ahmed and Biancofiore2021, Reference Ahmed and Biancofiore2023). Solving for larger
$\textit{De}$
did not change the trend given in figure 4, but increased the numerical complexity, particularly for the limiting cases of large and small
$a$
. Such limitations were not present for the
$\textit{De}$
-order model. However, despite the computational advantage offered by the linear model, a noticeable deviation from the VR model is observed at large
$\textit{De}$
.
Figure 5. The variation in the load-carrying capacity per unit width, predicted by the
$\textit{De}$
-order LIN model and the VR approach, for (a) the spanwise-varying channel and (b) the extruded channel versus the channel depth, using (
$a = \ell _x / \ell _z = 1$
).
The decreased accuracy versus increasing
$\textit{De}$
is not unexpected and has been similarly observed for the equivalent two-dimensional sliding cases. However, in thin-film lubrication, the accuracy of the linear perturbation approach also depends strongly on the channel surface gradients which do not explicitly appear in the Deborah number. We examine the variation in the load versus the channel depth, as shown in figure 5, for the two different channel configurations. The steepening of the channel implies an increase in the average shear rate distribution which directly contributes to an enhancement in the total stretch of the polymer additives. Regardless of whether the surface varies along the spanwise direction, a deviation in the estimated load is observed as the depth increases. In fact, for the extruded channel profile, the difference manifests strongly even for
$\textit{De} = 0.02$
.
In the most general case, the manifestation of viscoelasticity depends strongly on both the Deborah and Weissenberg numbers which reduce to the relation
$\textit{De} = \epsilon Wi$
in the thin-film limit, and therefore the Deborah number typically remains less than one (since
$\epsilon \ll 1$
). Interestingly, the first-order perturbation model in
$\textit{De}$
operates in the linearly viscoelastic limit, but as the channel surface steepens (increasing
$d$
) we observe a discrepancy between the VR and LIN methods even for small magnitudes of
$\textit{De}$
. We attribute this to the increase in the local Weissenberg number:
\begin{equation} {Wi^* = \frac {\lambda U_0}{(1 - d)H_0} = \frac {De}{\epsilon (1 - \textit{d})},} \end{equation}
where the reference film height
$(1 - d)H_0$
is the minimum of the channel height. For small values of
$d$
, the maximum and minimum channel heights are of the same order of magnitude and either choice is suitable. However, as
$d$
increases,
$Wi^*$
rises by nearly an order of magnitude. An increase in
$\textit{Wi}^*$
is a symptom of stronger viscoelastic effects also at low
$\textit{De}$
. Therefore, at large
$d$
we can expect a significant nonlinearity for smaller Deborah numbers with respect to the low-
$d$
case. Alternative scaling approaches do, however, extract
$\textit{Wi}$
but at the expense of the Deborah number. This is the case, for instance, when applying the shear scaling for all the stress components, i.e.
$\tau _{\textit{ij}}^* = ({\eta _0 U_0}/{H_0} )\tau _{\textit{ij}}$
, in the constitutive relation and the momentum equations (Akyildiz & Bellout Reference Akyildiz and Bellout2004; Venkatesh, Anand & Narsimhan Reference Venkatesh, Anand and Narsimhan2022; Abbaspur et al. Reference Abbaspur, Norouzi, Akbarzadeh, Vaziri, Sharghi, Kim and Kim2023).
3.3. Load-carrying capacity
In this section, we examine first the influence of the channel aspect ratio on the load-carrying capacity for both the spanwise-varying and the extruded channel types (see § 3.1). The load-carrying capacity is one of the main quantities for measuring the tribological performance of a lubricant. Another quantity to be checked, generally, is the surface friction arising from the polymer shear stress. However, we did not observe any meaningful dependence on either the aspect ratio or Deborah number. The independence of the friction on the fluid viscoelasticity is in agreement with the findings of previous works (Ahmed & Biancofiore Reference Ahmed and Biancofiore2021; Gamaniel et al. Reference Gamaniel, Dini and Biancofiore2021) and, therefore, not surprising. For the load, we delineate the contributions from the different components of the polymer extra stress tensor, paying particular attention to the dominant streamwise-normal stress component. The analysis is conducted only via the VR approach owing to the inaccuracy of the first-order model for larger
$\textit{De}$
, as demonstrated in § 3.2.
Figure 6. Load variation versus the channel aspect ratio for (a) the spanwise-varying channel and (b) the extruded channel, for three different values of the Deborah number using
$d=0.2$
.
3.3.1. Effect of the aspect ratio
We vary the channel aspect ratio and examine the influence of the polymeric stress on the load-carrying capacity, depicted in figure 6, for
$De=0.05$
,
$\textit{De} = 0.1$
and
$\textit{De} = 0.15$
. Increasing
$\textit{De}$
has the desired beneficial effect of increasing the load; however, the response to varying channel dimensions is quite different between the two configurations. We further observe that the Newtonian response is identically zero for all
$a$
, implying that the extra force is due to purely viscoelastic effects.
For the spanwise-varying case represented in figure 6(a), the load varies strongly versus increasing
$a$
. Close to the two limiting cases, (i) a wide channel (
$a \lt 1$
) and (ii) a slender channel (
$a \gt 1$
), the load is small and rises as we deviate from these limits, reaching a maximum when the width is of the order of the channel length. The variation is essentially governed by two factors: (i) the fraction of curved region and (ii) the fluid leakage along the spanwise boundary. When the width is large, the curved portion covers only a small portion of the total channel wall, which is insufficient to generate strong pressure gradients. As width and length start to be comparable (i.e. around
$a=1$
), the surface gradients increase proportionally generating strong pressure gradients driving the flow. This leads to an enhanced stretching of the polymer chains. However, upon further decrease in the width
$a \gt 1$
, the close proximity of the spanwise boundaries allows the fluid (forcibly drawn in due to the sliding action) to escape, leading to a decline in the pressure.
In the extruded channel we find a starkly contrasting trend of the load versus the channel aspect ratio as can be seen in figure 6(b). The maximum value is observed in the lower limit (
$a \ll 1$
) which tends to the two-dimensional channel. Unlike the spanwise-varying case, surface gradients are absent across the width but are present along the streamwise direction for all
$a$
. Hence, pressure gradients along
$x$
are always present, driving the flow and stretching the polymer additives over a larger portion of the channel, diminishing only as the spanwise boundaries are brought close to the bulk (in the upper limit
$a \gg 1$
). In both cases, the increasing channel aspect ratio promotes side leakage that can lead to a loss of useful polymeric stretch, thus diminishing the load-carrying capacity.
Figure 7. The variation in the components of the load-carrying capacity per unit width for (a) the spanwise-varying channel and (b) the extruded surfaces versus the channel aspect ratio (
$a = \ell _x / \ell _z$
) for
$\textit{De} = 0.1$
.
3.3.2. Streamwise-normal stress
The load-carrying capacity for finite
$\textit{De}$
is the result of a contribution from each component of the polymer stress,
$[\tau _{xx}, \tau _{xy}]$
and
$[\tau _{xx}, \tau _{xy}, \tau _{xz}, \tau _{yz}, \tau _{zz}]$
, for the two-dimensional and the three-dimensional cases, respectively. We decompose the individual contributions to the load by solving the VR equation but retaining only the desired stress components on the right-hand side of (2.11). The result is shown in figure 7, for
$d = 0.2$
,
$\textit{De} = 0.1$
and
$\beta = 0.8$
.
The primary contribution comes from the streamwise-normal stress (
$\tau _{xx}$
), similar to what was found by Ahmed & Biancofiore (Reference Ahmed and Biancofiore2021) for the two-dimensional case, and a secondary effect from the remaining stress components which, upon summation, appear to cancel out. The streamwise component dominates, not surprisingly, due to the boundary motion along
$x$
which pulls and stretches the polymers. Similar to the two-dimensional case, the shear induces a strong stretch along the flow (streamlines), since
$\tau _{xx}$
depends quadratically on the shear rate even for very small
$\textit{De}$
. During the converging sections of the channel, this causes a build-up of the normal stress (stretching) that relieves the pressure, since the net force required to propel the fluid is now both the pressure and streamwise-normal stress. The reverse occurs during the divergent segments where the shear rate is arguably lower causing the polymers to relax. This stretch–relaxation mechanism (also measured via the tension in the streamlines) was first reported by Hinch et al. (Reference Hinch, Boyko and Stone2024) and relates the pressure variation to the tension in the streamlines. A notable difference between this work and that of Hinch et al. (Reference Hinch, Boyko and Stone2024) is the prevailing presence of shear due to the Couette flow which prevents complete relaxation to zero stress, and the smallest value of the normal stress is that of the pure-shear-driven state in the absence of pressure gradients.
In this three-dimensional channel, the normal stress distribution along
$x$
varies strongly versus
$a$
, in both the spanwise-varying and extruded configurations. The variation is attributed not only to decreasing pressure gradients, hence less shearing in the channel, but also to insufficient room for relaxation due to a shortening of the channel width. This is reflected in figure 8, which shows the distribution of the film-averaged streamwise-normal stress along the spanwise direction for different values of the aspect ratio (
$\textit{De} = 0.1$
,
$\beta = 0.8$
and
$d = 0.2$
). When spanwise surface gradients are present, see figure 8(a), the channel height reaches a maximum beyond
$z \gt s = 0.1$
and the shear rate reduces, causing the polymers to relax (reaching the pure shear value
$\tau _{xx} = 2(1-\beta )De \dot {\gamma }_{xy}^2$
, where
$\dot {\gamma }_{xy} = 1$
). Furthermore, we continue to shrink the channel by increasing
$a$
(decreasing the channel width), losing the polymers across the spanwise boundary due to side leakage (which is neglected completely in the two-dimensional cases).
Figure 8. The film-averaged normal stress distribution along the spanwise direction (averaged along
$x$
) for (a) the spanwise-varying channel and (b) the extruded channel, for different aspect ratios (
$a = \ell _x / \ell _z$
), using
$\textit{De} = 0.1$
and
$d = 0.2$
.
However, for the extruded channel, shown in figure 8(b), the polymers are exposed to a greater shear rate (
$\dot {\gamma }_{xy}\approx {1}/{\min (\textit{h})}= {1}/{(1-\textit{d})}$
), preventing them from relaxing along the width. A sudden jump is observed at the exit due to a strong leakage of the flow along the spanwise boundary. It is evident that varying the aspect ratio impacts (i) the room available for the polymer to stretch/relax and (ii) the rate at which the polymers escape along the spanwise open boundaries dictated by the spanwise velocity
$w$
. The latter always contributes to a decline in the load-carrying capacity.
We focus on these two competing mechanisms for the two different channel configurations. We adopt a measure of the spanwise velocity connected with the spanwise pressure gradient (independent of the film height):
\begin{equation} w_p = \frac {a}{2} \bigg ( h^2 \frac {\partial p}{\partial z} \bigg ), \end{equation}
which is depicted across the channel in figure 9(a–c) for the spanwise-varying channel and in figure 9(d–f) for the extruded channel using three different values of the channel aspect ratio. We observe for the spanwise-varying channel an increase close to the bulk region (
$z=0$
) which declines versus
$a$
and a strong leakage rate at the boundary (
$z=1$
). For the extruded case, shown in figure 9(d–f), the bulk spanwise velocity is negligible, and rises as the width decreases.
Figure 9. The distribution of the spanwise velocity (given by (3.4)) across the channels for (a–c) the spanwise-varying channel and (d–f) the extruded channel. The analysed aspect ratios (
$a = \ell _x / \ell _z$
) are (a,d)
$a = 0.2$
, (b,e)
$a = 1$
and (c,f)
$a = 5$
. Note that
$d = 0.2$
.
Hence, we define two characteristic spanwise velocity scales: (i) the velocity in the bulk and (ii) the velocity at the boundary, which generate and diminish, respectively, the effect of the viscoelasticity on the pressure (and then on the load). For the spanwise-varying case, the surface gradients along
$z$
persist for a length
$s$
(defined as the span of the Gaussian surface; see § 3.1) and diminish beyond this limit. For small
$a$
, the spanwise velocity is effectively zero for
$z \gt s$
. Therefore, the reference spanwise velocity in this region is taken as
$W_s = \max |w_p(x,z=s)|$
. On the other hand, the velocity contributing to leakage along the spanwise boundary
$W_a = \max |w_p(x,z=1)|$
continues to increase versus
$a$
, draining the channel of useful stretched polymers. For the extruded channel, the
$xy$
cross-section is constant, covering the entirety of the width, leading to negligible velocity close to the channel bulk (
$z = 0$
). For this case, only the leakage rate is considered significant.
Owing to these two different and competing flow mechanisms, we have two different effective Deborah numbers along the width: (i) the spanwise Deborah number
\begin{equation} De_s = \frac {a}{s}W_s De \end{equation}
and (ii) the leakage Deborah number
\begin{equation} De_a = a W_a De, \end{equation}
depicted in figure 10 versus the channel aspect ratio for the spanwise-varying and the extruded channels. For the spanwise-varying case, the spanwise Deborah number tends to rise as the channel width decreases raising the polymeric contribution to the load, while
$De_a$
, acting against it, rises exponentially when the width is sufficiently small. A second increase in
$De_s$
is observed around
$a = 5$
because the Gaussian protrusion now touches the spanwise boundary and the bulk flow merges with the leakage region. At this point, the useful stretch is lost along the open boundary. Once
$De_s$
is commensurable with
$\textit{De}$
(as
$a\sim O(1)$
) we reach the maximum load-carrying capacity (see figure 6
a). For the extruded case given in figure 10(b), we only obtain a continual loss of polymers due to leakage, depicted by a rise in
$De_a$
, which eventually overtakes the streamwise
$\textit{De}$
. This leads to a continual decline in load as observed in figure 6(b).
Figure 10. The spanwise effective
$De_s$
and the leakage
$De_a$
versus the channel aspect ratio (
$a = \ell _x / \ell _z$
) for (a) the spanwise-varying case and (b) the extruded case, using
$d = 0.2$
and
$\textit{De} = 0.1$
.
3.3.3. Spanwise forcing versus spanwise diffusion
In § 3.3.2, we introduced the spanwise Deborah number
$De_s$
and the leakage Deborah number
$De_a$
which compete to enhance and diminish the load, respectively. Furthermore, when
$a$
is O(
$1$
), the net force in the spanwise-varying case approaches its maximum value, as we have seen in figure 6. As observed from (2.10) and (2.11) (for LIN and VR models, respectively), the pressure is governed by a spanwise diffusion term (i.e.
$\mathcal{D} = a^2 ({\partial }/{\partial z}) ( ({\partial h^3}/{12}) ({\partial p}/{\partial z}) $
) of order
$a^2$
and a forcing term scaling with different powers of
$a$
(i.e.
$\mathcal{F}$
and
$\mathcal{G}$
). When
$a = 1$
, we find all terms of the same order of magnitude in both models and, therefore, expect a balance of spanwise diffusion and pressure forcing due to (i) channel surface variation and (ii) stress gradients.
We can visualise the influence of these terms on the load-carrying capacity using the LIN model which also predicts a maximum load when
$a=O(1)$
. This is depicted in figure 11 in which the load versus the channel aspect ratio is compared between VR and LIN models. While LIN is qualitatively reproducing well the features of the VR model, it is not able to accurately predict the peak load, which is achieved around
$a = 1.5$
while it is close to
$a=2.0$
from the more accurate VR model. This implies a 50 % error in predicting the optimal width of a channel. However, the exact value of the optimal
$a$
is not relevant to understand why we have a maximum for
$a\sim 1$
; therefore, we focus on LIN since it is simpler to analyse.
We delineate the dependence of the pressure (and hence the load) on the diffusive and forcing terms by exploiting
$\mathcal{F}$
,
$\mathcal{G}_1$
and
$\mathcal{G}_2$
from (2.10):
\begin{align} & \frac {\partial }{\partial x} \left ( \frac { h^3}{12} \frac {\partial p}{\partial x} \right ) + \mathcal{D} = \mathcal{F}, \end{align}
\begin{align} & \frac {\partial }{\partial x} \left ( \frac { h^3}{12} \frac {\partial p}{\partial x} \right ) + \mathcal{D} = \mathcal{G}_1, \end{align}
\begin{align} & \frac {\partial }{\partial x} \left ( \frac { h^3}{12} \frac {\partial p}{\partial x} \right ) + \mathcal{D} = \mathcal{G}_2, \end{align}
\begin{align} & \frac {\partial }{\partial x} \left ( \frac { h^3}{12} \frac {\partial p}{\partial x} \right ) = \mathcal{F}+\mathcal{G}_1+\mathcal{G}_2. \end{align}
The pressure obtained (i) from solving (3.7a
) will reflect the streamwise influence, (ii) from (3.7b
) the spanwise order-
$a^2$
influence, (iii) from (3.7c
) the spanwise order-
$a^4$
influence and (iv) from (3.7d
) the influence due to the absence of spanwise diffusion. The explicit expressions of
$\mathcal{F}$
,
$\mathcal{G}_1$
and
$\mathcal{G}_2$
are defined in (2.10) and are not repeated here for brevity.
Figure 11. The net force versus the channel aspect ratio for the spanwise-varying geometry. (a) A comparison between LIN (dashed lines) and VR (continuous lines) models for different Deborah numbers and (b) considering different order of magnitude for
$a$
using the
$\textit{De}$
-order linearised model (2.10).
In figure 11, we compare the load-carrying capacity versus
$a$
due to the diffusive and forcing terms (3.7d
). Firstly, as expected, the load in the absence of pressure diffusion is significantly larger. When
$a \gt 1$
the spanwise diffusion is strong, considerably decreasing the load. More importantly, the influence of
$\mathcal{G}_1$
and
$\mathcal{G}_2$
is negligible when
$a\sim 1$
showing that the maximum around
$a=1$
is only due to the interaction between the spanwise diffusion and the streamwise forcing
$\mathcal{F}$
(which has a maximum around
$a=1$
). Moving towards the limit
$a \ll 1$
, all effects diminish including the streamwise effects, since the channel surface gradients have not fully developed. Similarly, when
$a\gg 1$
, the mass exchange from the boundary dominates and all terms diminish the load. The balance lies between these two limits and hence we observe a maximum for
$a\sim O(1)$
.
3.4. Boundary conditions
The boundary conditions for the pressure have a strong influence on the net load-carrying capacity within the channel. As pointed out in the literature, the pressure is not properly defined in the viscoelastic case (Tichy & Bou-Saïd Reference Tichy and Bou-Saïd2008), and setting it equal to the ambient pressure, or any value, is not appropriate (Sawyer & Tichy Reference Sawyer and Tichy1998). In a similar manner, excluding the stress altogether may also be equally inappropriate. Balancing the force at the open boundaries, we require the expressions for total force
\begin{equation} \boldsymbol{F}^* = \int _{s} \boldsymbol{\sigma }^* {\rm d}\boldsymbol{s}^* \end{equation}
to vanish at the streamwise (
${\rm d}\boldsymbol{s}^* = {\rm d}y^* \,{\rm d}z^* \boldsymbol{e}_x$
) and spanwise (
${\rm d}\boldsymbol{s}^* = {\rm d}y^* \,{\rm d}x^* \boldsymbol{e}_z$
) boundaries. Non-dimensionalising using (2.2) and (2.3), and noting
$\boldsymbol{\sigma }^* = \boldsymbol{\tau }^* - p^*\unicode{x1D644} + 2\eta _s\unicode{x1D63F}^{\kern1pt *}$
, we obtain the conditions, up to O(
$\epsilon ^2$
):
\begin{align} p=\overline {\tau _{xx}}(0,z), \quad p=\overline {\tau _{xx}}(1,z), \quad p=\overline {\tau _{zz}}(x,1). \end{align}
Note that the channels modelled in this work are symmetric about the spanwise direction. Hence, the pressure and stress condition along the axis of symmetry (
$z = 0$
) satisfies
$({\partial p}/{\partial z}) = ({\partial \tau _{\textit{ij}}}/{\partial z})=0$
.
Using (3.9), we compare the load-carrying capacity versus the channel aspect ratio, as shown in figure 12, for three sets of boundary conditions: (i) the Newtonian condition
$p=0$
(case A), (ii) the Newtonian condition along the streamwise open boundaries, i.e.
$p(0, z) = p(1, z) = 0$
, and normal stress conditions along the spanwise boundaries (case B) and (iii) normal stress conditions along all open boundaries using (3.9) (case C).
Figure 12. The load-carrying capacity versus the channel aspect ratio (
$a = \ell _x / \ell _z$
) for (a) the spanwise-varying channel and (b) the extruded channel, considering three different cases of boundary conditions: case A, Newtonian pressure condition along all open boundaries; case B, Newtonian pressure condition only along the streamwise open boundaries; and case C, the pressure balanced by the average normal stress along all open boundaries (fully viscoelastic), using
$\textit{De} = 0.1$
and
$d = 0.2$
.
It is evident from figure 12 that a finite increase in pressure, and consequently the load, occurs at the boundary owing to the extra stress of the polymer. For both surface configurations, when comparing the three cases, the addition of the spanwise-normal stress has a weak influence on the load (due to the absence of spanwise surface motion). As the channel width reduces further (
$a \gt 1$
), a net increase is observed, delaying an immediate drop to zero observed for the Newtonian conditions.
When all the open boundaries are influenced by the polymer normal stress components (i.e. case C), a large net increase in the force is obtained, dominated by the streamwise-normal stress component. This additional load persists over the entire range of
$a$
examined, diminishing only when the channel width becomes small in comparison with its length (
$a \gt 1$
). At this point, the net shear strain rate that stretches the polymer chains (both streamwise and spanwise) is reducing owing to a decreasing pressure gradient.
This additional gain is nearly one order of magnitude greater than that for the Newtonian conditions. Despite having two distinct surface profiles, for small values of
$a$
, we obtain the same values of the load and almost identical trends. As
$a$
decreases below one, the presence of any surface feature, no matter how sharp, is not felt by the fluid, since the terms involving spanwise variation scale with
$a^2$
. Thus, for small
$a$
, we have an effectively flat channel for which the extra stress is easily predicted via the linear model. This net increase in load was construed as spurious by Tanner (Reference Tanner1969), since classical lubrication theory for Newtonian fluids plainly predicts a zero pressure field for boundary-driven flows. Experiments may be necessary to shed light on the true conditions prevailing at the open boundaries under strongly sliding conditions. In addition to this, we observe for
$a \gt 1$
that the curved portion of the surface intersects with the open boundary, leading to strong re-entry of the fluid at
$z = 1$
. In this case, assuming that the fluid beyond
$z=1$
is at a fixed pressure condition, independent of the total stress (the blue curves in figure 12), may be inaccurate.
4. Conclusions
The flow of a thin viscoelastic lubricating film along a channel with boundary motion was studied numerically using the Oldroyd-B constitutive relation for a three-dimensional channel with a flat lower surface and an upper curved surface modelled as a Gaussian protrusion. Two special cases are considered for the curved surface: (i)
$h = h(x,z)$
and (ii) an extrusion of the
$xy$
cross-section,
$h=h(x)$
. For these two cases, we employ the VR approach as the numerical procedure for solving the governing equations and, also, compare the model’s performance with that of an arguably simpler and efficient
$\textit{De}$
-order perturbed model. The main observation is the possibility of enhancing the film’s load-bearing capability by exploiting the polymer’s elastic properties (as evidenced via experiments).
An interesting feature of three-dimensional channels is the ratio between the length and the width (aspect ratio) which strongly influences the forces developed in the channel. By varying the channel aspect ratio but keeping the apparent Deborah number constant (defined identically to that of a two-dimensional configuration), we can modify the strength of viscoelastic effects. The load versus
$a$
trend is similar to the two-dimensional case, i.e. monotonically increasing for small
$\textit{De}$
and then an eventual saturation for large
$\textit{De}$
. This was the case for both surface configurations.
A prediction of the load was also made via a
$\textit{De}$
-order perturbation model which served to validate the VR approach and demonstrate the viability of low-order models. Results compared favourably for low values of
$\textit{De}$
but exhibited a discrepancy with the VR model for large values. The reduced model in the three-dimensional case must be solved numerically, since a partial differential equation emerges for the leading order and the first order and, by extension, will also hold for the higher-order systems. While the
$\textit{De}$
-order case can still be tackled via classical FDMs (or any other discretisation method), the appearance of mixed- and higher-order derivatives requires careful treatment. In fact, we infer that for the
$De^2$
-order case, the stencil requirements for the FDM will become quite cumbersome and make the numerical treatment difficult.
The load enhancement was traced to the growth in the streamwise-normal stress in the channel because it aligns with the sliding direction. For the spanwise-varying channel, the increase is due to a rise in the bulk spanwise velocity, finite only in the vicinity of the protrusion and zero everywhere else. Defining as a spanwise Deborah number
$De_s$
, corresponding to this local spanwise velocity, we observe that
$De_s$
rises as
$a$
increases becoming commensurable to
$\textit{De}$
for
$a\sim O(1)$
(i.e. when the load is maximum). However, increasing
$a$
also narrows the width and exposes the useful stretched polymers to the open boundary. As the flow exits this boundary, the stretched polymers are lost and the load begins to drop. This leakage rate is measured by a competing leakage Deborah number
$De_a$
, which increases as well with
$a$
. Eventually
$De_a$
overpowers
$De_s$
and diminishes the load entirely. Note that for
$a \gt 5$
, the curved region covers the entirety of the channel, but the close proximity of the open boundaries causes the bulk flow to merge with the leakage, thereby nullifying the beneficial rise in the local effective
$\textit{De}$
. For the extruded channel, the bulk spanwise velocity remains zero, and as a consequence the spanwise
$\textit{De}$
is zero. On the other hand, the leakage Deborah number rises (similar to the spanwise-varying case), causing the load to always decline as
$a$
increases.
Finally, we examine the influence of including the polymer stress in the boundary conditions. The presence of surface sliding (Couette flow) leads to finitely large normal stresses even when the surface is flat. For both channel configurations, including the normal stress at the streamwise and spanwise boundaries leads to a dramatic enhancement of the film load-carrying capacity. To clarify whether this additional increase is indeed spurious will require an experimental effort that measures the excess pressure at the boundaries as they undergo sliding motion.
This work focused explicitly on channel surface configurations modelled neatly via a Cartesian coordinate system. However, many cases, e.g. tear film lubrication, synovial joint lubrication, to name a few, conform to a spherical coordinate system (if reduced-order modelling is desired). In such cases, the aspect ratio may be fixed (close to unity) and the problem prescribed entirely in three dimensions, especially in the absence of any geometrical symmetry. In these instances, the VR approach offers a means to predict the useful quantities such as surface shear stresses and pressure over the useful range of
$\textit{De}$
.
Funding
The authors would like to acknowledge the Turkish National Research Agency (TÜB
$\dot {\textrm{I}}$
TAK) for supporting this work under the project 221N576.
Declaration of interests
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Appendix A. Viscoelastic Reynolds equation in three dimensions
The VR equation for a two-dimensional channel can be extended to three-dimensional channels (in Cartesian configuration) following the procedure in Ahmed & Biancofiore (Reference Ahmed and Biancofiore2021):
\begin{align} &\,\,\frac {\partial }{\partial x} \left ( \frac {h^3}{12} \frac {\partial p}{\partial x} \right ) + a^2 \frac {\partial }{\partial z} \left ( \frac {h^3}{12} \frac {\partial p}{\partial z} \right ) = \frac {\partial \mathcal{F}}{\partial x} + \frac {\partial \mathcal{G}}{\partial z} , \end{align}
\begin{align} & \mathcal{F} = \frac {\beta }{2}h + \mathcal{A} \left (\frac {\partial \tau _{xx}}{\partial x}\right ) + \mathcal{B}(\tau _{xy}) + a\mathcal{A} \left ( \frac {\partial \tau _{xz}}{\partial z} \right )\!, \end{align}
\begin{align} &\,\,\, \mathcal{G} = a^2 \mathcal{A} \left (\frac {\partial \tau _{zz}}{\partial z}\right ) + a\mathcal{A} \left ( \frac {\partial \tau _{xz}}{\partial x}\right ) + a \mathcal{B}(\tau _{yz}), \end{align}
where
\begin{align} \mathcal{A}(\boldsymbol{\cdot }) &= \int _0^h \frac {y}{h} \int _{0}^{h} \int _{0}^{y^{\prime }} (\boldsymbol{\cdot }) \, {\rm d}y^{\prime }{\rm d}y{\rm d}y - \int _0^h \int _{0}^{y^{\prime \prime }} \int _{0}^{y^{\prime }} (\boldsymbol{\cdot }) \, {\rm d}y^{\prime }{\rm d}y^{\prime \prime }{\rm d}y, \end{align}
\begin{align} \mathcal{B}(\boldsymbol{\cdot }) &= \int _0^h \frac {y}{h} \int _{0}^{h} (\boldsymbol{\cdot }) \, {\rm d}y^{\prime }{\rm d}y - \int _0^h \int _{0}^{y^{\prime }} (\boldsymbol{\cdot }) \, {\rm d}y^{\prime }{\rm d}y \end{align}
are integral operators. These operators can be further reduced by (i) applying integration by parts and (ii) switching to the computational configuration;
$h(x,z) = H(X, Z)$
,
$y = YH(X)$
:
\begin{align} &\qquad\qquad \int _0^h \int _{0}^{y^{\prime }} q \, {\rm d}y^{\prime }{\rm d}y = H^2\int _0^1 (Y-1)q {\rm d}Y = H^2\overline {[(Y-1)q]}, \end{align}
\begin{align} &\qquad\qquad\qquad\qquad\quad \int _0^h \frac {y}{h}\int _{0}^{h} q \, {\rm d}y^{\prime }{\rm d}y = \frac {H^2}{2}\overline {[q]}, \end{align}
\begin{align} & \int _0^h \int _{0}^{y^{\prime }}\int _{0}^{y^{\prime \prime }} q \, {\rm d}y^{\prime \prime }{\rm d}y^{\prime }{\rm d}y = H^3\int _0^1 \int _{0}^{Y^{\prime }}\int _{0}^{Y^{\prime \prime }} q \, {\rm d}Y^{\prime \prime }{\rm d}Y^{\prime }{\rm d}Y = H^3\overline {\left[\left(\frac {Y}{2}-Y+\frac {1}{2}\right)q\right]}, \end{align}
\begin{align} &\qquad\qquad\qquad \int _0^h \frac {y}{h}\int _{0}^{y^{\prime }}\int _{0}^{y^{\prime \prime }} q \, {\rm d}y^{\prime \prime }{\rm d}y^{\prime }{\rm d}y = H^3\overline {\left[\left(-\frac {Y}{2}+\frac {1}{2}\right)q\right]}, \end{align}
for the quantity
$q$
, yielding (upon substitution into (A1a
) and (A2))
\begin{align} &\qquad\qquad\qquad\frac {\partial }{\partial X} \left ( \frac {H^3}{12} \frac {\partial p}{\partial X} \right ) + a^2 \frac {\partial }{\partial Z} \left ( \frac {H^3}{12} \frac {\partial p}{\partial Z} \right ) = \frac {\partial \mathcal{F}}{\partial X} + \frac {\partial \mathcal{G}}{\partial Z} , \\[-12pt] \nonumber \end{align}
\begin{align} & \mathcal{F} = \frac {\beta }{2}H + \frac {H^3}{2}\overline {\left[(Y-Y^2)\frac {\partial \tau _{xx}}{\partial x}\right]} + \frac {H^2}{2}\overline {[(2Y-1)\tau _{xy}]} + a\frac {H^3}{2}\overline {\left[(Y-Y^2)\frac {\partial \tau _{xz}}{\partial z}\right]}, \\[-12pt] \nonumber \end{align}
\begin{align} &\quad \mathcal{G} = a^2\frac {H^3}{2}\overline {\left[(Y-Y^2)\frac {\partial \tau _{zz}}{\partial z}\right]} + \frac {H^2}{2}\overline {[(2Y-1)\tau _{zy}]} + a\frac {H^3}{2}\overline {\left[(Y-Y^2)\frac {\partial \tau _{xz}}{\partial x}\right]}. \end{align}
Appendix B. Three-dimensional curvilinear domain
In this appendix, following the work of Boyko et al. (Reference Boyko, Hinch and Stone2024), we provide the curvilinear coordinates used in the numerical treatment of the governing system of equations. The coordinates are transformed through
\begin{align} X^* = x^* + \epsilon Q^*(x^*, y^*), \quad Y^* = H_0\frac {y^*}{H^*(x^*, z^*)} , \quad Z^* = z^* + \epsilon R^*(z^*, y^*), \end{align}
where the asterisk denotes dimensional quantities and
$Q^*$
and
$R^*$
are the functions accounting for the channel variation along the streamwise and the spanwise directions. Using (2.2), and following the procedure leading to an orthogonal coordinate system (see Appendix A in Boyko et al. (Reference Boyko, Hinch and Stone2024)), an orthogonal coordinate system is found by computing
\begin{align} \boldsymbol{\nabla }X & = \left[\epsilon \left(1 + \epsilon ^2 \frac {\partial Q}{\partial x}\right), \epsilon ^2 \frac {\partial Q}{\partial y}, 0\right] \!, \end{align}
\begin{align} \boldsymbol{\nabla }Y & = \left[-\epsilon \frac {y}{H^2} \frac {\partial H}{\partial x}, \frac {1}{H}, -a\epsilon \frac {y}{H^2} \frac {\partial H}{\partial z}\right]\!, \end{align}
\begin{align} \boldsymbol{\nabla }Z & = \left[0, \epsilon ^2 \frac {\partial R}{\partial y}, a\epsilon \left(1 + \epsilon ^2 \frac {\partial R}{\partial z}\right)\right]\!, \end{align}
where
$\boldsymbol{\nabla } = \epsilon ({\partial }/{\partial x})\boldsymbol{e}_x + ({\partial }/{\partial y})\boldsymbol{e}_y + a\epsilon ( {\partial }/{\partial z})\boldsymbol{e}_z$
is the dimensionless gradient operator. Setting
$\boldsymbol{\nabla }X \boldsymbol{\cdot }\boldsymbol{\nabla }Y = 0$
and
$\boldsymbol{\nabla }Z \boldsymbol{\cdot }\boldsymbol{\nabla }Y = 0$
, we obtain the desired expressions up to O(
$\epsilon ^4$
):
\begin{align} Q & = \frac {H}{2} \frac {\partial H}{\partial x} (Y^2 - 1), \end{align}
\begin{align} R & = a^2\frac {H}{2} \frac {\partial H}{\partial z} (Y^2 - 1), \end{align}
where
$Y = y/H$
. Finally, the curvilinear basis vectors are defined (in dimensional form) as
\begin{align} \boldsymbol{e}_X = \left .\frac {\partial \boldsymbol{x}^*}{\partial X^*} \right/ \left|\left|\frac {\partial \boldsymbol{x}^*}{\partial X^*}\right|\right| \!, \quad \boldsymbol{e}_Y = \left .\frac {\partial \boldsymbol{x}^*}{\partial Y^*} \right/ \left|\left|\frac {\partial \boldsymbol{x}^*}{\partial Y^*}\right|\right| \! , \quad \boldsymbol{e}_Z = \left .\frac {\partial \boldsymbol{x}^*}{\partial Z^*} \right/ \left|\left|\frac {\partial \boldsymbol{x}^*}{\partial Z^*}\right|\right| \!, \end{align}
where
$\boldsymbol{x}^* = x^* \boldsymbol{e}_x + y^* \boldsymbol{e}_y + z^* \boldsymbol{e}_z$
. Substituting (B3) into (B1
a
,
b
,
c
), we obtain
\begin{align} x & = X - \epsilon ^2 \frac {H}{2} \frac {\partial H}{\partial X} (Y^2 - 1), \end{align}
\begin{align} y & = YH, \end{align}
\begin{align} z & = Z - a^2\epsilon ^2 \frac {H}{2} \frac {\partial H}{\partial Z} (Y^2 - 1), \\[6pt] \nonumber \end{align}
and then computing the derivatives in (B4 a , b , c ) gives us the required basis vectors in the curvilinear frame (note that we use the scaling from (2.2)):
\begin{align} &\qquad\quad \boldsymbol{e}_X = \boldsymbol{e}_x + \epsilon Y\frac {\partial H}{\partial X}\boldsymbol{e}_y, \end{align}
\begin{align} &\,\, \boldsymbol{e}_Y = -\epsilon Y\frac {\partial H}{\partial X}\boldsymbol{e}_x + \boldsymbol{e}_y - a\epsilon Y\frac {\partial H}{\partial Z}\boldsymbol{e}_z, \end{align}
\begin{align} &\qquad\quad \boldsymbol{e}_Z = a\epsilon Y\frac {\partial H}{\partial Z}\boldsymbol{e}_y +\boldsymbol{e}_z. \\[12pt] \nonumber \end{align}
The components of the velocity related through
\begin{align} &\quad\boldsymbol{u}^* = \unicode{x1D648} \, \boldsymbol{U}^*, \end{align}
\begin{align} &\unicode{x1D648} = [\boldsymbol{e}_X \, \boldsymbol{e}_Y \, \boldsymbol{e}_Z], \end{align}
and are non-dimensionalised using (2.2):
\begin{align} U = u + O(\epsilon ^2), \quad V = v - Y\frac {\partial H}{\partial X} u - a Y \frac {\partial H}{\partial Z} w,\quad W=w+O(\epsilon ^2). \end{align}
Equations (B8 a , b , c ) are used in the numerical procedure, specifically in the flux-conserving discretisation of the stress advection.
Appendix C. Implicit discretisation
The system of equations presented in this work are solved via the classical FDM applied to the curvilinear (or computational) domain
$\boldsymbol{X}=(X, Y, Z)$
in which the grid spacing along each axis is constant. The numerical solution of the constitutive relation, a system of coupled hyperbolic equations, is achieved by first converting the advective component of the material derivative operator
$(\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{\tau }$
into
$ \boldsymbol{\nabla } \boldsymbol{\cdot }(\boldsymbol{u}\otimes \boldsymbol{\tau })$
, where
$\boldsymbol{u}\otimes \boldsymbol{\tau }$
is a third-order tensor. Note that by utilising the continuity equation, we can show that the two expressions are identical. However, from a numerical perspective, the former expression causes the FDM to suffer from numerical dissipation when the solution ceases to be smooth. In these cases, the attempt to increase the accuracy by refining the mesh tends to have little effect if higher-order schemes are used.
We pursue a fully implicit numerical approach which requires a sparse-matrix operator for each of the terms appearing in the constitutive relation (for each component of the extra stress tensor). First we convert the expression
$ \boldsymbol{\nabla } \boldsymbol{\cdot }(\boldsymbol{u}\otimes \boldsymbol{\tau })$
, written explicitly as
\begin{equation} \boldsymbol{\nabla } \boldsymbol{\cdot }(\boldsymbol{u}\otimes \boldsymbol{\tau }) = \bigg ( u\frac {\partial }{\partial x} + v\frac {\partial }{\partial y} + aw\frac {\partial }{\partial z} \bigg ) \boldsymbol{\tau }, \end{equation}
into the computational domain via (B8 a , b , c ). Carrying out the simplification and noting that
\begin{equation} \frac {\partial (HU)}{\partial X} + \frac {\partial V}{\partial Y} + a\frac {\partial (HW)}{\partial Z}=0, \end{equation}
we get
\begin{equation} \boldsymbol{\nabla } \boldsymbol{\cdot }(\boldsymbol{u}\otimes \boldsymbol{\tau }) = \frac {1}{H}\bigg ( \frac {\partial (HU\boldsymbol{\tau })}{\partial X} + \frac {\partial (V\boldsymbol{\tau })}{\partial Y} + a\frac {\partial (HW\boldsymbol{\tau })}{\partial Z} \bigg ). \end{equation}
Discretising (C3) via the classical FDM will result in a matrix-vector product
$[\boldsymbol{\nabla } \boldsymbol{\cdot }(\boldsymbol{u}\otimes \boldsymbol{\tau })] = \unicode{x1D63C}\boldsymbol{q}$
, where
$\unicode{x1D63C}$
is the sparse coefficient matrix and
$\boldsymbol{q}$
is the unknown vector containing the extra stress components
$\tau _{\textit{ij}}$
. The sparse-matrix operators needed for the implicit treatment are then obtained:
\begin{equation} \unicode{x1D63C} = \unicode{x1D63F}_{\unicode{x1D643}\unicode{x1D65E}} (\unicode{x1D63F}_{\partial \unicode{x1D653}} \unicode{x1D63F}_{\unicode{x1D643}\unicode{x1D650}} + \unicode{x1D63F}_{\partial \unicode{x1D654}} \unicode{x1D63F}_{\unicode{x1D651}} + a \unicode{x1D63F}_{\partial \unicode{x1D655}} \unicode{x1D63F}_{\unicode{x1D643}\unicode{x1D652}}), \end{equation}
where
$\unicode{x1D63F}_{\unicode{x1D643}\unicode{x1D65E}} = \textit{diag}(H^{-1})$
,
$\unicode{x1D63F}_{\unicode{x1D643}\unicode{x1D650}} = \textit{diag}(HU)$
,
$\unicode{x1D63F}_{\unicode{x1D651}} = \textit{diag}(V)$
and
$\unicode{x1D63F}_{\unicode{x1D643}\unicode{x1D652}} = \textit{diag}(HW)$
are sparse-diagonal matrices and
$\unicode{x1D63F}_{\partial \unicode{x1D653}}$
,
$\unicode{x1D63F}_{\partial \unicode{x1D654}}$
and
$\unicode{x1D63F}_{\partial \unicode{x1D655}}$
are the sparse-difference matrix operators.
In this setting, a fourth-order scheme can be used for the operators along
$Y$
and
$Z$
, and based on
$\textit{De}$
, a second-order flux-conserving scheme for
$X$
. The latter is more sensitive at high
$\textit{De}$
and leading to poor convergence (requiring a large number of iterations) if higher-order schemes are used. It was observed that increasing the nodes along
$X$
, and using a first-order hybrid scheme that is unconditionally numerically stable, was far more effective in terms of convergence and accuracy at higher
$\textit{De}$
. Despite these measures, no noticeable improvement in the HWNP was observed since it persists within the polynomial-based finite-difference, volume and element discretisation measures.
Appendix D. Coordinate transformation
The surfaces in lubricated channels are curved and the application of the classical FDM of constant mesh spacing is not directly applicable. We transform the derivatives appearing in the governing equations from the rectilinear domain
$\boldsymbol{\partial x}$
to the curvilinear domain
$\boldsymbol{\partial X}$
:
\begin{equation} \boldsymbol{\partial X} = \unicode{x1D63C} \boldsymbol{\partial x}, \end{equation}
where
\begin{align} & \boldsymbol{\partial X}^T = \bigg [ \frac {\partial }{\partial X} , \frac {\partial }{\partial Y} , \frac {\partial }{\partial Z} , \frac {\partial ^2 }{\partial X^2} , \frac {\partial ^2 }{\partial Y^2} , \frac {\partial ^2 }{\partial Z^2} , \frac {\partial ^2 }{\partial X \partial Y} , \frac {\partial ^2 }{\partial X \partial Z} , \frac {\partial ^2 }{\partial Y \partial Z} \bigg ], \end{align}
\begin{align} &\quad\,\, \boldsymbol{\partial x}^T = \bigg [ \frac {\partial }{\partial x} , \frac {\partial }{\partial y} , \frac {\partial }{\partial z} , \frac {\partial ^2 }{\partial x^2} , \frac {\partial ^2 }{\partial y^2} , \frac {\partial ^2 }{\partial z^2} , \frac {\partial ^2 }{\partial x \partial y} , \frac {\partial ^2 }{\partial x \partial z} , \frac {\partial ^2 }{\partial y \partial z} \bigg ], \\[12pt] \nonumber \end{align}
and
$\unicode{x1D63C}$
is the Jacobian of the transformation. The detailed contents of
$\boldsymbol {A}$
are omitted for the sake of brevity but are obtained via application of the chain rule and the total derivative. See Hirsch (Reference Hirsch2007) for the two-dimensional case which can be extended to the three-dimensional problem in a straightforward manner.
In this work, we only pursue the transformations given by (D1) to gain a computational advantage, namely that the operators
$\boldsymbol{\partial \!X}$
can be constructed independently of
$H$
and they enable discretisation via the classical finite-difference schemes which are easily available. The difference operators are then used to construct an implicit system of algebraic equations that results in a numerically efficient solution procedure (see Ahmed & Biancofiore (Reference Ahmed and Biancofiore2023) for more details).
Appendix E. Validation against direct numerical simulations
In this appendix, we present a comparison of the results obtained from LIN, VR and DNS to validate our approach. It should be noted that for a two-dimensional parabolic slider, whose height is given by
\begin{equation} h = 2x^2 -2x+1, \end{equation}
Ahmed & Biancofiore (Reference Ahmed and Biancofiore2021) found an excellent agreement between VR and DNS for low to moderate values of
$\textit{De}$
, until reaching the DNS limit (for
$De\gt 0.08$
, it was impossible to find a solution with DNS due to the HWNP). While LIN was shown to be accurate for low
$\textit{De}$
, the error was significantly increasing at moderate to large
$\textit{De}$
, clearly showing the higher degree of accuracy of VR with respect to LIN. Here, we extend the parabolic slider studied by Ahmed & Biancofiore (Reference Ahmed and Biancofiore2021) to three dimensions (see figure 13
a).
Figure 13. (a) The three-dimensional extruded parabolic slider analysed in Appendix E. Comparison of the pressure profiles between DNS, VR and LIN for the parabolic slider channel (b) along the spanwise direction at
$x = 0.33$
and (c) along the streamwise
$x$
direction at
$z = 0.5$
for
$a=1$
,
$\textit{De} = 0.04$
and
$\beta = 0.8$
.
As for the two-dimensional case, the DNS were performed via the robust RheoTool library implemented in OpenFOAM (Alves et al. Reference Alves, Oliveira and Pinho2021) that leverages the log-conformation representation of the constitutive relations for improved numerical stability (Fattal & Kupferman Reference Fattal and Kupferman2004). However, since the film thickness is small, we are constrained to Courant–Friedrichs–Lewy number
$ \lt 10^{-4}$
(since the library enforces a segregated explicit coupling between
$p$
,
$\boldsymbol{u}$
and
$\boldsymbol{\tau }$
). In fact, we used a mesh resolution of
$(N_x, N_y, N_z) = (96, 32, 96)$
which gave grid-independent results. However, the simulation took approximately 10 days to complete (on four Intel Xeon gold 6132 2.6 GHz processors and 128 gigabytes of RAM) due to the small time-step size
$\Delta t \approx 10^{-9}$
(adaptively modified to favour convergence and numerical stability).
We present a comparison of the film pressure for the three-dimensional parabolic slider channel with length-to-width ratio
$a = 1$
,
$\textit{De} = 0.04$
and
$\beta = 0.8$
. In figures 13(b) and 13(c), the pressure distributions along the spanwise
$z$
direction and streamwise
$x$
direction, respectively, are illustrated. The VR model agrees very well with the predictions of DNS along both the streamwise and the spanwise directions. However, the LIN model only shows a good agreement along the streamwise direction while deviating from DNS along the spanwise direction. Furthermore, the calculated load from DNS is
$F_{\ell , \textit{DNS}} = 0.00563$
, from the VR model is
$F_{\ell ,\textit{VR}} = 0.00584$
and from the LIN model is
$F_{\ell ,\textit{LIN}} = 0.00608$
. The errors are
$3.7\,\%$
and
$8\,\%$
for the VR and LIN models, respectively. We have shown in Ahmed & Biancofiore (Reference Ahmed and Biancofiore2021) that this deviation will grow versus
$\textit{De}$
in two-dimensional cases, and we do expect the same trend also for three-dimensional simulations. Unfortunately, for
$\textit{De} \gt 0.04$
, the simulations failed to converge due to the HWNP, exhibiting stress residuals that would not decrease below
$10^{-6}$
. It is possible that a fully implicit implementation could reduce this computation time significantly and perhaps is the direction to take for the DNS of viscoelastic thin films with boundary motion (Wittschieber, Demkowicz & Behr Reference Wittschieber, Demkowicz and Behr2022).
Open access
