2.1. Spheres widely separated

For $\delta _0\gg 1$, the governing equations for the fluid are the Stokes equation and continuity

(2.2a,b)\begin{equation} \boldsymbol{\nabla} p = \mu \nabla^{2} \boldsymbol{u}, \quad \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} = 0, \end{equation}

where $\boldsymbol{u}$ is the fluid velocity and p is the pressure.

Figure 1 shows the cylindrical coordinate system that will be used to analyse the flow in the gap. In cylindrical coordinates ($r,\theta ,z$), the dimensional gap separation $h(r)$ with minimum gap separation $h_0$ is given by

(2.3)\begin{equation} h(r) = h_0 + \left(a_S-\sqrt{a_S^{2}-r^{2}}\right) + \left(a_P-\sqrt{a_P^{2}-r^{2}}\right). \end{equation}

For $\varepsilon \ge O(1)$, the disturbance velocity field of a sphere can be written as a multipole expansion about its centre (Kim & Karrila Reference Kim and Karrila1991). The leading-order terms for the far field can be calculated by the method of reflections, but this becomes computationally expensive, as the number of terms required for a given accuracy increases with decreasing separation. Jeffrey & Onishi (Reference Jeffrey and Onishi1984) developed the twin multipole expansion method that expresses Lamb's solution expanded with respect to two different origins. The result is a tractable means for computing the HIs of two spheres of arbitrary size ratio down to separations as small as $\varepsilon \gtrsim 0.1$.

2.2. Spheres in near contact

In the limit $\varepsilon \ll 1$ and $\delta _0\gg 1$, the continuum flow in the gap obeys the Stokes equation and can be analysed using a matched asymptotic analysis. O'Neill & Stewartson (Reference O'Neill and Stewartson1967) identified an ‘inner region’, where they applied the standard lubrication approximation, and an ‘outer region’, where the full Stokes equations were used. The solutions are equated at a radial distance $r_{M}$. Figure 2 depicts the ‘inner’ and ‘outer’ regions. In the inner region, we are interested in the limit $r\ll a_P$ and $r\ll a_S$, allowing us to introduce a Taylor series expansion for the two square roots in (2.3), yielding

(2.4)\begin{align} h(r) = h_0 + a_S\left[1- \left(1-\frac{1}{2a_S^{2}}r^{2}-\frac{1}{8a_S^{4}}r^{4}\ldots\right) \right] \!+ a_P\left[1\!- \left(1\!-\frac{1}{2a_P^{2}}r^{2}\!-\frac{1}{8a_P^{4}}r^{4}\ldots\right) \right]. \end{align}

Rearranging (2.4) and introducing (2.1) yields

(2.5)\begin{equation} \frac{h(r)}{h_0} =1+ \frac{1}{2} \left( \frac{1}{ah_0} \right)r^{2} + \frac{1}{8} \left( \frac{1}{ah_0} \right)^{2} r^{4}\ldots. \end{equation}

The independent variables are made dimensionless as

(2.6a,b)\begin{equation} R\equiv\frac{r}{\sqrt{a h_0}} ,\quad Z\equiv\frac{{z}}{h_0} . \end{equation}

To leading order, the dimensionless gap of (2.5) is expressed as

(2.7)\begin{equation} H(R)\equiv\frac{h(r)}{h_0} = 1 + \frac{1}{2}R^{2}. \end{equation}

The equations defining the surfaces of the two spheres are then defined as

(2.8)\begin{equation} \left. \begin{gathered} Z_S = 1 + \tfrac{1}{2}\alpha R^{2}, \\ Z_P ={-}\tfrac{1}{2}(1-\alpha)R^{2}, \end{gathered} \right\} \end{equation}

where $\alpha$ is the size ratio

(2.9)\begin{equation} \alpha\equiv\frac{a}{a_S} , \end{equation}

and consequently $\alpha \rightarrow 1$ for $a_P/a_S\rightarrow \infty$ and $\alpha \rightarrow 0$ for $a_S/a_P\rightarrow \infty$.

Figure 2.

Sketch of the inner region defined as $0<r<r_{M}$ and outer region defined as $r>r_{M}$.

The velocity components and pressure can be written as (O'Neill & Stewartson Reference O'Neill and Stewartson1967)

(2.10)\begin{equation} \left. \begin{gathered} u_r(r,\theta,{z}) = U_S\; U(R,Z)\cos\theta, \\ u_\theta(r,\theta,{z}) = U_S\; V(R,Z)\sin\theta, \\ u_z(r,\theta,{z}) = U_S\; W(R,Z)\cos\theta, \\ p(r,\theta,{z}) = \frac{\mu U_S}{a} P(R,Z)\cos\theta, \end{gathered} \right\} \end{equation}

where U, V and W are the dimensionless velocity components and P is the dimensionless pressure. The resulting dimensionless governing equations become

(2.11)\begin{equation} \left. \begin{gathered} \varepsilon^{3/2} \frac{\partial P}{\partial R} ={-}\varepsilon \left[\frac{2(U+V)}{R^{2}} +\frac{\partial^{2} U}{\partial R^{2}} \right] + \frac{\partial^{2} U}{\partial Z^{2}},\\ -\varepsilon^{3/2}\frac{P}{R} ={-}\varepsilon\frac{ 2(U+V)}{R^{2}} + \frac{\partial^{2} V}{\partial Z^{2}},\\ \varepsilon^{3/2}\frac{\partial P}{\partial Z} = \varepsilon^{3/2} \left[ \frac{\partial^{2} W}{\partial R^{2}} +\frac{1}{R}\frac{\partial W}{\partial R} -\frac{W}{R^{2}} \right] + \varepsilon^{1/2} \frac{\partial^{2} W}{\partial Z^{2}} \end{gathered} \right\} \end{equation}

and the continuity equation is

(2.12)\begin{equation} \frac{\partial U}{\partial R} + \frac{U+V}{R} + \frac{1}{\varepsilon^{1/2}}\frac{\partial W}{\partial Z} = 0. \end{equation}

Equation (2.12) implies that both $U$ and $V$ are $O(1)$ and $W$ is order $\varepsilon ^{1/2}$. It follows from (2.11) that the pressure, $P$, scales as $\varepsilon ^{-3/2}$. A regular expansion in $\varepsilon$ is then assumed to take the form

(2.13)\begin{equation} \left. \begin{gathered} U(R,{Z}) = U_0(R,{Z}) + \varepsilon U_1(R,{Z}) + \ldots,\\ V(R,{Z}) = V_0(R,{Z}) + \varepsilon V_1(R,{Z}) + \ldots,\\ W(R,{Z}) = \varepsilon^{1/2} W_0(R,{Z}) + \varepsilon^{3/2} W_1(R,{Z}) + \ldots,\\ P(R,{Z}) = \varepsilon^{{-}3/2} P_0(R,{Z}) + \varepsilon^{{-}1/2} P_1(R,{Z}) + \ldots.\end{gathered} \right\} \end{equation}

To leading order, the governing equations (2.11) and (2.12) become

(2.14a–c)\begin{gather} \frac{\partial P_0}{\partial R} = \frac{\partial^{2} U_0}{\partial Z^{2}} ,\quad -\frac{P_0}{R} = \frac{\partial^{2} V_0}{\partial Z^{2}} , \quad \frac{\partial P_0}{\partial Z} = 0, \end{gather}

(2.15)\begin{gather} \frac{\partial U_0}{\partial R} + \frac{U_0+V_0}{R} + \frac{\partial W_0}{\partial Z} = 0. \end{gather}

3.1. Decomposition of the lubrication flow

