Chapter 10 Lubrication Theory

Let’s first consider the squeeze flow in a thrust bearing. Figure 10.3 shows the geometry and coordinate system where the bearing remains stationary and the plane underneath the fluid moves with horizontal velocity $-U$. The fluid film wedge generates a hydrodynamic pressure that supports an applied load on the bearing with a small tangential force for good lubrication. We would like to calculate the pressure in the fluid and hence the forces exerted on the bearing by the fluid. To solve the lubrication problem we: describe the geometry, solve for (almost) unidirectional flow, apply mass conservation to close the system, and calculate quantities of interest in this case the forces on the bearing.

The second squeeze flow problem we will consider is a cylinder approaching a wall. We will carry out a similar procedure to the previous problem. Figure 10.5 shows the geometry of the problem. A cylinder of radius $a$ falls vertically towards a horizontal plane with minimum gap $d(t)$ between the cylinder and the plane, where $d\ll a$ in the thin film limit. There are non-slip boundary conditions on the cylinder and on the plane.

Surface tension is a property of a free surface which appears due to the asymmetry of intermolecular forces acting in the interfacial region (typically of nanometric width). As one goes to smaller length scales, surface forces ($\sim L^2$) begin to dominate volume forces ($\sim L^3$). In particular, the surface tension force typically begins to dominate gravity when the length scale goes below millimetres-cm.

In the general case, for a free surface between a viscous liquid and a second viscous fluid (gas or liquid), if we assume that the velocity is continuous across the interface, there are 4 boundary conditions required for the 4 unknowns (3 velocities, which are not fixed, and 1 variable representing the free surface shape). These are provided by kinematic and dynamic boundary conditions.

E.g. For a two-dimensional flow, suppose the free surface is defined as $y=h(x,t)$. Then $f(x,y,t)=y-h(x,t)$. At the free surface we have $f=0$ and a unit normal vector $𝐧=(-{\partial }_{x}h,1)/\sqrt{1+{({\partial }_{x}h)}^2}$ which point out of fluid 1. For $f(x,y,t)=y-h(x,t)$ so that $\nabla f=(-{\partial }_{x}h,1)$ we find that

$$-{\partial }_{t}h+(u{\partial }_x+v{\partial }_y)(y-h(x,t))=0\text{so that}{\partial }_{t}h+u{\partial }_{x}h=v.$$

(10.34)

To calculate the forces from inside the interface we will be interested in surface tension. Surface tension acts like an elastic membrane on the surface of a liquid and forces it to attempt to minimise its area (crudely, this is because it is energetically favourable for molecules to be in the bulk of the liquid, rather than at its surface). This force is responsible for a whole host of effects ranging from the suspension of small objects (e.g. water striders) through to the breakup of liquid jets (e.g. from a tap in a kitchen), which occurs when the liquid can reduce its surface area by forming a series of spheres (drops) rather than one cylinder (the jet), see figure 10.6. Problems involving surface tension are known as capillary flows.

The surface tension force $\gamma$ (taken to be constant here), per unit length of line, is directed tangentially along the surface which is parameterised by coordinates $(s_1,s_2)$ with tangent vectors in these directions $({𝐭}_1,{𝐭}_2)$. Our control volume will have ‘lower edge’ at $(s_{10},s_{20})$ and now has a surface tension force tugging at all four ends in the tangential direction.

Consider a control volume of length and width along the interface equal to $\delta s$ whose height is negligible compared to $\delta s$, see figure 10.7. This gives a force balance of

$$\delta s^2{𝐓}^{+}\cdot 𝐧+\delta s^2{𝐓}^{-}\cdot (-𝐧)+\gamma \delta s\left({𝐭}_1(s_{10}+\delta s,s_2)-{𝐭}_1(s_{10},s_2)\right)+\gamma \delta s\left({𝐭}_2(s_1,s_{20}+\delta s)-{𝐭}_2(s_1,s_{20})\right)=\mathrm{𝟎}$$

at $f=0$. Dividing through by $\delta s^2$ and taking the limit $\delta s\to 0$ we find that

$\left({𝐓}^{+}-{𝐓}^{-}\right)\cdot 𝐧={\left[𝐓\cdot 𝐧\right]}_{-}^{+}=-\gamma \left({\displaystyle \frac{\partial {𝐭}_1}{\partial s_1}}+{\displaystyle \frac{\partial {𝐭}_2}{\partial s_2}}\right)\equiv \gamma \kappa 𝐧\text{at}f=0.$

(10.35)

where we have defined the curvature of the interface $\kappa$ using

$\kappa 𝐧=-\left({\displaystyle \frac{\partial {𝐭}_1}{\partial s_1}}+{\displaystyle \frac{\partial {𝐭}_2}{\partial s_2}}\right)\text{or}\kappa =[(𝑰-\mathrm{𝐧𝐧})\cdot \nabla ]\cdot 𝐧\equiv {\nabla }_s\cdot 𝐧.$

