2. Set-up

We consider the near-contact hydrodynamics of two spherical particles of radius $a$ , whose surfaces are decorated with small roughness asperities (bumps); see figure 1(a). The particle surfaces are separated by a nominal `macroscopic’ gap $D$ , while the `microscopic’ gap $d$ – defined as the minimum separation between the asperities – can be much smaller. Hydrodynamic features are highly localised to the small gaps between opposing asperities, so we focus on a single pair of asperities (one per particle). The bumps are treated as spherical caps with amplitude $A$ , radial width $w$ , a relative orientation angle $\phi$ , and a centre-to-centre distance $s$ in the plane of the gap. The amplitude is small relative to the particle radius, but can be comparable to the gap $D$ .

Figure 1.

(a) Two spherical particles of radius $a$ , with asperities of width $w$ and amplitude $A$ , rotate and translate relative to each other. The microscale separation $d$ between asperities can be much smaller than the nominal separation $D$ between the particles, leading to singular hydrodynamic forces in the approach to contact. (b) Pressure profile in the gap at $y=0$ , showing a sharp spike as the bumps approach contact (decreasing $\mathcal{S}$ ).

We allow the top particle (particle 1) to translate with velocity $V \boldsymbol{e}_x$ and rotate with angular velocity $\varOmega \boldsymbol{e}_y$ , while keeping the lower particle (particle 2) fixed; more general kinematics will be treated elsewhere. We consider small gaps and gentle bump profiles, i.e. $D \ll a$ and $A \ll w$ , providing necessary and sufficient conditions for lubrication theory, but will place no restriction on $A/D$ . Under these geometric conditions, all spherical surfaces can be approximated as paraboloidal. Thus, the rough-particle surfaces are encoded by $z = h_{1, 2}(\boldsymbol{x})$ , defined in (A1), where $\boldsymbol{x} = (x,y)$ represents locations in the plane of the gap; see figure 1(a).

Rather than follow the trajectory of the asperities through time, we characterise the system in configuration space. Defining the gap profile $h(\boldsymbol{x}) = h_1(\boldsymbol{x}) - h_2(\boldsymbol{x})$ and twice the mean gap surface $\overline {h}(\boldsymbol{x})=h_1(\boldsymbol{x})+h_2(\boldsymbol{x})$ , the pressure $p(\boldsymbol{x})$ is governed by the Reynolds equation (see Appendix A):

(2.1) \begin{align} \boldsymbol{\nabla }&\boldsymbol{\cdot }\left [\frac {h^3}{12\mu }\boldsymbol{\nabla }p+\frac {\overline {h}}{2}\left (V-a \varOmega \right )\boldsymbol{e}_x\right ] +\varOmega x=0, \end{align}

where $\mu$ is the viscosity of the fluid and $\boldsymbol{\nabla}$ is the gradient in the $xy$ -plane. We non-dimensionalise by scales characteristic of the smooth limit: distances across the gap by $D$ , horizontal distances by $L = \sqrt {2 a D}$ , and the pressure by $p_c = \mu VL/D^2$ . The rescaled problem involves the dimensionless amplitude $\mathcal{A} = A/D$ and width $\mathcal{W} = w/L$ of the asperities, their relative in-plane orientation $\phi$ and separation $\mathcal{S} = s/L$ , which together set the dimensionless gap $\delta = d/D$ . We numerically solve the rescaled form of (2.1) for the pressure distribution, and then integrate the stresses of the flow to obtain forces and torques exerted by the flow on the particles (see Appendix B).

To aid the discussion, we temporarily suppress rotation and focus on bumps oriented head-on ( $\phi = 0$ ); both constraints will be relaxed later. We first discuss the horizontal force $F_x$ , which we decompose into a combination of smooth $ F_{x,S}$ and rough $F_{x,R} = F_x - F_{x,S}$ contributions. Including Stokes drag, the smooth contribution is (Jeffrey & Onishi Reference Jeffrey and Onishi1984)

