Shear thickening in dense suspensions driven by particle interlocking

Abstract

The processing of dense suspensions is a crucial step in many industries including mining; the production of ceramics; and the manufacture of pharmaceuticals. It is widely reported that these suspensions exhibit nonlinear behaviours such as shear thinning and thickening, with particle surface contacts recently being accepted as a primary culprit in the latter. In light of this, the modelling community have started to explore the role of particle surface tribology, predominantly by incorporating Coulombic friction models borrowed from the field of dry granular matter. Full details of the interactions between particle surfaces remain unclear, however, and it is suggested that physical interlocking of particle asperities may be key. Here, we use particle-based simulations to explore explicitly the effect of interlocking on the rheology of dense suspensions of micron-sized solids in a Newtonian fluid. Our simplified model recovers shear thinning, thickening and jamming phenomena commonly seen in experiments.

1. Introduction

Dense suspensions of solid particles in a Newtonian fluid are found throughout nature, in everything from geophysical phenomena such as landslides and magma flow to processes like the transport of sediment in rivers (Denn, Morris & Bonn 2018; Guazzelli & Pouliquen 2018; Hsiao & Pradeep 2019; More & Ardekani 2020c; Ness, Seto & Mari 2022). Meanwhile, they are present in many industries such as the manufacture of pastes and ceramics; the production of vaccines and drug formulations; and the preparation of metal pastes for modern fuel cells (Denn et al. 2018; Guazzelli & Pouliquen 2018; More & Ardekani 2020c). Nevertheless, there is no unified model to describe their behaviour under flow, which typically departs substantially from classical laws and instead shows capricious rate dependence (Barnes 1989; Bender & Wagner 1996; Brown & Jaeger 2012; Seto et al. 2013; Mari et al. 2014; Ness & Sun 2015). Hence, the study of dense suspensions remains a diverse sphere of research where experimental and computational studies must be used in parallel to unveil the microscopic mechanisms underlying their bulk rheology and fluid mechanics (Denn et al. 2018; Duran 2012; More & Ardekani 2020a,c).

We focus here on suspensions of non-Brownian, non-attractive particles, which often shear thicken at large stresses. The new consensus (Denn et al. 2018; Ness et al. 2022) is that frictional contacts are responsible for this behaviour, although the mechanism by which these contacts occur is still debated (Jamali & Brady 2019). Experimental measurements of particle–particle-contact physics (Comtet et al. 2017; Hsu et al. 2018) reveal that static friction is present, although measured sliding coefficients are incongruous with the bulk rheology suggesting that other physics (rolling, twisting, adhesion, contact deformation leading to shear thinning) may also be at play (Lobry et al. 2019; Arshad et al. 2021).

Although numerical models which employ fictitious friction coefficients have led to substantial progress in establishing the constitutive behaviour of these suspensions (Mari et al. 2014; Ness, Xing & Eiser 2017b; Gillissen & Ness 2020; Gillissen et al. 2020) and have shed light on experimental findings (Lin et al. 2015; Guy et al. 2020; Singh et al. 2020), they necessarily miss more detailed aspects of the surface physics involved – especially the effect of physical interlocking between asperities. Other numerical works have made substantial progress in exposing the role of roughness height (More & Ardekani 2020a,b), but lack explicit asperity interlocking.

In particular, this picture and these models raise fundamental questions about the tribology of particle contacts: (i) What constitutes a contact between suspended particles? (ii) How do lubrication layers break down and are stresses sustained by fluid or solid? (iii) How do we get from surface chemistry to mechanical friction? (iv) How can we apply a ‘coarse-grained’ macroscopic friction coefficient to contacts that may comprise just a few interlocking asperities? Such questions present a new challenge to the field, and their answers will eventually inform the next generation of formulation and synthesis techniques for suspensions with engineered rheology.

Here, we present the first step in adapting our computational approach to this new paradigm. To do so, we omit the Coulombic friction coefficient typically used in dense suspension models (Seto et al. 2013; Mari et al. 2014) in favour of a more detailed particle surface representation comprising rigid (but frictionless) asperities and a short-range repulsive force. We then subject particle assemblies to simple shear flow and measure their bulk rheology, viz. the reduced suspension viscosity $\eta _r$ as a function of the solids fraction $\phi$ and shear rate $\dot {\gamma }$.

2. Simulation methodology

Our model is designed to invoke non-Newtonian suspension rheology (especially shear thickening) without the need for critical load models (Seto et al. 2013) nor imposed friction coefficients. Instead, we simulate gear-like aggregate particles, which are almost a two-dimensional (2-D) analogue to those developed experimentally by Hsu et al. (2018); figure 1.