(10.36)

To arrive at this expression we know that ${𝐭}_1\cdot 𝐧=0$ and ${𝐭}_2\cdot 𝐧=0$. Differentiating with respect to interface coordinates $(s_1,s_2)$ gives

${\displaystyle \frac{\partial {𝐭}_1}{\partial s_1}}\cdot 𝐧+{𝐭}_1\cdot {\displaystyle \frac{\partial 𝐧}{\partial s_1}}=0\mathrm{and}{\displaystyle \frac{\partial {𝐭}_2}{\partial s_2}}\cdot 𝐧+{𝐭}_2\cdot {\displaystyle \frac{\partial 𝐧}{\partial s_2}}=0$

(10.37)

$\Rightarrow \kappa =-\left({\displaystyle \frac{\partial {𝐭}_1}{\partial s_1}}+{\displaystyle \frac{\partial {𝐭}_2}{\partial s_2}}\right)\cdot 𝐧={𝐭}_1\cdot {\displaystyle \frac{\partial 𝐧}{\partial s_1}}+{𝐭}_2\cdot {\displaystyle \frac{\partial 𝐧}{\partial s_2}}={\nabla }_s\cdot 𝐧.$

(10.38)

This can be simplified further by noting that $𝐧\cdot (𝐧\cdot \nabla )𝐧=0$. Hence, the curvature can be written as

$\kappa ={\nabla }_s\cdot 𝐧=\nabla \cdot 𝐧.$

(10.39)

We have seen from equation (10.35) that the stresses tangential to the interface are continuous, whereas the normal stress has a jump caused by the surface tension force $\gamma \kappa$. Hence, the pressure generated by surface tension is of the order $\gamma /L$, where $L$ is a characteristic lengthscale in the problem which is approximately the radius of curvature of our interface. Therefore, capillary forces will be important when this is comparable to the pressure generated by gravitational forces, i.e. the hydrostatic pressure $\rho gL$ which occurs when

$Bo={\displaystyle \frac{\rho g{L}^2}{\gamma }}=1\text{i.e. when}{L}_{\gamma }=\sqrt{\gamma /(\rho g)},$

(10.40)

where $L_{\gamma }$ is known as the capillary length and $Bo$ is the Bond number. For water in air, where $\gamma =0.07$ N/m, we have $L_{\gamma }=3$ mm. Therefore, for large scale flows, such as sea waves, surface tension has no effect, but for small scale flows, such as those encountered in micro and nanofluidics, its importance far outweighs that of gravity.

In many cases, we are interested in knowing the shape of a liquid drop or bubble whose flow has negligible force in comparison to the surface tension forces. In the latter case, for low Reynolds number flows, where viscous forces dominate inertial ones, looking at the dynamic boundary condition we can see that the viscous force has magnitude $\mu U/L$ whilst the surface tension force is $\gamma /L$. The ratio of these forces gives the capillary number:

$Ca={\displaystyle \frac{\mu U}{\gamma }}$

(10.41)

and when $Ca\ll 1$ the flow within the drop has a negligible effect on the free surface, so that the shape can be calculated independently of the flow (a huge simplification).

Let’s consider a liquid drop sat on a horizontal surface. For simplicity, let’s consider a two-dimensional drop pinned at $x=\pm X$ where $h(\pm X,t)=0$, see figure 10.8. In the limit of small capillary number $Ca\ll 1$, we can neglect the stress in the drop and calculate the shape of the drop from the dynamic boundary condition:

$p^{-}-p^{+}=\gamma \nabla \cdot 𝐧\text{at}f=0.$

(10.42)

This is known as the Young-Laplace equation. We can parametrise the free surface in terms of a height function $y=h(x)$, so that $f=y-h(x,t)$, see figure 10.8. Then

$𝐧={\displaystyle \frac{(-h_x,1)}{\sqrt{(}1+h_x^2)}}\text{and}{\nabla }_{\cdot }𝐧={\displaystyle \frac{\partial n_x}{\partial x}}+{\displaystyle \frac{\partial n_y}{\partial y}}=-{\displaystyle \frac{h_{xx}}{{(1+h_x^2)}^{3/2}}}.$

(10.43)

From our $y$-momentum equation we find the pressure in the fluid is hydrostatic $p^{-}=p_0^{-}-\rho gy$. We assume the ambient fluid has a uniform pressure $p^{+}=p_0^{+}$. Substituting these into the dynamic boundary condition (10.42) gives

$-P-\rho gh=\gamma \nabla \cdot 𝐧\text{at}f=0.$

(10.44)

where $P=p_0^{+}-p_0^{-}$ is the (constant) pressure drop across the interface. Hence,