(2.2) \begin{align} F_{x, S} \approx 6\pi F_c \left [\frac {1}{6} \log \left (\frac {a}{D}\right ) + 1\right ]\!, \end{align}

where $F_c = \mu a V$ is the characteristic force scale on the particle. The logarithmic term is due to smooth-particle lubrication and depends only weakly on the gap (i.e. it is only approximately as large as Stokes drag even for $D/a$ as small as $10^{-3}$ ). The rough contribution to the force, which is our focus, depends on the inter-bump separation $\delta$ , which we control by varying the in-plane separation $\mathcal{S}$ between the bumps. The bump on particle 1 is centred for simplicity, but this is non-essential for our discussion.

From numerical solutions, we find that the pressure is highly localised in the gap between the asperities (figure 1 b). It is also much greater than the scale expected for smooth spheres. The pressure grows in magnitude and becomes further localised as the bumps are brought closer together by decreasing $\mathcal{S}$ , which decreases the gap width $\delta$ .

Figure 2.

Dimensionless horizontal force (rough contribution) versus gap height $\delta$ for a pair of opposing asperities. The force $F_{x, R}$ scales as $\delta ^{-1}$ for $\mathcal{A} \gt 1/2$ , and exceeds the smooth contribution (shaded region). Results are given at (a) fixed $\mathcal{W}=1/5$ for a range of $\mathcal{A}$ , (b) fixed $\mathcal{A}=0.6$ for a range of $\mathcal{W}$ . The force increases with both $\mathcal{A}$ and $\mathcal{W}$ .

These localised pressure features correlate with increased horizontal forces, unlike for smooth particles where the horizontal force is largely caused by shear. Increasing either the amplitude or width of the asperities increases $F_{x,R}$ (figure 2). When $\mathcal{A} \lt 0.5$ (or $A \lt D/2$ , cf. figure 1 a), all horizontal positions $\mathcal{S}$ are accessible, and the rough contribution to the force remains negligible relative to the smooth case (figure 2 a). This behaviour changes qualitatively as soon as $\mathcal{A}$ exceeds $0.5$ . Now, each asperity occupies more than half the nominal gap width, restricting the available configuration space (limiting the range of $\mathcal{S}$ ). The asperities are now in a state of impending contact under a tangential sliding of the particles (figure 1 a). In this regime, $\delta$ approaches zero as the bumps come horizontally closer, leading to a sharp increase in force. Even for a single pair of asperities, this rough hydrodynamic force can overwhelm the smooth-sphere sliding lubrication force (shaded regions in figure 2). The numerical results suggest that the force diverges as a power-law $\delta ^{-1}$ , which lies in stark contrast with the much weaker (logarithmic) smooth contribution. This scaling is consistent with a squeeze flow between the asperities due to tangential sliding of the particles, suggested in previous work (Wang et al. Reference Wang, Jamali and Brady2020). It is therefore of interest to understand how the prefactors of this scaling behaviour depend on the geometry of the particles and the asperities.

3. Theory

