2.1. Lubrication equations for transverse motions of permeable particles

A cylindrical coordinate system $(r,\theta,z)$ is used with $z$-coordinate coincident with the line of centres of the two particles, and with $z=0$ at the surface of particle 2 shown in figure 1. A Cartesian coordinate system $(x,y,z)$ is also defined with the same origin, $x$-coordinate aligned with $\theta =0$, and $y$ aligned with $\theta ={\rm \pi} /2$. The surfaces of the particles are approximately parabolic in the near-contact region, $r\ll a$, where $a=a_1 a_2/(a_1+a_2)$ is the reduced radius. The surface of particle 1 corresponds to $z=h_0+r^2/(2 a_1)$, and the surface of particle 2 corresponds to $z=-r^2/(2 a_2)$.

Figure 1.

Schematic showing two particles with radii and permeabilities $a_i$ and $k_i$ $(i=1,2)$, respectively, separated by a gap $h_0$; particle 1 has translational and angular velocities $U_1$ and $\omega _1$; particle 2 is stationary. The Cartesian and cylindrical coordinate systems are shown.

The leading-order lubrication equations for transverse motion of the particles are

(2.1a–d)\begin{equation} \frac{\partial^2 \bar v_r}{\partial \bar z^2}=\frac{\partial \bar p}{\partial \bar r},\quad \frac{\partial^2 \bar v_\theta}{\partial \bar z^2}=\frac{1}{\bar r}\frac{\partial \bar p}{\partial\theta},\quad \frac{\partial \bar p}{\partial \bar z}=0,\quad \frac{1}{\bar r}\frac{\partial }{\partial \bar r}\left(\bar r\bar v_r\right)+ \frac{1}{\bar r}\frac{\partial \bar v_\theta}{\partial\theta}+\frac{\partial \bar v_z}{\partial \bar z}=0. \end{equation}

The overbars denote dimensionless variables defined by

(2.2a–f)\begin{equation} \bar r=\frac{r}{L_1},\quad \bar z=\frac{z a_1}{L_1^2},\quad \bar v_r=\frac{v_r}{U_0},\quad \bar v_\theta=\frac{v_\theta}{U_0},\quad \bar w=\frac{w a_1}{U_0 L_1},\quad \bar p=\frac{p L_1^3}{\mu U_0 a_1^2}, \end{equation}

where $L_1=\sqrt {a_1 h_0}$ is the characteristic lateral length scale set by the geometry of the near-contact region, and $U_0$ is a characteristic velocity magnitude that depends on the boundary conditions of the problem. The lubrication approximation holds when $h_0\ll a$.

Only two transverse motions need consideration: (a) transverse motion of particle 1 with velocity in the $x$-direction (i.e. $\theta =0$) with magnitude $U_1$, and (b) rotation of particle 1 with angular velocity in the positive $y$-direction with magnitude $\omega _1$. Particle 2 is stationary in both problems. The resistances corresponding to the translation and rotation of particle 2 are obtained by relabelling and symmetry. Boundary conditions for these two problems are