$-P-\rho gh=-{\displaystyle \frac{\gamma h_{xx}}{{(1+h_x^2)}^{3/2}}}\text{at}y=h(x,t).$

(10.45)

This is a second order PDE for $h$ and hence we must specify a boundary condition on $h$ or its derivative at each end of the free surface. On top of this one either needs to specify the pressure drop $P$ or, more often, the volume of the liquid which can then be used to find $P$.

In the lubrication limit $h_x\ll 1$ for small bond numbers $Bo\ll 1$, this reduces to

$P=\gamma {\displaystyle \frac{d^{2}h}{d{x}^2}}$

(10.46)

with solution a parabola

$h={\displaystyle \frac{-P}{2\gamma }}\left(X^2-x^2\right),$

(10.47)

where we have the used the pinned boundary condition $h(\pm X,t)=0$. We see that for meaningful solutions we need as $P$ represents the pressure in the gas minus the pressure in the liquid (which we expect to be higher for a drop, due to surface tension). If we calculate the volume of the drop $V$ we can see how the pressure and volume of the drop are related:

$V=2{\displaystyle {\int }_0^X}{\displaystyle \frac{-P}{2\gamma }}\left(X^2-x^2\right)dx={\displaystyle \frac{-2P{X}^3}{3\gamma }}\Rightarrow -P={\displaystyle \frac{3\gamma V}{2X^3}}.$

(10.48)

In this solution we assumed the drop was pinned at points $x=\pm X$. This is known as a pinned contact line. An alternative boundary condition would be to prescribe the contact angle $\theta$, the angle at which the free surface meets the horizontal plane, see figure 10.8. This then gives a boundary condition on the gradient $h_x=\pm \tan \theta$ at $x=\pm X$. In this case the position of the contact line would be found as part of the solution.

The prescribed contact angle will depend on the liquid-solid-gas combination in question: if the solid is ‘wettable’ (aka hydrophilic for water) this angle will be small (i.e. close to zero), so the drop spreads out a lot (e.g. water on glass), whereas if the angle is high (close to ${180}^{\circ }$) then the solid is non-wettable (aka hydrophobic) and the drop beads up (e.g. water on a lotus leaf).

When the lengthscale of the liquid drop becomes large, gravitational forces begin to dominate over capillary forces (large Bond number) and viscous forces begin to dominate over surface tension (large capillary number). We can therefore neglect the jump in the normal stress at the free surface but now need to calculate the flow within the fluid. Let’s consider the shape of an axisymmetric drop on a horizontal surface under the influence of gravity, or equivalently, let’s consider the axisymmetric, gravitational spreading of a thin film layer on a horizontal surface, see figure 10.9.

Finally, integrating the momentum equation and applying the boundary conditions gives radial velocity profile

$u=-{\displaystyle \frac{\rho g}{2\mu }}{\displaystyle \frac{\partial h}{\partial r}}z(2h-z).$

(10.54)

We would like to solve the following PDE (non-linear diffusion equation) subject to the global constraint that the volume is constant:

${\displaystyle \frac{\partial h}{\partial t}}={\displaystyle \frac{g}{3\nu }}{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial }{\partial r}}\left(r{h}^3{\displaystyle \frac{\partial h}{\partial r}}\right),h(R,t)=0,\mathrm{and}V=2\pi {\displaystyle {\int }_0^{R(t)}}rhdr.$

(10.61)

Scaling the equation of motion and global mass conservation using characteristic height, radial and time scales $H,R$ and $T$ gives

${\displaystyle \frac{H}{T}}\sim {\displaystyle \frac{g{H}^4}{\nu R^2}}\mathrm{and}V\sim h{R}^2,$

(10.62)

which can be rearranged to

$H\sim {\left({\displaystyle \frac{\nu V}{gT}}\right)}^{1/4}\mathrm{and}R\sim {\left({\displaystyle \frac{g{V}^{3}T}{\nu }}\right)}^{1/8}.$

(10.63)

Since there are no more scales in the problem (assuming the initial condition was a sufficiently long ago to be irrelevant) it suggests we such look for a similarity solution of the form

$h(r,t)={\left(\nu V/gt\right)}^{1/4}H(\eta )\mathrm{where}\eta \equiv r/{\left(g{V}^{3}t/\nu \right)}^{1/8}$

(10.64)

with the edge of the film at ${\eta }_N=R/{\left(g{V}^{3}t/\nu \right)}^{1/8}$. These solutions are self-similar which means that by appropriately scaling the radius and height of the film by a power of $t$ the time evolution of the current falls onto a universal curve (we’ll plot the solution up later to show this more clearly). Let’s substitute this ansatz into the PDE and constraints. With some algebra with arrive at:

$-{\displaystyle \frac{1}{4}}H-{\displaystyle \frac{1}{8}}\eta H^{\prime }={\displaystyle \frac{1}{3\eta }}{(\eta H^3{H}^{\prime })}^{\prime },H({\eta }_N)=0,\mathrm{and}1=2\pi {\displaystyle {\int }_0^{{\eta }_N}}\eta Hd\eta ,$

(10.65)

where the prime denotes differentiation w.r.t $\eta$. The power of this type of solution is that we’ve turned a PDE into an ODE for $H(\eta )$ which is significantly easier to solve. We can now make progress with the ODE as follows:

	$-{\displaystyle \frac{1}{8}}\left(2\eta H+{\eta }^2{H}^{\prime }\right)$	$={\displaystyle \frac{1}{3}}{\left(\eta H^3{H}^{\prime }\right)}^{\prime }$		(10.66)
	$-{\displaystyle \frac{1}{8}}{\left({\eta }^{2}H\right)}^{\prime }$	$={\displaystyle \frac{1}{3}}{\left(\eta H^3{H}^{\prime }\right)}^{\prime }$		(10.67)
	$-{\displaystyle \frac{1}{8}}{\eta }^{2}H$	$={\displaystyle \frac{1}{3}}\eta H^3{H}^{\prime }+\underset{=0}{\underbrace{\mathrm{const}.}}$		(10.68)

since we want $H({\eta }_N)=0$. Integrating again we get

${\displaystyle \frac{1}{9}}{H}^3={\displaystyle \frac{1}{16}}({\eta }_N^2-{\eta }^2).$

(10.69)

Finally, applying our volume constraint we can determine ${\eta }_N$ and hence solution

$H={\left({\displaystyle \frac{9}{16}}({\eta }_N^2-{\eta }^2)\right)}^{1/3}\mathrm{where}{\eta }_N={\left({\displaystyle \frac{2^{10}}{3^5{\pi }^3}}\right)}^{1/8},$

(10.70)

or equivalently

$h={\left({\displaystyle \frac{\nu V}{gt}}\right)}^{1/4}{\left({\displaystyle \frac{9{\eta }_N^2}{16}}\left(1-{\displaystyle \frac{r^2}{R^2}}\right)\right)}^{1/3},R(t)={\eta }_N{\left({\displaystyle \frac{g{V}^{3}t}{\nu }}\right)}^{1/8}.$

(10.71)

Figure 10.11 plots the solution to demonstrate the self-similar behaviour.

To investigate the shape of the nose we let $h=A{(R-r)}^{\alpha }$. Substituting into the governing equation when $r\simeq R$ gives

	$A\alpha {(R-r)}^{\alpha -1}\dot{R}$	$={\displaystyle \frac{g}{3\nu }}{\displaystyle \frac{1}{R}}{\displaystyle \frac{\partial }{\partial r}}\left(R{A}^3{(R-r)}^{3\alpha }A\alpha (-){(R-r)}^{\alpha -1}\right)$		(10.72)
		$={\displaystyle \frac{g}{3\nu }}{A}^4\alpha (4\alpha -1){(R-r)}^{4\alpha -2}.$		(10.73)

Matching the powers of $R-r$ gives $\alpha =1/3$. The shape of the nose is then $h=A{(R-r)}^{1/3}$. The gradient $\partial h/\partial r$ is therefore singular at the nose, but not too singular that the flux boundary condition $Q(r\to R)\to 0$ is still satisfied.

Let’s now consider the shape of a two-dimensional drop spreading down an inclined plane of slope $\theta$. $x$ is the downslope coordinate and $y$ is the vertical coordinate perpendicular to the slope. The body force is then $𝐟=\rho g(\sin \theta ,-\cos \theta )$. The fluid is released from a fixed point such that $h(0,t)=0$ and the ambient has constant pressure $p_0$.

The Lubrication approximation is then

$0=-{\displaystyle \frac{\partial p}{\partial x}}+\mu {\displaystyle \frac{{\partial }^{2}u}{\partial y^2}}+\rho g\sin \theta ,0=-{\displaystyle \frac{\partial p}{\partial y}}-\rho g\cos \theta .$

(10.74)

Integrating the y-momentum equation gives hydrostatic pressure $p=p_0+\rho g\cos \theta (h-y)$. Calculating the velocity field, integrating across the depth to get the volume flux, and applying mass conservation, we arrive at evolution equation

${\displaystyle \frac{\partial h}{\partial t}}+{\displaystyle \frac{g\sin \theta }{3\nu }}{\displaystyle \frac{\partial h^3}{\partial x}}={\displaystyle \frac{g\cos \theta }{3\nu }}{\displaystyle \frac{\partial }{\partial x}}\left(h^3{\displaystyle \frac{\partial h}{\partial x}}\right).$

(10.75)

A similarity solution does not exist for the full governing equation. Instead we will consider one particular limit.