To understand the numerical results we zoom into a `local’ region characterising the gap between the asperities (see the side view of figure 1 a). In this region, the asperities themselves appear as spheres of radius $b= a \kappa /(1 + \kappa )$ , separated by a gap $d$ , where $\kappa = w^2/(2Aa) = \mathcal{W}^2/\mathcal{A}$ is a measure of the bump’s curvature relative to that of the particle. The hydrodynamics of the inter-bump gap are governed by a horizontal length scale $\ell = \sqrt {2 b d} \propto w\sqrt {d/A}$ , which becomes much smaller than both $w$ and $L$ as the bumps approach contact. We exploit this separation of scales to isolate the bumps from the rest of the geometry, and invoke lubrication theory once more, this time in the thin film between the two spherical bump surfaces in relative translation. The force on bump 1 in the local configuration takes the form

(3.1) \begin{align} \boldsymbol{F}=-\mu b \left [R_n\boldsymbol{p}\boldsymbol{p}+R_t\left (\boldsymbol{I}-\boldsymbol{p}\boldsymbol{p}\right )\right ]\boldsymbol{\cdot }(V \boldsymbol{e}_x), \end{align}

where $R_n$ and $R_t$ are dimensionless resistance coefficients for normal and tangential motion between bumps, respectively, $\boldsymbol{I}$ is the identity matrix, and $\boldsymbol{p} = \sin \theta \cos \phi \,\boldsymbol{e}_x+\sin \theta \sin \phi \,\boldsymbol{e}_y-\cos \theta \,\boldsymbol{e}_z$ is the unit vector which connects the centres of the two bumps and is normal to their surfaces (figure 1 a). Also, $\phi$ is the orientation angle introduced earlier, while $\theta$ is a polar angle that depends on the relative position of the bumps.

When the distance between asperities is small ( $d \ll b$ ), the resistances are

(3.2) \begin{align} R_n=\frac {3\pi }{2} \frac {b}{d}, \quad R_t=\frac {\pi }{2}\log \left (1+\frac {A}{d}\right )\!. \end{align}

The normal resistance $R_n$ is a well-known result (Cox Reference Cox1974), whereas the tangential resistance $R_t$ is slightly modified from the usual result as the outer cutoff of the logarithm is set by the bumps’ width $w$ rather than their radius of curvature $b$ (Appendix C). Importantly, $R_n$ scales as $d^{-1}$ and is associated with a squeeze flow between bumps, while $R_t$ , which is due to shear, grows weakly as $\log d^{-1}$ .

Substituting (3.2) into (3.1) leads to general expressions for the force components:

(3.3a) \begin{align} F_{x,R} &=-\mu b V \big[(R_n-R_t)\sin ^2\theta \cos ^2\phi +R_t \big], \end{align}

(3.3b) \begin{align} F_{y,R} &=-\mu b V(R_n-R_t)\sin ^2\theta \sin \phi \cos \phi , \end{align}

(3.3c) \begin{align} F_{z,R} &=\mu b V(R_n-R_t)\sin \theta \cos \theta \cos \phi . \end{align}

The polar angle $\theta$ is a function of the microscopic geometry ( $\mathcal{A}$ , $\mathcal{W}$ and $\mathcal{S}$ ) and thus depends implicitly on $\delta$ . Under a small-angle approximation consistent with the thin-film geometry, we find (Appendix A)

(3.4) \begin{align} \sin \theta \approx \sqrt {\frac {(2A-D+d)(1+\kappa )}{a \kappa (1+2\kappa ) }} . \end{align}

Substituting this result into (3.3a ) yields

(3.5) \begin{align} F_{x,R} &= -\mu a V {\mathcal{K}} \left [\frac {3\pi }{2} \left (\frac {2 \mathcal{A} - 1}{\delta } + 1\right )\cos ^2\phi +R_t\right ]\!. \end{align}

Here, we have defined the dimensionless geometric quantity

(3.6) \begin{align} \mathcal{K} &= \frac {\kappa }{(1+\kappa )(1 + 2 \kappa )} = \frac {(b/a)(1 - b/a)}{(1 + b/a)}, \end{align}

which depends on the ratio of the bump to the particle curvatures; note that ${\mathcal{K}} \sim \kappa = w^2/(2 A a)$ for small $\kappa$ . The tangential rough contribution $R_t$ is much smaller than even the smooth case (Appendix C) so we neglect it for the remainder of the discussion.

The analytic theory (3.5) fully reproduces the features of the numerical results in figure 2. When $\mathcal{A} \lt 1/2$ , the bumps can pass over each other, so the more singular $R_n$ contribution vanishes at closest approach ( $\theta = 0$ , where $\delta = 1 - 2\mathcal{A}$ ) and only the smaller logarithmic $R_t$ contribution survives. By contrast, when $\mathcal{A}\gt 1/2$ each bump is taller than half the nominal gap width, so a sliding motion of the particles decreases their separation between their surfaces down to zero. That is, sliding motion of the particles for $\mathcal{A} \gt 1/2$ necessitates a squeeze flow between asperities, generating a $\delta ^{-1}$ hydrodynamic force that resists contact. Despite its small prefactor, the $\delta ^{-1}$ singularity from a single bump-pair quickly overwhelms the combined `smooth’ Stokes drag and lubrication forces on the rest of the sphere (shaded region in figure 2). An estimate of this crossover for typical experiments in made in § 5.