Under the standard lubrication approximation, the local flow between the two spheres can be represented as a superposition of Poiseuille, Couette and plug flows, as depicted in figure 3. Poiseuille flow is characterized by the pressure gradient driving the flow, whereas Couette flow is characterized by the shear rate between the two surfaces. Plug flow assumes a uniform velocity across the channel. Its role is to allow arbitrary boundary conditions to be satisfied. Poiseuille and plug flows contribute to the mass flux through the gap while Poiseuille and Couette flows contribute to the shear stress on the sphere surfaces. The weights associated with the Couette and plug flows are determined by the boundary conditions.

Figure 3.

Sketch of the inner region in the xz plane decomposed into three planar flows: Couette, Poiseuille and plug flows. The Poiseuille flow depends upon the pressure gradient, the mass flux is determined by the plug and Poiseuille flows, and the shear stresses are obtained from the Poiseuille and Couette flows.

Recall this study is concerned with the forces and torques associated with tangential motions of the spheres. In the Cartesian coordinate system shown in figure 1, the tangential motion between the spheres arises from a combination of translation of sphere S in a direction perpendicular to the line of centres, given by ($U_S,0,0$), and its rotation along an axis perpendicular to the line of centres, defined as ($0,\varOmega a_S,0$). The resulting ${x}$ velocity of the surface of sphere S, $U_{wall}$ (see figure 3), can be approximated in the inner region as

(3.1)\begin{equation} U_{{wall}} = U_{{S}} (1 - \omega), \end{equation}

where

(3.2)\begin{equation} \omega = \frac{ \varOmega a_S}{U_{{S}}}. \end{equation}

Below, we show how the boundary conditions are used to obtain the weighting of Couette and plug flows. However, first we must account for the non-continuum Poiseuille and Couette flows in the gap.

3.1.1. Poiseuille flow

Cercignani & Daneri (Reference Cercignani and Daneri1963) analysed the pressure-driven flow of a rarefied gas between two parallel plates using the BGK model (Bhatnagar et al. Reference Bhatnagar, Gross and Krook1954). For consistency with the analysis for spheres given below, we designate the separation distance between the plates $h$ and the dimensionless separation, $\delta = h/\lambda$. The volumetric flux, defined as the average velocity integrated across the channel, is proportional to the average pressure gradient in all regimes, i.e. the pseudo-free molecular regime ($\delta \ll 1$), transition ($\delta = O(1$)) and continuum regime ($\delta \gg 1$). The general dimensionless expression is given by

(3.3)\begin{equation} \int_0^{1} U(Z)\,\textrm{d}Z ={-} \frac{g(\delta)}{2\delta} \frac{\textrm{d}P}{\textrm{d}X}, \end{equation}

where we apply the same scaling as described in § 2.2 and $g(\delta )$ is a dimensionless function that accounts for the transitions through the different regimes. The limit $\delta \ll 1$ corresponds to the pseudo-free molecular regime, where $g(\delta )$ has an analytical form. Likewise in the limit $\delta \gg 1$, the solution approaches the Hagen Poiseuille law with a slip coefficient, $k_1= 1.016$ (see § 3.4 for a definition and more complete explanation of the slip coefficient). In the transition regime, the values were numerically evaluated and tabulated by Cercignani & Daneri (Reference Cercignani and Daneri1963). The resulting expression for $g(\delta )$ is

