1 Introduction
The ‘art’ of lubrication by thin liquid layers has been known since ancient times (Dowson Reference Dowson1998), permitting motion between adjacent solid surfaces at low friction and wear. Lubrication is of paramount importance to the safe, reliable and controlled operation of many key elements in engineering applications ranging from very-large-scale (planes, wind turbines) to microfluidic devices. Synovial joints in mammals are the archetype of this mechanism in nature. From a theoretical point of view, the flow of a liquid within a narrow gap can be described by the lubrication approximation of Stokes equations, first developed by Reynolds more than a century ago (Reynolds Reference Reynolds1886). Since then, lubrication theory has been used to understand a wide range of phenomena like moving bubbles in a tube (Bretherton Reference Bretherton1961), motion of red blood cells in capillaries (Fitz-Gerald Reference Fitz-Gerald1969; Secomb et al. Reference Secomb, Skalak, Özkaya and Gross1986; Feng & Weinbaum Reference Feng and Weinbaum2000), biomechanics of articular cartilage (Hou et al. Reference Hou, Mow, Lai and Holmes1992; Mow, Ateshian & Spilker Reference Mow, Ateshian and Spilker1993) or the physics of the ‘Kugel fountain’ (Snoeijer & van der Weele Reference Snoeijer and van der Weele2014), to name a few examples.
Due to the reversibility of Stokes flow in a lubricating layer between rigid solids, the lift force
$$\begin{eqnarray}L=-\int \unicode[STIX]{x1D61B}_{zj}n_{j}\,\text{d}A\simeq \int p\,\text{d}A\end{eqnarray}$$
on the sliding or rotating bodies vanishes: in a non-cavitating liquid a lubricated contact could not support any load if the bodies were entirely rigid. Here, the
$\unicode[STIX]{x1D61B}_{zj}$
are the components of the total stress tensor in fluid, which in the lubrication approximation is dominated by the isotropic part (given by the fluid pressure
$p$
)
$n_{j}$
is the outward normal to the solid surface, and
$A$
is the area of the solid boundary. However, usually the counter-moving bodies are deformable. Importantly, the deformation has been found to break the reversibility of Stokes flow in the lubricating layer, which generates a lift force
$L>0$
between the bodies (Hooke & O’Donoghue Reference Hooke and O’Donoghue1972; Bissett Reference Bissett1989; Sekimoto & Leibler Reference Sekimoto and Leibler1993; Snoeijer, Eggers & Venner Reference Snoeijer, Eggers and Venner2013). Motivated in part by novel biological or bio-inspired engineering applications, this problem has been addressed on many occasions in the last decade. Subjects range from ‘soft lubrication’ (Martin et al.
Reference Martin, Clain, Buguin and Brochard-Wyart2002; Skotheim & Mahadevan Reference Skotheim and Mahadevan2004, Reference Skotheim and Mahadevan2005), motion of lubricated eyelid wipers (Jones et al.
Reference Jones, Fulford, Please, McElwain and Collins2008), sticking of particles on lubricated compliant substrates (Urzay Reference Urzay2010; Mani, Gopinath & Mahadevan Reference Mani, Gopinath and Mahadevan2012) to translating spinning particles near soft boundaries (Urzay, Llewellyn Smith & Glover Reference Urzay, Llewellyn Smith and Glover2007; Salez & Mahadevan Reference Salez and Mahadevan2015; Saintyves et al.
Reference Saintyves, Jules, Salez and Mahadevan2016).
In most studies, the deformation of elastic materials has been assumed to adapt instantaneously to the stresses at their boundaries. In practice, however, most soft materials like gels, elastomers or cartilage behave viscoelastically: due to dissipation, their relaxation behaviour is time-dependent, and deformation requires finite time to adapt to changes in the loading. For example, recent experiments show the viscoelastic nature of articular cartilage during osteoarthritis (Trickey, Lee & Guilak Reference Trickey, Lee and Guilak2000; Desrochers, Amrein & Matyas Reference Desrochers, Amrein and Matyas2012). Hence, the coupling between lubrication pressure and viscoelastic deformation is crucial in understanding these systems (Hooke & Huang Reference Hooke and Huang1997; Scaraggi & Persson Reference Scaraggi and Persson2014).
In this paper we investigate how the lubrication is affected by viscoelasticity of the deformable boundaries. Figure 1(a) shows four principal configurations in steady soft lubrication in a two-dimensional configuration: translation with constant velocity
$V$
(cases (i, ii)) or rotation with constant frequency
${\it\omega}$
(cases (iii, iv)) of cylinders near a wall. Either the wall (i, iii) or the cylinder (ii, iv) is assumed to be soft. Contrary to freely moving particles, here we impose translation or rotation on the cylinder to decouple the intricate particle dynamics from the viscoelastic effects. When the soft material is assumed to be perfectly elastic, adapting instantaneously to changes in loading, all four cases are equivalent. At constant separation distance, the lift force and the velocity are then related as
$L\sim V^{2}$
when the deformation is small (Skotheim & Mahadevan Reference Skotheim and Mahadevan2004) and as
$L\sim |V|^{3}$
when the deformation is large (Bissett Reference Bissett1989; Snoeijer et al.
Reference Snoeijer, Eggers and Venner2013).
(a) Four typical geometries of two-dimensional soft lubrication. The arrows show the direction of steady motion of the cylinder. If the deformed material is purely elastic in nature, all four have the same steady-state response when studied in a comoving frame attached to the cylinder. Viscoelasticity breaks the equivalence among the four cases (see the text for details). (b) Schematic of the problem studied here: a steadily moving rigid cylinder of radius
$R$
close to a viscoelastic wall (distance
${\it\Delta}$
), corresponding to case (i). The fluid has a dynamic viscosity
${\it\eta}$
.
$\mathscr{H}$
is the deformation of the viscoelastic wall. It should be noted that case (iv) can be recovered by setting
$V\rightarrow {\it\omega}R$
.
Importantly, viscoelasticity breaks the equivalence among the various cases in figure 1(a). In cases (ii, iii), the deformations are stationary relative to the material points. In stationary conditions, one only probes the long-time relaxation and the response is purely elastic. Contrary to this, in cases (i, iv), the material points in the deformed body are exposed to a dynamic loading and the deformations lag behind the steadily moving fluid pressure. This continuous reconfiguration of the viscoelastic material dissipates energy and probes the time-dependent rheology of the solid. The key question addressed here is how the viscoelastic rheology affects the generated lift force.
This paper is organised as follows. In § 2 we formulate the problem based on the linear viscoelastic response and discuss the solution strategy. Section 3 presents the deformation profiles for three canonical examples, namely Kelvin–Voigt, standard linear and power law rheology. We show how viscoelasticity affects the lift force in these cases, by a combination of numerical solution and asymptotic analysis. The paper concludes with a discussion in § 4.
2 Formulation
We focus on the geometry sketched in figure 1(b), where the cylinder is treated as rigid and the lower boundary as viscoelastic. We restrict the analysis to small deformations, in which case the analysis equally applies to case (iv), using the connection
$V\rightarrow {\it\omega}R$
. Below we first formulate the problem and introduce dimensionless variables, and subsequently explain the solution strategy to solve the viscoelastic lubrication problem.
2.1 Lubrication equation and viscoelastic deformation
A rigid cylinder of radius
$R$
moves with a constant velocity
$V$
within a fluid of dynamic viscosity
${\it\eta}$
. The minimum distance between the cylinder and the undeformed wall is
${\it\Delta}$
(see figure 1
b). In the limit
${\it\Delta}\ll R$
, the shape of the cylinder is described by a parabola in the first approximation and the gap profile is given by
$h_{0}(x)={\it\Delta}+x^{2}/2R={\it\Delta}(1+x^{2}/2R{\it\Delta})$
. Hence, the characteristic length of the contact zone becomes
$$\begin{eqnarray}\ell =\sqrt{2{\it\Delta}R}.\end{eqnarray}$$
The motion of the cylinder creates a lubrication pressure in the gap, and this in turn deforms the wall. This deformation is characterised by
$\mathscr{H}(x)$
. The deformed gap profile is
$h(x)=h_{0}(x)-\mathscr{H}(x)$
. The profile of the thin gap
$h(x)$
and the fluid pressure
$p(x)$
in the gap are the unknowns of this coupled problem. They are related by two equations: the steady-state hydrodynamic lubrication equation and the relation between load and deformation of a viscoelastic half-space.
We first compute the deformation by treating the lubrication pressure as a traction acting on a semi-infinite viscoelastic solid. The stress in an incompressible (Poisson ratio
${\it\nu}=1/2$
) linear viscoelastic material under dynamic strain is (Ferry Reference Ferry1961)
$$\begin{eqnarray}{\it\sigma}_{ij}(x,z,t)=\int _{-\infty }^{t}{\it\Psi}(t-t^{\prime })\dot{{\it\epsilon}}_{ij}(x,z,t^{\prime })\,\text{d}t^{\prime }-{\it\Pi}(x,z,t){\it\delta}_{ij},\end{eqnarray}$$
where
${\it\sigma}_{ij}$
are the components of the stress tensor,
${\it\epsilon}_{ij}$
are the components of the strain tensor,
${\it\Psi}(t)$
is the shear relaxation function and
${\it\Pi}$
is the isotropic part of the stress tensor. The dot represents a time derivative. In this two-dimensional setting, the deformation of the viscoelastic wall is assumed to be zero in the
$y$
direction – this leads to plane strain conditions, for which all the strain components normal to the page vanish. For the case of inertia-free dynamics, mechanical equilibrium is defined by
${\it\sigma}_{ij,j}=0$
. We apply a Fourier transform in time (defined as
$\widehat{f}({\it\omega})=\int _{-\infty }^{\infty }f(t)\text{e}^{-\text{i}{\it\omega}t}\,\text{d}t$
) to (2.2) and
${\it\sigma}_{ij,j}=0$
to obtain
$$\begin{eqnarray}\displaystyle & \displaystyle \widehat{{\it\sigma}}_{ij}(x,z,{\it\omega})={\it\mu}({\it\omega})\widehat{{\it\epsilon}}_{ij}(x,z,{\it\omega}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \widehat{{\it\sigma}}_{ij,j}(x,z,{\it\omega})=0 & \displaystyle\end{eqnarray}$$
${\it\mu}({\it\omega})$
is the complex shear modulus of the material and is given by 
$$\begin{eqnarray}{\it\mu}({\it\omega})=\text{i}{\it\omega}\int _{0}^{\infty }{\it\Psi}(t)\text{e}^{-\text{i}{\it\omega}t}\,\text{d}t=G^{\prime }({\it\omega})+\text{i}G^{\prime \prime }({\it\omega}),\end{eqnarray}$$
where
$G^{\prime }({\it\omega})$
and
$G^{\prime \prime }({\it\omega})$
are the storage and loss moduli. Equation (2.3) can be solved for the surface profile
$\mathscr{H}$
by a Green’s function approach. Recently, a similar formulation was used to study dynamic deformation of a viscoelastic substrate under a moving contact line (Karpitschka et al.
Reference Karpitschka, Das, van Gorcum, Perrin, Andreotti and Snoeijer2015). For an arbitrary dynamic traction
$\widehat{p}(x,{\it\omega})$
, applied at the top surface of the wall (
$z=0$
), the deformation is given by
$$\begin{eqnarray}\widehat{\mathscr{H}}(x,{\it\omega})=\int _{-\infty }^{\infty }\widehat{p}(x^{\prime },{\it\omega})\frac{\mathscr{K}(x-x^{\prime })}{{\it\mu}({\it\omega})}\,\text{d}x^{\prime },\end{eqnarray}$$
where
$\mathscr{K}(x)$
is the elastic Green’s function. Taking a Fourier transform of (2.5) in space (defined as
$\widetilde{f}(q)=\int _{-\infty }^{\infty }f(x)\text{e}^{-\text{i}qx}\,\text{d}x$
), we obtain
$$\begin{eqnarray}\widehat{\widetilde{\mathscr{H}}}(q,{\it\omega})=\widehat{\widetilde{p}}(q,{\it\omega})\frac{\widetilde{\mathscr{K}}(q)}{{\it\mu}({\it\omega})}.\end{eqnarray}$$
The Green’s function for an elastic half-space is
$\mathscr{K}(x)=\log |x|/2{\rm\pi}$
, or, in Fourier space,
$\widetilde{\mathscr{K}}(q)=-1/2|q|$
(Johnson Reference Johnson1987).
In the present problem the cylinder moves at a constant velocity, so that the dynamical loading has the form of a travelling pressure wave
$p(x-Vt)$
. This simple form of the temporal loading enables detailed analysis taking into account the full history-dependent response. Taking Fourier transforms of
$p(x-Vt)$
with respect to space and time, we reach
$\widehat{\widetilde{p}}(q,{\it\omega})=2{\rm\pi}\widetilde{p}(q){\it\delta}({\it\omega}+Vq)$
, where
${\it\delta}({\it\omega})$
is the Delta function. Then (2.6) simplifies to
$$\begin{eqnarray}\widehat{\widetilde{\mathscr{H}}}(q,{\it\omega})=-\frac{{\rm\pi}\widetilde{p}(q)}{|q|}\frac{{\it\delta}({\it\omega}+Vq)}{{\it\mu}({\it\omega})}.\end{eqnarray}$$
Taking a backward transform from
${\it\omega}$
to
$t$
(defined as
$f(t)=\int _{-\infty }^{\infty }\widehat{f}({\it\omega})\text{e}^{\text{i}{\it\omega}t}(\text{d}{\it\omega}/2{\rm\pi})$
), we obtain
$$\begin{eqnarray}\widetilde{\mathscr{H}}(q,t)=-\frac{\widetilde{p}(q)}{2|q|}\frac{\text{e}^{-\text{i}Vqt}}{{\it\mu}(-Vq)}.\end{eqnarray}$$
The only remaining time dependence is the phase factor
$\text{e}^{-\text{i}Vqt}$
, which describes the translational motion of the cylinder. Hence, in the comoving frame that travels with the cylinder, the profile becomes
$$\begin{eqnarray}\mathscr{H}(x)=\int _{-\infty }^{\infty }-\frac{\widetilde{p}(q)}{2|q|{\it\mu}(-Vq)}\text{e}^{\text{i}qx}\frac{\text{d}q}{2{\rm\pi}}.\end{eqnarray}$$
This is the first key equation of the problem, relating the lubrication pressure in the narrow gap to the deformation of the viscoelastic wall.
The second equation is obtained by the steady-state lubrication equation, describing the Stokes flow in the narrow gap. Using no-slip boundary conditions at both solids (horizontal motion due to viscoelastic deformation is negligible in the lubrication limit and the no-slip boundary condition holds on the soft surface), we obtain the horizontal velocity component
$$\begin{eqnarray}u(x,z)=V\left(1-\frac{z}{h}\right)+\frac{1}{2{\it\eta}}\frac{\text{d}p}{\text{d}x}(z^{2}-zh),\end{eqnarray}$$
where
${\it\eta}$
is the fluid viscosity and
$\text{d}p/\text{d}x$
is the pressure gradient in the horizontal direction. Integrating (2.10) over the gap we obtain the volume flux
$$\begin{eqnarray}Q=\int _{0}^{h(x)}u(x,z)\,\text{d}z,\end{eqnarray}$$
which for an incompressible flow at steady state is constant (
$\text{d}Q/\text{d}x=0$
). In the comoving frame, it thus leads to
$$\begin{eqnarray}\frac{\text{d}Q}{\text{d}x}=\frac{\text{d}}{\text{d}x}\left[\frac{1}{6{\it\eta}}h^{3}\frac{\text{d}p}{\text{d}x}+Vh\right]=0,\end{eqnarray}$$
where we recall that
$h(x)=h_{0}(x)-\mathscr{H}(x)$
is the thickness of the liquid layer. Equations (2.9) and (2.12) form a set of coupled equations that constitute the viscoelastic lubrication problem. For
${\it\mu}=\text{constant}$
, this is the same set of equations as for ‘classical’ 2D elastohydrodynamics (Bissett Reference Bissett1989; Venner & Lubrecht Reference Venner and Lubrecht2000; Snoeijer et al.
Reference Snoeijer, Eggers and Venner2013).
2.2 Non-dimensionalisation
We use the contact length
$\ell$
as the horizontal length scale and the gap height
${\it\Delta}$
as the vertical length scale. The lubrication pressure then scales as
$$\begin{eqnarray}P^{\ast }={\it\eta}\ell V/{\it\Delta}^{2}.\end{eqnarray}$$
The scale of the deformation induced by this lubrication pressure is given by
$$\begin{eqnarray}\mathscr{H}^{\ast }=P^{\ast }\ell /2G={\it\eta}VR/G{\it\Delta},\end{eqnarray}$$
where
$G$
is the static shear modulus of the viscoelastic material, defined as
$G^{\prime }({\it\omega}=0)$
. Hence, it is natural to introduce the first dimensionless parameter of the problem as
$$\begin{eqnarray}{\it\beta}\equiv \frac{\mathscr{H}^{\ast }}{{\it\Delta}}=\frac{{\it\eta}VR}{G{\it\Delta}^{2}},\end{eqnarray}$$
which is the ratio of the elastic deformation and the typical gap size. In this paper we solve the coupled equations in the limit where
$\mathscr{H}^{\ast }$
is small compared with
${\it\Delta}$
, i.e.
${\it\beta}\ll 1$
.
In contrast to the purely elastic case, the viscoelastic wall exhibits a relaxation time scale
${\it\tau}$
. This time scale needs to be compared with that of the dynamical loading due to the lubrication pressure. This pressure evolves on a time scale
${\it\tau}_{\text{p}}=\ell /V$
, which can be seen as the inverse shear rate at which the solid is excited. The ratio of these two time scales gives
$$\begin{eqnarray}\mathscr{T}\equiv \frac{{\it\tau}}{{\it\tau}_{\text{p}}}=\frac{{\it\tau}V}{\ell },\end{eqnarray}$$
the second dimensionless parameter in the problem. The parameter
$\mathscr{T}$
is the solid analogue of the Deborah number of a viscoelastic fluid. Equivalently,
$\mathscr{T}$
can be interpreted as a ratio of two length scales and compares the lateral extent of the viscoelastic deformation with the contact length. If the material relaxes much faster than the time scale of the changes of its load i.e.
$\mathscr{T}\ll 1$
, the material behaves purely elastically. If both time scales are comparable (
$\mathscr{T}\sim O(1)$
), viscoelasticity becomes important.
It turns out that
${\it\beta}$
and
$\mathscr{T}$
are the only two dimensionless groups in the problem. This is made explicit by introducing a set of non-dimensional variables,
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \overline{x}=\frac{x}{\ell },\quad \overline{z}=\frac{z}{{\it\Delta}},\quad \overline{h}=\frac{h}{{\it\Delta}},\quad \overline{\mathscr{H}}=\frac{\mathscr{H}}{\mathscr{H}^{\ast }},\quad \overline{p}=\frac{p}{P^{\ast }},\\ \displaystyle \overline{t}=\frac{t}{{\it\tau}},\quad \overline{{\it\mu}}=\frac{{\it\mu}}{G},\quad \overline{q}=q\ell ,\quad \overline{L}=\frac{L}{P^{\ast }\ell }.\end{array}\right\} & & \displaystyle\end{eqnarray}$$
In the remainder, we will only use dimensionless quantities and thus drop the overbars.
In dimensionless form (2.9) becomes
$$\begin{eqnarray}\mathscr{H}(x)=-\int _{-\infty }^{\infty }\frac{\widetilde{p}(q)}{|q|{\it\mu}(-\mathscr{T}q)}\text{e}^{\text{i}qx}\frac{\text{d}q}{2{\rm\pi}},\end{eqnarray}$$
which contains
$\mathscr{T}$
as a parameter. Likewise, the dimensionless lubrication equation becomes
$$\begin{eqnarray}\frac{\text{d}}{\text{d}x}\left[\frac{\text{d}p}{\text{d}x}h(x)^{3}+6h(x)\right]=0.\end{eqnarray}$$
The deformed gap profile
$h(x)$
couples (2.18) and (2.19) by
$$\begin{eqnarray}h(x)=h_{0}(x)-{\it\beta}\mathscr{H}(x),\end{eqnarray}$$
where
$h_{0}(x)=1+x^{2}$
. As anticipated, the problem contains only the two parameters
$\mathscr{T}$
and
${\it\beta}$
. In the case of a purely elastic deformation (
$\mathscr{T}=0$
), the solution depends only on
${\it\beta}$
(Hooke & O’Donoghue Reference Hooke and O’Donoghue1972; Bissett Reference Bissett1989; Snoeijer et al.
Reference Snoeijer, Eggers and Venner2013).
2.3 Solution strategy
We seek a perturbative solution of (2.19) and (2.20) for
${\it\beta}\ll 1$
, in the spirit of previous work on thin compressible elastic layers (Skotheim & Mahadevan Reference Skotheim and Mahadevan2004). In this limit, we can expand
$p(x)$
in
${\it\beta}$
as
$$\begin{eqnarray}p(x)=p_{0}(x)+{\it\beta}p_{1}(x)+O({\it\beta}^{2}).\end{eqnarray}$$
For the leading orders in
${\it\beta}$
we obtain from the lubrication equation (2.19)
$$\begin{eqnarray}\displaystyle & \displaystyle O(1):\frac{\text{d}}{\text{d}x}\left[\frac{\text{d}p_{0}}{\text{d}x}h_{0}^{3}+6h_{0}\right]=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle O({\it\beta}):\frac{\text{d}}{\text{d}x}\left[\frac{\text{d}p_{1}}{\text{d}x}h_{0}^{3}-3h_{0}^{2}\mathscr{H}\frac{\text{d}p_{0}}{\text{d}x}-6\mathscr{H}\right]=0. & \displaystyle\end{eqnarray}$$
$p_{0}(-\infty )=p_{0}(\infty )=0$
: 
$$\begin{eqnarray}p_{0}(x)=\frac{2x}{(1+x^{2})^{2}},\quad \widetilde{p}_{0}(q)=-\text{i}{\rm\pi}q\text{e}^{-|q|}.\end{eqnarray}$$
We note that the zeroth-order pressure
$p_{0}(x)$
is antisymmetric and does not contribute to the lift force. Substituting
$\widetilde{p}(q)$
by
$\widetilde{p}_{0}(q)$
in (2.18), we solve for the deformation
$$\begin{eqnarray}\mathscr{H}(x)=\int _{-\infty }^{\infty }\frac{{\rm\pi}q\text{e}^{-|q|}}{|q|}\left[\frac{\text{i}G^{\prime }(-\mathscr{T}q)+G^{\prime \prime }(-\mathscr{T}q)}{G^{\prime }(-\mathscr{T}q)^{2}+G^{\prime \prime }(-\mathscr{T}q)^{2}}\right]\text{e}^{\text{i}qx}\frac{\text{d}q}{2{\rm\pi}}.\end{eqnarray}$$
Subsequently we solve (2.22b
) for
$p_{1}(x)$
. Inspecting (2.24), one can see that the storage modulus (
$G^{\prime }$
) is associated with the imaginary part of
$\widetilde{\mathscr{H}}(q)$
which gives antisymmetric deformation in physical space. This antisymmetric part of
$\mathscr{H}(x)$
breaks the reversibility of the Stokes equation and generates a lift force. In contrast, the dissipation in the viscoelastic medium (
$G^{\prime \prime }$
) leads to the symmetric part of
$\mathscr{H}(x)$
. Given
$\mathscr{H}(x)$
, one can solve for
$p_{1}(x)$
, and calculate the lift force (per unit length) on the cylinder
$$\begin{eqnarray}L=\int _{-\infty }^{\infty }p(x)\,\text{d}x={\it\beta}\int _{-\infty }^{\infty }p_{1}(x)\,\text{d}x+O({\it\beta}^{2}).\end{eqnarray}$$
It is important to note that the functions
$p_{1}(x)$
and
$\mathscr{H}(x)$
depend only on
$\mathscr{T}$
, so that the lift scales as
$L\sim {\it\beta}$
. The proportionality factor, however, will have a subtle dependence on
$\mathscr{T}$
, and hence on the cylinder velocity. The primary goal of the analysis will be to identify this
$\mathscr{T}$
dependence for different rheological models.
3 Results
We consider viscoelastic lubrication for three different rheological models for the wall, each with one single characteristic time scale: the standard linear solid (SLS), the Kelvin–Voigt model (KV) and a power law (PL) gel. In the following, we first briefly discuss the elastic case and then introduce viscoelasticity through the three different models.
3.1 Elastic wall (
$\mathscr{T}=0$
)
For an elastic wall that adapts instantaneously to load changes, the scaled shear modulus reduces to
${\it\mu}=1$
. Solving (2.24), we obtain
$$\begin{eqnarray}\mathscr{H}(x)=-\frac{x}{(1+x^{2})}.\end{eqnarray}$$
This deformation is purely antisymmetric, just like
$p_{0}(x)$
. The first-order pressure
$p_{1}$
is obtained by solving (2.22b
) using boundary conditions
$p_{1}(-\infty )=p_{1}(\infty )=0$
,
$$\begin{eqnarray}p_{1}(x)=\frac{1-2x^{2}}{(1+x^{2})^{4}}.\end{eqnarray}$$
The lift force on the cylinder is obtained by integrating over
$p_{1}(x)$
and yields
$$\begin{eqnarray}L_{0}=\frac{3{\rm\pi}}{16}{\it\beta}.\end{eqnarray}$$
Our main interest is to study how the lift force differs from the purely elastic
$L_{0}$
for various viscoelastic models. The difference will appear when
$\mathscr{T}$
is of order unity (or larger), for which the viscoelastic solid is excited on time scales comparable to (or faster than) the relaxation time. For
$\mathscr{T}\ll 1$
the response will reduce to the elastic case, with
$L=L_{0}$
.
3.2 Standard linear solid model
The simplest rheological model for a viscoelastic solid is the superposition of an elastic material and a viscous fluid with frequency-independent viscosity. It is often represented by a spring and dashpot connected in parallel, which is the so-called KV solid. The KV model is, however, not applicable to many solids, and we start our analysis with a generalisation of it: the SLS model, which shows exponential stress relaxation (cf. (3.4a
)). In a spring–dashpot representation, a spring and a dashpot are connected in series to form one arm of the model, and a second spring is connected in parallel to this arm. The ratio of the two spring stiffnesses is given by
$c$
. Both the SLS and the KV model have equilibrium or long-time modulus, a signature of viscoelastic solid models. The extra spring in the SLS model provides an instantaneous modulus which the KV model lacks. For intermediate frequencies, viscous dissipation causes an exponential relaxation behaviour, characterised by a time scale
${\it\tau}$
. In terms of the dimensionless relaxation function
${\it\Psi}(t)$
and complex shear modulus
${\it\mu}({\it\omega})$
, the SLS model reads
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\Psi}(t)=1+c\text{e}^{-ct}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\mu}({\it\omega})=\frac{{\it\omega}^{2}+c^{2}+c{\it\omega}^{2}}{{\it\omega}^{2}+c^{2}}+\text{i}\frac{c^{2}{\it\omega}}{{\it\omega}^{2}+c^{2}}. & \displaystyle\end{eqnarray}$$
$G^{\prime }({\it\omega})=\text{Re}[{\it\mu}({\it\omega})]$
and loss modulus
$G^{\prime \prime }({\it\omega})=\text{Im}[{\it\mu}({\it\omega})]$
are plotted in figure 2(a) for
$c=100$
. One indeed observes two distinct values of
$G^{\prime }$
, at low and high frequency respectively.
(a) The storage and loss moduli of the SLS model for
$c=100$
. (b) The deformation of the SLS half-space (
$c=100$
) due to zeroth-order lubrication pressure
$p_{0}$
for
$\mathscr{T}=\{0.05,1,2,5,20,10^{4}\}$
(cf. (3.5)). Inset: the magnified deformation for
$\mathscr{T}=10^{4}$
shows that it regains the small-
$\mathscr{T}$
shape. (c) First-order pressure
$p_{1}$
for
$\mathscr{T}=\{0.05,1,2,5,20,10^{4}\}$
. Inset: the magnified pressure at
$\mathscr{T}=10^{4}$
. (d) Lift force on the cylinder as a function of
$\mathscr{T}$
for three different
$c$
values normalised by
$L_{0}$
. The red dashed line is the small-
$\mathscr{T}$
asymptotic (
$L/L_{0}=1$
). The blue dashed lines indicate the large-
$\mathscr{T}$
asymptotic given by (3.6). The grey curves represent the numerically calculated lift force using the full profile of (3.5). The green dashed curve gives the lift force for the KV model calculated using (3.9).
For the SLS model, the viscoelastic deformation is obtained in closed form by solving (2.24):
$$\begin{eqnarray}\displaystyle \mathscr{H}(x) & = & \displaystyle -\frac{x}{(1+c)(1+x^{2})}+\frac{c^{2}\displaystyle \exp \left(\frac{cx}{\mathscr{T}+c\mathscr{T}}\right)}{\mathscr{T}(1+c)^{2}}\nonumber\\ \displaystyle & & \displaystyle \times \,\text{Re}\left[\exp \left(\frac{\text{i}c}{\mathscr{T}+c\mathscr{T}}\right)\left(\text{Ei}\left[-\frac{c(\text{i}+x)}{\mathscr{T}+c\mathscr{T}}\right]+\text{i}{\rm\pi}\right)\right],\end{eqnarray}$$
where
$\text{Ei}$
is the exponential integral of a complex function
$z$
, defined as
$\text{Ei}(z)=-\int _{-z}^{\infty }(\text{e}^{-t}/t)\,\text{d}t$
. The deformation according to (3.5) is plotted in figure 2(b). At very small values of
$\mathscr{T}$
, the response is essentially elastic and
$\mathscr{H}$
is perfectly antisymmetric. Viscoelastic effects become apparent for increasing
$\mathscr{T}$
, for which the deformation decreases in amplitude and loses its perfectly antisymmetric form. However, at very high
$\mathscr{T}$
, the instantaneous elasticity of the SLS model becomes dominant. Hence, for
$\mathscr{T}\rightarrow \infty$
, one recovers the same profile as for
$\mathscr{T}=0$
, but with an amplitude reduced by a factor
$1/(1+c)$
(inset to figure 2
b).
Using (3.5), we solve (2.22b
) numerically for the first-order pressure
$p_{1}$
. Integrating
$p_{1}$
we obtain the lift force
$L$
on the cylinder. The resulting pressures are shown in figure 2(c), while we report the lift force
$L/L_{0}$
(i.e. normalised by the elastic case) in figure 2(d). As anticipated,
$L/L_{0}=1$
in the limit of small
$\mathscr{T}$
. The lift decreases upon increasing
$\mathscr{T}$
: the deformation
$\mathscr{H}$
drops in amplitude and its shape develops a symmetric component, both leading to a reduction in the generated lift force. At very large
$\mathscr{T}$
, where the SLS model responds as purely elastic with a higher modulus, the shape and pressure become independent of
$\mathscr{T}$
, and the lift force settles at a constant value. Since the effective modulus at short times is smaller by a factor
$1/(1+c)$
, the lift force in this limit becomes
$$\begin{eqnarray}L/L_{0}=(1+c)^{-1},\end{eqnarray}$$
which is indicated by the blue dashed lines in figure 2(d). All of the asymptotic results are obtained assuming
$\mathscr{T}>0$
.
Intriguingly, the numerical results in figure 2(d) suggest an intermediate asymptotic regime that emerges when
$c\gg 1$
, indicated by a dashed green curve. As we will show in the next section, this intermediate regime corresponds exactly to the KV model. It will be shown that for
$\mathscr{T}\gg 1$
, the lift force for the KV model reads
$L\simeq (2/3)|\mathscr{T}|^{-2}$
. Comparing with (3.6), we see that the intermediate asymptotics is described by
$$\begin{eqnarray}L/L_{0}={\textstyle \frac{2}{3}}|\mathscr{T}|^{-2},\quad \text{for}~1\ll \mathscr{T}\ll \sqrt{c}.\end{eqnarray}$$
This regime can indeed be observed when
$c\gg 1$
, which is naturally expected to be the case for solids that exhibit an instantaneous elasticity.
3.3 Kelvin–Voigt limit
In the limit of
$c\rightarrow \infty$
, the instantaneous relaxation of the SLS model is suppressed and one recovers the KV model. The relaxation function and the complex modulus of the KV model are given by
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\Psi}(t)=1+{\it\delta}(t), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\mu}({\it\omega})=1+\text{i}{\it\omega}. & \displaystyle\end{eqnarray}$$
${\it\omega}\ll 1$
, the KV model behaves as purely elastic, while for
${\it\omega}\gg 1$
it acts as a Newtonian fluid. The surface deformation of a KV half-space is obtained as
$c\rightarrow \infty$
in (3.5), 
$$\begin{eqnarray}\mathscr{H}(x;\mathscr{T})=\frac{\displaystyle \exp \left(\frac{x}{\mathscr{T}}\right)}{\mathscr{T}}\text{Re}\left[\exp \left(\frac{\text{i}}{\mathscr{T}}\right)\left(\text{Ei}\left[-\frac{(\text{i}+x)}{\mathscr{T}}\right]+\text{i}{\rm\pi}\right)\right].\end{eqnarray}$$
The green dashed curve in figure 2(d) gives the numerically calculated lift force, which is indeed the intermediate asymptotics of the SLS model.
We now calculate the asymptotic nature of the lift force for this material model. At large
$\mathscr{T}$
, (3.9) reduces to
$$\begin{eqnarray}\displaystyle \mathscr{H}(x;\mathscr{T}) & = & \displaystyle |\mathscr{T}|^{-1}\left({\it\gamma}+\frac{1}{2}\log (1+x^{2})\right)\nonumber\\ \displaystyle & & \displaystyle +\,|\mathscr{T}|^{-2}\left(x({\it\gamma}-1)+\frac{x}{2}\log (1+x^{2})+\tan ^{-1}x\right)+O(|\mathscr{T}|^{-3}),\end{eqnarray}$$
where
${\it\gamma}$
is Euler’s constant, which is equal to 0.577216. Integrating (2.22b
), we find
$$\begin{eqnarray}p_{1}=\int _{-\infty }^{x}3\left(\frac{1}{h_{0}}\frac{\text{d}p_{0}}{\text{d}x^{\prime }}+\frac{2}{h_{0}^{3}}\right)\mathscr{H}\,\text{d}x^{\prime }=|\mathscr{T}|^{-1}g_{1}(x)+|\mathscr{T}|^{-2}g_{2}(x).\end{eqnarray}$$
Here,
$g_{1}(x)$
is antisymmetric and makes no contribution to the lift force. The leading-order contribution to the lift force is thus
${\sim}|\mathscr{T}|^{-2}$
:
$$\begin{eqnarray}L=\int _{-\infty }^{\infty }p_{1}\,\text{d}x=|\mathscr{T}|^{-2}\int _{-\infty }^{\infty }g_{2}(x)\,\text{d}x=k|\mathscr{T}|^{-2}.\end{eqnarray}$$
Numerical integration over
$g_{2}(x)$
appears to give an exact ratio
$k={\rm\pi}/8$
up to eight decimals. As a result, the large-
$\mathscr{T}$
asymptotics for the lift force becomes
$$\begin{eqnarray}L/L_{0}={\textstyle \frac{2}{3}}|\mathscr{T}|^{-2},\end{eqnarray}$$
which is the scaling law anticipated before.
3.4 Power law rheology
Many crosslinked polymers like PDMS and polyurethane exhibit a PL relaxation behaviour, i.e. both
$G^{\prime }({\it\omega})$
and
$G^{\prime \prime }({\it\omega})$
scale as
${\it\omega}^{n}$
for large
${\it\omega}$
. In general, a large degree of polydispersity or a broad relaxation spectrum causes PL relaxation behaviour in polymeric systems (Ng & McKinley Reference Ng and McKinley2008). The value of the exponent
$n$
depends on the stoichiometry. For example,
$n=1/2$
for a stoichiometrically balanced PDMS, otherwise
$n$
varies between 1/2 and 1 (Chambon & Winter Reference Chambon and Winter1987). The rheology of such a PL gel can be modelled as
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\Psi}(t)=1+{\it\Gamma}(1-n)^{-1}\frac{1}{t^{n}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\mu}({\it\omega})=1+(\text{i}{\it\omega})^{n}, & \displaystyle\end{eqnarray}$$
${\it\Gamma}$
is the Gamma function. Equations (3.14) require
$0<n<1$
, to allow for integrability. Figure 3(a) shows the corresponding storage and loss moduli for
$n=3/4$
.
The case where
$n$
approaches unity is a singular limit, in the sense that the limit of large
${\it\omega}$
and
$n\rightarrow 1$
cannot be reversed. At a given frequency
${\it\omega}$
, the response of the PL gel approaches that of the KV rheology in the limit
$n\rightarrow 1$
, see (3.8b
). However, for any value of
$n<1$
, the high-frequency asymptotics of the PL gel is
$G^{\prime }\sim G^{\prime \prime }\sim {\it\omega}^{n}$
, while for the KV model one has
$G^{\prime }\sim {\it\omega}^{0}$
. For a given material of
$n<1$
, we therefore anticipate the lift in the high-velocity regime (
$\mathscr{T}\gg 1$
) to be different from the KV behaviour. Again, we will find that the KV model serves as an intermediate regime for the PL solid.
The Fourier transform of the deformation of the PL solid reads
$$\begin{eqnarray}\widetilde{\mathscr{H}}(q)=\frac{\text{i}{\rm\pi}\text{e}^{-|q|}q}{|q|(1+(-\text{i}\mathscr{T}q)^{n})}.\end{eqnarray}$$
The inversion to the physical space has to be performed numerically. The results are shown in figure 3(b). The inset shows a zoom of a profile of
$\mathscr{H}(x)$
for large
$\mathscr{T}$
. Unlike the SLS model, a PL solid does not exhibit an elastic response at large
$\mathscr{T}$
, and the profile is not perfectly antisymmetric. Figure 3(c) shows the corresponding pressure profiles (for
$n=0.75$
), and the lift forces are shown in figure 3(d) (grey curves). For small
$\mathscr{T}$
, the lift force is similar to the elastic case. For large
$\mathscr{T}$
,
$L$
decreases algebraically with an exponent that depends on
$n$
. The green curve corresponds to the KV model.
(a) The storage and loss moduli of the PL model for
$n=3/4$
. (b) The deformation of the PL half-space (with
$n=3/4$
) obtained by numerical backward transform of (3.15) for
$\mathscr{T}=\{0.05,1,2,5,20,10^{4}\}$
. Inset: the magnified deformation for
$\mathscr{T}=10^{4}$
shows that unlike the SLS model, the PL model does not exhibit elastic deformation for large
$\mathscr{T}$
. (c) First-order pressure
$p_{1}$
for
$\mathscr{T}=\{0.05,1,2,5,20,10^{4}\}$
. Inset: the magnified pressure profile at
$\mathscr{T}=10^{4}$
. (d) Lift force on the cylinder as a function of
$\mathscr{T}$
for four different
$n$
values. The red dashed line is the small-
$\mathscr{T}$
approximation given by (3.3). The blue dashed lines indicate the large-
$\mathscr{T}$
asymptotic given by (3.19). The grey curves represent the numerically calculated lift force using the full profile of (3.15). The green dashed curve gives the lift force for the KV model. Inset: the data points show the numerically calculated prefactor
$f(n)$
as
$n$
varies between 0 and 1.
We now extract the large-
$\mathscr{T}$
behaviour of
$L$
, and in particular its dependence on the rheological exponent
$n$
. We expand (3.15), and to leading order obtain
$$\begin{eqnarray}\widetilde{\mathscr{H}}(q)=-\frac{\text{e}^{-|q|}{\rm\pi}(-\text{i}q)^{1-n}|\mathscr{T}|^{-n}}{|q|},\end{eqnarray}$$
which can be inverted to physical space as
$$\begin{eqnarray}\mathscr{H}(x;\mathscr{T})=-{\it\Gamma}(1-n)\text{Re}\left[\frac{\text{i}^{-n}(1+\text{i}x)^{n}}{-\text{i}+x}\right]|\mathscr{T}|^{-n}.\end{eqnarray}$$
Using (3.11) we can write
$$\begin{eqnarray}p_{1}=|\mathscr{T}|^{-n}\int _{-\infty }^{x}-3{\it\Gamma}(1-n)\text{Re}\left[\frac{\text{i}^{-n}(1+\text{i}x^{\prime })^{n}}{-\text{i}+x^{\prime }}\right]\left(\frac{1}{h_{0}}\frac{\text{d}p_{0}}{\text{d}x^{\prime }}+\frac{2}{h_{0}^{3}}\right)\,\text{d}x^{\prime }=|\mathscr{T}|^{-n}f_{1}(x;n)\end{eqnarray}$$
and the lift force
$$\begin{eqnarray}L/L_{0}=\frac{16}{3{\rm\pi}}\int _{-\infty }^{\infty }|\mathscr{T}|^{-n}f_{1}(x;n)\,\text{d}x=|\mathscr{T}|^{-n}f(n).\end{eqnarray}$$
Hence, we find that
$L\sim |\mathscr{T}|^{-n}$
. The prefactor
$f(n)$
is evaluated numerically and is shown in the inset of figure 3(d). The blue dashed lines in figure 3(d) show the large-
$\mathscr{T}$
asymptotes, in agreement with the full numerical evaluation.
The KV model again serves as an intermediate regime (green curve), which now emerges as
$n\rightarrow 1$
. Comparing (3.19) and (3.13), we find the scaling
$$\begin{eqnarray}L/L_{0}={\textstyle \frac{2}{3}}|\mathscr{T}|^{-2},\quad \text{for}~1\ll \mathscr{T}\ll (1-n)^{-1},\end{eqnarray}$$
where we used that
$f\sim (1-n)$
near
$n=1$
. The upper bound reflects that for
$0<1-n\ll 1$
, the PL gel does not converge to the KV model at very high frequencies.
4 Discussion
We have analysed how the mechanics of lubricated contacts is affected by the viscoelastic properties of the solid. We focused on two-dimensional cylindrical contacts in the limit of small deformations, and considered several different rheological models for the lubricated solid. Here, we briefly summarise the key findings, presented in physical units, and discuss how the results are generalised to the arbitrary form of
${\it\mu}({\it\omega})=G^{\prime }({\it\omega})+\text{i}G^{\prime \prime }({\it\omega})$
.
At low lubrication velocities, all rheologies with a non-vanishing static modulus
$G=G^{\prime }({\it\omega}=0)$
give rise to a purely elastic response. In this case, the lift force (per unit length) becomes
$$\begin{eqnarray}L_{0}=\frac{3{\rm\pi}{\it\eta}^{2}R^{2}}{8G{\it\Delta}^{3}}|V|^{2}.\end{eqnarray}$$
For the Kelvin–Voigt solid, the large-velocity asymptote reads
$$\begin{eqnarray}L=\frac{{\rm\pi}{\it\eta}^{2}R^{3}}{2G{\it\tau}^{2}{\it\Delta}^{2}},\end{eqnarray}$$
which interestingly corresponds to a lift that is independent of velocity. For the PL model we find the scaling law
$$\begin{eqnarray}L\sim \frac{{\it\eta}^{2}R^{2}\ell ^{n}}{G{\it\Delta}^{3}{\it\tau}^{n}}|V|^{2-n},\end{eqnarray}$$
with a prefactor that depends on the rheological exponent
$n$
.
The effect of viscoelasticity is twofold: (i) the resulting deformation is reduced in amplitude with respect to the purely elastic case, and (ii) the deformation profile develops a symmetric part. The former gives a reduction in the lift force, but the details of this reduction depend on the rheological model. It is of interest to generalise these findings to arbitrary rheology. One can identify the relevant scale of the solid deformation upon inspection of (2.24), bearing in mind that only the antisymmetric deformation contributes to the lift. From this, we derive that the lift scales as
$$\begin{eqnarray}L\sim \frac{{\it\eta}^{2}R^{2}}{{\it\Delta}^{3}}\frac{G^{\prime }(|V|/\ell )}{G^{\prime }(|V|/\ell )^{2}+G^{\prime \prime }(|V|/\ell )^{2}}|V|^{2},\end{eqnarray}$$
which is indeed consistent with all of the scaling laws mentioned above. Clearly, this implies a reduction of the lift force with respect to the elastic response (4.1), whenever the rheology contains a significant contribution of the loss modulus. This expression also highlights the importance of both the storage and a loss modulus to determine the hydrodynamic lift in viscoelastic lubrication: neither a vanishing nor an infinite
$G^{\prime }$
will lead to lift.
The presented formulation may be considered as a rheological tool. Indeed, lubrication has recently been exploited for in situ atomic force microscopy measurements of elastic properties of thin films at the microscale (Leroy & Charlaix Reference Leroy and Charlaix2011; Leroy et al.
Reference Leroy, Steinberger, Cottin-Bizonne, Restagno, Léger and Charlaix2012; Wang, Dhong & Frechette Reference Wang, Dhong and Frechette2015). Similar experiments can be performed to measure the lift force on a particle moving steadily in parallel to a soft boundary. The large-
$V$
limit of the velocity–lift data thus obtained may show the viscoelastic signature of the underlying substrate according to (4.4). In principle, the entire range of velocity–lift data provides access to the full rheological spectrum of very soft layers. From a theoretical point of view, interesting future directions would be to consider the effect of viscoelasticity on freely moving particles and to release the condition of small deformation.
Acknowledgements
We thank B. Andreotti and T. Salez for discussions, and two anonymous reviewers for comments. S.K. acknowledges financial support from NWO through VIDI grant no. 11304. A.P. and J.S. acknowledge financial support from ERC (the European Research Council) Consolidator grant no. 616918.