Figure 1. Schematic of aggregate particles (coloured by particle). Shown are suspensions with asperity sizes $k_r=$ (a) $0.05$; (b) $0.10$; (c) $0.20$. Green arrows indicate the imposed velocity profile $u_x(y) = \dot {\gamma }y$. Blue dashed circle represents contact with $P=0$ (see (3.1)).

We follow an established dense suspension simulation technique (Ness 2021) consisting of a discrete element method (Cundall & Strack 1979) complemented by short-range hydrodynamics (Ball & Melrose 1997). The basic approach is to: (i) initialise a system of non-overlapping particles at a desired solids fraction $\phi$; and (ii) evaluate the trajectory of each particle over a series of timesteps by numerically solving Newton's second law under a prescribed background fluid velocity gradient $\boldsymbol {\nabla }\boldsymbol {u}^\infty$ and a set of pairwise interactions. When needed, the bulk stress tensor $\mathbb {\Sigma }$ is calculated from the interaction forces and particle positions, thus generating rheology data comprising the stress $\mathbb {\Sigma }$ as a function of deformation rate $\mathbb {E}$ (with $\mathbb {E}\equiv \frac {1}{2}(\boldsymbol {\nabla }\boldsymbol {u}^\infty + {\boldsymbol {\nabla }\boldsymbol {u}^\infty }^\mathrm {T}$)) and $\phi$.

We consider a dense ($\phi$ close to the jamming fraction $\phi _m$), non-Brownian ($k_bT$ omitted) suspension of aggregate particles (described below) under a uniform flow with an imposed deformation rate. The particle properties that set the length, mass and time scales are the characteristic radius $a$ (length), the density $\varrho$ (mass/length$^3$) (taken to be equal to the fluid density) and the normal stiffness $k_n$ (mass/time$^2$) (with tangential counterpart $k_t$). We also define a fluid viscosity $\eta _{f}$ $({\rm mass}/({\rm length}\times {\rm time}))$ and a particle–particle friction coefficient $\mu$ (dimensionless) (set to 0 for all aggregate suspensions but used to simulate smooth frictional particles for comparison). The background flow is characterised by a velocity field $\boldsymbol {u}^\infty$ (length/time) and its gradient (a tensor, taken to be spatially uniform) $\boldsymbol {\nabla }{\boldsymbol {u}}^\infty$ (1/time), the time $t$ for which it is applied, and a stress tensor $\mathbb {\Sigma }$ (${\rm mass}/({\rm time}^2\times {\rm length})$). Below we address a scalar velocity gradient as $\dot {\gamma }$ (${\equiv }\partial u_x/\partial y$) and a scalar stress as $\varSigma _{xy}$ (the $xy$ component of $\mathbb {\Sigma }$). The non-dimensional parameters necessary to fully define a simple suspension under given flow conditions are then: $\dot {\gamma } \sqrt {\varrho a^3/k_n}$, $\varrho \dot {\gamma }a^2/\eta _{f}$, $\mu$, $\phi$, $\eta _r \equiv \varSigma _{xy}/\eta _{f}\dot {\gamma }$ and $\dot {\gamma }t$. Setting the first two quantities $\ll 1$ ensures hard particles and inertia-free conditions. The model then produces rate-independent rheology, that is $\eta _r=\eta _r(\phi )$ at large $\dot {\gamma }t$ (Boyer, Guazzelli & Pouliquen 2011).

Particles experience four types of force and torque: drag, pairwise lubrication, pairwise repulsion and pairwise contact. The full form of these is reported by several authors (Trulsson, Andreotti & Claudin 2012; Mari et al. 2014; Cheal & Ness 2018; Ge & Brandt 2020) and we give a simplified description here.

The Stokes drag on particle $i$ (radius $a_i$) is proportional to the difference between its velocity $\boldsymbol {u}_i$ and the fluid streaming velocity at its centre $\boldsymbol {u}^\infty (\boldsymbol {x}_i)$

(2.1)\begin{equation} \boldsymbol{F}^{d}_i ={-}6{\rm \pi}\eta_{f}a_i(\boldsymbol{u}_i - \boldsymbol{u}^\infty(\boldsymbol{x}_i)) {.} \end{equation}

This force induces flow in the simulation, causing particles to conform to the streaming velocity set by $\boldsymbol {u}^\infty$. Similarly, a torque acts to cause the particles to rotate with angular velocity set by $\frac {1}{2}(\boldsymbol {\nabla }\times \boldsymbol {u}^\infty )$. Neighbouring particles $i$ and $j$ with centre-to-centre vector $\boldsymbol {r}_{i,j}$ (with corresponding unit vector $\boldsymbol {n}_{i,j}$) experience lubrication forces (see Kim & Karrila 1991; Jeffrey 1992) dependent on the gap $h$ between them and their relative velocity. The leading term of the force (written below for equal-sized spheres) scales with $1/h$ and the normal part of the velocity difference