(3.4)\begin{equation} g(\delta)= \left\{ \begin{array}{@{}ll} -\dfrac{1}{\sqrt{\rm \pi}}(-\ln\delta +0.8), & (\delta<0.01),\\ V(\delta), & (0.01<\delta<10),\\ \dfrac{\delta}{6}+ k_1, & (\delta>10), \end{array} \right. \end{equation}

where $\delta < 0.01$ corresponds to the pseudo-free molecular regime, $\delta > 10$ corresponds to the continuum flow regime and $V(\delta )$ represents a curve fit of the tabulated results from Cercignani & Daneri (Reference Cercignani and Daneri1963). Figure 4 shows a plot of $g(\delta )$. The markers ($\bullet$) are the tabulated values, the solid line is the cubic spline fit, and the limiting solutions for the pseudo-free molecular and continuum regimes are the dashed and dash-dot lines, respectively.

Figure 4.

The $g(\delta )$ function (3.4) showing the tabulated results of Cercignani & Daneri (Reference Cercignani and Daneri1963) ($\bullet$), the cubic spline fit used to interpolate between the points (line) and the limiting solutions for the pseudo-free molecular (dashed line) and slip (dash-dot line) regimes.

3.1.2. Couette flow

In the case of planar Couette flow, the variable of interest is the shear stress in terms of the velocity difference between the top and bottom walls. Linearity and symmetry imply that the normal stress differences in low-Mach-number non-continuum planar Couette flow are zero, just as they are for low-Reynolds-number continuum Couette flow. We consider the motion of a rarefied gas between two parallel plates separated by a distance, $h$, with the top plate moving with a velocity $U_{wall}/2$ and the bottom plate moving with a velocity $-U_{wall}/2$, creating a velocity difference across the channel of $U_{wall}$. In the free molecular limit (i.e. $\delta = h/\lambda \rightarrow 0$), the linearized Boltzmann equation can be solved analytically for small Mach number (Kogan Reference Kogan1969). The velocity of the rarefied gas in this limit is constant over the gap and equal to the arithmetic mean of the wall velocities such that the dimensional free molecular shear stress $\tau _{{xz}}^{{fm}}$ in the ${xz}$ plane is given by

(3.5)\begin{equation} \tau_{{{xz}}}^{{fm}} ={-}\frac{\mu U_{wall}}{\sqrt{\rm \pi} \lambda }, \end{equation}

where $\mu = (1/2) \rho c \lambda$ is the molecular viscosity, $c = \sqrt {2p/\rho }$ is the mean molecular speed, $\rho$ is the density and $p$ is the thermodynamic pressure. Sharipov (Reference Sharipov2016) used the method of successive approximations based on the linearized Boltzmann equation to obtain the first order correction to the stress for small values of $\delta$. If we non-dimensionalize the stress with the free molecular stress in (3.5), the resulting expression is

(3.6)\begin{equation} \frac{\tau_{{{xz}}}}{\tau_{{{xz}}}^{{fm}}} = 1 - \frac{\sqrt{\rm \pi}}{2} \delta. \end{equation}

In the continuum limit (i.e. $\delta \gg 1$), the solution approaches the Couette solution with slip boundary conditions

(3.7)\begin{equation} \frac{\tau_{{{xz}}}}{\tau_{{{xz}}}^{{fm}}} = \frac{\sqrt{\rm \pi}}{\delta + 2 k_1} , \end{equation}

where once again $k_1$ is the slip coefficient.

In the transition region, Cercignani & Pagani (Reference Cercignani and Pagani1966) solved the linearized Boltzmann equation using a variational principle with cutoff points approaching the limits. We can define a comprehensive expression for the stress across all regimes by introducing the dimensionless function $\eta (\delta )$, defined as

(3.8)\begin{equation} \eta(\delta) = \frac{ \tau_{{{xz}}}}{\tau_{{{xz}}}^{{fm}}}, \end{equation}

which takes the form

(3.9) \begin{equation} \eta(\delta)= \left\{ \begin{array}{@{}ll} 1 - \dfrac{\sqrt{\rm \pi}}{2}\delta, & (\delta<0.1),\\ Y(\delta), & (0.1<\delta<10),\\ \dfrac{\sqrt{\rm \pi}}{\delta+2k_1}, & (\delta>10), \end{array} \right. \end{equation}

where $Y(\delta )$ is a cubic spline fit of the tabulated values obtained by Cercignani & Pagani (Reference Cercignani and Pagani1966). The plot of $\eta (\delta )$ as a function of $\delta$ is shown in figure 5. The markers ($\bullet$) are the tabulated values, the solid line is the cubic spline fit, and the limiting solutions for the free molecular and continuum regimes are the dashed and dash-dot lines, respectively.

Figure 5.

The $\eta (\delta )$ function (3.9) showing the tabulated results of Cercignani & Pagani (Reference Cercignani and Pagani1966) ($\bullet$), the cubic spline fit used to interpolate between the points (line) and the limiting solutions for the free molecular (dashed line) and slip (dash-dot line) regimes.

Lastly, we note that there are no non-continuum effects to plug flow. So the above solutions for Poiseuille and Couette flows are all that we need to define the inner solution to the lubrication analysis. The next section discusses how the pressure field in the gap between the spheres is determined.

3.2. Inner pressure field

Sundararajakumar & Koch (Reference Sundararajakumar and Koch1996) extended the results of Cercignani & Daneri (Reference Cercignani and Daneri1963) to include the effects of curvature so that the analysis could be applied to the inner region of the lubrication flow between the spheres. The approach taken here is identical except now applied to the translation of sphere $S$ in the ${x}$-direction and rotation along the ${y}$-axis, i.e. perpendicular to the line of centres. The flow in the gap is represented locally as a superposition of Couette, Poiseuille and plug flows. However, because of curvature, the resulting velocity is two-dimensional. Integrating the continuity equation (2.15) across the gap yields

(3.10)\begin{equation} \int_{Z_P}^{Z_S} \frac{\partial U_0}{\partial R} \,\textrm{d}Z + \frac{1}{R}\int_{Z_P}^{Z_S}(U_0 + V_0) \,\textrm{d}Z + \int_{Z_P}^{Z_S} \frac{\partial W_0}{\partial Z} \,\textrm{d}Z = 0. \end{equation}

The final term can be evaluated in terms of the kinematic boundary condition at the two sphere surfaces, which can be written as

(3.11a,b)\begin{equation} W_0(R,Z_P) = 0, \quad W_0(R,Z_S) = [(U_0(R,Z_S)-1)]\alpha R . \end{equation}

We rewrite the first term in (3.10) using Leibniz's rule

(3.12)\begin{equation} \int_{Z_P}^{Z_S} \frac{\partial U_0}{\partial R} \,\textrm{d}Z = \frac{\textrm{d}}{\textrm{d}R}\left(\int_{Z_P}^{Z_S} U_0 \,\textrm{d}Z\right) - U_0(R,Z_S)\frac{\textrm{d} Z_S}{\textrm{d}R} + U_0(R,Z_P)\frac{\textrm{d} Z_P}{\textrm{d}R}. \end{equation}

Introducing the variable change

(3.13)\begin{equation} \tilde{x} = \tfrac{1}{2} R^{2} \quad\rightarrow\quad {\textrm{d}\tilde{x}} = R\,\textrm{d}R, \end{equation}

(3.10) is then written as

(3.14)\begin{equation} \frac{\textrm{d}}{\textrm{d}\tilde{x}}\left(\int_{Z_P}^{Z_S} U_0\, \textrm{d}Z\right) + \frac{1}{2\tilde{x}} \int_{Z_P}^{Z_S} (U_0 + V_0)\,\textrm{d}Z - \alpha = 0 . \end{equation}

Note that the integrated continuity result shown in (3.14), expressed in terms of the mass flux through the lubrication gap, is solely due to the kinematic boundary condition at the solid surfaces and is not altered by the slip of the gas at the surfaces.

Recall that we decompose the velocities $U_0$ and $V_0$ into a sum of Poiseuille, Couette and plug flows, denoted by superscripts $p, c, g$ respectively, yielding

(3.15)\begin{equation} \left. \begin{gathered} U_0 = U^{p}_0 + U^{c}_0 + U^{g}_0, \\ V_0 = V^{p}_0 + V^{c}_0 + V^{g}_0. \end{gathered} \right\} \end{equation}

The pressure field is related to the flux through the gap. To leading order, the pressure is uniform across the gap (i.e. $\partial P_0/\partial Z = 0$). Generalizing (3.3) to the two-dimensional flow within the gap, the relations for the integral of velocities $U^{p}_0$ and $V^{p}_0$ are (Cercignani & Daneri Reference Cercignani and Daneri1963)

(3.16)\begin{gather} \int_{Z_P}^{Z_S} U^{p}_0\, \textrm{d}Z ={-}\frac{1}{2} \frac{\delta^{2}}{\delta_0^{3}} g(\delta) \frac{\partial P_0}{\partial R}, \end{gather}

(3.17)\begin{gather} \int_{Z_P}^{Z_S} V^{p}_0\, \textrm{d}Z = \frac{1}{2} \frac{\delta^{2}}{\delta_0^{3}} g(\delta) \frac{P_0}{R}. \end{gather}

Couette flow does not contribute to the flux

(3.18)\begin{equation} \int_{Z_P}^{Z_S} U^{c}_0\, \textrm{d}Z = \int_{Z_P}^{Z_S} V^{c}_0\, \textrm{d}Z = 0, \end{equation}

and for plug flow the integrals of the velocities are

(3.19)\begin{equation} \int_{Z_P}^{Z_S} U^{g}_0\, \textrm{d}Z ={-} \int_{Z_P}^{Z_S} V^{g}_0\, \textrm{d}Z = \tfrac{1}{2}(1-\omega). \end{equation}

Substituting these expressions into (3.14) yields a differential equation for the non-continuum pressure. To simplify the result, we introduce the notation

(3.20)\begin{equation} \tilde q= \frac{P_0}{R} . \end{equation}

Notice that only the integrals of the velocities due to Poiseuille flow contain the pressure, while (3.19) can be combined with the ‘$\alpha$’ in (3.14) to give $(2\alpha - 1 + \omega )$. Factoring that quantity out of (3.14) yields the following

(3.21)\begin{equation} 2\tilde{x}\delta^{2} g q'' + \left[2\tilde{x} \frac{\textrm{d}(\delta^{2} g)}{\textrm{d}\tilde{x}} + 4\delta^{2} g \right] q'+ \frac{\textrm{d}(\delta^{2} g)}{\textrm{d}\tilde{x}} q ={-}\delta_0^{3} , \end{equation}

where

(3.22)\begin{equation} \left. \begin{gathered} q' \equiv \dfrac{\textrm{d}q}{\textrm{d}\tilde{x}} = \frac{1}{R^{2}}\frac{\partial P_0}{\partial R} -\frac{P_0}{R^{3}} ,\\ q'' \equiv \dfrac{\textrm{d}^{2} q}{{\textrm{d}\tilde{x}}^{2}} =\frac{1}{R^{3}}\frac{\partial^{2} P_0}{\partial R^{2}} - \frac{3}{R^{4}}\frac{\partial P_0}{\partial R} + \frac{3P_0}{R^{5}} , \end{gathered} \right\} \end{equation}

and $q$ is $\tilde {q}$ rescaled as follows to separate translational and rotational motions

(3.23)\begin{equation} \tilde{q} = q (2 \alpha - 1 + \omega). \end{equation}

In the limit of pure translation, $\omega \rightarrow 0$ and therefore $\tilde {q} = q(2 \alpha -1)$. In the limit of pure rotation, $\omega \rightarrow \infty$ , $\tilde {q} = q \omega$ and from (2.10) the pressure is rescaled as $p = \mu \varOmega P \cos \theta$. Note that the definition of $q$ precisely eliminates the singularity (at $r = 0$) in the differential equation. The numerical solution for $q$ of interest is finite at $\tilde{x} = 0$ and asymptotes to the continuum pressure $q_\infty$ as $\tilde{x} \rightarrow \infty$ (O'Neill & Majumdar Reference O'Neill and Majumdar1970)

(3.24)\begin{equation} q_\infty = \tfrac{1}{5 H^{2}} . \end{equation}

The ratio of the non-continuum pressure to the continuum pressure as a function of $\delta _0$ is shown in figure 6.

Figure 6.

Graph of non-continuum pressure normalized by the continuum pressure on the axis between the centres plotted against rarefied parameter $\delta _0$.

3.3. Force and torque from inner region

The force exerted by the fluid in the inner region of the surface of the particles is comprised of the force due to pressure (Poiseuille flow) and the force arising from the shear stresses (Poiseuille and Couette flows).

The force due to pressure is calculated by integrating minus the pressure times the unit outward normal on the surface, $\boldsymbol {n}_P$ and $\boldsymbol {n}_S$ for spheres $P$ and $S$, respectively, over the sphere surface. The unit normals can be approximated as

(3.25)\begin{equation} \left. \begin{gathered} \boldsymbol{n}_P = \boldsymbol{e}_{{z}} + \frac{r}{a_P}\sin\theta\; \boldsymbol{e}_{{y}} + \frac{r}{a_P}\cos\theta\ \boldsymbol{e}_{{x}} ,\\ \boldsymbol{n}_S ={-} \boldsymbol{e}_{{z}} + \frac{r}{a_S}\sin\theta\; \boldsymbol{e}_{{y}} + \frac{r}{a_S}\cos\theta\ \boldsymbol{e}_{{x}}, \end{gathered} \right\} \end{equation}

