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 Reference Kausch, Fesko and Tschoegl1971). 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 (Reference Reichhardt and Reichhardt2014) and O'Hern et al. (Reference O'Hern, Silbert, Liu and Nagel2003). 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. (Reference Hsu, Ramakrishna, Zanini, Spencer and Isa2018) as well as several works recounted by Hsiao & Pradeep (Reference Hsiao and Pradeep2019). Similar to Cheal & Ness (Reference Cheal and Ness2018) 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. Reference Nordstrom, Verneuil, Arratia, Basu, Zhang, Yodh, Gollub and Durian2010). 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 (Reference Egres and Wagner2005), Boersma, Laven & Stein (Reference Boersma, Laven and Stein1990) and Hsu et al. (Reference Hsu, Ramakrishna, Zanini, Spencer and Isa2018) 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. (Reference Pan, de Cagny, Weber and Bonn2015) and Hsu et al. (Reference Hsu, Ramakrishna, Zanini, Spencer and Isa2018), and a comparison with Tanner & Dai (Reference Tanner and Dai2016) 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 (Reference Tanner and Dai2016) 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 Reference Wyart and Cates2014), experiment (Bender & Wagner Reference Bender and Wagner1996; Lootens et al. Reference Lootens, Van Damme, Hémar and Hébraud2005; Brown & Jaeger Reference Brown and Jaeger2012; Guy, Hermes & Poon Reference Guy, Hermes and Poon2015) and other simulations (Mari et al. Reference Mari, Seto, Morris and Denn2014; More & Ardekani Reference More and Ardekani2020b,Reference More and Ardekanic). 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. Reference Ness, Seto and Mari2022). 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. (Reference Mari, Seto, Morris and Denn2014) 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 Reference Sun and Sundaresan2011) 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. Reference Meyer, Song, Jin, Wang and Makse2010; Connelly et al. Reference Connelly, Gortler, Solomonides and Yampolskaya2020).

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 (Reference Krieger and Dougherty1959) and Wyart & Cates (Reference Wyart and Cates2014). 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.