Frictional effects on shear-induced diffusion in suspensions of non-Brownian particles

Abstract

The shear-induced diffusivity of non-Brownian spheres in monodisperse suspensions undergoing viscous flow was calculated using simulations that account for particle roughness and friction as independent parameters. The diffusivity increases significantly as the friction coefficient is increased, and the effect is largest on rougher particles. Roughness reduces the transverse diffusivities relative to smoother particles for sufficiently concentrated suspensions of frictionless and low-friction particles. However, the diffusivity of roughened particles is larger than smoother ones at high values of the friction coefficient. The increase of the diffusivity with friction is associated with a significant broadening of the variance of the rotational velocities. The most prevalent observation, when correlating the microstructure to changes in diffusivity for frictionless particles, is that less diffusive systems, with larger roughness, form layers along the flow direction. These results confirm previous experimental and simulation results that roughness can decrease diffusivity at large concentrations using a more detailed model. Also, comparisons of the simulation results with previously published experimental measurements indicate that friction improves the alignment of the results with experiments.

1. Introduction

In shearing flows of concentrated suspensions, particle collisions are frequent and significantly alter the dynamics and rheology of the suspensions. Many of these changes cannot be solely attributed to viscous interactions; allowing for solid–solid contacts during collisions, even when the suspending fluid is viscous and inertial effects are minimal, resolves many qualitative discrepancies between models and experimental observations. For instance, contact forces create asymmetry in the pair-distribution functions, as experimentally verified for relatively dilute suspensions (Rampall, Smart & Leighton 1997; Blanc, Peters & Lemaire 2011a), and research indicates that the loss of reversibility of particle positions in oscillatory shear is due to solid contacts rather than chaos from hydrodynamic interactions (Metzger, Pham & Butler 2013a; Pham, Metzger & Butler 2015). Likewise, the abrupt changes in rheology measured upon reversing the direction of a fully developed shearing flow of concentrated suspensions (Kolli, Pollauf & Gadala-Maria 2002; Blanc, Peters & Lemaire 2011b) can be explained by the sudden release of contact forces at the moment of reversal (Bricker & Butler 2007; Peters et al. 2016). These, and other examples, are described in detail in a recent and thorough review by Lemaire et al. (2023).

Simulations incorporate contacts by applying a repulsive force between particle pairs along their common normal vector when the separation distance drops below a specified threshold. This range over which this contact force acts can be correlated with particle roughness (Rampall et al. 1997). Modifying the contact range, or roughness of the particle, can significantly alter the dynamics, even though the range is a small fraction ($10^{-2}$ to $10^{-4}$) of the particle diameter (Lemaire et al. 2023). For example, the contact force causes a net displacement between colliding particle pairs in a shearing flow (Da Cunha & Hinch 1996), leading to an effective shear-induced diffusivity after multiple random collisions. Consequently, larger roughness results in higher diffusivity for small concentrations. However, measurements and simulations indicate that roughened particles diffuse less than smoother ones in concentrated suspensions (Zhang et al. 2023). Corresponding simulations indicated that the rougher particles form layers in the flow direction and that this layered structure results in a reduced diffusivity.

The presence of solid–solid contact suggests that frictional forces may also play a role in the suspension dynamics, and models incorporating friction have been developed to study the rheology of suspensions (Seto et al. 2013; Gallier et al. 2014; Mari et al. 2014). Including friction increases the viscosity, and discontinuous shear thickening (Brown & Jaeger 2014) can be predicted by also adding electrostatic interactions to friction (Seto et al. 2013) or using a load-dependent friction (Mari et al. 2014). These methods have also been used to simulate the pressure-imposed rheology of frictional particles (Athani et al. 2022). However, few results are available regarding the impact of friction on shear-induced diffusion. While simulations predict that friction significantly increases shear-induced diffusivity (Gallier 2014; Singh & Saitoh 2023), these results have not been fully explored, especially in light of recent work which suggests roughness decreases diffusivity at high concentrations (Zhang et al. 2023).

To address the combined effects of friction and roughness on diffusivity, we have developed a simulation that integrates a frictional contact model with lubrication interactions, similar to the approach of Mari et al. (2014). Given the lack of a clear relation between roughness and friction, both were varied independently to explore their effects on shear-induced diffusivities. Our simulations confirm that in sufficiently concentrated suspensions, rough particles exhibit lower transverse diffusivity compared with smooth particles when friction is absent. As the friction coefficient increases, diffusivity rises significantly, with the effect being more pronounced for rougher particles. This larger increase is substantial enough that, at sufficiently high friction coefficients, the diffusivity of roughened particles surpasses that of smoother ones. Comparing the results with experimental data shows improved predictions when friction is included. These findings provide insights into the suspension dynamics and also have practical implications for mixing and heat transfer within flowing viscous suspensions (Lopez & Graham 2008; Wang et al. 2009; Metzger, Rahli & Yin 2013b; Souzy et al. 2015).

2. Methods

A monodisperse suspension of neutrally buoyant non-Brownian spheres of radius $a$ suspended in a Newtonian fluid is considered. Both the fluid and particle inertia are neglected under the assumption that the Reynolds number, $Re=\rho a^2\dot {\gamma }/\eta$, is low. Here, $\rho$ is the density of the fluid and the particles, $\dot {\gamma }$ is the shear rate and $\eta$ is the viscosity of the fluid. Similarly, Brownian diffusion also is ignored since the particles of interest are relatively large (i.e. ${>}10\ {\mathrm {\mu }}{\rm m}$). The equations that follow are presented non-dimensionally using the particle radius $a$ as the characteristic length, the inverse of the shear rate $\dot {\gamma }^{-1}$ as the characteristic time and $6{\rm \pi} \eta \dot {\gamma } a^2$ as the characteristic force.

The suspension motion is driven by a simple shear flow $\boldsymbol {U}^\infty$

(2.1)\begin{equation} \boldsymbol{U}^\infty(\boldsymbol{r})=\boldsymbol{\varOmega}^\infty\times\boldsymbol{r} + \boldsymbol{E}^\infty\boldsymbol{\cdot}\boldsymbol{r}, \end{equation}