where $\boldsymbol {e}_{{x}}$, $\boldsymbol {e}_{{y}}$ and $\boldsymbol {e}_{{z}}$ are unit normals pointing in the ${x}$, ${y}$ and ${z}$ directions, respectively. Recall the pressure to leading order in $\varepsilon$ can be expressed as

(3.26)\begin{equation} p(r,\theta,{{z}}) = \varepsilon^{{-}3/2} \frac{\mu U_S}{a} P_0(R) \cos\theta , \end{equation}

where $P_0(R)$ is determined by integrating (3.21) and using the relationships $P_0 = \tilde {q} \sqrt {2\tilde{x}}$ and $\tilde{x} = R^{2}/2$. We are interested in the contribution to the force acting in the ${x}$-direction (i.e. tangential), which is obtained by integrating minus the pressure times the ${x}$ component of the unit normal over the surface of both spheres. Because the pressure acts normal to the sphere surfaces it produces no torque on the spheres. The result is

(3.27)\begin{equation} - \int_0^{r_M} \int_0^{2{\rm \pi}} p \boldsymbol{n}_P \boldsymbol{\cdot} \boldsymbol{e}_{x}\, \mathrm{d} \theta\ r\, \mathrm{d}r ={-}6{\rm \pi}\mu U_S a (2\alpha-2) \int_0^{\tilde{x}_M} \tilde{x} \tilde{q}\, \mathrm{d}\tilde{x}. \end{equation}

A momentum balance on the inner region is used to construct the components of the shear stress due to Poiseuille flow in the $r{z}$ and $\theta {z}$ planes, as shown in figure 7. The resulting expressions are

(3.28)\begin{gather} \tau^{p}_{r{{z}}} = \frac{1}{2} h(r) \frac{\partial p}{\partial r} , \end{gather}

(3.29)\begin{gather}\tau^{p}_{\theta {{z}}} = \frac{1}{2}\frac{h(r)}{r}\frac{\partial p}{\partial\theta} , \end{gather}

where the superscript ‘$p$’ indicates the shear stress due to Poiseuille flow and the subscripts are the coordinate directions. Resolving the stress acting in the $x$ direction yields

(3.30)\begin{equation} \tau^{p}_{{x{z}}} = \tau^{p}_{r{{z}}} \cos\theta - \tau^{p}_{\theta {{z}}} \sin\theta . \end{equation}

Substituting (3.26) and integrating the result over the sphere surfaces allows us to calculate the force and torque due to the shear stress from Poiseuille flow.

Figure 7.

Momentum balance in the $r{z}$ plane (a) and in the $\theta {z}$ plane (b) due to Poiseuille flow only in the local gap in cylindrical polar coordinates.

Finally, we consider the shear stress due to Couette flow in the gap, which can be written as (cf. (3.5), (3.8) and (3.9))

(3.31)\begin{equation} \tau^{c}_{{x{z}}} ={-}\frac{\mu U_{wall}}{\sqrt{\rm \pi} h_0}\delta_0\eta(\delta) , \end{equation}

where the superscript ‘$c$’ indicates the shear stress due to Couette flow. Integrating this quantity over the surface of the two spheres enables the calculation of the force and torque due to the shear stress from Couette flow.

We obtain the total force acting on each sphere in the tangential ($x$) direction by summing the contribution from the pressure (Poiseuille flow) and shear stresses (Poiseuille and Couette flows). Likewise, we can compute the total torque by summing the contributions due to the Poiseuille and Couette stresses. The tangential force and torque acting on a sphere due to its translation and rotation are non-dimensionalized as follows

(3.32a,b)\begin{gather} \hat{f}^{{\perp}}_{t,i} = \frac{F^{{\perp}}_{t,i}}{-6{\rm \pi}\mu a U_S},\quad \hat{f}^{{\perp}}_{r,i} = \frac{F^{{\perp}}_{r,i}}{-6{\rm \pi}\mu a (a_S \varOmega)} , \end{gather}

(3.33a,b)\begin{gather} \hat{t}^{{\perp}}_{t,i} = \frac{T^{{\perp}}_{t,i}}{-8{\rm \pi}\mu a a_i U_S},\quad \hat{t}^{{\perp}}_{r,i} = \frac{T^{{\perp}}_{r,i}}{-8{\rm \pi}\mu a a_i (a_S \varOmega)} , \end{gather}

where the subscript $i$ refers to sphere $P$ or $S$ and the hat above the variable indicates the contribution from the inner region.

Integrating the radial direction up to the non-dimensional radial distance $R_{M}$, which lies in the range $\delta _0^{-1/2}\ll R_M\ll \varepsilon ^{-1/2}$, yields the following expressions for the dimensionless non-continuum inner forces and torques

(3.34)\begin{gather} \left. \begin{gathered} \hat{f}^{{\perp} {NC}}_{t,S} = \frac{1}{3\sqrt{\rm \pi}}\int_{\delta_0}^{\delta_{M}} \eta(\delta)\,\mathrm{d} \delta + (2\alpha-1)(\tilde{x}Hq)_{M} + (2\alpha - 1)^{2}\int_0^{\tilde{x}_{M}} \tilde{x}q\, \mathrm{d}\tilde{x},\\ \hat{f}^{{\perp} {NC}}_{t,P} ={-}\frac{1}{3\sqrt{\rm \pi}}\int_{\delta_0}^{\delta_{M}} \eta(\delta)\,\mathrm{d} \delta + (2\alpha-1)(\tilde{x}Hq)_{M} - (2\alpha - 1)^{2}\int_0^{\tilde{x}_{M}} \tilde{x}q\, \mathrm{d}\tilde{x},\\ \hat{f}^{{\perp} {NC}}_{r,S} ={-}\frac{1}{3\sqrt{\rm \pi}}\int_{\delta_0}^{\delta_{M}} \eta(\delta)\,\mathrm{d} \delta + (\tilde{x}Hq)_{M} + (2\alpha - 1)\int_0^{\tilde{x}_{M}} \tilde{x}q \,\mathrm{d}\tilde{x},\\ \hat{f}^{{\perp} {NC}}_{r,P} = \frac{1}{3\sqrt{\rm \pi}}\int_{\delta_0}^{\delta_{M}} \eta(\delta)\,\mathrm{d} \delta + (\tilde{x}Hq)_{M} - (2\alpha - 1)\int_0^{\tilde{x}_{M}} \tilde{x}q\, \mathrm{d}\tilde{x}, \end{gathered} \right\} \end{gather}