Figure 3.

Horizontal force $F_{x,R}$ rescaled by the theoretically predicted scale $F_c \mathcal{K} (2 \mathcal{A} - 1)$ , plotted against $\delta$ . Symbols are data from figure 2, which collapse into the theoretical prediction $3\pi /(2\delta)$ for small $\delta$ (solid line) for all $\mathcal{A}$ and $\mathcal{W}$ , without fitting parameters. The collapse persists for larger $\delta$ , and is captured by adding a subdominant empirical constant $c = 15$ to the theoretical prediction (dashed curve).

We quantitatively test the theory in the singular regime $\mathcal{A} \gt 1/2$ . As the bumps approach contact, (3.5) reduces to

(3.7) \begin{align} F_{x,R} \sim - F_c \frac {3\pi {\mathcal{K}}}{2} \frac {(2\mathcal{A}-1)}{\delta } \cos ^2\phi , \quad \mathcal{A}\gt 1/2. \end{align}

We use this result to rescale the numerical results of figure 2 by the scaling factor $(2 \mathcal{A} - 1) {\mathcal{K}}$ . The rescaled data collapse onto the theoretical prediction $3\pi /(2\delta)$ for small $\delta$ , without adjustable parameters (figure 3). Unexpectedly, this collapse persists even when $\delta$ is not small, and is described rather well by a universal curve $-(3 \pi /(2 \delta) + c)$ , where we find the subdominant term $c \approx 15$ from a fit.

