1. Introduction
The random and irreversible displacements of particles driven by shear flow, or shear-induced diffusion, are a fundamental feature of non-Brownian suspensions that emerges once the concentration exceeds the dilute regime. Beyond enhancing mixing and dispersion in industrial and natural flows (Lopez & Graham Reference Lopez and Graham2008; Wang et al. Reference Wang, Koch, Yin and Cohen2009; Metzger, Rahli & Yin Reference Metzger, Rahli and Yin2013; Souzy et al. Reference Souzy, Yin, Villermaux, Abid and Metzger2015), shear diffusion reflects the underlying particle-scale interactions and mechanisms that govern suspension transport under shear. Consequently, many have characterised (Eckstein, Bailey & Shapiro Reference Eckstein, Bailey and Shapiro1977; Leighton & Acrivos Reference Leighton and Acrivos1987; Breedveld et al. Reference Breedveld, van den Ende, Tripathi and Acrivos1998, Reference Breedveld, van den Ende, Bosscher, Jongschaap and Mellema2001, Reference Breedveld, Van Den Ende, Bosscher, Jongschaap and Mellema2002; Metzger et al. Reference Metzger, Rahli and Yin2013; Zhang et al. Reference Zhang, Pham, Metzger, Kopelevich and Butler2023) and modelled (da Cunha & Hinch Reference da Cunha and Hinch1996; Foss & Brady Reference Foss and Brady1999; Marchioro & Acrivos Reference Marchioro and Acrivos2001; Sierou & Brady Reference Sierou and Brady2004; Lopez & Graham Reference Lopez and Graham2007; Gallier Reference Gallier2014; Zhang, Kopelevich & Butler Reference Zhang, Kopelevich and Butler2024) the diffusivity.
Many of these studies emphasise the importance of particle–particle contact interactions in determining shear-induced diffusivity, contributing to the experimental and computational evidence that contacts significantly influence transport in semi-dilute and dense suspensions (Lemaire et al. Reference Lemaire, Blanc, Claudet, Gallier, Lobry and Peters2023). For example, experiments (Zhang et al. Reference Zhang, Pham, Metzger, Kopelevich and Butler2023) indicate that roughening particles reduce the diffusivity at high concentrations, in contrast to expectations based on dilute suspension models (da Cunha & Hinch Reference da Cunha and Hinch1996). Simulations also indicate that friction between particles can significantly enhance diffusivity at high volume fractions and improve agreement with experimental data (Gallier Reference Gallier2014). A more recent study (Zhang et al. Reference Zhang, Kopelevich and Butler2024) independently varied particle roughness and friction, confirming that friction increases diffusivity at high concentrations and improves agreement with measurements. That work also showed that roughness reduces diffusivity under low-friction conditions by promoting flow-aligned layering, but increases it under high-friction conditions by enhancing rotational fluctuations.
The simulations by Zhang et al. (Reference Zhang, Kopelevich and Butler2024) considered an idealised case of monodisperse particles in a fully periodic domain. In contrast, real suspensions are bounded and typically polydisperse, and both factors can significantly alter suspension dynamics and affect measurements of shear-induced diffusion. The detailed particle dynamics are qualitatively altered in the vicinity of bounding walls (Zurita-Gotor, Blawzdziewicz & Wajnryb Reference Zurita-Gotor, Blawzdziewicz and Wajnryb2007) and confinement can suppress shear-induced diffusivity (Leighton & Acrivos Reference Leighton and Acrivos1987; Yeo & Maxey Reference Yeo and Maxey2010a ). At large concentrations, particles near walls form layered structures, as observed in experiments (Blanc et al. Reference Blanc, Lemaire, Meunier and Peters2013; Metzger et al. Reference Metzger, Rahli and Yin2013) and in simulations with and without frictional contacts (Singh & Nott Reference Singh and Nott2000; Yeo & Maxey Reference Yeo and Maxey2010b ; Gallier et al. Reference Gallier, Lemaire, Lobry and Peters2016). This layering modifies the local dynamics, leading to diffusivities near the wall that can differ markedly from those in the bulk.
While monodisperse suspensions can develop ordered structures under shear, polydispersity tends to suppress such organisation. In rheological studies, polydisperse suspensions have been found to exhibit lower viscosities and weaker normal stress differences compared with monodisperse systems at the same volume fraction (Pednekar, Chun & Morris Reference Pednekar, Chun and Morris2018). To prevent crystallisation, bidisperse suspensions are commonly used when simulating rheology of frictional suspensions (Gallier et al. Reference Gallier, Lemaire, Peters and Lobry2014; Mari et al. Reference Mari, Seto, Morris and Denn2014). Singh & Saitoh (Reference Singh and Saitoh2023) used bidisperse particles in two-dimensional simulations to examine the rheology of frictional suspensions and examined correlating the stress with particle diffusivity. In addition to affecting rheology, size segregation occurs when shearing polydisperse suspensions; it occurs radially in torsional flows (Krishnan, Beimfohr & Leighton Reference Krishnan, Beimfohr and Leighton1996) and smaller particles tend to migrate towards bounding walls in planar shear flows (Howard, Maxey & Gallier Reference Howard, Maxey and Gallier2022). This migration is accompanied by differences in particle mobility, with smaller particles exhibiting higher diffusivities than larger ones (Abbas et al. Reference Abbas, Climent, Simonin and Maxey2006).
While studies have examined the influence of polydispersity and confinement on shear induced diffusion, their influence in frictional and rough particle suspensions remains largely unexplored. Building on the study of Zhang et al. (Reference Zhang, Kopelevich and Butler2024), which systematically characterised the roles of surface roughness and interparticle friction for monodisperse particles in unbounded domains, the present work extends that framework to confined and polydisperse suspensions. By holding
$\phi = 0.45$
fixed, the variations in diffusion with confinement, roughness and friction, with and without polydispersity, are isolated for the first time. The model and simulations are described in the next section. Results are reported and discussed in § 3 for two friction coefficients and roughness values, spanning a broad range of confinement and polydispersity. The conclusions in § 4 then summarise the main results concerning the effects of polydispersity and confinement on shear-induced diffusion.
2. Methods
The simulation method extends the approach developed by Zhang et al. (Reference Zhang, Kopelevich and Butler2024) for monodisperse suspensions with periodic boundaries to handle confined geometries and polydisperse particle distributions. As in prior studies of concentrated suspension rheology (Gallier Reference Gallier2014; Mari et al. Reference Mari, Seto, Morris and Denn2014), the model includes short-range lubrication forces, forces to maintain excluded volume and interparticle friction. Particle motion is computed in the Stokesian limit of zero inertia, where the sum of hydrodynamic and contact forces and torques on each particle
$\alpha$
must vanish,
\begin{equation} \biggl (\begin{array}{l} {\boldsymbol{F}}_\alpha ^{(H)} \\ \boldsymbol{T}_\alpha ^{(H)} \end{array}\biggl ) + \biggl (\begin{array}{l} {\boldsymbol{F}}_\alpha ^{(C)} \\ \boldsymbol{T}_\alpha ^{(C)} \end{array}\biggl ) = \boldsymbol{0}. \end{equation}
The hydrodynamic force and torque on particle
$\alpha$
arise from Stokes drag and pairwise lubrication interactions,
\begin{equation} {\boldsymbol{F}}_\alpha ^{(H)} = a_{\alpha }({\boldsymbol{U}}^\infty (\boldsymbol{r}_\alpha ) - \boldsymbol{U}_\alpha ) +\sum _{\beta \ne \alpha }\!{\boldsymbol{F}}_{{\alpha \beta }}^{(L)} \end{equation}
and
\begin{equation} \boldsymbol{T}_\alpha ^{(H)} = (4/3)a^3_{\alpha }(\boldsymbol{\varOmega }^\infty - \boldsymbol{\varOmega }_\alpha ) + \sum _{\beta \ne \alpha }\boldsymbol{T}_{{\alpha \beta }}^{(L)}. \end{equation}
All lengths are non-dimensionalised by the mean particle radius
$\overline {a}$
, such that the dimensionless radius of particle
$\alpha$
is
$a_\alpha$
. Time is scaled by the inverse shear rate
$\dot {\gamma }^{-1}$
and forces by
$6\pi \eta \dot {\gamma } \overline {a}^2$
, where
$\eta$
is the viscosity of the suspending Newtonian fluid. The translational and angular velocities of particle
$\alpha$
are
${\boldsymbol{U}}_\alpha$
and
$\boldsymbol{\varOmega }_\alpha$
, and its position is
$\boldsymbol{r}_\alpha$
. The imposed shear flow is in the
$x$
-direction with gradient in the
$y$
-direction:
${\boldsymbol{U}}^\infty (\boldsymbol{r}) = y \boldsymbol{\delta }_x$
and angular velocity
$\boldsymbol{\varOmega }^\infty = -\boldsymbol{\delta }_z / 2$
.
Lubrication forces and torques,
${\boldsymbol{F}}_{{\alpha \beta }}^{(L)}$
and
$\boldsymbol{T}_{{\alpha \beta }}^{(L)}$
, are computed pairwise as detailed in the Appendix. These short-range hydrodynamic interactions (HI) are applied only when the normalised distance between particle surfaces,
$h_{{\alpha \beta }}=2s_{{\alpha \beta }}/(a_{\alpha }+a_{\beta })$
, is less than 0.5. Here,
$s_{{\alpha \beta }} = r_{{\alpha \beta }}-a_{\alpha }-a_{\beta }$
is the surface separation distance, where
$r_{{\alpha \beta }}$
is the length of the vector
$\boldsymbol{r}_{{\alpha \beta }} = \boldsymbol{r}_\beta - \boldsymbol{r}_\alpha$
. Long-range HI are neglected, as the focus is on concentrated suspensions at a volume fraction
$\phi = 0.45$
, where lubrication dominates. Zhang et al. (Reference Zhang, Kopelevich and Butler2024) showed that simulations without long-range HI at this concentration closely match results from Stokesian Dynamics simulations (Sierou & Brady Reference Sierou and Brady2004) that include them. However, long-range interactions may exert a greater influence when frictional contacts are present (Gallier, Peters & Lobry Reference Gallier, Peters and Lobry2018).
The contact forces include contributions that enforce excluded volume constraints between particles and with any confining walls, as well as interparticle friction,
\begin{equation} {\boldsymbol{F}}_\alpha ^{(C)} = \sum _{\beta \ne \alpha } \left ( {\boldsymbol{F}}_{{\alpha \beta }}^{(\textit{EV})} + {\boldsymbol{F}}_{{\alpha \beta }}^{(F)} \right ) + {\boldsymbol{F}}_\alpha ^{(W)}. \end{equation}
Only friction contributes to the contact portion of the particle torques,
\begin{equation} \boldsymbol{T}_\alpha ^{(C)} = \sum _{\beta \ne \alpha } \boldsymbol{T}_{{\alpha \beta }}^{(F)}. \end{equation}
A Hertzian contact force,
${\boldsymbol{F}}_{{\alpha \beta }}^{(\textit{EV})} = -k_n\lvert s_{\alpha \beta }-2\epsilon \rvert ^{3/2}\boldsymbol{n}_{{\alpha \beta }}$
, is applied between particles
$\alpha$
and
$\beta$
when
$s_{{\alpha \beta }} \leqslant 2\epsilon$
. Here,
$\epsilon$
is the particle roughness, which is assumed to be the same for all particles. The force repels the particles along their common normal,
$\boldsymbol{n}_{{\alpha \beta }} = \boldsymbol{r}_{{\alpha \beta }}/r_{{\alpha \beta }}$
. The stiffness coefficient is set to
$k_n = (\epsilon /10)^{-3/2}$
, following the rationale provided by Zhang et al. (Reference Zhang, Kopelevich and Butler2024).
When particles are in contact (
$s_{{\alpha \beta }} \leqslant 2\epsilon$
), friction can be included using the Cundall–Strack spring algorithm (Cundall & Strack Reference Cundall and Strack1979; Luding Reference Luding2008). The frictional force and torque are given by
\begin{equation} {\boldsymbol{F}}_{{\alpha \beta }}^{F} = k_t\boldsymbol{\xi }_{{\alpha \beta }} \quad \text{and}\quad \boldsymbol{T}_{{\alpha \beta }}^{F} = a_{\alpha }\boldsymbol{n}_{{\alpha \beta }}\times {\boldsymbol{F}}_{{\alpha \beta }}^{F}, \end{equation}
where
$\boldsymbol{\xi }_{{\alpha \beta }}$
is the tangential spring vector, which tracks the integrated tangential displacement during contact, and
$k_t = 20/(7\epsilon )$
is the spring stiffness. At the onset of contact, the spring length is initialised to zero. It evolves according to the relative tangential velocity of the particle surfaces
$\Delta {\boldsymbol{U}}_{{\alpha \beta }}$
, such that the trial update at strain
$\gamma + \Delta \gamma$
is
\begin{equation} \boldsymbol{\xi }^{(t)}_{{\alpha \beta }} (\gamma + \Delta \gamma ) = \boldsymbol{\xi }_{{\alpha \beta }} (\gamma ) + \Delta {\boldsymbol{U}}_{{\alpha \beta }} \Delta \gamma . \end{equation}
Then, the magnitude of
$\boldsymbol{\xi }_{{\alpha \beta }} (\gamma + \Delta \gamma )$
is set to the smaller magnitude of
$\boldsymbol{\xi }^{(t)}_{{\alpha \beta }} (\gamma + \Delta \gamma )$
and
$\mu {\boldsymbol{F}}_{{\alpha \beta }}^{(\textit{EV})} / k_t$
, and its direction is projected onto the plane perpendicular to
$\boldsymbol{n}_{{\alpha \beta }}$
. This algorithm closely approximates Coulomb friction with friction coefficient
$\mu$
; additional details on its implementation and parameter choices are given by Zhang et al. (Reference Zhang, Kopelevich and Butler2024).
If bounding walls are included at
$y=\pm H/2$
, an additional Hertzian contact force is applied to each particle in contact with them,
\begin{equation} {\boldsymbol{F}}^{(W)}_{\alpha } = \left \{ \begin{array}{ll} k_n \left | (H/2 - a_\alpha - \lvert y_\alpha \rvert ) - \epsilon _{\alpha ,w} \right |^{3/2} \boldsymbol{n}_{w} & \mbox{if } \quad H/2 - a_\alpha - \lvert y_\alpha \rvert \le \epsilon _{\alpha ,w}, \\ 0 & \mbox{if } \quad H/2 - a_\alpha - \lvert y_\alpha \rvert \gt \epsilon _{\textbf {}\alpha ,w}. \end{array} \right . \end{equation}
Here,
$\epsilon _{\alpha ,w} = \epsilon a_\alpha /\bar {a}$
is the roughness used for particle–wall interactions. The wall normal vector
$\boldsymbol{n}_{w}$
points into the fluid and
$y_\alpha$
is the position of particle
$\alpha$
in the gradient direction. Lubrication interactions and friction between the walls and particles are not included.
Simulations were conducted for both bounded and unbounded systems. For simulations with boundaries,
$H$
was set to 10, 20 and 40 to test the effects of confinement. Periodic boundaries of dimension 20 were used in the flow (
$x$
) and vorticity (
$z$
) directions, regardless of the value of
$H$
. Unbounded systems were simulated using a fully periodic cubic domain of size
$20\times 20\times 20$
. Depending on the level of confinement and polydispersity, the number of particles in a simulation ranged from 430 to 1719.
Simulations were performed with monodisperse, bidisperse and polydisperse particle size distributions. In monodisperse systems, all particles had the same radius,
$a_\alpha = 1$
for all
$\alpha$
. In bidisperse systems, two particle sizes with a 1.4 size ratio were used, with each size contributing equally to the total volume fraction and the mean particle radius maintained at 1, as in prior studies (O’Hern et al. Reference O’Hern, Silbert, Liu and Nagel2003; Mari et al. Reference Mari, Seto, Morris and Denn2014; Gallier et al. Reference Gallier, Peters and Lobry2018). For polydisperse systems, the particle radii were selected from a normal distribution with mean 1 and standard deviation
$\sigma$
. To exclude extreme sizes, radii smaller than 0.1 or larger than 2.5 were discarded.
Particle velocities were obtained by solving (2.1) and particle positions were updated using the fourth-order Adams–Bashforth method (Hairer, Nørsett & Wanner Reference Hairer, Nørsett and Wanner1993; Butcher Reference Butcher2016). In these simulations, where strain rather than time governs evolution, a strain increment of
$2\boldsymbol{\times }10^{-4}$
was used for particles with roughness
$\epsilon = 2\boldsymbol{\times }10^{-3}$
and
$5\boldsymbol{\times }10^{-4}$
for particles with
$\epsilon = 2\boldsymbol{\times }10^{-2}$
. Particles with
$\epsilon = 2\boldsymbol{\times }10^{-2}$
are referred to as rough and those with
$\epsilon = 2\boldsymbol{\times }10^{-3}$
as smooth. Convergence of the step size was confirmed by verifying that excluded volume was maintained and that the computed diffusivities were insensitive to further reductions. For smooth particles, where the effective roughness range is smaller, a smaller strain increment was required to accurately resolve the short-range contact forces that act over this reduced distance.
The system was initialised in a non-overlapping configuration by first reducing particle sizes to 80 % of their target radius at a volume fraction of
$0.45$
, then increasing them to their full size while simulating shear flow using the simplified model from Zhang et al. (Reference Zhang, Pham, Metzger, Kopelevich and Butler2023) until all overlaps were resolved. After that, simulations with the more detailed model described previously were performed until steady state was reached. Steady state was identified by monitoring the number of contacts per particle, with most periodic systems stabilising at a strain of 20 or less. However, the monodisperse periodic system with rough frictionless particles (
$\epsilon = 2\boldsymbol{\times }10^{-2}$
) required a strain of 100. Systems with bounding walls take longer to reach steady state: in monodisperse suspensions of smooth (
$\epsilon = 2\boldsymbol{\times }10^{-3}$
) and rough particles, steady state required approximately 30 and 100 strains, respectively. In polydisperse bounded suspensions, it required up to 200 strains, due to the strains needed for the particles to demix.
After reaching steady state, systems were simulated for a minimum of 50 additional strains to ensure sufficient statistical sampling of their dynamics and microstructure. To account for variability due to initial configurations, standard errors were estimated from eight independent simulation runs. Diffusivities were calculated by integrating the velocity autocorrelation function,
$C_{\textit{yy}}(\gamma ) = \langle U_{\!y}(\gamma ')U_{\!y}(\gamma '+\gamma )\rangle$
for strains
$\gamma '$
larger than that needed to reach steady state. These values closely agreed with those obtained from the slope of the mean square displacement, except in highly polydisperse systems where the displacements of the smaller, more mobile particles complicated the analysis. In such cases, the autocorrelation-based diffusivities were less noisy and more reliable. Likewise, position-dependent diffusivities,
$D(y)$
, were computed from the local slope of the mean square displacement profile, averaged over particles located near position
$y$
, rather than using velocity autocorrelations.
3. Results and discussion
Section 3.1 presents results for suspensions of monodisperse particles confined between parallel walls, highlighting the effects of confinement on shear-induced diffusion. Section 3.2 explores how polydispersity modifies these effects. As noted previously, all simulations were performed at a volume fraction of
$\phi =0.45$
, where long-range HI are expected to be negligible in the presence of frictional contacts.
3.1. Effect of confinement
Figure 1(a) shows velocity autocorrelation functions in the gradient direction from simulations with walls separated by a distance of
$H = 20$
. The corresponding transverse diffusivities, extracted from these and from simulations at other confinement widths, are summarised in figure 1(b). Fully periodic systems are included for comparison and are denoted by
$H = \infty$
. Only the gradient (transverse) component is presented; vorticity-direction diffusivities follow similar trends as shown in figure S1 and discussed in § S1 of the supplementary material available at https://doi.org/10.1017/jfm.2025.11082.
Figure 1. (a) Representative normalised velocity autocorrelation functions
$C_{\textit{yy}}$
for suspensions with a confinement of
$H=20$
. (b) Dimensionless diffusion coefficient as a function of confinement
$H$
for frictionless (
$\mu =0$
) and frictional (
$\mu =0.5$
) particles that are relatively smooth (
$\epsilon =2\boldsymbol{\times }10^{-3}$
) and rough (
$\epsilon =2\boldsymbol{\times }10^{-2}$
). The data represent an average over all particles and the standard error is smaller than the markers. Simulations using periodic boundary conditions are included for comparison and are denoted by
$H = \infty$
.
For all combinations of particle roughness
$\epsilon$
and friction coefficient
$\mu$
, the transverse diffusivity decreases with increasing confinement (i.e. decreasing
$H$
). At
$H=40$
, the diffusivities of rough frictionless particles (
$\epsilon = 2 \boldsymbol{\times }10^{-2}$
,
$\mu = 0.0$
) are nearly identical to those in the fully periodic system. For the other three cases in figure 1(b), the diffusivities continue to increase between
$H=40$
and
$H=\infty$
, indicating that a wall separation of 40 is not fully sufficient to approximate unbounded conditions.
Consistent with earlier results for frictionless particles (
$\mu = 0$
) (Zhang et al. Reference Zhang, Pham, Metzger, Kopelevich and Butler2023, Reference Zhang, Kopelevich and Butler2024), the smooth particles (
$\epsilon =2\boldsymbol{\times }10^{-3}$
) exhibit higher diffusivity than the rough ones under periodic boundary conditions. This trend persists as confinement increases, though the difference between them narrows slightly as
$H$
approaches 10. As also noted previously (Gallier Reference Gallier2014; Zhang et al. Reference Zhang, Kopelevich and Butler2024), the presence of friction enhances diffusivity at high volume fractions, such as
$\phi = 0.45$
. Frictional particles consistently exhibit higher diffusivity than frictionless ones, but their sensitivity to confinement is more pronounced: at
$H = 10$
, their diffusivity drops to nearly the same level as that of frictionless, smooth particles.
Confinement also influences the relative diffusivity of smooth and rough frictional particles. While the differences appear small on the scale of figure 1(b), they are statistically significant (i.e. larger than the standard error). For weak confinement (
$H \geqslant 40$
), rough particles are slightly more diffusive than smooth ones, but this trend reverses at strong confinement (
$H = 10$
), where smooth particles become more diffusive.
Figure 2. Snapshots of rough (
$\epsilon = 2\boldsymbol{\times }10^{-2}$
), frictionless particle configurations in confined systems with (a)
$H = 10$
, (b)
$H = 20$
, (c)
$H = 40$
and (d) the unconfined (fully periodic) system. Particle radii have been reduced and all particles projected onto a common plane to better reveal structural features.
Figure 3. Effect of the channel height
$H$
on the volume fraction profiles
$\bar {\phi }(y)$
of suspensions of (a) rough (
$\epsilon = 2\boldsymbol{\times }10^{-2}$
) and (b) smooth (
$\epsilon = 2\boldsymbol{\times }10^{-3}$
) frictionless particles. The profiles are shifted so that the origin corresponds to the box wall and only half of the channel height is shown for the case of
$H=40$
.
Figure 2 shows projections of particle positions onto the gradient–vorticity (
$y$
–
$z$
) plane for rough, frictionless particles at three levels of confinement (
$H = 10$
, 20 and 40), as well as for the unconfined system. In the latter case, frictionless rough particles exhibit an intrinsic layering structure, as reported in earlier work (Zhang et al. Reference Zhang, Kopelevich and Butler2024). The dense zone near
$z = 0$
and the hexagonal-like features visible in figure 2 reflect this natural layering behaviour rather than any confinement effect. As discussed previously, such layering suppresses diffusion because particles remain largely trapped within their layer since attempts to move between layers are frequently reflected by adjacent ones, making inter-layer transitions rare. Consequently, the diffusivity of rough, frictionless particles in the unbounded system is lower than that of the other cases (see figure 1
b).
Introducing confinement further modifies this intrinsic ordering: for the strongest confinement (
$H = 10$
), the particles form five distinct layers, with those adjacent to the walls being most sharply defined and the central layer more diffuse. These layered structures for the bounded suspensions are captured by the volume fraction profiles
$\bar {\phi }(y)$
shown in figure 3. The profile
$\bar {\phi }(y)$
represents the averaged volume fraction across the gradient direction. For each simulation,
$\phi (y)$
is first computed at individual strain steps beyond the initial transient, and then averaged over strain and an ensemble of simulations initialised with different particle configurations. The quantity
$\phi (y_0)\Delta y$
gives the fraction of particle volume located between planes
$y = y_0$
and
$y = y_0 + \Delta y$
. In systems exhibiting layering,
$\bar {\phi }(y)$
displays oscillations with characteristic wavelength
$\lambda \approx 2$
, corresponding roughly to a particle diameter.
In the
$H=10$
case, five peaks are visible in
$\bar {\phi }(y)$
, corresponding to the layers observed in figure 2. The reduced amplitude of the central peaks reflects weaker structural order away from the walls. As
$H$
increases to 20 and 40, the layering near the walls remains strong, but the structure becomes increasingly disordered towards the centre: although oscillations in
$\bar {\phi }(y)$
are still visible at larger
$H$
, their amplitudes decrease towards the channel centre. Similar trends are observed for smooth frictionless particles (figure 3
b), though comparing with figure 3(a) shows that rough frictionless particles maintain more persistent layering, particularly away from the bounding walls.
To quantify the layering more systematically, spectral intensities
$\sigma ^2_\phi (q)$
of
$\bar {\phi }(y)$
were computed and are plotted in figure 4 as a function of wavelength
$\lambda =2 \pi /q$
, where
$q$
is the wavenumber. The spectral intensity was obtained by Fourier transforming each frame of data
$\phi (y)$
to compute
$\hat {\phi }(q)$
, then calculating the power
$|\hat {\phi }(q)|^2$
for each frame and averaging over all frames. For an individual frame, the amplitude of the density profile oscillations with the wavelength
$\lambda$
is
$A_\phi (\lambda ) = 2|\hat {\phi }(2\pi /\lambda )|$
. The spectral intensity therefore corresponds to one quarter of the mean-squared amplitude of the oscillations. In this representation, the reduction in structural ordering with increasing
$H$
is evident from the decreasing peak values of
$\sigma ^2_\phi$
. For all
$H$
, rough frictionless particles (figure 4
a) exhibit higher maximum intensities than smooth frictionless particles (figure 4
b). The wavelength
$\lambda _m$
at which the peak occurs corresponds to the interlayer spacing. It depends on
$H$
, shifting below the expected value of 2 as
$H$
increases. This shift likely reflects the tendency of the structures to form hexagonal arrangements under weaker confinement (
$H \geqslant 20$
), see figure 2(b–d), which imposes a preferred interlayer spacing of
$\lambda \approx 1.8.$
In confined systems, the layers located next to the walls are separated by a distance of at least one particle radius from the wall. This leads to a non-integer number of layers that can fit in a weakly confined system with the preferred interlayer spacing of
$\lambda \approx 1.8$
. To accommodate this, the systems develop disordered regions in the centre of the box (see figure 2
b,c). In contrast, as confinement strengthens, the walls increasingly enforce an integer number of layers, driving the spacing towards
$\lambda = 2$
.
Figure 4. Power spectra of the volume fraction profiles
$\phi (y)$
for (a) rough (
$\epsilon = 2\boldsymbol{\times }10^{-2}$
) and (b) smooth (
$\epsilon = 2\boldsymbol{\times }10^{-3})$
frictionless particles under three confinements and fully periodic conditions (
$H = \infty$
). Spectral intensity
$\sigma ^2_\phi$
is plotted as a function of wavelength
$\lambda$
to quantify the degree and spacing of particle layers.
To compare alignment across systems more directly, figure 5 presents the peak amplitude of the dominant density oscillation, defined as
$\bar {A}_\phi = 2 \sqrt { \sigma ^2_\phi (\lambda _m) }$
, where
$\lambda _m$
is the wavelength near 2 at which the power spectrum exhibits a maximum. Higher values of
$\bar {A}_\phi$
for rough frictionless particles, compared with smooth frictionless and frictional particles (regardless of roughness), indicate stronger and more persistent alignment of the particles in the flow direction and across the gradient direction. Errors in
$\bar {A}_\phi$
are estimated from the values of
$A_\phi (\lambda _{\textit{m}})$
of individual samples, where
$\lambda _{\textit{m}}$
corresponds to the location of the maximum of the average spectral intensity. These errors are smaller than the symbols shown in the figure. Previous studies (Zhang et al. Reference Zhang, Pham, Metzger, Kopelevich and Butler2023, Reference Zhang, Kopelevich and Butler2024) found that the strong alignment in rough frictionless systems reduces diffusivity by trapping particles along streamlines and suppressing cross-streamline collisions. That effect is confirmed in figure 1, where rough particles have lower diffusivity than smooth ones for
$\mu =0$
and all considered values of
$H$
.
As confinement increases, alignment intensifies for all particle types, leading to suppression of transverse diffusivity. The increased structural order reduces the frequency of collisions between adjacent layers, and the presence of bounding walls imposes additional constraints on particle motion. These effects are evident in figure 6, which shows a pronounced decrease in transverse diffusivity near the walls.
Figure 5. Effect of the confinement
$H$
on
$\bar {A}_\phi$
, the mean amplitude of the oscillations of the density profile with wavenumber
$\lambda \approx 2$
, which indicates the extent of particle layering.
$\bar {A}_\phi$
was estimated from the power spectrum of the volume fraction profile.
Figure 6. Position-dependence of diffusivities of (a) rough (
$\epsilon = 2\boldsymbol{\times }10^{-2}$
) and (b) smooth (
$\epsilon = 2\boldsymbol{\times }10^{-3}$
) frictionless particles in confined systems. The black dotted horizontal lines show diffusivities in unconfined systems. The position-dependent diffusivities were obtained by fitting local mean-squared displacements to straight lines.
Figure 7. Position dependence of the amplitude
$A_\phi (y)$
of oscillations of the volume fraction profiles
$\bar {\phi }(y)$
for (a) rough (
$\epsilon = 2\boldsymbol{\times }10^{-2}$
) and (b) smooth (
$\epsilon = 2\boldsymbol{\times }10^{-3}$
) frictionless particles in the confined systems. The black dotted horizontal lines show the root-mean-squared average amplitudes
$\bar {A}_\phi$
in the unconfined (periodic) systems.
Figure 6 also shows qualitative differences in the position-dependent diffusivity depending on particle roughness. For smooth frictionless particles at
$H=40$
, the diffusivity near the centre closely matches that of the periodic system, with significant reductions only near the walls. In contrast, for rough frictionless particles, the diffusivity at the centre exceeds that of the corresponding periodic system. This difference can be understood by examining the local structure, quantified in figure 7. The local amplitude
$A_\phi (y)$
of density oscillations, derived from
$\bar {\phi }(y)$
, is computed as half of the difference between neighbouring maxima and minima. For smooth particles (figure 7
b), the structure throughout much of the domain at
$H=40$
closely resembles that of the periodic system. For rough particles (figure 7
a), the structure at the centre is less aligned than in the unconfined case. This weaker alignment results in more frequent collisions and higher diffusivity. Conversely, in the periodic system, stronger layering reduces collisions and suppresses diffusion.
These structural differences clarify an otherwise misleading observation in figure 1(b), where the average diffusivity of rough frictionless particles appears nearly identical for
$H=40$
and
$H=\infty$
. The spatially resolved data in figures 6 and 7 reveal that this agreement is coincidental: although the mean diffusivities are similar, the underlying structure and local dynamics are quite different. In contrast, for smooth frictionless particles, the average diffusivities at
$H=40$
and
$\infty$
differ appreciably, yet the local structure and diffusivity profiles are quite similar across most of the domain. The difference in both average diffusivity and
$\bar {A}_\phi$
arises primarily from the distinct particle dynamics near the confining walls.
To further investigate correlation between local particle ordering and diffusivity, figure 8 shows dependence of the local diffusivity
$D(y)$
on the local alignment amplitude
$\bar {A}_\phi (y)$
. For rough frictionless particles, the alignment alone does not predict their diffusivity in confined systems and the
$D{-}\bar {A}_\phi$
relations depend on the channel width
$H$
(see figure 8
a). This is likely related to the tendency of these particles to spontaneously form ordered structures even in the absence of walls, which causes subtle differences (not captured by
$\bar {A}_\phi$
) in particle alignment due to competition between the wall-induced and spontaneous particle alignment. However, for the diffusivity of smooth frictionless particles, which form layered structures only in the presence of walls, the correlation between the local alignment and diffusivity does not depend on the channel width
$H$
.
Another factor that may affect the local particle diffusivity is the local particle concentration
$\bar \phi (y)$
. In addition to the oscillations with the wavelength of
$\lambda \approx 2$
, caused by the particle layering,
$\bar \phi (y)$
exhibits a slower variation on a length scale comparable with the channel width. This is evident from figure S2, which shows moving averages
$\bar \phi _m(y)$
of the volume fraction profiles. The moving averages are computed with the window size of 5 particle radii, which substantially reduces the magnitude of the oscillations caused by particle layering. It is observed that in the centre of the channel, the local concentration is slightly lower than near the walls (
$\bar \phi _m(y)$
is between 0.42 and 0.44 in the channel centre versus
$\bar \phi _m = 0.5$
near the walls). A lower concentration would normally imply reduced diffusivity, yet the diffusivity is higher in the centre (see figure 6). Thus, variations in local volume fraction cannot explain the observed trend, and particle alignment and layering are responsible.
Figure 8. Dependence of local diffusivity
$D(y)$
on the local alignment amplitude
$A_\phi (y)$
for (a) rough and (b) smooth frictionless particles.
The transverse diffusivities for rough and smooth frictional particles, also shown in figure 1(b), exhibit trends distinct from the frictionless case. Friction increases diffusivity, particularly in weakly confined or unbounded systems. For
$H=\infty$
, the alignment (as quantified by
$\bar {A}_\phi$
shown in figure 5) for smooth and rough frictional particles is comparable to, though slightly lower than, that of smooth frictionless particles, yet their diffusivity is higher by a factor of approximately 2.5. This is consistent with the idea that weaker alignment permits more inter-layer collisions and mixing, but the increase is not fully explained by the the small difference in
$\bar {A}_\phi$
. Also, this relationship does not hold at moderate confinement. At
$H=40$
, frictional particles are more aligned than smooth frictionless ones (higher
$\bar {A}_\phi$
) but still show higher diffusivity. This suggests that alignment alone is insufficient to explain differences in diffusivity between frictional and frictionless particles.
Previous work (Zhang et al. Reference Zhang, Kopelevich and Butler2024) on systems with periodic boundaries identified a strong correlation between rotational fluctuations and elevated shear-induced diffusion of frictional particles. Figure 9 confirms that friction greatly increases the variance
$\sigma _\varOmega ^2$
of the rotational velocity
$\boldsymbol{\varOmega }$
. This increase in the variance is accompanied by enhanced velocity fluctuations, which is a direct result of the frictional contacts. The increased velocity fluctuations result in a larger diffusivity as well as a weaker layering. The microstructure and position dependence of diffusivity of both rough and smooth frictional particles are qualitatively similar to those of smooth frictionless particles (see figures S3–S8) which exhibit a weaker layering than the rough frictionless particles. However, there is a notable difference under the strongest confinement
$(H = 10)$
: the frictional particles exhibit an approximately constant local diffusivity
$D(y)$
near the channel centre, independent of the local degree of alignment
$\bar {A}_\phi (y)$
(see figures S6 and S8). This is likely caused by limits on diffusivity imposed by the confinement; the diffusivity in the centre of the channel with
$H = 10$
does not exceed 0.03 in all considered systems.
As confinement increases,
$\sigma _\varOmega ^2$
decreases, though it remains significantly higher for frictional particles, differing by an order of magnitude from frictionless ones at
$H=10$
. This decline in
$\sigma _\varOmega ^2$
, along with increased
$\bar {A}_\phi$
, is consistent with the reduction in diffusivity observed at smaller wall spacings.
Figure 9. Dependence of the variance
$\sigma _\varOmega ^2$
of the rotational velocity
$\boldsymbol{\varOmega }$
on the confinement
$H$
for smooth and rough particles with (
$\mu =0.5$
) and without (
$\mu =0.0$
) friction.
Figure 10. Dimensionless diffusion coefficient as a function of polydispersity
$\sigma$
for frictionless (
$\mu = 0$
) and frictional (
$\mu = 0.5$
) particles, shown for two values of roughness in (a) unbounded systems and (b) confined systems with wall separation
$H = 20$
. In panel (a), additional data for bidisperse particles are shown in blue (smooth,
$\epsilon = 2 \boldsymbol{\times }10^{-3}$
) and green (rough,
$\epsilon = 2 \boldsymbol{\times }10^{-2}$
).
3.2. Effects of particle size polydispersity
To investigate the influence of particle size polydispersity on shear-induced diffusion, simulations were conducted for rough and smooth particles with both frictionless and frictional interactions. Each simulation used a specified size distribution, characterised by the standard deviation of particle radii,
$\sigma$
. Monodisperse suspensions correspond to
$\sigma = 0$
, while bidisperse suspensions, used in a subset of unbounded cases, have an effective
$\sigma$
slightly greater than 0.2. Figure 10 shows dimensionless transverse diffusivities,
$D$
, as a function of
$\sigma$
. Diffusivities in the vorticity direction follow similar trends (see figure S9). Results are shown for two geometries: an unbounded, fully periodic domain (figure 10
a) and a bounded system with wall separation
$H=20$
(figure 10
b). The bidisperse data for unbounded systems are overlaid in panel (a) for reference.
In unbounded systems, transverse diffusivity increases with polydispersity for all particle types. This trend is particularly pronounced for rough, frictionless particles, whose diffusivity is markedly lower than that of all other cases at
$\sigma = 0$
. This is associated with particle layering at
$\sigma = 0$
and a substantial disruption of particle layers due to polydispersity. Even a relatively small polydispersity with
$\sigma = 0.1$
substantially reduces the degree of layering, with
$\bar {A}_\phi$
reducing from 0.066 at
$\sigma = 0$
to 0.031 at
$\sigma = 0.1$
. A snapshot of the latter system shown in figure 11(a) demonstrates that layers present at
$\sigma = 0$
(see figure 2
d) essentially disappear at
$\sigma = 0.1$
. Although the diffusivity of rough, frictionless particles increases with
$\sigma$
and approaches that of the smooth, frictionless particles, it remains consistently lower across the entire range. Increasing the friction coefficient from
$\mu = 0$
to 0.5 raises diffusivity in both smooth and rough systems across all
$\sigma$
. Notably, the bidisperse data align closely with those of the polydisperse suspensions at a matching
$\sigma \approx 0.2$
, indicating that
$\sigma$
serves as a reliable descriptor of the impact of size variation on transport, at least in this case.
Figure 11. Snapshots of polydisperse systems of rough frictionless particles with
$\sigma = 0.1$
: (a) unbounded system; (b) bounded system (
$H = 20$
). Particle radii have been reduced and all particles projected onto a common plane to better reveal structural features.
In confined systems, transverse diffusivity also increases with polydispersity (figure 10
b). For a given combination of roughness, friction and polydispersity, the diffusivity in the confined geometry is consistently lower than in the unbounded case. Also as in the unbounded system, rough, frictionless particles exhibit the lowest diffusivity across the entire range of
$\sigma$
, while increasing the friction coefficient enhances diffusivity for both rough and smooth particles.
The observed trends in diffusivity with increasing polydispersity reflect changes in the suspension microstructure that closely parallel the mechanisms described in § 3.1 for monodisperse suspensions. Although microstructural features become more difficult to quantify at high
$\sigma$
, consistent trends emerge. First, confined systems exhibit lower diffusivities than unbounded systems under the same conditions of friction, roughness and polydispersity. This reduction arises from increased structural ordering near the walls (see figures 11
b and S10) and the geometric constraint imposed on cross-stream displacements.
Further evidence of layering near the walls of confined systems is shown in figure 12(a), which presents power spectra of the transverse concentration profiles
$\bar {\phi }(y)$
for the frictionless, rough particles. The position of the peak in the spectra shifts to wavelengths
$\lambda \gt 2$
with increasing
$\sigma$
, and the peak amplitude,
$\bar {A}_\phi ^2/4$
, decreases. These changes indicate that layering becomes weaker and less well-defined as polydispersity increases. The trend in
$\bar {A}_\phi$
is summarised in figure 12(b), which shows a monotonic decrease with
$\sigma$
for all particle types. This breakdown of layering enables more cross-stream motion and collisions, leading to the observed increase in transverse diffusivity.
Figure 12. Effect of polydispersity
$\sigma$
on the microstructure of confined suspensions (
$H=20$
): (a) spectral intensities
$\sigma ^2_\phi$
of the volume fraction profile
$\bar {\phi }$
of rough (
$\epsilon = 2\boldsymbol{\times }10^{-2}$
), frictionless (
$\mu = 0$
) particles at various
$\sigma$
; (b) dependence of the magnitude
$\bar {A}_\phi$
of the oscillations of
$\bar {\phi }$
on polydispersity for smoth and rough particles with and without friction.
While layering strength provides one explanation for differences in diffusivity, it does not account for all observed trends. For example, smooth, frictional and frictionless particles have similar values of
$\bar {A}_\phi$
at high
$\sigma$
, but exhibit different diffusivities. These differences correlate with the rotational velocity fluctuations shown in figure 13. The variance in angular velocity,
$\sigma ^2_\varOmega$
, is higher for frictional particles, consistent with enhanced transverse diffusion. Thus, both layering disruption and rotational dynamics contribute to the increase in diffusivity with polydispersity, and their relative importance depends on friction and surface roughness.
Figure 13. Dependence of the variance
$\sigma _\varOmega ^2$
of the rotational velocity
$\boldsymbol{\varOmega }$
on polydispersity
$\sigma$
of particles in a confined system (
$H = 20$
) for smooth and rough particles with and without friction.
Figure 14. Spatial variation of particle size in confined shear flow (
$H = 20$
) at
$\sigma = 0.4$
: (a) local mean particle radius
$\bar {a}$
and (b) local standard deviation
$\sigma _l$
.
Under confinement, polydisperse suspensions exhibit size-based demixing, in contrast to unbounded systems where no such effect is observed. On the time scale of several hundred strain units, larger particles migrate towards the channel centre while smaller ones accumulate near the walls. This segregation is reflected in the profiles of the local mean particle radius,
$\bar {a}(y)$
, and local standard deviation,
$\sigma _l(y)$
, for the case
$\sigma = 0.4$
shown in figure 14(a) and 14(b). The mean radius increases towards the centre, indicating enrichment of larger particles, while
$\sigma _l(y)$
is lower than the global
$\sigma$
, consistent with each region containing a narrower local size distribution than the imposed size distribution. The influence of particle surface properties is modest:
$\bar {a}(y)$
is largely insensitive to roughness or friction, whereas
$\sigma _l(y)$
tends to be slightly higher near the walls for frictional particles.
The dependence on
$\sigma$
is shown in figures S11 and S12 for additional cases and summarised in figure 15, which plots the normalised demixing magnitude
$(\max \bar {a}-1)/\sigma$
versus
$\sigma$
for all four combinations of roughness and friction. Differences among the conditions are most evident at low
$\sigma$
, where microstructural ordering enhances the influence of surface properties. At higher
$\sigma$
, the curves converge, indicating that strong polydispersity suppresses these structural distinctions and diminishes the effects of roughness and friction on demixing.
Figure 15. Comparison of the extent of demixing as a function of polydispersity
$\sigma$
for simulation results with a gap of
$H=20$
.
To assess how the size segregation affects the local dynamics, the position-dependent transverse diffusivity,
$D(y;\sigma )$
, calculated in confined systems is compared with a corresponding prediction,
$D_{\textit{pred}}(y;\sigma )$
, based on unconfined systems. The prediction is constructed using the local mean particle radius
$\bar {a}(y)$
and local standard deviation
$\sigma _l(y)$
, according to
\begin{equation} D_{\textit{pred}}(y;\sigma ) = \bar {a}(y)^2 D_u\! \left (\frac {\sigma _l(y)}{\bar {a}(y)}\right )\!, \end{equation}
where
$D_u(\sigma )$
is the transverse diffusivity in unbounded suspensions with global polydispersity
$\sigma$
. Here,
$D_u(\sigma _l/\bar {a})$
was obtained by linear interpolation of the data shown in figure 10(a). The local polydispersity
$\sigma _l(y)$
is rescaled by
$\bar {a}(y)$
to account for the fact that
$D_u(\sigma )$
was computed using particles with a global mean radius of one. Likewise, the factor
$\bar {a}(y)^2$
multiplies
$D_u(\sigma )$
to maintain a consistent non-dimensionalisation with respect to the global mean particle size.
Figure 16. Position dependence of diffusivities of polydisperse suspensions in confined system for (a) frictionless and (b) frictional rough particles as well as (c) frictionless and (d) frictional smooth particles. The dashed lines show predictions
$D_{\textit{pred}}(y; \sigma )$
of the diffusivity based on the diffusivity
$D_u$
in the unconfined system and the local particle size distribution (see (3.1)).
Figure 16 compares the position-dependent diffusivities from the simulations and predictions for all four combinations of roughness and friction. In contrast to the monodisperse case, where
$D_u$
is constant and
$D_{\textit{pred}}$
is uniform (see figure 6), the polydisperse prediction varies with position due to the spatial dependence of the local particle size distribution. In all cases, the diffusivity
$D(y;\sigma )$
is highest at the channel centre and decreases towards the walls, reflecting the influence of layering. In the wall region,
$D_{\textit{pred}}$
consistently overestimates
$D(y;\sigma )$
, since it does not incorporate the suppression of cross-stream motion caused by wall-induced structural ordering.
In the channel centre, the agreement depends on the particle properties. For rough, frictionless particles with polydispersity
$\sigma \ge 0.2$
,
$D_{\textit{pred}}$
closely matches
$D(y;\sigma )$
, suggesting that layering is negligible in the central region and that local particle size statistics dominate transport. The negligible layering in the channel centre for
$\sigma \ge 0.2$
is confirmed by the concentration profiles shown in figure S10(a). At
$\sigma = 0.1$
, however, layering in the channel centre persists and the predicted value of diffusivity underestimates the direct simulation result. This mirrors the behaviour observed in the monodisperse system (see figure 6) and may reflect reduction of layering in the channel centre of the confined geometry in comparison with the unconfined suspension.
For smooth frictionless particles with
$\sigma \ge 0.3$
, and rough and smooth frictional particles,
$D_{\textit{pred}}$
systematically overestimates
$D(y;\sigma )$
across the entire channel (see figure 16
b–d). These particles do not form layers in the channel centre (see figure S10
b–d), i.e. their microstructure is expected to be similar to that of an unconfined suspension. Similarly, to the monodisperse suspension, the confined polydisperse suspensions exhibit a slight deviation of the local volume fraction in the channel centre from the bulk value of 0.45, as evident from the moving averages of the concentration profiles shown in figure S13. However, this deviation is very small and is not nearly enough to account for the discrepancy between the predicted and measured diffusivities. This suggests that other factors, beyond layering, particle size distribution and changes in local concentration, contribute to the diffusivity of frictional and smooth frictionless particles in confined polydisperse suspensions.
4. Conclusions
This study examined the effects of confinement, polydispersity and interparticle interactions on shear-induced diffusion in suspensions of non-Brownian spheres at a concentration of 45 % by volume. Particle-resolved simulations with short-range HI, surface roughness and contact friction predict that confinement reduces transverse diffusivity compared with unbounded suspensions, whereas polydispersity increases it. These trends hold across the studied combinations of roughness and friction and are broadly consistent with expectations. However, the changes in diffusivity due to confinement and polydispersity can vary with system parameters and, in some cases, alter the relative importance of roughness and friction.
The suppression of diffusivity under confinement arises from the formation of layered microstructures near the confining walls, which restrict cross-stream mobility, and from a reduction in rotational velocity fluctuations. These effects mirror trends observed in unbounded systems with varying friction and roughness, as reported by Zhang et al. (Reference Zhang, Kopelevich and Butler2024). As discussed in § 3.1 and in detail in our earlier work (Zhang et al. Reference Zhang, Kopelevich and Butler2024), friction reduces relative tangential velocities during collisions, thereby enhancing rotational fluctuations. These fluctuations are transmitted to translational motion, leading to increased diffusivity. In those systems, friction enhances diffusion by promoting rotational velocity fluctuations and disrupting strong layering, while in the absence of friction, increasing surface roughness suppresses diffusion by promoting particle layers. This layering reduces collisions and cross-stream motion.
Polydispersity disrupts layering and enhances diffusion in both confined and unconfined systems. In particular, high levels of polydispersity largely eliminate the difference in diffusivity between rough and smooth particles in the frictionless case, regardless of confinement. Polydispersity also induces size segregation in confined suspensions. Larger particles concentrate near the channel centreline, while smaller ones migrate towards the walls.
These findings highlight the complex interplay between confinement, particle size distribution and interparticle interactions in determining shear-induced diffusion. While accurate predictions require knowledge of all four parameters, including channel width (
$H$
), polydispersity (
$\sigma$
), roughness (
$\epsilon$
) and friction coefficient (
$\mu$
), the results suggest that diffusivity becomes relatively insensitive to roughness when either strong friction or sufficient polydispersity is present, and especially when both are. Nonetheless, the diffusivity of rough, frictionless particles remains at least slightly lower than that of smooth ones across all confinement and polydispersity levels examined in this suspension with volume fraction 0.45.
Acknowledgements
Acknowledgement is made to the donors of the American Chemical Society Petroleum Research Fund for partial support of this research through grant #61501-ND9. D.K. acknowledges partial support of the Harry and Bertha Bernstein Professorship at the University of Florida.
Supplementary material
Supplementary material is available at https://doi.org/10.1017/jfm.2025.11082.
Declaration of interests
The authors report no conflict of interest.
Appendix. Lubrication interactions
The lubrication force and torque are given by
\begin{equation} \left (\begin{array}{l} {\boldsymbol{F}}_{{\alpha \beta }}^{(L)} \\[5pt] {\boldsymbol{F}}_{{\beta \alpha }}^{(L)} \\[5pt] \boldsymbol{T}_{{\alpha \beta }}^{(L)} \\[5pt] \boldsymbol{T}_{{\beta \alpha }}^{(L)} \end{array}\right ) = \boldsymbol{R}_{{\alpha \beta }}^{(L)}\boldsymbol{\cdot }\left (\begin{array}{c} {\boldsymbol{U}}^\infty (\boldsymbol{r}_\alpha ) - {\boldsymbol{U}}_\alpha \\[5pt] {\boldsymbol{U}}^\infty (\boldsymbol{r}_\beta ) - {\boldsymbol{U}}_\beta \\[5pt] \boldsymbol{\varOmega }^\infty - \boldsymbol{\varOmega }_\alpha \\[5pt] \boldsymbol{\varOmega }^\infty - \boldsymbol{\varOmega }_\beta \end{array}\right ) + \boldsymbol{\widehat {R}}_{{\alpha \beta }}^{(L)}:\left (\begin{array}{l} {\boldsymbol{E}}^\infty \\[5pt] {\boldsymbol{E}}^\infty \end{array}\right )\!, \end{equation}
where
${\boldsymbol{E}}^\infty$
is the rate of strain tensor, and
$\boldsymbol{R}_{{\alpha \beta }}^{(L)}$
and
$\boldsymbol{\widehat {R}}_{{\alpha \beta }}^{(L)}$
are the resistance matrices (Jeffrey & Onishi Reference Jeffrey and Onishi1984; Jeffrey Reference Jeffrey1992; Kim & Karrila Reference Kim and Karrila2005),
\begin{equation} \boldsymbol{R}_{{\alpha \beta }}^{(L)} = \left (\begin{array}{llll} {\boldsymbol{A}}^{(\alpha \alpha )}& {\boldsymbol{A}}^{({\alpha \beta })} & \boldsymbol{\widetilde {B}}^{(\alpha \alpha )}& \boldsymbol{\widetilde {B}}^{({\alpha \beta })}\\ {\boldsymbol{A}}^{({\beta \alpha })} & {\boldsymbol{A}}^{({\beta \beta })} & \boldsymbol{\widetilde {B}}^{({\beta \alpha })} & \boldsymbol{\widetilde {B}}^{({\beta \beta })}\\ {\boldsymbol{B}}^{(\alpha \alpha )}& {\boldsymbol{B}}^{({\alpha \beta })} & \boldsymbol{C}^{(\alpha \alpha )} & \boldsymbol{C}^{({\alpha \beta })}\\ {\boldsymbol{B}}^{({\beta \alpha })} & {\boldsymbol{B}}^{({\beta \beta })} & \boldsymbol{C}^{({\beta \alpha })} & \boldsymbol{C}^{({\beta \beta })} \end{array}\right ) \end{equation}
and
\begin{equation} \boldsymbol{\widehat {R}}_{{\alpha \beta }}^{(L)} = \left (\begin{array}{ll} \boldsymbol{\widetilde {G}}^{(\alpha \alpha )} & \boldsymbol{\widetilde {G}}^{({\alpha \beta })} \\ \boldsymbol{\widetilde {G}}^{({\beta \alpha })} & \boldsymbol{\widetilde {G}}^{({\beta \beta })} \\ \boldsymbol{\widetilde {H}}^{(\alpha \alpha )} & \boldsymbol{\widetilde {H}}^{({\alpha \beta })} \\ \boldsymbol{\widetilde {H}}^{({\beta \alpha })} & \boldsymbol{\widetilde {H}}^{({\beta \beta })} \end{array}\right )\!. \end{equation}
Here,
$\boldsymbol{A}$
,
$\boldsymbol{B}$
and
$\boldsymbol{C}$
are second-order tensors, and
$\boldsymbol{G}$
and
$\boldsymbol{H}$
are third-order tensors. These tensors obey the following symmetry relations:
\begin{align} A_{\textit{ij}}^{(\alpha \beta )} &= A_{\textit{ji}}^{(\beta \alpha )}, \end{align}
\begin{align} B_{\textit{ij}}^{(\alpha \beta )} &= \widetilde {B}_{\textit{ji}}^{(\beta \alpha )}, \end{align}
\begin{align} C_{\textit{ij}}^{(\alpha \beta )} &= C_{\textit{ji}}^{(\beta \alpha )}. \\[9pt] \nonumber \end{align}
The tensor elements are
\begin{align} A_{\textit{ij}}^{(\alpha \beta )} &= X_{\alpha \beta }^A n_i n_{\!j}+Y_{\alpha \beta }^A (I_{\textit{ij}}-n_i n_{\!j}), \end{align}
\begin{align} B_{\textit{ij}}^{(\alpha \beta )} &= Y_{\alpha \beta }^B\epsilon _{\textit{ijk}}n_k, \end{align}
\begin{align} C_{\textit{ij}}^{(\alpha \beta )} &= X_{\alpha \beta }^C n_i n_{\!j}+Y_{\alpha \beta }^C (I_{\textit{ij}}-n_i n_{\!j}), \end{align}
\begin{align} \widetilde {G}_{\textit{ijk}}^{({\alpha \beta })} &= X_{{\beta \alpha }}^G (n_{\!j} n_k-\frac {1}{3}I_{jk}) n_i+Y_{{\beta \alpha }}^G (n_{\!j} I_{ik} + n_k I_{\textit{ij}}-2n_i n_{\!j} n_k), \end{align}
\begin{align} \widetilde {H}_{\textit{ijk}}^{({\alpha \beta })} &= Y_{{\beta \alpha }}^H(\epsilon _{\textit{jil}}n_l n_k +\epsilon _{kil}n_l n_{\!j}). \end{align}
Here,
$n_i$
is the
$i$
th component of the unit vector pointing from the centre of particle
$\alpha$
to the centre of particle
$\beta$
,
$\boldsymbol{I}$
is the identity matrix and
$\epsilon _{\textit{ijk}}$
is the Levi–Civita symbol. Here,
$X$
and
$Y$
are scalars which depend on the gap
$h_{{\alpha \beta }}$
between particles
$\alpha$
and
$\beta$
. Following Mari et al. (Reference Mari, Seto, Morris and Denn2014), in the current work, we keep only the leading order terms of these scalars,
\begin{align} X_{\alpha \beta } = g^{X_{\alpha \beta }}h_{{\alpha \beta }}^{-1}, \end{align}
\begin{align} Y_{\alpha \beta } = g^{Y_{\alpha \beta }}\ln {h_{{\alpha \beta }}^{-1}}. \end{align}
Here, the coefficients
$g^{X_{\alpha \beta }}$
and
$g^{Y_{\alpha \beta }}$
are functions of the radius ratio of two particles,
$\lambda =a_{\beta }/a_{\alpha }$
. The coefficients
$g^{X_{\alpha \beta }}$
and
$g^{Y_{\alpha \beta }}$
for each element in the resistance matrices are the same as those used elsewhere (Jeffrey & Onishi Reference Jeffrey and Onishi1984; Jeffrey Reference Jeffrey1992; Kim & Karrila Reference Kim and Karrila2005; Mari et al. Reference Mari, Seto, Morris and Denn2014).
Open access