(2.2)\begin{equation} \boldsymbol{F}^l_{i,j} = \frac{3}{2}{\rm \pi} a_i^2 \eta_{f}\frac{1}{h}(\boldsymbol{u}_j - \boldsymbol{u}_i)_n {.} \end{equation}

This force opposes relative motion between particle pairs. We also include the $\log (1/h)$ terms and the tangential components as described in Cheal & Ness (2018). A lower limit on $h$ (typically ${O}(10^{-3}a_i)$) prevents divergence at contact. A torque also acts to resist relative rotation between $i$ and $j$ (see Cheal & Ness 2018). Neighbouring particles also experience a repulsive force given by the simple form

(2.3)\begin{equation} \boldsymbol{F}^r_{i,j} = F^{rep}\exp\left(\frac{(a_i+a_j)-|\boldsymbol{r}_{i,j}|}{\rho}\right) \boldsymbol{n}_{i,j} {,} \end{equation}

with $F^{rep}$ the force at contact and $\rho$ setting the rate of decay of the force. An overlapping particle pair $i$ and $j$ experience a pseudo-hard-sphere contact force dependent upon the scalar overlap $\delta$ and the tangential displacement accumulated over the duration of the contact $\boldsymbol {\xi }$

(2.4)\begin{equation} \boldsymbol{F}^c_{i,j} = k_n\delta \boldsymbol{n}_{i,j} - k_t\boldsymbol{\xi} {.} \end{equation}

The friction coefficient $\mu$ sets an upper bound on $\boldsymbol {\xi }$ through $|\boldsymbol {\xi }|\leqslant \mu k_n\delta /k_t$ (when $\mu =0$ there is no tangential contact force).

The stress contribution from drag forces is proportional to $\mathbb {E}$. The $\alpha,\beta$ component of the stress due to lubrication, repulsion and contact is found, respectively, by summing $(F^{l,\alpha }_{i,j}r^\beta _{i,j} + F^{l,\beta }_{i,j}r^\alpha _{i,j})/2$, $F^{r,\alpha }_{i,j}r^\beta _{i,j}$ and $F^{c,\alpha }_{i,j}r^\beta _{i,j}$ over all interacting pairs. The forces are summed on each particle and the trajectories are then updated using a numerical scheme with timestep chosen to be small compared with $\sqrt {\varrho a^3/k_n}$ and $\varrho a^2/\eta _{f}$.

Here we address gear-like aggregate particles comprising central hosts with affixed surface asperities, figure 1. Simplifying the contact of two rough surfaces to an interaction between a number of spherically tipped asperities is a longstanding approach in tribology that follows Greenwood & Williamson (1966). For tractability in a particle-based code aimed at obtaining rheological predictions, we interpret this approach in its simplest form, by fixing spherical asperities to particle surfaces (hemispherical asperities would yield analogous but quantitatively different results since a lower degree of interlocking would occur). This approach offers an explicit description of particle interlocking enabling study of microphysics at a more resolved level compared with those models that apply Coulombic friction.

We consider a 1 : 1 (by number) mixture of host particles with radii $a_h$ and $1.4a_h$ (following conventional practice (Cheal & Ness 2018)). Aggregates are constructed by fixing small asperity particles (radii $k_ra_h$ and $1.4k_ra_h$) to host surfaces at equal intervals, and we explore roughnesses $k_r = (0.05, 0.10, 0.20)$. The quantity $nk_r$ is set to unity, thus specifying the number $n$ of asperities on each host. The drag, lubrication and contact forces (with $\mu =0$) described above are summed over all particle constituents (host plus asperities) of each aggregate at each timestep before the resultant total force is applied uniformly across all particle constituents taking into account the differing masses of host and asperity. Repulsive forces ($\boldsymbol {F}^r$) act on hosts only. Thus, the aggregates behave as rigid bodies (no relative translation or rotation between host and asperities). For comparison, we also simulate smooth particles under the same shearing conditions. These are essentially aggregates with $k_r=0$ and with a variable surface friction coefficient $\mu$ set either to 0 (frictionless) or 100 (frictional). For simplicity we use a 2-D domain, simulating a monolayer of spherical or spherical-aggregate particles. We consider systems of ${O}(10^2)$ aggregate particles, having verified that larger systems produce statistically equivalent results.

We apply simple shear flow with rate $\dot {\gamma }$ (green lines, figure 1) by tilting the periodic simulation box (at fixed volume) and mapping particle motion across the upper and lower $x$ planes using Lees–Edwards conditions. Each simulation is run to $\dot {\gamma t}=4$ and bulk properties are time averaged across the steady state (typically observed for $\dot {\gamma }t\gtrsim 1$), and across 10 independent realisations. The repulsive force magnitude $F^{rep}$ can be used to define a dimensionless shear rate $\dot {\gamma }^*=6{\rm \pi} \eta _{f}a^2\dot {\gamma }/F^{rep}$, quantifying rate dependence in our model. The model is implemented in LAMMPS (Plimpton 1995).