where the rate-of-strain tensor $\boldsymbol {E}^\infty$ and vorticity $\boldsymbol {\varOmega }^\infty$ have non-zero elements of ${E}_{12}^\infty ={E}_{21}^\infty =1/2$ and ${\varOmega }_3^\infty =-1/2$. The position, translational and rotational velocities of any particle $\alpha$ are denoted as $\boldsymbol {r}_\alpha$, $\boldsymbol {U}_\alpha$ and $\boldsymbol {\varOmega }_\alpha$, respectively. In the absence of inertia, the forces acting on particles are balanced and the equation of motion is

(2.2)\begin{equation} \left(\begin{array}{@{}c@{}} \boldsymbol{F}^{(H)} \\ \boldsymbol{T}^{(H)} \end{array}\right) + \left(\begin{array}{@{}c@{}} \boldsymbol{F}^{(C)} \\ \boldsymbol{T}^{(C)} \end{array}\right) = \boldsymbol{0}, \end{equation}

where $\boldsymbol {F}$ and $\boldsymbol {T}$ are the force and torque vectors, respectively. The superscripts $H$ and $C$ indicate hydrodynamic and contact forces.

2.1. Hydrodynamic forces

Long-range hydrodynamic interactions are not included in the current model since the goal is to study concentrated suspensions where short-range lubrication forces dominate. Therefore, only the Stokes drag and lubrication interactions were included, and the hydrodynamic force and torque acting on each particle $\alpha$ are

(2.3a,b)\begin{equation} \boldsymbol{F}_\alpha^{(H)} = \boldsymbol{F}_\alpha^{(S)}+\sum_{\beta\ne\alpha}\boldsymbol{F}_{\alpha\beta}^{(L)} \quad {\rm and}\quad\boldsymbol{T}_\alpha^{(H)} = \boldsymbol{T}_\alpha^{(S)}+\sum_{\beta\ne\alpha}\boldsymbol{T}_{\alpha\beta}^{(L)}, \end{equation}

where

(2.4a,b)\begin{equation} \boldsymbol{F}_\alpha^{(S)} = \boldsymbol{U}^\infty(\boldsymbol{r}_\alpha) - \boldsymbol{U}_\alpha \quad\mathrm{and}\quad \boldsymbol{T}_\alpha^{(S)} = (4/3)(\boldsymbol{\varOmega}^\infty - \boldsymbol{\varOmega}_\alpha) \end{equation}

are the Stokes drag force and torque and $\boldsymbol {F}_{\alpha \beta }^{(L)}$ and $\boldsymbol {T}_{\alpha \beta }^{(L)}$ are the short-range lubrication force and torque between particles $\alpha$ and $\beta$. The lubrication forces are applied only when the gap between particle surfaces, $h_{\alpha \beta } = r_{\alpha \beta }-2$, satisfies the constraint $h_{\alpha \beta } \leqslant 0.5$. Here, $r_{\alpha \beta }$ is the length of the vector $\boldsymbol {r}_{\alpha \beta } = \boldsymbol {r}_\beta - \boldsymbol {r}_\alpha$ connecting the particles $\alpha$ and $\beta$. It is shown in figure S1 of the supplementary material available at https://doi.org/10.1017/jfm.2024.1121 that the cutoff lubrication range has a relatively small effect on diffusivity. The short-range lubrication force and torque are given by

(2.5)\begin{equation} \left(\begin{array}{@{}c@{}} \boldsymbol{F}_{\alpha\beta}^{(L)} \\ \boldsymbol{F}_{\beta\alpha}^{(L)} \\ \boldsymbol{T}_{\alpha\beta}^{(L)} \\ \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 \\ \boldsymbol{U}^\infty(\boldsymbol{r}_\beta) - \boldsymbol{U}_\beta \\ \boldsymbol{\varOmega}^\infty - \boldsymbol{\varOmega}_\alpha \\ \boldsymbol{\varOmega}^\infty - \boldsymbol{\varOmega}_\beta \end{array}\right) + \boldsymbol{\hat{R}}_{\alpha\beta}^{(L)}:\left(\begin{array}{@{}c@{}} \boldsymbol{E}^\infty \\ \boldsymbol{E}^\infty \end{array}\right), \end{equation}

where $\boldsymbol {R}_{\alpha \beta }^{(L)}$ and $\boldsymbol {\hat {R}}_{\alpha \beta }^{(L)}$ are resistance matrices (Jeffrey & Onishi 1984; Kim & Karrila 1991; Jeffrey 1992) Their elements are summarized in Appendix A.

2.2. Contact and friction forces

Contact forces are applied when $h_{\alpha \beta }\leqslant 2\epsilon$, where $\epsilon$ is the particle roughness. The normal contact force is expressed using a Hertzian contact model (Pöschel & Schwager 2005; Radjai & Dubois 2011). This model accounts for the deformation of contacting spherical asperities on the surface and the force is given by