(3.35)\begin{gather} \left. \begin{gathered} \hat{t}^{{\perp}{NC}}_{t,S} ={-} \frac{1}{4\sqrt{\rm \pi}}\int_{\delta_0}^{\delta_{M}} \eta(\delta)\,\mathrm{d} \delta - \frac{3}{4}(2\alpha-1)(\tilde{x}Hq)_{M} + \frac{3}{4}(2\alpha - 1)\int_0^{\tilde{x}_{M}} \tilde{x}q \,\mathrm{d}\tilde{x},\\ \hat{t}^{{\perp} {NC}}_{t,P} ={-} \frac{1}{4\sqrt{\rm \pi}}\int_{\delta_0}^{\delta_{M}} \eta(\delta)\,\mathrm{d} \delta + \frac{3}{4}(2\alpha-1) (\tilde{x}Hq)_{M} - \frac{3}{4}(2\alpha - 1)\int_0^{\tilde{x}_{M}} \tilde{x}q\, \mathrm{d}\tilde{x}, \\ \hat{t}^{{\perp} {NC}}_{r,S} = \frac{1}{4\sqrt{\rm \pi}}\int_{\delta_0}^{\delta_{M}} \eta(\delta)\,\mathrm{d} \delta - \frac{3}{4}(\tilde{x}Hq)_{M} + \frac{3}{4}\int_0^{\tilde{x}_{M}} \tilde{x}q\, \mathrm{d}\tilde{x},\\ \hat{t}^{{\perp} {NC}}_{r,P} = \frac{1}{4\sqrt{\rm \pi}}\int_{\delta_0}^{\delta_{M}} \eta(\delta)\,\mathrm{d} \delta + \frac{3}{4}(\tilde{x}Hq)_{M} - \frac{3}{4}\int_0^{\tilde{x}_{M}} \tilde{x}q\, \mathrm{d}\tilde{x} , \end{gathered} \right\} \end{gather}

where $\tilde{x}_{M} = R_{M}^{2}/2$, $\delta _{M}\equiv \delta _0(1+\tilde{x}_{M})$ and $(\tilde{x} H q)_{M}\equiv \tilde{x}_{M} (1 + \tilde{x}_{M}) q(\tilde{x}_{M})$. The superscript ‘$NC$’ is added to indicate the non-continuum inner solution. To complete the analysis, we supplement the inner contribution for the forces and torques with the outer contributions computed for continuum flow using a multipole solution to obtain an expression that is valid for all $\varepsilon$ in § 3.5. However, we first demonstrate that the non-continuum inner expressions derived herein are consistent with the classical analysis for a continuum fluid with ‘slip’ boundary conditions in § 3.4.

3.4. Relationship to classical slip analysis

The earliest work to introduce the first effects of the discrete nature of the gas was carried out by Hocking (Reference Hocking1973). For $\delta _0\gg 1$, he modelled the flow using the continuum momentum equation, but replaced the no-slip boundary condition with the Maxwell slip boundary condition, as shown below for translational and rotational motions

(3.36)\begin{equation} \left. \begin{gathered} U_0(R, Z_P) - \beta\left. \frac{\partial U_0}{\partial Z} \right|_{Z_P}= V_0(R, Z_P) - \beta \left.\frac{\partial V_0}{\partial Z} \right|_{Z_P} = 0 ,\\ U_0(R, Z_S) + \beta \left.\frac{\partial U_0}{\partial Z} \right|_{Z_S} = V_0(R, Z_S) + \beta\left. \frac{\partial V_0}{\partial Z} \right|_{Z_S} = 1 - \omega, \end{gathered} \right\} \end{equation}

where $\beta$ is the slip coefficient defined as

(3.37)\begin{equation} \beta \equiv \frac{k_1}{\delta_0} , \end{equation}

and recall $k_1 = 1.016$ as predicted by the BGK theory (Loyalka & Ferziger Reference Loyalka and Ferziger1967).

We now show that this earlier analysis is consistent with the present work. The differential equation for the non-continuum pressure is shown in (3.21). If we replace $g(\delta )$ by the form it takes in the slip regime, $g(\delta )=\delta /6 + k_1$, we obtain the following differential equation for the pressure, originally derived by Hocking (Reference Hocking1973)

(3.38)\begin{gather} 2\tilde{x} H^{2}(H + 6\beta) q_S'' + H[2\tilde{x}(3H + 12\beta) + 4H(H + 6\beta)] q_S' + H(3H + 12\beta) q_S ={-}1 , \end{gather}

where $q_S$ indicates the continuum pressure with slip boundary conditions, scaled similarly as $\tilde q$ in (3.23). Likewise, substituting $\eta _S(\delta ) = \sqrt {{\rm \pi} }/(\delta + 2k_1)$ into the expressions for the force and torque given in (3.34) and (3.35), the inner resistivities reduce to the form derived by Hocking (Reference Hocking1973) (note that Hocking presented combined effects due to translational and rotational motions and defined the resistivity for torque with the opposite sign convention to this paper)

(3.39)\begin{gather} \left. \begin{gathered} \hat{f}^{{\perp} {S}}_{t,S} = \frac{1}{3}\ln\left(\frac{1+{\tilde{x}_{M}}+2\beta}{1+2\beta} \right) + (2\alpha-1) (\tilde{x}Hq_S)_{M} + (2\alpha-1)^{2}\int_0^{\tilde{x}_{M}}\tilde{x}q_S \,\mathrm{d}\tilde{x} ,\\ \hat{f}^{{\perp} {S}}_{t,P} ={-}\frac{1}{3}\ln\left(\frac{1+{\tilde{x}_{M}}+2\beta}{1+2\beta} \right) + (2\alpha-1) (\tilde{x}Hq_S)_{M} - (2\alpha-1)^{2}\int_0^{\tilde{x}_{M}}\tilde{x}q_S \,\mathrm{d}\tilde{x} ,\\ \hat{f}^{{\perp} {S}}_{r,S} ={-}\frac{1}{3}\ln\left(\frac{1+{\tilde{x}_{M}}+2\beta}{1+2\beta} \right) + (\tilde{x}Hq_S)_{M} + (2\alpha-1)\int_0^{\tilde{x}_{M}}\tilde{x}q_S\, \mathrm{d}\tilde{x} ,\\ \hat{f}^{{\perp} {S}}_{r,P} = \frac{1}{3}\ln\left(\frac{1+{\tilde{x}_{M}}+2\beta}{1+2\beta} \right) + (\tilde{x}Hq_S)_{M} - (2\alpha-1)\int_0^{\tilde{x}_{M}}\tilde{x}q_S\, \mathrm{d}\tilde{x} , \end{gathered} \right\} \end{gather}