3. Simulation results

The interpretation of the bulk rheology resulting from our model is qualitatively consistent with much of the recent literature (Ness et al. 2022). We give a brief account of the $\phi$, $\dot {\gamma }^*$ and $\rho$ dependence of $\eta _r$ before focussing on the microstructural behaviour.

3.1. The effect of solid loading

To study the role of solid loading we consider $\phi = 0.55\unicode{x2013}0.86$ at intervals of $\Delta \phi =0.01$ (where $\phi$ takes into account the total aggregate areas (i.e. the combined areas of hosts plus asperities)), ensuring that $\phi _m$ is captured in all cases (Kausch, Fesko & Tschoegl 1971). Here, $\phi _m$ represents the value of $\phi$ at the inflexion in $\eta _r(\phi )$ (see below and figure 2). We focus on the particle-contact regime by setting the magnitude of the repulsive force $F^{rep}$ negligibly small (equivalent to $\dot {\gamma }^*\rightarrow \infty$, and thus recovering rate independence). In all cases $\eta _r$ increases according to a power law as $\phi \rightarrow \phi _m$, with aggregate particles having smaller $\phi _m$ than smooth ones, figure 2. For frictionless smooth particles $\phi _m\approx 0.84$, which is in quantitative agreement with Reichhardt & Reichhardt (2014) and O'Hern et al. (2003). Frictional smooth particles and each of the aggregate suspensions exhibit $\phi _m$ well below this value.

Figure 2. Divergence of the relative viscosity $\eta _r$ with area fraction $\phi$ at $\dot {\gamma }^*=\infty$ for each of the suspension types explored. Shown in dashed lines are fits to $\eta _r = \nu (1-\phi /\phi _m)^{-\lambda }$. The maximum solid fractions for smooth particles $\mu =(0,100)$ and aggregates with $k_r=(0.05,0.10,0.20)$ are, respectively, $\phi _m=(0.84,0.77,0.74,0.71,0.69)$ (see table 1).

Generally, jamming occurs at lower $\phi$ as the size of the tangential forces between particles increases in smooth suspensions (i.e. increasing $\mu$); while for aggregates $\phi _m$ decreases systematically as surface roughness $k_r$ is increased. Our chosen values of $k_r$ demonstrate a transition from the $k_r\to 0$ limit (smooth particles with large $\mu$) to saturation (there is not much variation between $k_r=0.10$ and $k_r=0.20$). The behaviour may deviate at larger $k_r$ as aggregates resemble more fractal, floc-like structures that are beyond the scope of this work.

This agrees broadly with the experimental work of Hsu et al. (2018) as well as several works recounted by Hsiao & Pradeep (2019). Similar to Cheal & Ness (2018) the divergence is in the contact contribution to the stress, although interestingly the value of $\phi$ above which contacts dominate is roughly independent of $k_r$ for aggregates despite their differing $\phi _m$. Due to the finite stiffness of our simulated particles we observe an overturn in the viscosity at $\phi \geqslant \phi _m$, ordinarily only observed experimentally for soft particles (e.g. Nordstrom et al. 2010). Focussing on the behaviour below $\phi _m$, a reasonable approximation of the rheology is obtained via the phenomenological expression $\eta _r=\nu (1-\phi /\phi _m)^{-\lambda }$, with values of $\phi _m$ given in the figure 2 caption, and full fitting parameters given in table 1.

Table 1. Fitting parameters for the relative suspension viscosity as a function of the solid fraction, with functional form $\eta _r = \nu (1-\phi /\phi _m)^{-\lambda }$.

3.2. Shear rate dependence

We explored shear rates $\dot {\gamma }^* = {O}(10^{-4}\unicode{x2013}10^3)$ at modest repulsive force decay rates $\rho =3k_ra_h$ ($\rho =0.1(a_i+a_j)$ for smooth particles), and a range of $\phi <\phi _m$, figure 3. In all cases the model predicts a sequence of shear thinning at $\dot {\gamma }^*\lesssim 1$, thickening at $1\lesssim \dot {\gamma }^*\lesssim 10^2$ and a quasi-Newtonian plateau at $\dot {\gamma }^*\gtrsim 10^2$. The sharp transition from shear thinning to dramatic shear thickening observed by Egres & Wagner (2005), Boersma, Laven & Stein (1990) and Hsu et al. (2018) is recovered close to $\phi _m$ in each case. Shear thickening in our aggregate model occurs as the short-ranged repulsive force gives way to asperity–asperity interactions with increasing $\dot {\gamma }^*$. For smooth spheres we observe very modest shear thickening when $\mu =0$ as the transition is from repulsive interactions to frictionless contacts as the stress is increased. For $\mu =100$ we observe strong shear thickening, similar to the aggregate cases.