(2.6) \begin{equation} \boldsymbol{F}^{(C,n)}_{\alpha\beta} = \left\{\begin{array}{@{}ll@{}} {-}k_n\lvert h_{\alpha\beta}-2\epsilon\rvert^{3/2}\boldsymbol{n}_{\alpha\beta} & \mbox{if } h_{\alpha\beta} \leqslant 2\epsilon, \\ 0 & \mbox{if } h_{\alpha\beta} > 2\epsilon. \end{array} \right. \end{equation}

Here, $k_n$ is the normal stiffness and $\boldsymbol {n}_{\alpha \beta } = \boldsymbol {r}_{\alpha \beta }/r_{\alpha \beta }$ is the unit vector pointing from particle $\alpha$ to particle $\beta$. In principle, the normal stiffness $k_n$ can be estimated from the particle's mechanical properties. However, in this work, the approach of Gallier et al. (2014) is adopted and $k_n$ is chosen so that the repulsive force balances with the characteristic force (i.e. $6{\rm \pi} \eta \dot {\gamma } a^2$ when written dimensionally). Making the balance when the separation distance is $h_{\alpha \beta } = 2(\epsilon -\delta )$ gives

(2.7)\begin{equation} k_n(2\delta)^{3/2} = 1, \end{equation}

where $\delta$ is set to $0.05\epsilon$. This choice for $\delta$ has a negligible effect on diffusivity if $\delta$ is sufficiently small ($\delta \leqslant 0.15\epsilon$), as demonstrated in §§ 3.1 and S1. Also note that the value of $\delta$ is limited from above: a large $\delta$ (and, hence, a small $k_n$) leads to overlap between particles, as detailed in § S1 of the supplementary material.

In addition to the normal contact force, particles in contact ($h_{\alpha \beta }\leqslant 2 \epsilon$) experience friction in the tangential direction. According to Coulomb's law of friction, when two particles are in contact, the magnitude of the friction force $\boldsymbol {F}^{(C,t)}_{\alpha \beta }$ opposing their relative motion is

(2.8)\begin{equation} F^{(C,t)}_{\alpha\beta} = \left\{\begin{array}{@{}ll@{}} \mu F^{(C,n)}_{\alpha\beta}, & \mathrm{if}\ |\Delta\boldsymbol{U}^{(t)}_{\alpha\beta}| \ne 0 \mathrm{\ (dynamic\ regime)},\\ 0\dots\mu F^{(C,n)}_{\alpha\beta}, & \mathrm{if}\ |\Delta\boldsymbol{U}^{(t)}_{\alpha\beta}| = 0\ \textrm{(static}\ \textrm{regime).} \end{array} \right. \end{equation}

Here, the friction coefficient is $\mu$ and

(2.9)\begin{equation} \Delta\boldsymbol{U}^{(t)}_{\alpha\beta} = (\boldsymbol{I}-\boldsymbol{n}_{\alpha\beta}\boldsymbol{n}_{\alpha\beta})\boldsymbol{\cdot}[\boldsymbol{U}_\beta-\boldsymbol{U}_\alpha-(\boldsymbol{\varOmega}_\alpha + \boldsymbol{\varOmega}_\beta)\times\boldsymbol{n}_{\alpha\beta}] \end{equation}

is the relative tangential velocity of the particle surfaces and $\boldsymbol {I}$ is the identity tensor. In (2.8), the notation in the last condition indicates that the static frictional force can take on multiple values and is not known a priori. Also note that the static and dynamic friction coefficients are set to the same value within this model.

The main challenge in a numerical implementation of this friction law is to obtain the static tangential force $F^{(C,t)}_{\alpha \beta }$ such that $\Delta \boldsymbol {U}^{(t)}_{\alpha \beta } = 0$ for a given configuration. To address this challenge, we use the Cundall–Strack tangential spring method (Cundall & Strack 1979; Luding 2008), which is commonly used for simulating granular materials (Walton 1983; Tsuji, Tanaka & Ishida 1992; Lee 1994) and has been incorporated into models of suspensions (Gallier et al. 2014; Mari et al. 2014). This method introduces a fictitious spring so that the tangential force and the torque are given by

(2.10a,b)\begin{equation} \boldsymbol{F}^{(C,t)}_{\alpha\beta}=k_t\boldsymbol{\xi}_{\alpha\beta} \quad \text{and}\quad \boldsymbol{T}^{(C,t)}_{\alpha\beta}=\boldsymbol{n}_{\alpha\beta}\times\boldsymbol{F}^{(C,t)}_{\alpha\beta}, \end{equation}

respectively. Here, $k_t$ is the tangential stiffness and $\boldsymbol {\xi }_{\alpha \beta }$ is the tangential spring vector corresponding to the colliding particles ($\alpha$ and $\beta$). Upon initial contact, the spring length is set to zero and then incremented by the relative displacement of the particle surfaces

(2.11)\begin{equation} \tilde{\boldsymbol{\xi}}_{\alpha\beta}(\gamma+\Delta\gamma) = \boldsymbol{\xi}_{\alpha\beta}(\gamma) + \Delta\boldsymbol{U}^{(t)}_{\alpha\beta}(\gamma)\Delta\gamma, \end{equation}

where $\gamma$ is the strain, $\Delta \gamma$ is the step size and $\tilde {\boldsymbol {\xi }}_{\alpha \beta }$ refers to an intermediate spring vector during a simulation step. The particles are in the static friction regime when $k_t\tilde {\xi }_{\alpha \beta } < \mu F ^{(C,n)}_{\alpha \beta }$ and are in the dynamic regime otherwise. In the latter case, the spring length is reset to ensure that (2.10a,b) yields the correct dynamic friction. Hence, the spring length is

(2.12)\begin{equation} \xi_{\alpha\beta}(\gamma+\varDelta\gamma) = \min\left(\tilde{\xi}_{\alpha\beta}(\gamma+\Delta\gamma),\frac{\mu}{k_t}F^{(C,n)}_{\alpha\beta}(\gamma)\right). \end{equation}

After the spring length is determined, the spring is re-oriented in the tangential direction, i.e.

(2.13)\begin{equation} \boldsymbol{\xi}_{\alpha\beta} = \boldsymbol{t}_{\alpha\beta}\xi_{\alpha\beta}, \end{equation}

where

(2.14)\begin{equation} \boldsymbol{t}_{\alpha\beta} = \frac{(\boldsymbol{I} - \boldsymbol{n}_{\alpha\beta}\boldsymbol{n}_{\alpha\beta})\boldsymbol{\cdot}\tilde{\boldsymbol{\xi}}_{\alpha\beta}} {|(\boldsymbol{I} - \boldsymbol{n}_{\alpha\beta}\boldsymbol{n}_{\alpha\beta})\boldsymbol{\cdot}\tilde{\boldsymbol{\xi}}_{\alpha\beta}|}. \end{equation}

When the contact ends, the spring vector $\boldsymbol {\xi }_{\alpha \beta }$ is reset to zero.

Following the literature (Shäfer, Dippel & Wolf 1996; Silbert et al. 2001; Gallier et al. 2014), the tangential stiffness is set to

(2.15)\begin{equation} k_t=mk_n\sqrt{2\delta}, \end{equation}

where the prefactor $m$ is $2/7$ unless stated otherwise. Since $k_t$ is an artificial parameter that determines a threshold for switching between the static and dynamic regimes, its specific value should not affect the simulation result. For a sufficiently large tangential stiffness $k_t$, (2.12) ensures that the system is in a dynamic regime when the particles are sliding past each other and is in a static regime when their relative surface velocity $\Delta \boldsymbol {U}^{(t)}_{\alpha \beta }$ is close to zero. This is illustrated by an example of two spheres colliding in § 3.2. Also, we verify in §§ 3.2 and S1 of the supplementary material that $k_t$ has a relatively small effect on diffusivity for $m\geqslant 2/7$.

2.3. Numerical implementation

The particle velocities were obtained from their instantaneous positions and contact forces by solving the force balance equations (2.2) using a conjugate gradient solver (Press et al. 1992) that takes advantage of sparsity of the resistance matrix. The particle positions were updated using the fourth-order Adams–Bashforth method (Hairer, Nørsett & Wanner 1993; Butcher 2016). The step for the numerical integration was varied between $2\times 10^{-4}$ for particles with roughness ${\epsilon < 2\times 10^{-2}}$ and $5\times 10^{-4}$ for particles with roughness $2\times 10^{-2}$. Selection of the timestep was confirmed by analysis of convergence of the measured diffusivity as a function of the step size.

All reported simulation results were performed in a periodic cubic cell with sides of length $20$ and Lees–Edwards boundary conditions (Lees & Edwards 1972) were employed for shearing of the periodic cell. Simulations where the box height was doubled produced the same diffusivity within the error. Initially, the particles were placed in the simulation box sequentially with their positions sampled from the uniform random distribution. When a particle overlapped the others, the particle was reassigned to a new random position until no overlap existed. For simulations at the largest considered concentration, $\phi = 0.45$, the particle radius was reduced to 80 % of its original size to accelerate the generation of initial configurations. Then the particles were restored to their original size and the system was sheared forwards until all overlaps disappeared.

After obtaining a non-overlapping configuration, the system was sheared until it reached steady state. Approach to steady state was determined by monitoring the number of contacts in the system. Once the number of contacts per particle approached a constant, the system was considered to be steady. Most systems reached steady state by the strain of 20, whereas rough (roughness $\epsilon = 2\times 10^{-2}$) frictionless particles at the volume fraction $\phi = 0.45$ required a strain of 80 to reach steady state. Data collection for evaluating the reported microstructures and diffusivities began only after reaching steady state.

Each reported diffusivity represents the average across all particles for eight different simulations. Variations in the computed diffusivities for each set of conditions arises due to the assignment of the initial particle positions, hence the standard error was computed from the diffusivity values of each of the separate eight runs. The errors are approximately the size of the symbols representing the diffusivities, save for the results for small values of $k_t$ (see figure S3 of the supplementary material).

3. Results: pairwise collisions

3.1. Collisions of frictionless particles

Representative trajectories of pairwise collisions of frictionless particles are shown in figure 1(a). In order to visualize the trajectories near their apex, figure 1(b) shows $h/2\epsilon$ vs $\theta$, where $h$ is the distance between particle surfaces, $\epsilon$ is the particle roughness and $\theta$ is the angle between the vector connecting the particle centres and the negative direction of the $x$-axis, as shown in figure 1(a).

Figure 1. Pairwise collisions of frictionless particles: (a) trajectories of smooth ($\epsilon = 2\times 10^{-3}$) and rough ($\epsilon = 2\times 10^{-2}$) particles are shown by dashed blue and solid red lines, respectively; (b) dependence of the normalized gap $h/2\epsilon$ between the particle surfaces on the angle $\theta$ between the vector connecting the particle centres and the negative direction of the $x$-axis for trajectories experiencing a contact. Trajectories corresponding to different values of $k_n$ (determined by $\delta$, see (2.7)) are shown. The border between reversible (no contact) and irreversible (contact) segments of the trajectories is shown by the dashed magenta line at $h/2 \epsilon = 1$.

If particles do not come into contact during a collision (i.e. $h$ is always greater than $2\epsilon$, as for trajectory 1 in figure 1a), their trajectories are reversible and result in a zero offset in the transverse direction. If particles come into the contact (e.g. trajectories 2, 3 and 4 in figure 1), their motion becomes irreversible, which results in a non-zero offset in the transverse direction. The offset depends on the particle roughness (with rougher particles having a larger offset), but not on the initial conditions, as long as a contact between particles takes place, as illustrated by trajectories 2 and 3 in figure 1(a). After the particles contact, they follow the same trajectory determined by the particle roughness.

Moreover, during contact, the trajectory normalized by the particle roughness is nearly independent of particle roughness, as evident from figure 1(b). Of course, particles of different roughness come into contact ($h \leqslant 2\epsilon$) at different $\theta$. However, they escape the contact region of $h \leqslant 2\epsilon$ at approximately the same $\theta$, $\theta \approx 90^\circ.$ The offset is then determined by a reversible trajectory with the initial conditions at $h = 2\epsilon$, $\theta \approx 90^\circ$.

The trajectory during contact depends on the value of the elastic force constant $k_n$, as shown in figure 1(b). Note that the closest approach $h_{min}$ of two particles is closely related to the parameter $\delta$ used in the definition (2.7) of $k_n$, $h_{min} \approx 2(\epsilon - \delta )$. Despite the fact that the particle trajectories for different $k_n$ differ during their contact, the value of the offset at the end of a collision is nearly independent of $k_n$, since it is determined by the point at which the contact is lost, i.e. $h = 2\epsilon$, $\theta \approx 90^\circ$, which is very similar in all considered cases. Our calculations show that, for $\delta$ varying between 0.01 and 0.15, the value of the offset in pairwise collisions changes by less than 0.5 %. Consequently, the effect of the normal elasticity $k_n$ on diffusivity in dense suspensions is also relatively small, as shown in figure S2 of the supplementary material.

3.2. Collisions of particles with friction

Friction substantially reduces the relative tangential velocity $\Delta \boldsymbol {U}^{(t)}_{\alpha \beta }$ of particle surfaces during contact and, for a substantial segment of the trajectory during contact, this velocity is nearly zero, as shown in figure 2(a). In the absence of friction, particle surfaces slide past each other, whereas friction leads to particles rolling on each other. However, friction has a relatively small effect on the translational motion of the centre of mass of the particles during a pairwise collision: the trajectories are nearly identical (see figure S4 of the supplementary material) and reduction of the translational velocity by friction is typically less than $7\,\%$ (see figure 2b). Consequently, friction has a negligible effect on the particle offset due to a collision and a marginal effect on the time of the collision. In all considered cases, the duration of particle contact is extended by less than $5\,\%$ in the presence of friction.

Figure 2. Pairwise collisions of smooth particles ($\epsilon = 2\times 10^{-3}$) with and without friction: (a) relative tangential surface velocities $|\Delta \boldsymbol {U}^{(t)}_{\alpha \beta }|$, (b) relative translational velocities, (c) relative rotational velocities $\varOmega _r= \varOmega _{\alpha,z}+\varOmega _{\beta,z}-2\varOmega _{z}^\infty$ and (d) normalized tangential spring length, $k_t\xi _{\alpha \beta }$, and the dynamic friction force, $\mu F ^{(C,n)}_{\alpha \beta }$. The contact takes place between points $A$ and $D$. The friction coefficient of the frictional particles is $\mu = 0.5$. Panel (d) shows $k_t\xi _{\alpha \beta }$ for several values of $k_t$ ($k_t$ is determined by the prefactor $m$, see (2.15)).

Reduction of the relative surface velocity $\Delta \boldsymbol {U}^{(t)}_{\alpha \beta }$ by friction is accomplished largely by changing the rotational velocity of the particles. Friction substantially increases the magnitude of the relative rotational velocity, as seen in figure 2(c). The typical increase of the rotational velocity for the considered particles is approximately $30\,\%$. As will be shown in § 4.2, this increase in rotational velocity plays an important role in the dynamics of dense suspensions.

Figure 2(d) shows forces associated with friction during a collision, illustrating the Cundall–Strack method for modelling friction for different values of the tangential spring stiffness $k_t$. The particles are in a static friction regime immediately after contact, even though their relative surface velocity $\Delta \boldsymbol {U}^{(t)}_{\alpha \beta }$ is large. However, this static regime is very short lived: the tangential spring $\boldsymbol {\xi }_{\alpha \beta }$ quickly becomes large enough for the system to transition into a dynamic friction regime (point $A$ in figure 2).

While the particles are in the dynamic friction regime, the tangential spring length is reset to $(\mu /k_t)F^{(C,n)}_{\alpha \beta }$ at every step (see (2.12)), so that the spring forgets the history associated with the relative surface velocity $\Delta \boldsymbol {U}^{(t)}_{\alpha \beta }$. This continues as long as the magnitude of $\Delta \boldsymbol {U}^{(t)}_{\alpha \beta }$ is large enough so that the increment $\Delta \boldsymbol {U}^{(t)}_{\alpha \beta }\Delta \gamma$ of the spring length exceeds the change in the normalized dynamic friction force, $(\mu /k_t)F^{(C,n)}_{\alpha \beta }$, during a numerical step.

The friction slows down the relative motion of particles and the particle velocity eventually approaches zero, leading (2.12) to switch the friction force to the static regime at point B in figure 2. This transition takes place at approximately the same time for a wide range of the spring force parameter $k_t$. Note, however, that if $k_t$ is too small ($m \leqslant 0.01$, see (2.15)), the relative tangential velocity substantially deviates from zero even in the static regime. In general, the tangential velocity, albeit small, does not become exactly zero in the static friction regime and the magnitude of $\Delta \boldsymbol {U}^{(t)}_{\alpha \beta }$ in this regime is inversely proportional to $k_t$.

The friction force remains static until the normal force $F^{(C,n)}_{\alpha \beta }$ decreases to the value stored in the spring length (i.e. until $\mu F^{(C,n)}_{\alpha \beta } \leqslant k_t\tilde {\xi }_{\alpha \beta }$), at which point the friction transitions to the dynamic regime (point $C$ in figure 2). While in the static regime, the spring elongates according to (2.11). This introduces a correction to the value $k_t\tilde {\xi }_{ab}/\mu$ to which the normal force needs to decrease for the transition from the static to the dynamic regime. This correction ($k_t\xi _{\alpha \beta }$) is relatively independent of the spring constant $k_t$, since $\Delta \boldsymbol {U}^{(t)}_{\alpha \beta }$, and hence $\xi _{\alpha \beta }$, is inversely proportional to $k_t$. Therefore, the specific value of $k_t$ has a relatively small effect on pairwise collisions. As $m$ ranges from $10^{-2}$ to 2, the duration of particle contact is reduced by less than $1\,\%$ and the offset remains unchanged. It is shown in figure S3 of the supplementary material that the effect of the $k_t$ value on particle diffusivity in a dense suspension is also small, as long as $k_t$ is sufficiently large ($m\geqslant 2/7$).

4. Results: many-body systems

Representative mean-squared displacements $\langle (\Delta y)^2 \rangle = \langle (y(\gamma )-y(0))^2 \rangle$ in the transverse direction of particles in suspensions are shown in figure 3(a); note that $\gamma =0$ corresponds to the strain at which the suspension has already achieved steady state and the angle brackets ($\langle {\cdot } \rangle$) indicate an average over all particles and simulations. After an initial transient, the mean-squared displacement scales linearly with strain, thus confirming that the particles undergo a diffusive motion. Their diffusivity is obtained by fitting the mean-squared displacement to a linear relationship

(4.1)\begin{equation} \langle (\Delta y(\gamma))^2 \rangle - \langle (\Delta y(\gamma_0) )^2 \rangle = 2 D_{yy}\left( \gamma-\gamma_0 \right), \end{equation}

over strains for which $\langle (\Delta y)^2 \rangle$ scales linearly with strain. The additional strain ($\gamma _0$) required to reach the diffusive regime after reaching steady state ($\gamma =0$) varied between 30 and 70, depending on the system. Further details on diffusivity calculations are provided in § S3 of the supplementary material.

Figure 3. (a) Representative mean-squared displacements for suspensions with volume fraction $\phi = 0.45$; (b) diffusion coefficient as a function of volume fraction for frictionless ($\mu =0$) and frictional ($\mu =0.5$) particles for two values of roughness $\epsilon$.

Figure 3(b) shows the dependence of the diffusion coefficient $D_{yy}$ in the gradient direction on suspension concentration $\phi$ for particles with different roughness and friction coefficient. In the absence of friction, the diffusivity exhibits a non-monotonic dependence on $\phi$, with $D_{yy}$ growing at small $\phi$ and decreasing at large $\phi$. This result is in qualitative agreement with the simulations of Zhang et al. (2023), which utilized less detailed models for lubrication and contact forces. In contrast, frictional particles exhibit a monotonic increase of diffusivity with $\phi$.

The remainder of the results focus on suspensions with $\phi = 0.45$, where the differences between frictional and frictionless diffusivities are largest. Also note that the diffusivities predicted for the range of volume fractions studied here (0.25 to 0.45) closely match those predicted by Stokesian dynamics simulations (Sierou & Brady 2004) that include long-range hydrodynamic interactions, at least in the absence of friction (see figure S5 of the supplementary material). However, long-range hydrodynamic interactions may have a greater effect when friction is considered, although at higher concentrations such as $\phi = 0.45$, its influence is expected to be minimal, further justifying the focus on this volume fraction.

4.1. Effect of particle roughness

Figure 3(b) shows that, at sufficiently high concentrations, the diffusivity of rough, frictionless particles is lower than that of the smooth particles. Using simulations that neglected tangential and rotational components of the lubrication force, Zhang et al. (2023) previously demonstrated that the reduction in diffusivity for rough particles is due to the layering of rough particles. Here, results indicate that the same layering mechanism holds for the more accurate lubrication and normal contact models considered in the current simulations. Typical snapshots of suspensions of frictionless rough and smooth particles are shown in figure 4. Layering is clearly visible for the rough particles in figure 4(a).

Figure 4. Snapshots of steady-state system configurations of frictionless particles at the volume fraction ${\phi = 0.45}$ with (a) $\epsilon = 2\times 10^{-2}$ and (b) $\epsilon = 2\times 10^{-3}$. The plots show projections of the particle positions onto the $x\unicode{x2013}y$ plane. For clarity, particles are represented by circles with a radius smaller than the actual particle radius.

Layering is also evident from the pair-correlation function averaged over the azimuthal angle

(4.2)\begin{equation} \bar{g}(r,\theta) = \frac{\displaystyle\int_0^{2{\rm \pi}} g(r,\theta,\psi)\,{\rm d}\psi} {2{\rm \pi} r^2 \sin\theta}, \end{equation}

where $r$ is the length and $\theta$ and $\psi$ are the polar and azimuthal angles of a vector connecting particle centres, $g(r,\theta,\psi )$ is the pair-correlation function and the denominator in (4.2) corresponds to the volume over which $g(r,\theta,\psi )$ is averaged. The averaged pair-correlation function $\bar {g}(r,\theta )$ is shown in figure 5. Note the strong layering for rough particles without friction (see figure 5a). The pair-correlation function suggests that particles also exhibit some layering at other conditions, but it is much weaker. In particular, adding friction to the rough particle model substantially reduces layering (see figure 5c), which is consistent with the larger diffusivity of those particles than particles without friction, see figure 3(b).

Figure 5. Logarithm of the averaged pair-correlation functions $\bar {g}(r,\theta )$ of particles at the volume fraction ${\phi = 0.45}$ with (a) $\epsilon = 2\times 10^{-2}$, $\mu = 0$, (b) $\epsilon = 2\times 10^{-3}$, $\mu = 0$, (c) $\epsilon = 2\times 10^{-2}$, $\mu = 0.25$ and (d) $\epsilon = 2\times 10^{-3}$, $\mu = 0.25$. The function $\bar {g}$ is normalized so that $\bar {g}(r,\theta ) \to 1$ as $r\to \infty$.

To better quantify the layering, we computed the suspension structure factor

(4.3)\begin{equation} S(\boldsymbol{q}) = \frac{1}{N} \sum_{\alpha,\beta = 1}^N\langle \exp({{\rm i}\boldsymbol{q}\boldsymbol{\cdot}(\boldsymbol{r}_\alpha - \boldsymbol{r}_\beta)}) \rangle, \end{equation}

where $N$ is the total number of particles in the system and $\boldsymbol {q} = (q_x, q_y, q_z)$ is the wavevector. Up to an additive constant and a multiplicative factor, the structure factor corresponds to the Fourier transform of the pair-correlation function $g(\boldsymbol {r})$. To focus on layering in the gradient ($y$) direction, we consider the structure factor corresponding to the wavevector $\boldsymbol {q} = (0, q_y, 0)$. These structure factors for several systems are shown in figure 6(a). Clearly, there is a large peak at the wavelength of $\lambda \approx 2$ for the rough frictionless particles, which indicates layering with spacing between layers corresponding to the particle diameter. The peak corresponding to layering is observed in all systems plotted in figure 6(a), but its magnitude is at least an order of magnitude higher for rough frictionless particles than for other considered particles, confirming that the rough frictionless particles are much more ordered. Additionally, the strong layering for the rough frictionless particles remains the same even when the box size is doubled in the gradient direction, indicating that the layering does not depend on box size when it exceeds 20.

Figure 6. Microstructure characteristics of suspensions at the volume fraction $\phi = 0.45$: (a) the structure factors $S$ as a function of the wavenumber $\lambda =2{\rm \pi} /q_y$ in the gradient direction and (b) mean number of contacts $N_c$ per particle per timestep.

In earlier work, Zhang et al. (2023) assumed that layering lowers diffusivity by reducing the number of particle contacts. However, a direct calculation shows that the layered structures of frictionless rough particles have more contacts per particle than the relatively disordered structures of frictionless smooth particles, see figure 6(b).

The smaller diffusivity of rough particles is explained by the fact that the particles in layered structures tend to stay in their layer. Typical trajectories are shown in figure 7. It is evident that, for rough frictionless particles, most trajectories stay within the range $-1 \leqslant y-y(0) \leqslant 1$, where $\gamma = 0$ is the strain at which the system reached steady state. Also, when a particle moves away from its layer, it is often reflected back by another layer. This is also evident from the oscillations of the velocity autocorrelation function $C_{yy}(\tau ) = \langle U_y(t)U_y(t+\tau )\rangle$ for these particles shown in figure 8, which suggest collisions in which particles are reflected back rather than translating into another layer. Hence, transitions between layers of frictionless rough particles are rare events. In contrast, smooth frictionless particles, as well as particles with friction, exhibit a relatively uniform (in time) diffusion (see figure 7) and their velocity autocorrelation functions do not exhibit oscillations (see figure 8).

Figure 7. Representative trajectories of particles in a suspension with the volume fraction $\phi = 0.45$ with (a) $\epsilon = 2\times 10^{-2}$, $\mu = 0$, (b) $\epsilon = 2\times 10^{-3}$, $\mu = 0$, (c) $\epsilon = 2\times 10^{-2}$, $\mu = 0.25$, (d) $\epsilon = 2\times 10^{-3}$, $\mu = 0.25$. Here, $\gamma = 0$ is the strain at which the system reached steady state. The dashed lines in panel (a) indicate the locations of layers.

Figure 8. Representative normalized velocity autocorrelation functions $C_{yy}$ for suspensions at the volume fraction $\phi = 0.45$.

4.2. Effect of friction

Figure 3(b) shows that increasing the friction coefficient $\mu$ from 0 to 0.5 substantially increases the shear-induced diffusivity of both smooth and rough particles. Effects of further increasing the friction coefficient are shown in figure 9(a), which indicates that the diffusivity $D_{yy}$ continues to increase as $\mu$ increases. The rate of increase of $D_{yy}$ with $\mu$ is larger for rougher particles, and the diffusivity of rough particles surpasses that of smooth particles for $\mu \geqslant 0.25$.

Figure 9. (a) Effect of the friction coefficient $\mu$ on the diffusion coefficient $D_{yy}$. (b) Comparison of diffusivity of particles with roughness $\epsilon = 3.8\times 10^{-3}$ obtained by the simulations with the experimental result of Metzger et al. (2013b). At $\mu = 0.18$, results of simulations for several values of roughness are shown. Data for suspensions with volume fraction $\phi = 0.45$ are shown.

Figure 9(b) indicates that the addition of friction to the model improves the agreement with the experiments of Metzger et al. (2013b). The experiments were performed for a suspension of particles similar to that of Zhang et al. (2023) (roughness $\epsilon =3.8\times 10^{-3}$) at a volume fraction of $\phi = 0.45$. Simulations of frictionless particles at these conditions predict a diffusivity lower than that measured experimentally. However, since friction increases the diffusivity, at a sufficiently large friction coefficient ($\mu \geqslant 0.25$), the predictions for diffusivity exceed the experimental result. Linear interpolation of the simulation results suggests that the simulations are in agreement with experiment at $\mu \approx 0.2$. Also note that the diffusivity values show relatively little variation with changes in roughness for this value of the friction coefficient, as shown in figure 9(b).

The increase of diffusivity with increasing particle friction is somewhat counterintuitive, as one typically associates increased friction with reduced mobility and slower transport processes. This phenomenon can be partly attributed to decreased particle layering as $\mu$ increases. Layering in suspensions is quantified by the structure factor's peak ($S_{max}$) at $\lambda \approx 2$ (see figure 6a). Figure 10 shows $S_{max}$ decreasing with rising $\mu$, indicating reduced suspension layering. For rough particles, $S_{max}$ decreases by more than an order of magnitude as $\mu$ increases from 0 to 0.25. This disruption of layers in the suspension induced by friction facilitates diffusivity, as confirmed by sample trajectories of rough particles with friction shown in figure 7(c). Unlike the rough frictionless particles (shown in figure 7a), the rough frictional particles do not remain in the same layer for large strains before jumping to another layer. Further evidence of the changed dynamic of rough particles is provided by the velocity autocorrelation functions (see figure 8), which do not exhibit oscillations at $\mu = 0.25$, in contrast with the frictionless rough particles.

Figure 10. Effect of the friction coefficient $\mu$ on $S_{max}$, the peak value of the structure factor $S$ at $\lambda \approx 2$, which indicates the extent of particle layering. Data for suspensions with volume fraction $\phi =0.45$ are shown.

The disruption of layering by friction within suspensions of monodisperse particles also leads to higher viscosities (Goyal et al. 2022) as well as diffusivities. However, friction-induced disorder provides only a partial explanation for the diffusivity increase with friction, since the magnitude $S_{\max }$ of the structure factor peak reaches an asymptote at large $\mu$ (see figure 10), whereas the diffusivity increase continues (see figure 9a). Moreover, smooth particles are disordered even at $\mu = 0$.

Further insights into the mechanism of increasing diffusivity with friction are provided by an examination of the rotational velocity fluctuations of particles. The magnitude of fluctuations of the rotational velocity $\boldsymbol {\varOmega }$ of frictional particles is larger than that of frictionless particles (see figure 11a), and the variance $\sigma ^2_\varOmega$ of $\boldsymbol {\varOmega }$ grows with friction (see figure 11b).

Figure 11. (a) Probability distribution of the relative rotational velocity $\varOmega _z - \varOmega _z^\infty$ in the $z$-direction for frictional ($\mu =0.25$) and frictionless ($\mu =0.0$) particles with two values of roughness, $\epsilon$. (b) Dependence of the variance $\sigma _\varOmega ^2$ of the rotational velocity $\boldsymbol {\varOmega }$ on the friction coefficient $\mu$ for smooth and rough particles. Data for suspensions with volume fraction $\phi =0.45$ are shown.

The increase of $\sigma ^2_\varOmega$ with $\mu$ is caused by reduction of the relative tangential velocity $\Delta \boldsymbol {U}^{(t)}_{\alpha \beta }$ by friction. During pairwise collisions, friction reduces the relative tangential velocity $\Delta \boldsymbol {U}^{(t)}_{\alpha \beta }$ of particles, which is accomplished by increasing the magnitude of their relative rotational velocity (see § 3.2 and figure 2). A similar reduction in $\Delta \boldsymbol {U}^{(t)}_{\alpha \beta }$ upon contact is observed for frictional particles in dense suspensions. Figure 12 illustrates this phenomenon by showing the dependence of the mean magnitude of $\Delta \boldsymbol {U}^{(t)}_{\alpha \beta }$ on relative particle positions during and near contact. Achieving this reduction requires increasing the relative rotational velocities of particles, similarly to the pairwise collisions, which in turn leads to larger fluctuations of the rotational velocity (see figure 11).

Figure 12. Mean relative tangential velocity $\Delta \boldsymbol {U}^{(t)}_{\alpha \beta }$ near contact at the volume fraction $\phi = 0.45$ with(a) $\epsilon = 2\times 10^{-2}$, $\mu = 0$, (b) $\epsilon = 2\times 10^{-3}$, $\mu = 0$, (c) $\epsilon = 2\times 10^{-2}$, $\mu = 0.25$ and (d) $\epsilon = 2\times 10^{-3}$, $\mu = 0.25$.

The larger rotational velocity fluctuations are then transmitted to the translational velocity via perturbations of the fluid and the network of contacts formed by the particles. As a result, the variance $\sigma ^2_U$ of the translational velocity exhibits trends similar to those of $\sigma ^2_\varOmega$, as evident from figure 13. The growth of $\sigma ^2_U$ with friction gives a corresponding increase of diffusivity, as well as the destruction of layers by friction.

Figure 13. Effect of the friction coefficient $\mu$ on the variance $\sigma ^2_U$ of the transverse velocity of particles in suspensions with volume fraction $\phi =0.45$.

The rate of increase of $\sigma ^2_\varOmega$ with $\mu$ is larger for rough particles, which is consistent with the larger rate of increase of the diffusivity of these particles. A likely reason for the larger growth of $\sigma ^2_\varOmega$ with $\mu$ of the rough particles is the larger number of contacts per particle (see figure 6b), which leads to larger fluctuations of the rotational velocity.

5. Discussion and conclusions

Simulations were conducted to evaluate the effects of particle roughness and friction on the shear-induced diffusivity of suspensions. The results predict that concentrated suspensions of rough particles exhibit lower diffusivity compared with smooth particles when friction is either not considered or is sufficiently low. These simulations confirm the experimental and computational findings of Zhang et al. (2023), while incorporating additional components of hydrodynamic interactions and an improved contact model. The observed reduction in diffusivity is attributable to enhanced flow-aligned layering of rough particles. This well-organized structure retains particles on the same streamline for large strains, with infrequent transitions between layers.

At high concentrations, introducing friction significantly disrupts the layering of rough particles and increases their diffusivity, irrespective of roughness. Friction has a negligible impact on pairwise collisions, indicating that the disruption of the microstructure by friction is an inherently many-body effect. When examining the dynamics, friction is also observed to greatly increase the variation of the rotational velocities of the particles which, in turn, causes an increase in the centre of mass fluctuations and diffusivity, as argued in § 4.2.

Notably, the addition of friction enhances the alignment of simulated diffusivities with experimental measurements, as originally determined by Gallier (2014). In the current simulations, however, the friction coefficient and the roughness (the range $\epsilon$ over which the Hertzian contact acts) have been independently varied. Ideally, understanding the relationship between them would enable more detailed comparisons with experimental data. Tribology literature (Spikes 1997) suggests that, in the mixed lubrication regime, where both lubrication and asperity contact occur simultaneously, increasing roughness generally results in a higher coefficient of friction. However, friction also depends on other factors, such as the texture or pattern of the roughness, which can change that trend (Spikes & Olver 2003).

If rougher particles indeed experience higher frictional forces, the results presented here suggest a potential cancellation of effects, which may limit changes in diffusion due to varying particle roughness. At high concentrations, roughening the particles reduces diffusivity when holding the friction at a low value. However, roughening probably increases the friction coefficient, which increases the diffusivity. Hence, experimental observations made by Zhang et al. (2023) of a lower diffusivity for roughened particles suggests that the friction coefficient (of the smooth and rough particles) was very low.

Roughness can significantly alter tangential lubrication forces even in the absence of solid–solid contact, effectively changing the frictional behaviour (Spikes 1997). Recent studies focused on suspensions (Jamali & Brady 2019; Yariv et al. 2024) have also suggested that lubrication interactions between asperities can mimic the effects of friction. The strength of tangential lubrication interactions between smooth surfaces scales as $\log (1/h)$, offering little resistance to relative motion. However, the approach of asperities during tangential motion provides a large, additional resistance. This concept has been used to predict discontinuous shear thickening (Wang, Jamali & Brady 2020) and might also be a useful framework for predicting increased diffusivities, similar to the effects of friction explored in this study.

In this paper, all results were computed using periodic boundary conditions and spheres of uniform size. Deviating from these conditions can alter the predicted diffusivities. Confining the suspension between bounding walls impacts the microstructure and diffusivity of the suspensions. For example, Zhang et al. (2023) found that confining roughened particles between walls separated by less than 10 diameters results in increased layering and lower diffusivity. Regardless of the roughness, the simulations predicted that the diffusivity remained constant when the wall separation exceeded 10 diameters. When friction is considered, bounding walls would likely further influence the structure and diffusivity predictions. This is an important consideration for comparing with experimental measurements, as these are often made within confined systems. Real systems also typically contain some degree of polydispersity in particle size. Zhang et al. (2023) found that polydispersity reduces the differences in diffusivity between frictionless rough and smooth particles. Specifically, increasing the polydispersity was observed to increase the diffusivity of rougher particles. However, the effect of polydispersity on diffusivity in suspensions with friction remains unknown and warrants further investigation.