(3.40)\begin{gather} \left. \begin{gathered} \hat{t}^{{\perp}{S}}_{t,S} ={-} \frac{1}{4}\ln\left(\frac{1\!+\!{\tilde{x}_{M}}\!+\!2\beta}{1\!+\!2\beta} \right) \!-\! \frac{3}{4}(2\alpha\!-\!1)(\tilde{x}Hq_S)_{M} \!+\! \frac{3}{4}(2\alpha-1)\int_0^{\tilde{x}_{M}}\tilde{x}q_S \,\mathrm{d}\tilde{x} ,\\ \hat{t}^{{\perp}{S}}_{t,P} ={-} \frac{1}{4}\ln\left(\frac{1\!+\!{\tilde{x}_{M}}\!+\!2\beta}{1\!+\!2\beta} \right) \!+\! \frac{3}{4}(2\alpha-1)(\tilde{x}Hq_S)_{M} \!-\! \frac{3}{4}(2\alpha-1)\int_0^{\tilde{x}_{M}}\tilde{x}q_S\,\mathrm{d}\tilde{x},\\ \hat{t}^{{\perp}{S}}_{r,S} \!=\! \frac{1}{4}\ln\left(\frac{1\!+\!{\tilde{x}_{M}}\!+\!2\beta}{1\!+\!2\beta} \right) \!-\! \frac{3}{4}(\tilde{x}Hq_S)_{M} \!+\! \frac{3}{4}\int_0^{\tilde{x}_{M}}\tilde{x}q_S\,\mathrm{d}\tilde{x} ,\\ \hat{t}^{{\perp}{S}}_{r,P} \!=\! \frac{1}{4}\ln\left(\frac{1\!+\!{\tilde{x}_{M}}\!+\!2\beta}{1\!+\!2\beta} \right) \!+\! \frac{3}{4}(\tilde{x}Hq_S)_{M} \!-\! \frac{3}{4}\int_0^{\tilde{x}_{M}}\tilde{x}q_S\,\mathrm{d}\tilde{x} , \end{gathered} \right\} \end{gather}

where the superscript ‘$S$’ designates the continuum equations with slip boundary conditions. The above expressions differ from Hocking (Reference Hocking1973) in that the term $(\tilde{x}Hq_S)_{M}$ was approximated by $\tilde{x}^{2}_{M} q_S(\tilde{x}_{M})$; however, assuming $\tilde{x}_{M}\gg 1$, the difference is inconsequential.

The inner solution for the continuum equations with ‘no-slip’ boundary conditions are recovered by setting $k_1=0$ and therefore $\beta =0$. The equation for the pressure reduces to

(3.41)\begin{equation} 2\tilde{x}H^{3} q_C'' + 2H^{2} (2H + 3\tilde{x}) q_C' + 3H^{2} q ={-}1 , \end{equation}

where $q_C$ designates the continuum pressure with no-slip boundary conditions, scaled similarly as $\tilde q$ in (3.23). The result is consistent with the equation derived by O'Neill & Majumdar (Reference O'Neill and Majumdar1970) when written in terms of $P_0$ and $R$, and using their definition for the size ratio. A particular solution for this differential equation is

(3.42)\begin{equation} q_C = \frac{1}{5 H^{2}}. \end{equation}

Setting $k_1=\beta =0$ and introducing (3.42) into (3.39) and (3.40) yields the inner contribution to the resistivity functions for the continuum equations with no-slip boundary conditions

(3.43)\begin{gather} \left. \begin{gathered} \hat{f}^{{\perp}{C}}_{t,S} = \left[\tfrac{1}{3} + \tfrac{1}{5} (2\alpha-1)^{2}\right] \ln(1+\tilde{x}_{M}) + \tfrac{2}{5}(1-\alpha)(2\alpha-1) ,\\ \hat{f}^{{\perp}{C}}_{t,P} = \left[-\tfrac{1}{3} - \tfrac{1}{5} (2\alpha-1)^{2}\right] \ln(1+\tilde{x}_{M}) + \tfrac{2}{5}\alpha(2\alpha-1) ,\\ \hat{f}^{{\perp}{C}}_{r,S} = \left[-\tfrac{1}{3} + \tfrac{1}{5} (2\alpha-1)\right] \ln(1+\tilde{x}_{M}) + \tfrac{2}{5}(1-\alpha) ,\\ \hat{f}^{{\perp}{C}}_{r,P} = \left[\tfrac{1}{3} - \tfrac{1}{5} (2\alpha-1)\right] \ln(1+\tilde{x}_{M}) + \tfrac{2}{5}\alpha , \end{gathered} \right\} \end{gather}

(3.44)\begin{gather} \left. \begin{gathered} \hat{t}^{{\perp}{C}}_{t,S} = \left[-\tfrac{1}{4} + \tfrac{3}{20}(2\alpha-1)\right] \ln(1+\tilde{x}_{M}) - \tfrac{3}{10}(2\alpha-1) , \\ \hat{t}^{{\perp}{C}}_{t,P} = \left[-\tfrac{1}{4} - \tfrac{3}{20}(2\alpha-1)\right] \ln(1+\tilde{x}_{M}) - \tfrac{3}{10}(2\alpha-1) ,\\ \hat{t}^{{\perp}{C}}_{r,S} = \tfrac{2}{5} \ln(1+\tilde{x}_{M}) - \tfrac{3}{10} , \\ \hat{t}^{{\perp}{C}}_{r,P} = \tfrac{1}{10} \ln(1+\tilde{x}_{M}) + \tfrac{3}{10} . \end{gathered} \right\} \end{gather}

The superscript ‘$C$’ indicates these are solutions to the continuum equations with no-slip boundary conditions. O'Neill & Stewartson (Reference O'Neill and Stewartson1967) derived results for the pure translation of a sphere parallel to a plane, which corresponds to setting $\alpha =1$ yielding $\hat {f}^{\perp {C}}_{t,S} = 8(\ln (1 + \tilde{x}_{M}))/15$. Consistent with Jeffrey & Onishi (Reference Jeffrey and Onishi1984) and O'Neill & Stewartson (Reference O'Neill and Stewartson1967), the $O(1)$ term vanishes for $\alpha = 1$.

3.5. Uniformly valid solution

As discussed by O'Neill & Stewartson (Reference O'Neill and Stewartson1967), for a lubrication flow, the inner region corresponds to the domain where the pressure and velocity gradients are large and the lubrication approximation is valid. The inner solution will transition from the free molecular regime ($\delta _0\ll 1$) to the continuum regime ($\delta _0\gg 1$) depending upon the value of $\delta _0$. As we are restricting this investigation to the limit $Kn\ll 1$ (recall the particle Knudsen number is defined as $Kn\equiv \lambda /a$), the outer flow, which scales with the sphere radius, will always be in the continuum regime, regardless of the regime of the inner flow. Hence the outer region is governed by the Stokes equation (2.2a,b). Complete solutions for the forces and torques are obtained by summing the contributions from the inner and outer regions.

To simplify the nomenclature, we will temporarily drop references to a particular resistivity, and define $\hat {f}$ as the ‘inner’ contribution to a given resistivity and $\check {f}$ as the continuum ‘outer’ contribution to the same resistivity. The inner and outer solutions for a continuum fluid take the form

(3.45)\begin{equation} \left. \begin{gathered} \hat{f}^{C} = \hat{A}^{C}(\alpha) + B(\alpha)\ln\left(\frac{r_{M}}{h_0}\right) ,\\ \check{f} = \check{A}(\alpha) + B(\alpha)\ln\left(\frac{a}{r_{M}}\right), \end{gathered} \right\} \end{equation}