Figure 3. Rheology of smooth and aggregate particles, showing the relative viscosity $\eta _r$ as a function of shear rate $\dot {\gamma }^*$ (a,c,e,g,i) and stress $\eta _r\dot {\gamma }^*$ (b,d,f,h,j) for several $\phi$ (see legends). Shown are: (a,b) smooth frictionless ($\mu =0$) particles; (c,d) smooth frictional ($\mu =100$) particles; (e,f) aggregates with $k_r=0.05$; (g,h) aggregates with $k_r=0.10$; (i,j) aggregates with $k_r=0.20$. For (a–d) $\rho =0.1(a_i+a_j)$; for (e–j) $\rho =3k_ra_h$. Insets show sketches of smooth and aggregate particles represented in each panel.

When plotting $\eta _r(\dot {\gamma }^*)$ (figure 3), we see that aggregates exhibit $\eta _r$ significantly higher than those of frictionless and frictional smooth particles when comparing equal $\phi$. Meanwhile, as the aggregate $k_r$ is increased $\eta _r$ generally increases, although comparing $k_r=0.10$ with $k_r=0.20$ the effect is rather subtle suggesting, interestingly, that these two asperity sizes produce comparable $\phi _m$, consistent with figure 2. For example, at $\phi =0.68$, we see $\eta _r$ increase from ${O}(10^1)$ to ${O}(10^2)$ for aggregates as their surface roughness is increased from 0.05 to 0.20, compared with $\eta _r\approx 20$ for smooth particles with $\mu =100$. This is consistent with Pan et al. (2015) and Hsu et al. (2018), and a comparison with Tanner & Dai (2016) reveals that, at comparable proximity to jamming (comparing data at $\phi /\phi _m\approx 0.85$), moving from smooth particles to aggregates with $k_r=0.05$ leads to viscosity increases of 1.4 (experiment) and 2 (simulation). The particles used by Tanner & Dai (2016) were produced via grinding and have highly irregularly shaped asperities so a precise quantitative match is untenable. Notwithstanding the difficulty in comparing 2-D simulation with 3-D experiment, our aggregate model successfully captures the broad effects of varying surface roughness on suspension viscosity. A more quantitative comparison of simulation and experiment would require extensive 3-D computation involving aggregate particles with geometry closer to those of the experiment.

Looking at the stress dependence (figure 3), we see that for all suspensions thickening starts at a critical dimensionless stress $(\eta _r\dot {\gamma }^*)^{\dagger} \approx 1$ and reaches a plateau at a second $(\eta _r\dot {\gamma }^*)^\ddagger \approx 10^2$. These results being independent of $\phi$ and the particle roughness (both $\mu$ and $k_r$) is consistent with theory (Wyart & Cates 2014), experiment (Bender & Wagner 1996; Lootens et al. 2005; Brown & Jaeger 2012; Guy, Hermes & Poon 2015) and other simulations (Mari et al. 2014; More & Ardekani 2020b,c). The generality of this result could, therefore, be used to determine the critical shear rates required for the onset and plateau of shear thickening in a given system.

3.3. Range of repulsive force

Generically, the presence of a repulsive force leads to shear thinning (Ness et al. 2022). When a suspension is at very low $(\eta _r\dot {\gamma }^*)$, the relatively large value of $F^{rep}$ ensures particles will not come into contact with each other. Each particle is thus surrounded by a repulsive force field into which other particles cannot penetrate. This leads to them appearing artificially larger, with their radius becoming $a_{eff}=a_h+\varDelta$. As such we may define a $\phi _{eff}$ (larger than $\phi$) applicable to the shear thinning regime. Given a viscosity relation such as $\eta _r\sim (1-\phi /\phi _m)^{-\lambda }$, increases to $\phi$ will naturally lead to increased $\eta _r$.

The rate of decay $\rho$ associated with the repulsive force $F^{rep}$ plays a key role, most importantly affecting the rate of shear thinning. For frictionless smooth particles, increasing $\rho$ systematically lowers $\eta _r$ in the shear thinning region ($(\eta _r\dot {\gamma }^*)<(\eta _r\dot {\gamma }^*)^{\dagger}$) as particles are kept at further distance from one another, forcing a more isotropic structure. Increasing the stress shifts the system from being controlled by repulsive forces to being controlled by frictionless contact forces.