(2.3)\begin{gather} \bar z = \bar z_1(\bar r): \left\{\begin{array}{@{}l} \bar v_r=(\bar U_1-\bar \omega_1)\cos\theta,\quad \bar v_\theta={-}(\bar U_1-\bar \omega_1)\sin\theta \\ \bar v_z={-}\bar\omega_1 \bar r\cos\theta +\bar j_1[\bar p](\bar r,\theta) \end{array}\right., \end{gather}

(2.4)\begin{gather}\bar z= \bar z_2(\bar r):\quad \bar v_r=\bar v_\theta=0,\quad \bar v_z={-}\bar j_2[\bar p](\bar r,\theta), \end{gather}

where $\bar z_1$ and $\bar z_2$ define the surfaces of particles 1 and 2, respectively,

(2.5a,b)\begin{equation} \bar z_1=1+\frac{1}{2}\bar r^2,\quad \bar z_2={-}\frac{1}{2\kappa}\bar r^2. \end{equation}

Here, $\bar U_1=U_1/U_0$ and $\bar \omega _1=\omega _1 a_1/U_0$ are the dimensionless translational and rotational velocities of particle 1. Recall that $\kappa =a_2/a_1$.

The quantities $j_1, j_2$ are the fluxes of fluid into the surfaces of the permeable particles. Given that the intraparticle pressure fields satisfy Laplace's equation, are equal to the lubrication pressure at the particle surfaces, and decay to zero inside the particles, as shown in Appendix A, it follows that the intraparticle fluxes are linear functionals of the lubrication pressure. The intraparticle fluxes are normal to the particle surfaces but, as discussed in Appendix A, act in the $z$-direction to leading order, and thus only enter the boundary condition for the velocity in the $z$-direction, as indicated in (2.3)–(2.4). The fluxes are made dimensionless by the characteristic velocity in the $z$-direction, $U_0 L_1/a_1$,

(2.6)\begin{equation} \bar j_i=\frac{j_i a_1}{U_0 L_1}, \quad i=1,2. \end{equation}

Boundary conditions (2.3)–(2.4) impose the following $\theta$-dependence:

(2.7a)\begin{gather} \bar v_r=\bar U(\bar r,\bar z)\cos\theta,\quad \bar v_\theta=\bar V(\bar r,\bar z)\sin\theta,\quad \bar v_z=\bar W(\bar r,\bar z)\cos\theta, \end{gather}

(2.7b)\begin{gather}\bar p=\bar P(\bar r)\cos\theta,\quad \bar j_i=\bar J_i[\bar P](\bar r)\cos\theta,\quad i=1,2. \end{gather}

Inserting these forms into the lubrication equations (2.1a–d) and boundary conditions (2.3)–(2.4) yields,

(2.8a,b)\begin{gather} \frac{\partial^2 \bar U}{\partial \bar z^2}=\bar P',\quad \frac{\partial^2 \bar V}{\partial \bar z^2}={-}\frac{\bar P}{\bar r}, \end{gather}

(2.9)\begin{gather} \frac{\partial \bar U}{\partial \bar r}+\frac{1}{\bar r}\left(\bar U+\bar V\right)+\frac{\partial \bar W}{\partial \bar z}=0, \end{gather}

(2.10)\begin{gather}\bar z=\bar z_1(\bar r): \left\{\begin{array}{@{}l} \bar U=(\bar U_1-\bar \omega_1),\quad \bar V={-}(\bar U_1-\bar \omega_1),\\ \bar W ={-}\bar \omega_1 \bar r +\bar J_1[\bar P](\bar r). \end{array}\right. \end{gather}

(2.11)\begin{gather}\bar z= \bar z_2(\bar r):\quad \bar U=\bar V=0,\quad \bar W={-}\bar J_2[P](\bar r). \end{gather}

Note that $\bar P$ depends only on $\bar r$, and the prime in (2.8a,b) denotes a derivative.

Integrating the (2.8a,b) with boundary conditions (2.10)–(2.11) for $\bar U$ and $\bar V$ yields

(2.12a)\begin{gather} \bar U=\frac{1}{2}\bar P'(\bar r)(\bar z-\bar z_1(\bar r))(\bar z-\bar z_2(\bar r))+\frac{\bar U_1-\bar \omega_1}{\bar h(\bar r)}(\bar z-\bar z_2(\bar r)), \end{gather}

(2.12b)\begin{gather}\bar V={-}\frac{1}{2\bar r}\bar P(\bar r)(\bar z-\bar z_1(\bar r))(\bar z-\bar z_2(\bar r))-\frac{\bar U_1-\bar \omega_1}{\bar h(\bar r)}(\bar z-\bar z_2(\bar r)), \end{gather}

where

(2.13)\begin{equation} \bar h(\bar r)=\bar z_{1}-\bar z_{2}=1+\tfrac{1}{2}(1+\kappa^{{-}1})\bar r^2. \end{equation}

Inserting these results into the continuity equation (2.9) and integrating using the boundary conditions for $\bar W$, yields the Reynolds lubrication equation

(2.14a)\begin{equation} \frac{1}{12\bar r}\left(\bar r \bar P'\bar h^3\right)' -\frac{1}{12\bar r^2}\bar P \bar h^3-2\bar J\left[\bar P\right]={-}\frac{1}{2}C\left(1+\kappa^{{-}1}\right)\bar r, \end{equation}

which satisfies the homogeneous boundary conditions

(2.14b)\begin{equation} \bar P(0)=\bar P(\infty)=0. \end{equation}

Here, primes are used to denote differentiation with respect to $\bar r$ and $C$ is the constant

(2.15)\begin{equation} C=\left[\left(\bar U_1+\bar \omega_1\right)-\kappa^{{-}1}\left(\bar U_1-\bar \omega_1\right)\right] \left(1+\kappa^{{-}1}\right)^{{-}1}. \end{equation}

The quantity $2\bar J=\bar J_1+\bar J_2$ is the combined flux into both particles. For hard spheres, i.e. $\bar J=0$, the solution of (2.14) is

(2.16)\begin{equation} \bar P_{\infty}(\kappa,\bar r)=C\frac{6 }{5} \frac{\bar r}{\bar h^{2}}. \end{equation}

2.2. Forces, torques and stresslets

The forces, torques and stresslets on each of the spheres are calculated using the relations (O'Neill & Majumdar Reference O'Neill and Majumdar1970b; Corless & Jeffrey Reference Corless and Jeffrey1988b),

(2.17)\begin{gather} \frac{F_{x_1}}{{\rm \pi} \mu U_0 a_1}= \int^{\bar R_{0}}_{0} \left[-\bar P r + \left(\frac{\partial \bar V}{\partial \bar z}- \frac{\partial \bar U}{\partial \bar z}\right)_{\bar z=\bar z_{1}}\right] r \,{\rm d} r, \end{gather}

(2.18)\begin{gather}\frac{F_{x_2}}{{\rm \pi} \mu U_0 a_1}=\int^{\bar R_{0}}_{0} \left[-\bar P r \kappa^{{-}1} + \left(\frac{\partial \bar U}{\partial \bar z}-\frac{\partial \bar V}{\partial \bar z}\right)_{\bar z=\bar z_{2}}\right] r \,{\rm d} r, \end{gather}

(2.19)\begin{gather}\frac{T_{y_1}}{{\rm \pi} \mu U_0 a^{2}_{1}}=\int^{\bar R_{0}}_{0} \left(\frac{\partial \bar U}{\partial \bar z}-\frac{\partial \bar V}{\partial \bar z}\right)_{\bar z=\bar z_{1}} r \,{\rm d} r, \end{gather}

(2.20)\begin{gather}\frac{T_{y_2}}{{\rm \pi} \mu U_0 a^{2}_{1}}=\kappa\int^{\bar R_{0}}_{0} \left(\frac{\partial \bar U}{\partial \bar z}-\frac{\partial \bar V}{\partial \bar z}\right)_{\bar z=\bar z_{2}} r\,{\rm d} r, \end{gather}

(2.21)\begin{gather}\frac{S_{xz,1}}{{\rm \pi} \mu U_0 a^{2}_{1}}= \int^{\bar R_{0}}_{0} \left[\bar P r + \frac{1}{2}\left(\frac{\partial \bar U}{\partial \bar z}-\frac{\partial \bar V}{\partial \bar z}\right)_{\bar z=\bar z_{1}}\right] r\,{\rm d} r, \end{gather}

(2.22)\begin{gather}\frac{S_{xz,2}}{{\rm \pi} \mu U_0 a^{2}_{1}}= \int^{\bar R_{0}}_{0} \left[-\bar P r + \frac{\kappa}{2}\left(\frac{\partial \bar U}{\partial \bar z}-\frac{\partial \bar V}{\partial \bar z}\right)_{\bar z=\bar z_{2}}\right]r \,{\rm d} r, \end{gather}

where the upper limit of integration is $\bar R_{0}=O(\sqrt {a_1/h_0})$, a precise definition is not required (O'Neill & Stewartson Reference O'Neill and Stewartson1967; Kim & Karrila Reference Kim and Karrila2005).

Inserting (2.12a,b) into (2.17)–(2.22) and integrating by parts to separate the pressure contributions yields

(2.23)\begin{gather} \frac{F_{x_1}}{{\rm \pi} \mu U_0 a_1}=\frac{1}{2} C \left(\kappa^{{-}1}-1\right) (I_1-I_K )-2 B I_2, \end{gather}

(2.24)\begin{gather}\frac{F_{x_2}}{{\rm \pi} \mu U_0 a_1}={-}\frac{1}{2} C \left(\kappa^{{-}1}-1\right) (I_1-I_K) + 2 B I_2, \end{gather}

(2.25)\begin{gather}\frac{T_{y_1}}{{\rm \pi} \mu U_0 a^{2}_{1}}={-}\frac{1}{2} C\left(1+\kappa^{{-}1}\right)(I_1-I_K) + 2 B I_2, \end{gather}

(2.26)\begin{gather}\frac{T_{y_2}}{{\rm \pi} \mu U_0 a^{2}_{1}}= \frac{1}{2}C \left(1+\kappa\right)(I_1-I_K) + 2 \kappa B I_2, \end{gather}

(2.27)\begin{gather}\frac{S_{xz,1}}{{\rm \pi} \mu U_0 a^{2}_{1}}=\frac{1}{4}C\left(3-\kappa^{{-}1}\right) (I_1-I_K) + B I_2, \end{gather}

(2.28)\begin{gather}\frac{S_{xz,2}}{{\rm \pi} \mu U_0 a^{2}_{1}}={-}\frac{1}{4}C(3-\kappa) (I_1-I_K) + \kappa B I_2, \end{gather}

where $B$ is given by

(2.29)\begin{equation} B=\bar U_1-\bar \omega_1 \end{equation}

and $C$ is defined by (2.15). Here, we define the integrals

(2.30a,b)\begin{equation} I_1(\xi,\kappa,\bar R_{0})=\int^{\bar R_{0}}_{0} \bar{\mathcal{P}}_\infty r^2 \,{\rm d} r,\quad I_2(\xi,\kappa,\bar R_{0})=\int^{\bar R_{0}}_{0} \frac{r}{\bar h} \,{\rm d} r, \end{equation}

and

(2.31)\begin{equation} I_K(\xi,\kappa,K)=\int^{\infty}_{0} (\bar{\mathcal{P}}_\infty-\bar{\mathcal{P}}) r^2 \,{\rm d} r, \end{equation}

where $\xi =h_0/\bar a$ is the gap normalized by the average radius, and $\kappa =a_2/a_1$ is the size ratio. Here, $\mathcal {P}$ is the rescaled pressure

(2.32)\begin{equation} P=C\mathcal{P}, \end{equation}

which is introduced for convenience in the solution of the Reynolds lubrication equation that follows below.

This rearrangement of (2.17)–(2.22) is useful because it isolates the influence of particle permeability. The integrals $I_1$ and $I_2$ describe the dynamics of hard spheres; only the permeability integral, $I_K$, depends on the permeability. The upper limit for $I_K$ can be extended to infinity because $\bar {\mathcal {P} }_\infty -\bar {\mathcal {P}}$ decays sufficiently fast for $\bar r\to \infty$, indicating that permeability does not affect flow in the matching region.

3.1. Numerical method

Integral (3.10), was discretized on a set of $N$ points using a uniform mesh $\hat r_i$, ($i=1,\ldots, N$) on the interval $0 \leq \hat r \leq \hat r_N$. An $N\times N$ system of equations was thus generated for the values $\hat w_i$ at each of the points $\hat r_i$ using a piecewise linear representation of $\hat w$. The non-singular integrals $\hat I_A$ and $\hat I_B$, respectively, yield upper and lower triangular matrices with elements evaluated by trapezoid-rule integration on the intervals between points. The matrix obtained by discretizing the flux integral (3.8) was obtained by analytically evaluating the log-singular portion of the Green's function (A20) on all intervals and evaluating the non-analytic remainder with adaptive Gaussian quadratures. This method yields $O(1/N^2)$ convergence, and the permeability function $g(q)$ was obtained to $O(1/N^3)$ accuracy by extrapolating for $N\to \infty$ using results from calculations with different numbers of points. The contributions to integrals from the region $\hat r_N\leq \hat r\leq \infty$ were approximately incorporated using $\hat w\approx \hat w_\infty$ for $\hat r>\hat r_N$ to accelerate the convergence for $\hat r_N\to \infty$. The system of equations were iteratively solved using a Gauss–Seidel scheme.

The numerical values provided in the Supplementary Material available at https://doi.org/10.1017/jfm.2022.171 were obtained from calculations with $N\approx 400$ and $\hat r_N\approx 20$ and are accurate to 3 digits.

3.2. Numerical and asymptotic results for the permeability function $g(q)$

Results for the permeability function $g(q)$ are shown graphically in figure 2 and provided in tabular form in the Supplementary Material. The limiting behaviours for small and large $q$, derived in Appendix B, are

(3.17)\begin{equation} g(q)={-}\tfrac{12}{5}\log q+b_1+b_2 q + O(q^2),\quad q\ll 1, \end{equation}

with $b_1\doteq -0.48$ and $b_2\doteq -1.5$, and

(3.18)\begin{equation} g(q)=c_1 q^{{-}5/2}+O\left(q^{{-}5}\right),\quad q\gg 1, \end{equation}

with $c_1\doteq 2.12$. The results show that $g(q)>0$, indicating that permeability reduces the lubrication pressure between particles undergoing transverse near-contact motion, i.e. $\hat P_\infty > \hat P$ according to (3.13).

Figure 2.

Transverse and axisymmetric permeability functions $g(q)$ and $f(q)$, defined by (3.13) and (D5), respectively; numerical solutions (solid lines), asymptotic forms (3.17)–(3.18) and (D6)–(D7) (dashed lines). Tables of numerical values for $f(q)$ and $g(q)$ are provided in the Supplementary Material.

4.1. Transverse resistance functions

Rearranging (2.23)–(2.28) to isolate the type of motion yields the forces, torques and stresslets in terms of resistance functions. The notation and general definition of resistance functions used in the literature is followed (Jeffrey Reference Jeffrey1992; Kim & Karrila Reference Kim and Karrila2005). The results are

(4.1)\begin{gather} -\frac{F_{x_1}}{ {\rm \pi}\mu U_0 a_1}= 6 Y^{A}_{11} \bar U_1+4Y^{B}_{11} \bar \omega_1, \end{gather}

(4.2)\begin{gather}-\frac{F_{x_2}}{ {\rm \pi}\mu U_0 a_1}= 3(1+\kappa) Y^{A}_{21} \bar U_1+ (1+\kappa)^2 Y^{B}_{21} \bar \omega_1, \end{gather}

(4.3)\begin{gather}- \frac{T_{y_1}}{ {\rm \pi}\mu U_0 a^2_1}= 4 Y^{B}_{11} \bar U_1+8 Y^{C}_{11} \bar \omega_1, \end{gather}

(4.4)\begin{gather}-\frac{T_{y_2}}{{\rm \pi} \mu U_0 a^2_1}= (1+\kappa)^2 Y^{B}_{21} \bar U_1+ (1+\kappa)^3 Y^{C}_{21} \bar \omega_1, \end{gather}

(4.5)\begin{gather}- \frac{S_{xz,1}}{{\rm \pi} \mu U_0 a^{2}_{1}}={-}4 Y^{G}_{11} \bar U_1+ 8 Y^{H}_{11} \bar \omega_1, \end{gather}

and

(4.6)\begin{equation} - \frac{S_{xz,2}}{{\rm \pi} \mu U_0 a^{2}_{1}}={-} (1+\kappa)^2 Y^{G}_{21} \bar U_1+ (1+\kappa)^3 Y^{H}_{21}\bar \omega_1. \end{equation}

Here, $Y^{R}_{\alpha \beta }(\xi,\kappa,q)$ are the transverse resistance functions, where the superscript $R$ refers to one of the resistance tensors $A, B, C, G$ or $H$ in the resistance matrix (C1), and subscripts $\alpha$ and $\beta$ refer to the particle labels 1 or 2. Using symmetry relations and the Lorentz reciprocal theorem, the grand resistance matrix can be derived from this set of resistance functions, as shown in Appendix C.

The transverse resistance functions are

(4.7)\begin{equation} Y^{R}_{\alpha \beta}(\xi,\kappa,q)= Y^{R,0}_{\alpha \beta}(\xi,\kappa)-g(q) \varUpsilon^{R}_{\alpha \beta}(\kappa)+O\left(\xi \log \xi^{{-}1}\right), \end{equation}

where $Y^{R,0}_{\alpha \beta }$ are the hard-sphere resistance functions (C5)–(C14), $g(q)$ is the transverse permeability function (3.13) shown in figure 2 and $\varUpsilon ^{R}_{\alpha \beta }$ are the size-ratio-dependent coefficients given by

(4.8a,b)\begin{gather} \varUpsilon^{A}_{11}=\frac{\kappa}{12}\frac{(1-\kappa)^2 }{(\kappa +1)^3},\quad \varUpsilon^{A}_{21}={-}\frac{\kappa}{6}\frac{(1-\kappa)^{2}}{(\kappa +1)^{4}}, \end{gather}

(4.9a,b)\begin{gather}\varUpsilon^{B}_{11}={-}\frac{\kappa}{8}\frac{(1-\kappa )}{(\kappa +1)^2},\quad \varUpsilon^{B}_{21}=\frac{\kappa^2}{2}\frac{1-\kappa}{(\kappa +1)^4}, \end{gather}

(4.10a,b)\begin{gather}\varUpsilon^{C}_{11}=\frac{\kappa}{16 (1+\kappa)},\quad \varUpsilon^{C}_{21}={-}\frac{ \kappa^2}{2 (1+\kappa)^4}, \end{gather}

(4.11a,b)\begin{gather}\varUpsilon^{G}_{11}={-}\frac{(3\kappa-1)(1-\kappa) \kappa}{16 (1+\kappa)^3},\quad \varUpsilon^{G}_{21}=\frac{(3-\kappa)(1-\kappa) \kappa^2}{4 (1+\kappa)^5}, \end{gather}

(4.12a,b)\begin{gather}\varUpsilon^{H}_{11}={-}\frac{\kappa}{32}\frac{3\kappa-1}{(1+\kappa)^2},\quad \varUpsilon^{H}_{21}=\frac{\kappa^2}{4}\frac{3-\kappa}{(1+\kappa)^5}. \end{gather}

Recall that $\xi =h_0/\bar a$ is the gap normalized by the average radius, $\kappa =a_2/a_1$ is the size ratio and $q=K^{-2/5}h_0/a$ is the permeability parameter. The result (4.7) indicates that particle permeability additively affects the transverse resistance functions.

4.2. Axisymmetric resistance functions

The leading-order axisymmetric resistance functions are presented here using the results of our recent analysis (Reboucas & Loewenberg Reference Reboucas and Loewenberg2021a). Following the presentation in § 4.1, the forces, torques and stresslets are (Jeffrey Reference Jeffrey1992; Kim & Karrila Reference Kim and Karrila2005)

(4.13a,b)\begin{gather} -\frac{F_{z_1}}{ {\rm \pi}\mu U_0 a_1}= 6 X^{A}_{11} \bar U_1,\quad -\frac{F_{z_2}}{ {\rm \pi}\mu U_0 a_1}= 3(1+\kappa) X^{A}_{21} \bar U_1, \end{gather}

(4.14a,b)\begin{gather}- \frac{T_{z_1}}{ {\rm \pi}\mu U_0 a^2_1}= 8 \, X^{C}_{11} \bar \omega_1,\quad -\frac{T_{z_2}}{{\rm \pi} \mu U_0 a^2_1}= (1+\kappa)^3 \, X^{C}_{21} \bar \omega_1, \end{gather}

(4.15a,b)\begin{gather}- \frac{S_{zz,1}}{{\rm \pi} \mu U_0 a^{2}_{1}}={-}4X^{G}_{11} \bar U_1,\quad - \frac{S_{zz,2}}{{\rm \pi} \mu U_0 a^{2}_{1}}={-} (1+\kappa)^2X^{G}_{21} \bar U_1, \end{gather}

where $\bar U_1$ and $\bar \omega _1$ are translational and rotational velocities of particle 1 along the line of centres, in the $z$-direction. Here, $X^{R}_{\alpha \beta }(\xi,\kappa,q)$ are the axisymmetric resistance functions, the superscript $R$ refers to the resistance tensor $A, C$ or $G$ and subscripts $\alpha$ and $\beta$ are particle labels 1 or 2. The remaining resistance functions for the grand resistance matrix (C1) for axisymmetric motion are derived from this set of functions, as shown in Appendix C.

The leading-order, axisymmetric resistance functions are

(4.16a,b)\begin{equation} X^{A}_{\alpha \beta}(\xi,\kappa,q)= X^{A,0}_{\alpha \beta}(\xi,\kappa)f(q),\quad X^{G}_{\alpha \beta}(\xi,\kappa,q)= X^{G,0}_{\alpha \beta}(\xi,\kappa)f(q), \end{equation}

where $X^{A,0}_{\alpha \beta }$ and $X^{G,0}_{\alpha \beta }$ are the hard-sphere resistance functions (C15)–(C18), and $f(q)$ is the axisymmetric permeability function (D5) shown in figure 2. This function was recently analysed (Reboucas & Loewenberg Reference Reboucas and Loewenberg2021a), and its primary features are summarized in Appendix D.

Particle permeability is seen to have a multiplicative effect on the axisymmetric resistance functions, in contrast to its additive effect on the transverse functions. Axisymmetric rotation does not generate a lubrication pressure thus, $X^{C}_{\alpha \beta }(\xi,\kappa,q)=X^{C,0}_{\alpha \beta }(\xi,\kappa )$.

4.3. Matching to the outer region

Away from the near-contact region, hard-sphere resistance functions describe pairwise hydrodynamic interactions for permeable spheres to $O(K)$. According to (3.18) and (4.7) and (D7) and (4.16a,b), the resistance functions for permeable spheres presented in § 4.1–4.2 are equal to the corresponding hard-sphere functions in the overlapping region

(4.17)\begin{equation} K^{2/5} \ll \epsilon \ll 1, \end{equation}

for $K\to 0$, where $\epsilon =h_0/a$ is the gap normalized by the reduced radius $a=a_1a_2/(a_1+a_2)$. Accordingly, the resistance functions for permeable spheres match asymptotically to the hard-sphere functions. This provides a uniformly valid approximation for the resistance functions

(4.18)\begin{equation} Z^R_{\alpha\beta}(\xi,\kappa,q)= Z^{R,0}_{\alpha\beta}(\xi,\kappa)+Z^R_{L\alpha\beta} (\xi,\kappa,q)-Z^{R,0}_{L\alpha\beta}(\xi,\kappa), \end{equation}

where $Z^R_{\alpha \beta }$ is a pairwise resistance function for permeable spheres, $Z^{R,0}_{\alpha \beta }$ is the corresponding hard-sphere resistance function, $Z^{R,0}_{L\alpha \beta }$ is the lubrication approximation for the hard-sphere function and $Z^R_{L\alpha \beta }$ is the lubrication resistance function for permeable spheres, as developed in our study.

4.4. Non-singular particle motions

The qualitative effects of particle permeability on the lubrication resistances are discussed below. Permeability relieves the lubrication pressure in the near-contact region, as discussed at the end of § 3.2 and at the end of Appendix D. We show here that permeability qualitatively alters near-contact motion of particles and gives rise to additional non-singular contact motions that are inaccessible to hard spheres. For hard spheres, only rigid-body motions, i.e. pair translation and dumbbell rotation, and relative rotations about the symmetry axis are non-singular at contact. As a result of permeability, non-shearing motions of particles in contact also become non-singular. These include rolling without slip and axisymmetric approach.

4.4.1. Transverse motions

Particle permeability lessens the magnitude of the transverse lubrication function $Y^{R}_{\alpha \beta }$ if the contributions $Y^{R,0}_{\alpha \beta }$, $\varUpsilon ^{R}_{\alpha \beta }$ in (4.7) have the same algebraic sign; the latter is determined by formulas (4.8a,b)–(4.10a,b) and (C5)–(C10). It is thus seen that permeable particles of unequal size have lower translational resistances, $Y^{A}_{11}$, $Y^{A}_{21}$; translational resistances for equal-size particles are unaffected by permeability. Permeability reduces the rotational resistance $Y^C_{11}$ of a particle spinning close to a stationary particle but it enhances the rotational coupling between the particles, increasing $Y^C_{21}$.

To generally explain the role of permeability in transverse particle motions, it is convenient to examine the forces, torques and stresslets given by (2.23)–(2.28). The effect of permeability depends on whether the contribution from the lubrication pressure, $I_1-I_K$, reinforces or opposes the contribution from shearing motion of the particle surfaces, $I_2$; permeability reduces the magnitude of the former contribution, but has no effect on the latter. Accordingly, the net effect of particle permeability depends on the coefficients of $I_1-I_K$ and $I_2$ in (2.23)–(2.28).

The forces, torques and stresslets for the special case of particle motion prescribed by

(4.19a,b)\begin{equation} (\kappa -1)\bar U_1 + (\kappa +1)\bar \omega_1=0,\quad \bar U_2=\bar \omega_2=0 \end{equation}

are unaffected by particle permeability because $C=0$ in this case, according to (2.15), and thus the pressure contribution, associated with $I_1-I_K$, vanishes in (2.23)–(2.28). This result explains why the translational resistances for equal-size particles are unaffected by permeability.

Conversely, the shearing contribution vanishes in (2.23)–(2.28) for rolling motion of the particles without slip,

(4.20a,b)\begin{equation} \bar U_1=\bar \omega_1=1,\quad \bar U_2=\bar \omega_2=0, \end{equation}

because $B=0$, according to (2.29). The forces, torques, and stresslets in this case are due entirely to the pressure contribution $I_1-I_K$, and thus have reduced magnitudes for permeable particles. Moreover, rolling motion of permeable particles is non-singular at contact, as shown by combining (3.1), (3.17) and (C4a) to yield,

(4.21)\begin{equation} \lim_{\xi\to 0}\left(I_1-I_K\right)= \frac{24}{25}\left(1+\kappa^{{-}1}\right)^{{-}2}\log K^{{-}1}+C_0, \end{equation}

where $C_0$ depends only on size ratio. Other rolling motions without slip (e.g. pure rotation; $-\kappa \bar \omega _2=\bar \omega _1=1$, $\bar U_1=\bar U_2=0$) are also non-singular, and can be generated by a superposition of (4.20a,b) with non-singular rigid-body translation and rotation.

4.4.2. Axisymmetric motions

Given the form of the axisymmetric resistances (C15)–(C18), and formula (D6), it follows that axisymmetric resistances (4.16a,b) have the limiting forms,

(4.22a,b)\begin{equation} \lim_{\xi\to 0}X^{A}_{\alpha \beta}(\xi,\kappa,q)= \chi^{A}_{\alpha \beta}(\kappa)K^{{-}2/5},\quad \lim_{\xi\to 0}X^{G}_{\alpha \beta}(\xi,\kappa,q)= \chi^{G}_{\alpha \beta}(\kappa)K^{{-}2/5}, \end{equation}

where the functions $\chi ^{A}_{\alpha \beta }(\kappa )$ and $\chi ^{G}_{\alpha \beta }(\kappa )$ depend only on the size ratio. The result demonstrates that the axisymmetric resistances for permeable particles are non-singular at contact in contrast to the $\xi ^{-1}$ singular resistances for hard spheres.

5.1. Transverse mobilities

The transverse mobilities $H$, $M$ and $B$ have the near-contact form

(5.4)\begin{equation} \varLambda(\xi,\kappa, K) =\frac{\lambda_1+\lambda_2\log\xi^{{-}1} + \lambda_3\left[\log\xi^{{-}1}-\frac{5}{12}g(q)\right]\log\xi^{{-}1} + \lambda_6 g(q)}{\lambda_4+\lambda_5\log\xi^{{-}1} + \left[\log\xi^{{-}1}-\frac{5}{12}g(q)\right]\log\xi^{{-}1} +\lambda_7 g(q)}, \end{equation}

where $\xi =h_0/\bar a$ is the gap normalized by the average radius, and $q=K^{-2/5}h_0/a$ is the permeability parameter. The coefficients $\lambda _i$ ($i=1\text {--}7$) depend only on the size ratio; numerical values for several size ratios are listed in tables 1–3.

Table 1.

Coefficients of transverse mobility $H$ given by (5.4).

Table 2.

Coefficients of transverse mobility $M$ ($\gamma =1$) given by (5.4).

Table 3.

Coefficients of transverse mobility $B$ given by (5.4).

The results shown in figures 3–5 demonstrate that permeability quantitatively affects the transverse mobilities for $\xi \lesssim O(K^{2/5})$, corresponding to $q=O(1)$, and has the largest effect for extreme size ratios. Particle permeability has no effect for equal-size particles, as seen in figures 3 and 5, because no lubrication pressure is generated by the motion, as discussed in § 4.4.1.

Figure 3.

Transverse mobility functions $H$ from exact solution for hard spheres using code from A. Z. Zinchenko (solid lines) and lubrication solutions for hard spheres, $K=0$, and permeable spheres, $K=10^{-7}$ (dashed lines) vs gap, $\xi$, for size ratios indicated. Lubrication solution for hard spheres appears as a continuation of the exact solution; the lubrication solution for permeable spheres deviates for $\xi \lesssim K^{2/5}$.

Figure 4.

Same as figure 3, except for the transverse mobility function $M$.

Figure 5.

Same as figure 3, except for the transverse mobility function $B$.

Figures 3–4 show that mobility functions $H$ and $M$ for permeable particles are larger than for hard spheres, whereas the opposite is true for mobility function $B$, as shown in figure 5. These observations indicate that particle permeability diminishes the strength of hydrodynamic pair interactions in all cases, given that $H=M=1$ and $B=0$ in the absence of hydrodynamic interactions.

Inserting (3.17) into (5.4) and taking the limit as $\xi \to 0$ yields the contact values of the transverse mobilities

(5.5)\begin{equation} \varLambda_c(\kappa, K) =\lambda_3\frac{\log K^{{-}1}+\lambda'_1}{\log K^{{-}1}+\lambda'_2},\quad \xi=0, \end{equation}

where $\varLambda _c$ is the contact value of the transverse mobility function $H$, $M$ or $B$, and the coefficients $\lambda '_1$ and $\lambda '_2$ are given by

(5.6a)\begin{gather} \lambda'_1=6\lambda_6\lambda_3^{{-}1}+\tfrac{5}{2}\lambda_2\lambda_3^{{-}1}-\tfrac{5}{2} \log\nu-\tfrac{25}{24}b_1, \end{gather}

(5.6b)\begin{gather}\lambda_2'=6\lambda_7+\tfrac{5}{2}\lambda_5-\tfrac{5}{2}\log\nu -\tfrac{25}{24}b_1, \end{gather}

where $\nu =2\kappa (1+\kappa )^{-2}$ and $b_1\doteq -0.48$. For convenience, values for $\lambda _1'$ and $\lambda _2'$ are listed in table 4. Contact values for the mobilities of hard spheres are given by $\lambda _3$; even small permeabilities can significantly alter contact mobilities because its effect decays only logarithmically.

Table 4.

Coefficients for the contact values of transverse mobilities given by (5.5).

According to formula (5.5), $\varLambda _c$ increases with $K$ for $\lambda '_1>\lambda '_2$, decreases for $\lambda '_1<\lambda '_2$, the effect is largest for disparate values of these parameters, and $\varLambda _c$ is independent of $K$ for $\lambda '_1=\lambda '_2$. Consistent with the discussion above and the results shown in figures 3–5, the values in table 4 thus indicate that contact values of mobilities $H$ and $M$ increase with $K$, $B$ decreases with $K$, the effect is largest for extreme size ratios and $\varLambda _c$ is independent of $K$ for equal-size particles.

5.2. Axisymmetric mobilities

The axisymmetric mobility $G$ has the near-contact form

(5.7)\begin{equation} G=\frac{(1+\kappa)^2}{2 \kappa} \frac{\xi}{f(q)}=\frac{\epsilon}{f(q)},\quad \xi\ll 1, \end{equation}

where $\xi =h_0/\bar a$ and $\epsilon =h_0/a$ are the gaps normalized by the average and reduced radius, respectively, $\kappa$ is the size ratio and $q=K^{-2/5}h_0/a$ is the permeability parameter. Inserting (D6)–(D7) into (5.7) yields

(5.8)\begin{equation} G =\left(b_0+b_1 q+O(q^2)\right)K^{2/5},\quad q\ll 1, \end{equation}

and

(5.9)\begin{equation} G=\left(1+c_0 q^{{-}5/2}+O(q^{{-}5})\right)\epsilon,\quad q\gg 1, \end{equation}

where $b_0\doteq 1.332$, $b_1\doteq 0.397$ and $c_0 \doteq 1.8402$. The result for hard spheres is recovered from the latter for $q\to \infty$. Equation (5.8) indicates that the particles undergo a constant approach velocity at contact under the action of a constant force, consistent with the non-singular axisymmetric resistances for permeable particles given by (4.22a,b).

In the near-contact region, $L$ and $1-A$ are proportional to $G$, motivating the definition of modified mobility functions $\hat L$ and $\hat A_1$

(5.10a,b)\begin{equation} \hat L(\epsilon,\kappa,q)=\frac{L(\epsilon,\kappa,q)}{\hat F_c^{(g)}(\kappa)},\quad \hat A_1(\epsilon,\kappa,q)=\frac{1-A(\epsilon,\kappa,q)}{\hat F_c^{(E)}(\kappa)}, \end{equation}

with the property

(5.11)\begin{equation} \lim_{\epsilon\to 0}\hat L(\epsilon,\kappa,q)= \lim_{\epsilon\to 0}\hat A_1(\epsilon,\kappa,q)= \lim_{\epsilon\to 0}G(\epsilon,q), \end{equation}

where $\hat F_c^{(g)}$ and $\hat F_c^{(E)}$ are dimensionless contact forces, obtained by the procedure described in Appendix E and given by (E3)–(E4). These forces arise between two spheres in point contact, migrating along their line of centres in gravity-driven motion, and in axisymmetric straining flow, respectively. The mobility function $L$ also depends on the ratio of the particle density differences, $\gamma$, defined by (5.3). Numerical values of the contact forces are provided in table 6 for several size ratios.

Axisymmetric mobilities $G, \hat L$ and $\hat A_1$ are shown in figure 6 as a function of gap where the common behaviour of the modified axisymmetric mobilities (5.11) is demonstrated; permeability qualitatively changes the near-contact motion for $h_0/a \leq O(K^{2/5})$, corresponding to $q\leq O(1)$. Figure 7 shows the universal near-contact behaviour predicted by combining (5.7) and (5.11)

(5.12)\begin{equation} G K^{{-}2/5}=\hat A_1 K^{{-}2/5}=\hat L K^{{-}2/5}=\frac{q}{f(q)},\quad \epsilon\to 0. \end{equation}

Figure 6.

Axisymmetric mobility functions $\hat A_1$ and $\hat L$ ($\gamma =1$), defined by (5.10a,b), and $G$ for hard spheres vs gap, $\epsilon =h_0/a$, for size ratios $\kappa =0.5$ (solid lines), and $\kappa =0.75$ (dash-dotted lines) using code from A. Z. Zinchenko (curves for $\kappa =0.5$, and $\kappa =0.75$ are indistinguishable for G); lubrication solution for hard spheres, $K=0$, and permeable spheres $K=10^{-7}$ (dashed lines).

Figure 7.

Universal near-contact behaviour (5.12) vs permeability parameter, $q$ (solid line); small- and large-$q$ asymptotes given by (5.8)–(5.9) (dashed lines).

5.3. Mobility of a particle moving parallel to a wall

In this section, the motion of a sphere in close contact with a wall is considered, corresponding to the limit $\kappa \to \infty$. Two problems are studied: (I) motion of a particle in a quiescent fluid under the action of an imposed force $\boldsymbol {F}=F_0\boldsymbol {e}_x$ parallel to the wall and (II) motion of a force-free particle in an imposed shear flow parallel to the wall, $\boldsymbol {U}_{\infty }=E_\infty z \boldsymbol {e}_x$. The particle is torque free in both problems.

The resistance formulation for problems (I) and (II), derived from the general resistance formulation (C1) for the particle–wall configuration ($\kappa \to \infty$), is

(5.13)\begin{gather} -\frac{F_{x_1}}{ {\rm \pi}\mu U_0 a_1}= 6 Y^{A,\infty}_{11} \bar U_1+4 Y^{B,\infty}_{11} \bar \omega_1 - (6 Y^{A,\infty}_{11}+2Y^{B,\infty}_{11}-4 Y^{G,\infty}_{11})\bar E_{\infty}, \end{gather}

(5.14)\begin{gather}- \frac{T_{y_1}}{ {\rm \pi}\mu U_0 a^2_1}= 4 Y^{B,\infty}_{11} \bar U_1+8 Y^{C,\infty}_{11} \bar \omega_1 -(4 Y^{B,\infty}_{11}+4Y^{C,\infty}_{11}+8 Y^{H,\infty}_{11})\bar E_{\infty}, \end{gather}

where $a_1$ is the radius of the sphere, $\bar U_1=U_1/U_0$ and $\bar \omega _1=\omega _1 a_1/U_0$ are the dimensionless translational and rotational velocities, $\bar E_{\infty }=E_{\infty }a_1/U_0$ is the dimensionless shear rate and $Y^{R,\infty }_{11}$ are the resistance functions (4.7) corresponding to the limit $\kappa \to \infty$

(5.15)\begin{gather} Y^{A,\infty}_{11}={-}\frac{g(q)}{12}+\frac{8}{15}\log \epsilon^{{-}1} + A^{Y,\infty}_{11} + O(\epsilon \log \epsilon^{{-}1}), \end{gather}

(5.16)\begin{gather}Y^{B,\infty}_{11}={-}\frac{g(q)}{8}-\frac{1}{5}\log \epsilon^{{-}1} + B^{Y,\infty}_{11} + O(\epsilon \log \epsilon^{{-}1}), \end{gather}

(5.17)\begin{gather}Y^{C,\infty}_{11}={-}\frac{g(q)}{16}+\frac{2}{5}\log \epsilon^{{-}1} + C^{Y,\infty}_{11} + O(\epsilon \log \epsilon^{{-}1}). \end{gather}

(5.18)\begin{gather}Y^{G,\infty}_{11}={-}\frac{3}{16} g(q)+\frac{7}{10}\log \epsilon^{{-}1} + G^{Y,\infty}_{11} + O(\epsilon \log \epsilon^{{-}1}). \end{gather}

(5.19)\begin{gather}Y^{H,\infty}_{11}= \frac{3}{32} g(q)-\frac{1}{10}\log \epsilon^{{-}1} + H^{Y,\infty}_{11} + O(\epsilon \log \epsilon^{{-}1}). \end{gather}

Here, $\epsilon =h_0/a_1$ is the dimensionless gap, $q=K^{-2/5}h_0/a_1$ is the permeability parameter, $g(q)$ is the transverse mobility function and $R^{Y,\infty }_{11}$ are the matching constants to the outer solution for this geometry (Goldman et al. Reference Goldman, Cox and Brenner1967a,Reference Goldman, Cox and Brennerb; Corless & Jeffrey Reference Corless and Jeffrey1988b).

The translational and rotational mobilities of the particle are determined by solving (5.13) and (5.14), separately, for problems (I) and (II). The results are

(5.20)\begin{equation} \bar{V}_1 =\frac{\lambda_0+\lambda_1\log\epsilon^{{-}1} +\lambda_2 g(q)}{\lambda_3+\lambda_4\log\epsilon^{{-}1} + \left[\log\epsilon^{{-}1}-\frac{5}{12}g(q)\right]\log\epsilon^{{-}1} +\lambda_5 g(q)}, \end{equation}

where $\bar V_1=\bar U_1$ or $\bar V_1=\bar \omega _1$, and $U_0=(6 {\rm \pi}\mu a_1)^{-1}F_0$ for problem (I), and $U_0=E_{\infty }a_1$ for problem (II). The numerical coefficients appearing in (5.20) are listed in table 5 for both problems, and the results are plotted in figures 8 and 9. The results for hard spheres (Goldman et al. Reference Goldman, Cox and Brenner1967a,Reference Goldman, Cox and Brennerb; Corless & Jeffrey Reference Corless and Jeffrey1988b) are recovered for $g=0$, corresponding to $q\to \infty$.

Table 5.

Coefficients for particle velocity parallel to wall (5.20).

Figure 8.

Translational and rotational velocities of a particle moving parallel to a plane wall under the action of a force; mobility formula (5.20) (solid lines); simplified mobility formula for $q\ll 1$ obtained by inserting (3.17) into (5.20) (dashed lines); permeabilities as labelled.

Figure 9.

Same as 8, except for a particle moving parallel to a plane wall in shear flow.

Inserting (3.17) into (5.20) yields a simplified analytical mobility formula depicted by the dashed lines in figures 8–9. Evaluating the result at contact yields

(5.21)\begin{equation} \bar U_1=\bar\omega_1=\frac{d_1}{\log K^{{-}1}+d_2}, \end{equation}

indicating that permeable particles roll without slipping. Here, $d_1=3.125$, $d_1=7.280$ for problems (I), (II), respectively, and $d_2=6.666$ for both problems. This finding illustrates the non-singular rolling motion accessible to permeable particles at contact, as discussed in § 4.4. Hard spheres, by contrast, become stationary at contact and exhibit a slipping approach to contact with $\bar U_1>\bar \omega _1$ (Goldman et al. Reference Goldman, Cox and Brenner1967a,Reference Goldman, Cox and Brennerb), as seen in figures 8–9.