where the superscript ‘$C$’ denotes the continuum inner solution, $r_{M}$ is the boundary between the inner and outer regions and $\hat {A}^{C}(\alpha )$, $\check {A}(\alpha )$ and $B(\alpha )$ are coefficients that depend on the size ratio, $\alpha$. Note that the functional form for the inner solution in (3.45) is taken from Jeffrey & Onishi (Reference Jeffrey and Onishi1984) for no-slip boundary conditions in the limit $\tilde{x}_{M}\gg 1$. The matched asymptotic solution for the continuum model of lubrication flow then takes the form (Jeffrey & Onishi Reference Jeffrey and Onishi1984; O'Neill & Majumdar Reference O'Neill and Majumdar1970)

(3.46)\begin{equation} f^{C} = \hat{f}^{C} + \check{f} = \hat{A}^{C}(\alpha) + \check{A}(\alpha) + B(\alpha)\ln\left(\frac{1}{\varepsilon}\right). \end{equation}

Notice that the dependence on $r_{M}$ cancels out of the expression, so long as it lies within $(a h_0)^{1/2} \ll r_{M} \ll a$, where the approximations in (3.45) are valid.

Equations (3.34) and (3.35) have the following asymptotic form

(3.47)\begin{equation} \hat{f}^{NC} = \hat{A}^{NC}(\alpha, \delta_0) + B(\alpha)\ln\left(\frac{r_{M}}{h_0}\right). \end{equation}

The matched asymptotic expression for the uniformly valid resistivity for lubrication flow is then given by

(3.48)\begin{equation} f^{UV} = \hat{f}^{NC} + \check{f} = \hat{A}^{NC}(\alpha, \delta_0) + \check{A}(\alpha) + B(\alpha)\ln\left(\frac{1}{\varepsilon}\right) , \end{equation}

where the superscript ‘$UV$’ refers to the uniformly valid lubrication flow solution corresponding to $\varepsilon \ll 1$ and $\delta _0 = O(1)$. Alternatively, the uniformly valid resistivity could be constructed as follows

(3.49)\begin{equation} f^{UV} = f^{C} + \Delta f , \end{equation}

where

(3.50)\begin{equation} \Delta f \equiv \hat{f}^{NC} - \hat{f}^{C} . \end{equation}

Then, $\Delta f$ is found by subtracting the inner lubrication solutions for a continuum, (3.43) and (3.44), from (3.34) and (3.35). For example, evaluating $\Delta f^{\perp }_{t,S}\equiv \hat {f}^{\perp \mathrm {NC}}_{t,S} - \hat {f}^{\perp {C}}_{t,S}$ yields

(3.51) \begin{align} \Delta f^{{\perp}}_{t,S} &= \frac{1}{3\sqrt{\rm \pi}}\int_{\delta_0}^{\delta_{M}} \Delta\eta(\delta)\,\mathrm{d} \delta + (2\alpha - 1)^{2}\int_0^{\tilde{x}_{M}} \tilde{x}\Delta q \,\mathrm{d}\tilde{x} \nonumber\\ &\quad + \frac{1}{3\sqrt{\rm \pi}}\int_{\delta_0}^{\delta_{M}} \eta_C(\delta)\,\mathrm{d} \delta + (2\alpha-1)(\tilde{x}Hq)_{M} + (2\alpha - 1)^{2}\int_0^{\tilde{x}_{M}} \tilde{x}q_C\, \mathrm{d}\tilde{x} - \hat{f}^{{\perp}{C}}_{t,S} , \end{align}

where $\Delta \eta (\delta )\equiv \eta (\delta ) - \eta _C(\delta )$ and $\Delta q\equiv q - q_C$. In light of the analysis presented in § 3.4, the expression on the second line of (3.51) can be shown to vanish, yielding

(3.52)\begin{equation} \Delta f^{{\perp}}_{t,S} = \frac{1}{3\sqrt{\rm \pi}}\int_{\delta_0}^{\delta_{M}} \Delta\eta(\delta)\,\mathrm{d} \delta + (2\alpha - 1)^{2}\int_0^{\tilde{x}_{M}} \tilde{x}\Delta q \,\mathrm{d}\tilde{x} . \end{equation}

It is convenient to define the integrals as

(3.53)\begin{gather} \tilde{W}(\delta_0) \equiv \int_{\delta_0}^{\delta_{M}} \Delta\eta(\delta)\,\mathrm{d} \delta, \end{gather}

(3.54)\begin{gather}\tilde{Q}(\delta_0) \equiv \int_0^{\tilde{x}_{M}} \tilde{x} \Delta q\, \mathrm{d}\tilde{x} . \end{gather}

Figure 8 shows the behaviour of $\tilde {W}(\delta _0)$ and $\tilde {Q}(\delta _0)$ as a function of $\delta _0$. This behaviour can be well approximated by the following fitted functions

(3.55)\begin{equation} \tilde Q = \frac{b_1+b_2 \ln\left(\dfrac{1}{\delta_0}+b_3\right)}{1+b_4\delta_0} , \end{equation}

where $b_1 = -0.1580, b_2 = -1/5, b_3=0.65632, b_4 = 0.16330$ and

(3.56)\begin{equation} \tilde{W} = \frac{b_5+b_6 \ln\left(\dfrac{1}{\delta_0}+b_7\right)}{1+b_8\delta_0} , \end{equation}

where $b_5 = -1.6448, b_6 = -\sqrt {\rm \pi}, b_7=0.59098, b_8 = 0.19167$. The relative error is less than $2\,\%$ over the entire range.

Figure 8.

Graph of generalized pressure force function $\tilde Q$ (a) and stress function $\tilde W$ (b) as functions of $\delta _0$ to describe the tangential resistivities for any size ratio and at any gap separation.

Applying this procedure to all the expressions in (3.34) and (3.35) yields

(3.57)\begin{gather} \left. \begin{gathered} \Delta f^{{\perp}}_{t,S} ={-} \Delta f^{{\perp}}_{t,P} = \tfrac{1}{3\sqrt{\rm \pi}}\tilde W + (2\alpha - 1)^{2} \tilde Q ,\\ \Delta f^{{\perp}}_{r,S} ={-}\Delta f^{{\perp}}_{r,P} ={-}\tfrac{1}{3\sqrt{\rm \pi}}\tilde W + (2\alpha - 1) \tilde Q , \end{gathered} \right\} \end{gather}

(3.58)\begin{gather} \left. \begin{gathered} \Delta t^{{\perp}}_{t,S} ={-} \tfrac{1}{4\sqrt{\rm \pi}}\tilde W + \tfrac{3}{4} (2\alpha - 1) \tilde Q ,\\ \Delta t^{{\perp}}_{t,P} ={-} \tfrac{1}{4\sqrt{\rm \pi}}\tilde W - \tfrac{3}{4} (2\alpha - 1) \tilde Q ,\\ \Delta t^{{\perp}}_{r,S} = \tfrac{1}{4\sqrt{\rm \pi}}\tilde W + \tfrac{3}{4} \tilde Q ,\\ \Delta t^{{\perp}}_{r,P} = \tfrac{1}{4\sqrt{\rm \pi}}\tilde W - \tfrac{3}{4} \tilde Q . \end{gathered} \right\} \end{gather}

Figure 9 shows $\Delta f_{t,S}^{\perp }$ as a function of $\delta _0$. Notice that the function approaches zero for large $\delta _0$ since the non-continuum inner solution approaches the continuum inner solution in that limit. Jeffrey & Onishi (Reference Jeffrey and Onishi1984) derived generalized expressions, and developed FORTRAN codes, for the continuum resistivities that smoothly transition from twin multipole expansions at large separations (O'Neill & Majumdar Reference O'Neill and Majumdar1970) to the matched asymptotic solutions in the lubrication regime shown in (3.46). These expressions can be extended to incorporate non-continuum lubrication flows, for $Kn\ll 1$, yielding uniformly valid resistivity functions, valid for arbitrary $\varepsilon$, as follows