The behaviour for frictional smooth particles is rather more subtle. For $\rho \geqslant 0.1(a_i+a_j)$ there remains a systematic (but very weak) decrease in $\eta _r$ with increasing $\rho$ during shear thinning, with all cases eventually reaching the same plateau. For $\rho =0.01(a_i+a_j)$, however, the behaviour is qualitatively different: here, we see shear thinning at much lower $(\eta _r\dot {\gamma }^*)$, before a low $\eta _r$ plateau at $(\eta _r\dot {\gamma }^*)\approx 1$ that gives way to thickening at $(\eta _r\dot {\gamma }^*)>1$. In this case the particles behave almost as frictionless hard spheres until the stress reaches high enough values to induce frictional contact. This is because the added effective volume argument made above is less relevant (since in principle $(\phi _{eff}-\phi )\sim \rho ^3$) but nonetheless $F^{rep}$ is sufficient to prevent direct particle contact at low stress. Note that figures 4(a) and 4(b) are measured for comparable $\phi /\phi _m$ (which means they are at different $\phi$). The shear thickening behaviour in figure 4(b) is of the same type as that observed in the ‘electrostatic repulsion model’ of Mari et al. (2014) i.e. a short-ranged repulsive force gives way to contacts with Coulombic friction.

Figure 4. Rheology of smooth and aggregate particles, showing the relative viscosity $\eta _r$ as a function of repulsive force decay rate $\rho$ (see legend). In each case the solid fraction $\phi$ is chosen close to frictional jamming, so that $\phi /\phi _m\approx 0.93$. Shown are: (a) smooth frictionless ($\mu =0$) particles; (b) smooth frictional ($\mu =100$) particles; (c) aggregates with $k_r=0.05$; (d) aggregates with $k_r=0.10$; (e) aggregates with $k_r=0.20$. Insets show sketches of smooth and aggregate particles represented in each panel.

For aggregate suspensions the effect of $\rho$ is largely independent of the asperity size. Interestingly, for small $\rho$ (crucially this means $\rho \ll k_r a_h$) there is a regime in which the aggregate suspensions are no longer shear thinning. Here, the repulsive force decays so rapidly that it becomes too weak to prevent asperity interlocking: the suspension thus enters a ‘permanently thickened’ state in which asperities are always protruding. For larger $\rho$, the repulsive force extends beyond the perimeter of the asperities so that changing $\dot {\gamma }^*$ affects the potential for interlocking and therefore controls the rheological response. Specifically, aggregate suspensions show shear thinning followed by shear thickening as the stress is increased (qualitatively similar to frictional smooth particles).

Overall, where a shear thinning region does exist (that is, at low stresses where the coordination number $Z=0$) $\eta _r$ will be independent of the surface particle roughness (i.e. $\mu$ or $k_r$) at a given $\phi$ and $\rho$. Meanwhile, the plateau in $\eta _r$ above $(\eta _r\dot {\gamma }^*)^\ddagger$ is independent of $\rho$ but sensitive to $\mu$, $k_r$ and $\phi$.

3.4. Micromechanical mechanisms

We now consider the microstructural behaviour, focussing on the average number of contacts or interlocking particle pairs present at various $\phi$ and $\dot {\gamma }^*$. We use the average mechanical coordination number $Z$ (following Sun & Sundaresan 2011) calculated as $Z={2N_c}/{N}$, where $N_c$ is the number of particle–particle contacts and $N$ is the number of ‘whole’ particles (i.e. smooth particles or aggregates). To calculate $N_c$ under a given set of conditions, we first evaluated the number of particles separated by a distance less than a minimum distance required for contact. Importantly, smooth and aggregate particles will have different criteria for this. For smooth particles, contacts are identified based on intersecting surfaces (i.e. the scalar separation of the particle centres is less than the sum of their radii). For aggregates, contact detection is more involved as it becomes necessary to identify the point at which sufficient asperity interdigitation occurs to induce particle–particle interlocking (ultimately leading to shear thickening). To do this, we investigated a range of asperity interdigitation thresholds and determined which was most appropriate. We define a contact between aggregates $i$ and $j$ as occurring when the scalar separation between host particle centres $r$ is

(3.1)\begin{equation} r\leqslant (a_h^i + a_h^j) + (2-P)(k_r^ia_h^i + k_r^ja_h^j), \end{equation}

with $P$ the fractional interdigitation of the asperities for a pair of contacting particles (see figure 5a). For $P=0$ a contact occurs when one aggregate enters the outer perimeter set by the asperities of another (blue dashed circle in figure 1); for $P=1$ a contact occurs when the asperity of one aggregate intersects with the host of another. We thus define a set of new mechanical coordination numbers (subscripts indicate $P$ values): $Z_0$, $Z_{0.125}$, $Z_{0.25}$, $Z_{0.5}$, $Z_{0.75}$. We seek a suitable definition of $Z$ sufficiently stringent that the coordination is low (${\approx }0$) at low stress and approaching the isostatic value (${\approx }2\unicode{x2013}3$) at high stress (when $\phi$ is close to $\phi _m$). Presented in figure 5(b) are $Z$ data for aggregates with $k_r=0.20$, $\rho =3k_ra_h$ and $\phi =0.62$.