Thus, the introduction of roughness with small but finite amplitude $A$ such that $D/2 \leqslant A \ll a$ leads to a singular perturbation of the hydrodynamic resistance. While the smooth contribution scales as $\log (\delta ^{-1})$ , the rough contribution scales more singularly as $\delta ^{-1}$ , albeit with a prefactor that is smaller by a factor of $(b/a)$ . The rough contribution is essentially `switched on’ when $A \gt D/2$ ( $\mathcal{A} \gt 1/2$ ), which is the regime of impending contact between asperities under a sliding of the particles. This threshold is observed clearly in the numerical results, and emerges naturally in the theory. In a situation with multiple bumps, the threshold $A \gt D/2$ corresponds to configurations where particles can interlock, which has been identified as important to DST in past work (Hsu et al. Reference Hsu, Ramakrishna, Zanini, Spencer and Isa2018). However, while the interactions of interlocking particles are normally thought to be driven by contact friction, we see that similar interactions can arise due to hydrodynamics alone in the approach to contact. This can occur only when both surfaces are rough; if one of the surfaces is smooth, only the subdominant $R_t$ term survives since there is no geometric mechanism to generate a squeeze flow. We refer the reader to Yariv et al. (Reference Yariv, Brandão, Wood, Szafraniec, Higgins, Bazazi, Pearce and Stone2024) for a two-dimensional analogue with a single rough surface, and Kargar-Estahbanati & Rallabandi (Reference Kargar-Estahbanati and Rallabandi2025) for a treatment that also includes elastic deformation.

The theory provides a natural framework for arbitrary orientations $\phi$ between asperities. The predicted $\phi$ dependences in (3.3) are confirmed by our numerical solutions (Appendix C). At fixed $\phi$ , all three force components depend on $R_n$ , diverging as $\delta ^{-1}$ for $\mathcal{A} \gt 1/2$ . However, not all three components are equally important on average. Any physically realisable system of rough particles will involve an ensemble of pairwise bump interactions with different orientations at any instant of time. If bumps are distributed isotropically on the particles, all orientations $\phi \in [0, 2\pi )$ are realised with equal probability. Then, averaging (3.3) over $\phi$ leads to the average force per pair of opposing bumps,

(3.8) \begin{align} \langle F_{x,R} \rangle =- F_c \frac {3\pi {\mathcal{K}}}{4} \frac {(2\mathcal{A}-1)}{\delta }, \quad \langle F_{y,R} \rangle = \langle F_{z,R} \rangle =0, \end{align}

valid for $\mathcal{A} \gt 1/2$ (observe that the average force is half of the maximum). Forces in the direction of motion ( $x$ ) survive the averaging since $F_x$ always opposes particle motion regardless of orientation, while the other components average out to zero, due to orientational symmetry at the asperity scale. Thus, the configurational averaging retrieves the macroscopic Stokesian symmetry of the two-sphere system. We note that different averages may be more or less appropriate depending on the constraints of the problem (e.g. prescribed force versus prescribed motion). We expect different averaging methods to lead to similar results up to $O(1)$ prefactors.

4. Coupling of rotation and translation

We now generalise to include rotation $\varOmega \boldsymbol{e}_y$ in addition to translation $V \boldsymbol{e}_x$ , focusing on how these degrees of freedom are coupled by roughness. We also keep track of the torque on the particle, only the $y$ -component of which ( $T_y$ ) is relevant in the average sense for isotropically distributed roughness. For the same reason, we neglect vertical translation, which will not, on average, couple to horizontal translation or rotation due to symmetry. Under this configurationally averaged setting, the mobility relation connecting a tangential force and torque due to tangential sliding and rotation of particle 1 is

(4.1) \begin{align} \begin{bmatrix} F_x \\ T_y \\ \end{bmatrix} = -\mu a \begin{bmatrix} \mathcal{R}_{F V} & a \mathcal{R}_{F \varOmega } \\ a \mathcal{R}_{T V} & a^2 \mathcal{R}_{T \varOmega } \\ \end{bmatrix} \boldsymbol{\cdot }\begin{bmatrix} V \\ \varOmega \\ \end{bmatrix}\!, \end{align}

where the $\mathcal{R}_{ij}$ are dimensionless resistance coefficients (note that $\mathcal{R}_{F\varOmega } = \mathcal{R}_{TV}$ due to symmetry).

The $\mathcal{R}_{FV}$ coefficient relates the force to the particle velocity, and is identical to the quantity $F_{x,R}/F_c$ in figures 2 and 3, where we had set $\varOmega = 0$ . To understand the rough contributions to the other coefficients we revisit the local theory. Including particle rotation merely modifies the approach velocity of the asperities in (3.1) from $V \boldsymbol{e}_x$ to $(V - a \varOmega ) \boldsymbol{e}_x$ , which carries through to all subsequent theoretical results for the force. The rough contribution to the torque on particle 1 is $\boldsymbol{T}_R = \int \boldsymbol{r} \times \boldsymbol{n} \boldsymbol{\cdot }\boldsymbol{\sigma }_R \text{d}S$ , where $\boldsymbol{\sigma }_R$ is the stress due to roughness, $\boldsymbol{n}$ is the normal vector,  and $\boldsymbol{r}$ is the moment arm connecting the centre of the particle to a point on its surface S. As shown in figure 1(b) and predicted by the theory, stresses are localised to a small region of length $\ell \ll L$ surrounding the minimum separation of the bumps. Within this region, the moment arm is effectively the constant vector $\boldsymbol{r}_c$ connecting the centre of particle 1 to the point of minimum separation. The torque is therefore $\boldsymbol{T}_R \approx \boldsymbol{r}_c \times \boldsymbol{F}_R$ , where $\boldsymbol{F}_R = \int \boldsymbol{n} \boldsymbol{\cdot }\boldsymbol{\sigma }_R \text{d}S$ is the rough contribution to the force discussed previously. This cross product does not vanish even though $\boldsymbol{F}_R$ is largely a normal force (directed along $\boldsymbol{p}$ ) since $\boldsymbol{p}$ is misaligned with the moment arm $\boldsymbol{r}_c$ (this is not the case for smooth spheres). Then, the torque on the particle is simply a consequence of the rough force acting at an effective point of contact. We thus conclude that all four resistance coefficients are identical up to sign:

(4.2) \begin{align} \mathcal{R}_{FV} &= -\mathcal{R}_{F\varOmega } = -\mathcal{R}_{TV} = \mathcal{R}_{T\varOmega } = \frac {3 \pi {\mathcal{K}}}{2} \frac {(2 \mathcal{A} - 1)}{\delta } \cos ^2\phi , \quad (\mathcal{A} \gt 1/2). \end{align}

Figure 4.

(a) Resistance coefficient $\mathcal{R}_{TV}$ (showing rough contributions only) versus $\delta$ , for fixed $\mathcal{W}=1/5$ and various $\mathcal{A}$ . The data are virtually identical to those in figure 2(a). Shaded region indicates a typical range for smooth particles for $D/a \in (10^{-3}, 10^{-1})$ . (b) Ratio of resistance coefficients $\mathcal{R}_{TV}/\mathcal{R}_{T \varOmega } = a \varOmega /V$ versus $\delta$ . For all $\mathcal{A}$ and all $\delta$ , the ratio is always close to 1, indicating strong rotation–translation coupling.

Our numerical calculations quantitatively corroborate these predictions. The torque due to translation, $\mathcal{R}_{TV}$ , is virtually identical to $\mathcal{R}_{FV}$ (compare figure 4 a with figure 2 a), diverging as $\delta ^{-1}$ for $\mathcal{A}\gt 1/2$ and exceeding the torque on a smooth particle. If particle 1 is torque-free, its rotation and translation are coupled according to

(4.3) \begin{align} \frac {\varOmega a}{V}=-\frac {\mathcal{R}_{T V}}{\mathcal{R}_{T \varOmega }}. \end{align}

For small gaps, the torque from roughness dominates that from the rest of the two-particle geometry, and this ratio approaches unity according to the theory (4.2), corresponding to $\varOmega a/V \approx 1$ . This prediction is in agreement with the numerical results (figure 4 b), and corresponds precisely to vanishing relative surface velocity $V - a \varOmega$ . This hydrodynamic ‘pure-rolling’ condition recovers the kinematic constraint of classical dry rolling friction. It is interesting to note that the rolling condition persists across the range of $\delta$ , even when the resistance coefficients are not particularly large (note the short range of the vertical axis in figure 4 b). This is due primarily to the localisation of the force on the sub-asperity scale, which allows it to behave similarly to a point-like contact force. We contrast this with the case of smooth particles, where the coupling is negligible (Cox Reference Cox1974; Hsiao et al. Reference Hsiao, Saha-Dalal, Larson and Solomon2017b ), i.e. torque-free smooth particles experience a very weak (logarithmic) coupling that decays rapidly with $D$ . The roughness-induced torque dominates over the smooth contribution close to contact (i.e. when $A \gt D/2$ ), and so we expect rotation and translation to be strongly coupled in this regime. This coupling is insensitive to the geometric details of the roughness, and stems from the generation of large forces and localisation of the stresses on sub-asperity length scales in the interlocking regime.

The constraining of rotational modes to translational ones has been identified as a central piece in the rheology of dense suspensions. Past work has interpreted it as a kinematic constraint that occurs on contact, akin to the rotation of interlocked gears (see e.g. Singh et al. Reference Singh, Ness, Seto, de Pablo and Jaeger2020). The present work provides an alternative interpretation of this feature. If the particles did not roll when slid against each other, there would be relative motion between asperities in near-contact. This would generate squeeze flows, producing $d^{-1}$ singular forces and torques. The need to relieve this impending hydrodynamic singularity as asperities approach (or pull apart) underpins the hydrodynamic interpretation of the rolling constraint.