(3.59)\begin{equation} \left. \begin{gathered} f^{{\perp}{UV}}_{t,S} = f^{{\perp}{JO}}_{t,S}(\varepsilon,\alpha) + \Delta f^{{\perp}}_{t,S}(\delta_0,\alpha) , \\ f^{{\perp}{UV}}_{t,P} = f^{{\perp}{JO}}_{t,P}(\varepsilon,\alpha) + \Delta f^{{\perp}}_{t,P}(\delta_0,\alpha) ,\\ f^{{\perp}{UV}}_{r,S} = f^{{\perp}{JO}}_{r,S}(\varepsilon,\alpha) + \Delta f^{{\perp}}_{r,S}(\delta_0,\alpha) , \\ f^{{\perp}{UV}}_{r,P} = f^{{\perp}{JO}}_{r,P}(\varepsilon,\alpha) + \Delta f^{{\perp}}_{r,P}(\delta_0,\alpha) , \\ t^{{\perp}{UV}}_{t,S} = t^{{\perp}{JO}}_{t,S}(\varepsilon,\alpha) + \Delta t^{{\perp}}_{t,S}(\delta_0,\alpha) , \\ t^{{\perp}{UV}}_{t,P} = t^{{\perp}{JO}}_{t,P}(\varepsilon,\alpha) + \Delta t^{{\perp}}_{t,P}(\delta_0,\alpha) ,\\ t^{{\perp}{UV}}_{r,S} = t^{{\perp}{JO}}_{r,S}(\varepsilon,\alpha) + \Delta t^{{\perp}}_{r,S}(\delta_0,\alpha) ,\\ t^{{\perp}{UV}}_{r,P} = t^{{\perp}{JO}}_{r,P}(\varepsilon,\alpha) + \Delta t^{{\perp}}_{r,P}(\delta_0,\alpha), \end{gathered} \right\} \end{equation}

where the superscript ‘$JO$’ refers to the generalized resistivity functions taken from Jeffrey & Onishi (Reference Jeffrey and Onishi1984). Further details on the three-dimensional formulation of these equations are provided in Appendix A.

Figure 9.

Example of $\Delta f^{\perp }_t = \hat {f}^{NC}_t - \hat {f}^{C}_t$ obtained from (3.50) for $\alpha =1$. Non-continuum effects, corresponding to $\Delta f^{\perp }_t < 0$, are readily apparent for small $\delta _0$, whereas the function approaches zero in the continuum limit, corresponding to $\delta _0\gg 1$.

Figure 10 shows $f^{\perp }_{t,S}$, $t^{\perp }_{t,S}$ and $t^{\perp }_{r,S}$ computed in various ways vs $\varepsilon$ at the indicated Knudsen numbers. In the large $\varepsilon$ limit, all three solutions approach the continuum solution, $f^{C}_{t}$ (solid line). The open circles were computed from (3.59) using (3.57), (3.58) and the curve fits for $\tilde Q$ and $\tilde W$ given in (3.55) and (3.56). The dotted green line shows the non-continuum lubrication approximation (3.49), which is valid in the limit of small $\varepsilon$ where the overall resistivity is dominated by the non-continuum inner resistivity. In this limit, each resistivity can be shown analytically to approach a constant. In the limit as $\delta _0 \rightarrow 0$, (3.55) and (3.56) reduce to

(3.60)\begin{equation} \left. \begin{gathered} \tilde Q \rightarrow b_1+b_2 \ln\left(\frac{1}{\delta_0}\right),\\ \tilde W \rightarrow b_5+b_6 \ln\left(\frac{1}{\delta_0}\right), \end{gathered} \right\} \end{equation}

which when substituted into (3.59) yields the asymptote as a function of $Kn$

(3.61)\begin{align} \left. \begin{gathered} f^{{\perp}{UV}}_{t,S} \rightarrow \left[ -\frac{4}{15} \alpha(3 \alpha^{2}-3\alpha +2) \right] \ln\left(Kn\right) + \alpha \left[ b_1(2\alpha-1)^{2} +\frac{b_5}{3\sqrt{\rm \pi}} \right]+ A^{Y}_{22}(\alpha),\\ t^{{\perp}{UV}}_{t,S} \rightarrow \left[ \frac{1}{5} \alpha(3 \alpha -4) \right] \ln\left(Kn\right) - 2\alpha \left[ \frac{3}{4}b_1(2\alpha-1) - \frac{b_5}{4\sqrt{\rm \pi}} \right]+ B^{Y}_{22}(\alpha),\\ t^{{\perp}{UV}}_{r,S} \rightarrow \left[ -\frac{2}{5} \alpha \right] \ln\left(Kn\right) + \alpha \left[ \frac{3}{4}b_1 +\frac{b_5}{4\sqrt{\rm \pi}} \right]+ C^{Y}_{22}(\alpha), \end{gathered} \right\} \end{align}

where $A^{Y}_{22}(0.5) = 0.9983$, $B^{Y}_{22}(0.5) = -0.2390$ and $C^{Y}_{22}(0.5) = 0.7028$ are the values for the $O(1)$ terms, respectively, for $\alpha = 0.5$ (Jeffrey & Onishi Reference Jeffrey and Onishi1984). There are only a few tabulated values provided in Jeffrey & Onishi (Reference Jeffrey and Onishi1984) so we modified their FORTRAN code to output the $O(1)$ terms for any size ratio. Figure 10 shows the asymptotic values of the resistivities in the limit as $\varepsilon \rightarrow 0$ for the different Knudsen numbers that were calculated.

Figure 10.

The force on sphere $S$ due to translation (a), $f^{\perp }_{t,S}$, torque due to translation (b), $t^{\perp }_{t,S}$, and torque due to rotation (c), $t^{\perp }_{r,S}$, is plotted as a function of $\varepsilon$ for Knudsen numbers 0.01, 0.003 and 0.001 and $\alpha = 0.5$. The open circles use (3.59), the dotted green line uses (3.49), the solid green line shows the generalized expression from Jeffrey & Onishi (Reference Jeffrey and Onishi1984) and the magenta dashed line is using (3.46).

It is well known that the tangential force due to particle rotation and the torque due to tangential relative translation satisfy the reciprocal relation

(3.62)\begin{equation} \frac{f^{{\perp}{JO}}_{r,S}}{t^{{\perp}{JO}}_{t,S}} = \frac{4}{3} . \end{equation}

This equation was derived by Lorentz (Reference Lorentz1896) using the reciprocal theorem for the continuum Stokes flow equations and a no-slip boundary condition. However, a derivation based on Onsager's reciprocal relations (Onsager Reference Onsager1931a,Reference Onsagerb) for near equilibrium fluids can be used to show that the same relation holds for particles in any fluid exhibiting a linear response due to a small deviation from equilibrium (Landau & Lifshitz Reference Landau and Lifshitz1980; Ying & Peters Reference Ying and Peters1989; Zhang, Xu & Qian Reference Zhang, Xu and Qian2015) including a non-continuum gas governed by the linearized Boltzmann equation and a gas described using the Stokes equations with a Maxwell slip boundary condition. Consistent with this expectation, we can see from (3.57) and (3.58) that

(3.63)\begin{equation} \frac{\Delta f^{{\perp}}_{r,S}}{\Delta t^{{\perp}}_{t,S}} = \frac{4}{3} , \end{equation}

implying

(3.64)\begin{equation} \frac{f^{{\perp}{UV}}_{r,S}}{t^{{\perp}{UV}}_{t,S}} = \frac{4}{3} \end{equation}

at all separations.