Figure 5. Mechanical coordination number $Z$ for smooth and aggregate particle suspensions. (a) Sketches of the relative particle positions for various values of the interdigitation fraction $P$; (b) The coordination number measured by different criteria of interlocking, for particles with $k_r=0.20$, $\phi =0.62$ and $\rho =3k_ra_h$; (c) The coordination number ($Z, Z_{0.5}$) measured at different $\phi$ for smooth and aggregate particles (isostatic points depicted by arrows).

For contacts defined with weak interdigitation ($Z_0, Z_{0.125}, Z_{0.25}$) the coordination numbers obtained at low stress are too high. That is, $Z$ measured by these metrics indicate substantial particle interlocking, inconsistent with the bulk shear thinning behaviour observed. The large value of $Z$ means that particles do come into close contact at all stresses, but for small $(\eta _r\dot {\gamma }^*)$ the repulsive force is sufficiently high to prevent substantive interlocking.

In contrast, strong interdigitation ($Z_{0.5}, Z_{0.75}$) only occurs at shear stresses above $(\eta _r\dot {\gamma }^*)^\ddagger$, appearing to align with the onset of shear thickening. While an overlap of 75 % only results in small increases in coordination number following shear thickening (settling at values of 0.5, figure 5b), there is a convincing increase in $Z_{0.5}$. This suggests that interlocking extending to around 50 % of the asperity size is a good microstructural indicator of the move to shear thickening (indeed the $Z_{0.5}$ behaviour mirrors that of the number of mechanical contacts per aggregate, while the effective area fraction of particles with radii $a_h + 0.5k_r a_h$ differs from that of our aggregate particles only by a factor of $0.25k_r^2$). Using the above definition, we proceed to address the contact number as a function of $\phi$ and $\dot {\gamma }^*$.

We find that increasing $\phi$ results in greater ($Z_0, Z_{0.5}$), figure 5(c), as expected from the bulk rheology (figure 2). Meanwhile, at a given $\phi$, ($Z_0, Z_{0.5}$) increases in line with the particle surface roughness, or asperity size. This latter point is easily explained by realising that particles of equivalent host size but different asperity size will exhibit different degrees of reach. In other words, particles with larger surface asperities will be more likely to come into contact since (for a given $\phi$) they sweep a larger overall perimeter than particles with smaller asperities. It follows from this that the size of the asperities on a particle will have a direct effect on $\phi _m$. To investigate this, we sought to establish the critical coordination number $Z_C$ for each suspension i.e. the mechanical coordination number at the jamming solids fraction $\phi _m$ (Meyer et al. 2010; Connelly et al. 2020).

The theoretical critical coordination number for a 2-D system with particle surface roughness between zero and infinity is a well-known result. For frictionless ($\mu =0$) discs $Z_C=4$, whilst discs with large sliding friction coefficient $\mu \gg 0$ have $Z_C = 3$. The constraints necessary to reduce $Z_C$ may be generated by either: (i) tangential forces resulting from an imposed friction coefficient, $\mu$, or; (ii) the presence of particle asperities which lead to particle–particle interlocking. It is thus expected that the aggregates would exhibit $Z_C\approx \leqslant 3$ since, when interlocking, their asperities should fully constrain the sliding motion of two contacting particles (analogous to large $\mu$) and may constrain further degrees of freedom such as rolling too. We report in figure 5(c) the values of $Z$ corresponding to the values of $\phi _m$ reported above, showing that indeed $Z$ for smooth particles matches our prediction. Interestingly we found $Z_{0.5}<3$ for all of the aggregate particles, with $Z_{0.5}$ increasing as $k_r$ increases. Surprisingly, this suggests that aggregates with large asperities behave most like frictional (large $\mu$) smooth particles, whereas those with small asperities have $Z_C<3$ and thus likely have additional constraints to motion other than sliding. If this trend continues down to smaller asperity sizes it would suggest that sliding friction alone is not a good approximation for the interaction between spheroidal particles with small asperities.

We also considered the average number of mechanical contacts associated with each suspension at various shear stresses, figure 6. Considering first the smooth particle suspensions, we see that both $\mu =0$ and $\mu >0$ experience no particle contacts at low stress, with $Z$ then sharply increasing, reaching a plateau similar to the trends seen in $\eta _r$ (figure 3). The same trend occurs for the aggregates: we find $Z_{0.5}\approx 0$ when below the critical stress (i.e. while the suspension is in the shear thinning regime) before rising sharply to values of between 1.5 and 3 as we enter shear thickening. The increase in particle coordination occurs over the same range of $\eta _r\dot {\gamma }^*$ as the shear thickening reported in figure 3. In light of these results, we conclude that the observed shear thickening behaviour for aggregates is a consequence of a sudden increase in the number of interlocking particle contacts (to around $50\,\%$ of the asperity size) above a critical shear stress.

Figure 6. Mechanical coordination number as a function of dimensionless shear stress $(\eta _r\dot {\gamma }^*)$ at various solid fractions $\phi$. (a) Smooth particles with $\mu =0$; (b) smooth particles with $\mu =100$; (c) aggregates with $k_r=0.05$; (d) aggregates with $k_r=0.10$; (e) aggregates with $k_r=0.20$. For (a,b) $\rho =0.1(a_i+a_j)$; for (c–e) $\rho =3k_ra_h$ and we use $Z_{0.5}$. Insets show sketches of smooth and aggregate particles represented in each panel.

On the whole our rheology data suggest that aggregate particle rheology can be described analogously to that of smooth particles, using an appropriate $\phi _m$ along with models for the viscosity such as Krieger & Dougherty (1959) and Wyart & Cates (2014). Relating the maximum packing to just a sliding friction coefficient is not straightforward, however, as our measured values of $Z$ suggest that roughness introduces other constraints such as rolling and twisting.

4. Concluding remarks

The importance of particle interlocking on shear thickening in suspensions of rough particles has been demonstrated through the use of a simple model comprising hydrodynamic contributions, normal contact forces and medium-range repulsion. Our model is both minimal in its physics and tractable in its computational expense, yet shows agreement with relevant experimental findings.

We successfully predict the relative increase in viscosity for suspensions of increasing surface roughness (Tanner & Dai 2016; Hsu et al. 2018) by explicitly modelling the contacts that arise between surface asperities on aggregate particles, a feat which is not possible in models that apply fictitious friction coefficients to smooth particles. Our model is complementary to the latter, allowing explicit interaction of asperities and offering a more resolved description of frictional contacts. We also demonstrated that particle interlocking is important in describing the evolution of $\eta _r$ with stress and jamming with increasing $\phi$, and were able to recover isostatic points in close agreement with theory.

We also invoked the idea of an effective solid fraction to explain the roughness-independent rheology in the shear thinning regime. In particular, we concluded that in the shear thinning regime particles are held in a near-homogenous, non-contacting state and their surface morphology becomes irrelevant to their interactions. As the imposed stress is increased, however, we essentially shrink the ‘protective bubble’ around the particles, exposing their asperities to interlocking. Here, surface details become important to the bulk rheology and allow shear thickening. We do not observe shear thinning in aggregates when the repulsive force is short ranged since the suspension is in a permanently thickened, interlocking state. Moreover, if $\phi$ is increased, sterics mean the particles will constantly be forced to contact one another regardless of the repulsive forces; leading to greater levels of interlocking and higher viscosities.

From this collection of observations we identify 3 key practical outcomes: (i) all things being equal, shear thickening will occur at lower shear rates for higher $\phi$ and/or higher $k_r$; (ii) as the length over which the repulsive force existing between particles is increased, the shear rate required for thickening will increase; (iii) for aggregate suspensions with comparable $\phi /\phi _m$ and the same range of repulsive force decay (with respect to the size of their asperities), the surface roughness will determine the point at which shear thickening takes hold. In particular, as surface roughness is increased the shear rate at which shear thickening occurs will decrease and the transition to thickening will take place at lower $\phi$. (This trend was also observed during shear stress controlled simulations (More & Ardekani 2020a).)

As is true for most dense suspension models (Ness & Sun 2016), we are unable to fully reproduce experimental normal stress behaviour. Our model overpredicts $N_1$ during shear thinning, yielding highly positive values in contrast to experiments which suggest $N_1\approx 0$. Hence, further investigation is required to understand the mechanism leading to this discrepancy. Our model quickly becomes intractable for $k_r<0.05$ or for systems much larger than ${O}(10^2)$ aggregate particles; thus, there is a clear need to develop new methods for dealing with particles of complex geometry.

We focused here on a simplistic geometry and a single value of $nk_r$ (the surface density of asperities). A complete understanding of the link between tribology of rough particles and rheology requires systematic exploration of a vast parameter space, evaluating the role of geometry on sliding and rolling friction and particle interlocking and possibly requiring a revisiting of our contact definition once $n$ becomes large. Similarly, since many processes involving dense suspensions (such as injection moulding, the spraying of agrochemicals and fibre spinning) involve mixed shear and extensional flows (Galindo-Rosales, Alves & Oliveira 2013; Andrade et al. 2020), future studies of aggregate particles should address more generalised deformation types as has been started for smooth particles (Ness et al. 2017a; Cheal & Ness 2018). Our model and its predictions demonstrate the new tribological insight that can be gained by moving beyond the current paradigm of imposed friction coefficients and seeking to build rheological understanding using models with more detailed contact physics. This sets the scene for further fundamental advances in achieving a robust fluid mechanics of dense suspensions, the outcomes of which may be used to allow more informed design and processing of suspensions in engineering applications.