Force chains and networks: wet suspensions through dry granular eyes

Recent advances in shear-thickening suspension rheology suggest a relation between (wet) suspension flow below jamming and (dry) granular physics. To probe this connection, we simulated the contact force networks in suspensions of non-Brownian spheres using the discrete element method (DEM), varying the particle friction coefficient and volume fraction. We find that force networks in these suspensions show quantitative similarities to those in jammed dry grains. As suspensions approach the jamming point, the extrapolated volume fraction and coordination number at jamming are similar to critical values obtained for isotropically compressed spheres. Similarly, the shape of the distribution of contact forces in flowing suspensions is remarkably similar to that found in granular packings, suggesting potential refinements for analytical mean field models for the rheology of shear thickening suspensions.

vided the particles are well-stabilised, such suspensions typically shear thicken under flow, where the viscosity η increases with increasing stress (or shear rate) [1,2]. There is a growing consensus [2,3,4], based on recent experiments and simulations [5,6,7,8,9,10], that shear thickening results from the formation of frictional contacts. Friction constrains sliding motion, shifting the jamming volume fraction φ J from random close packing, φ RCP , for frictionless particles at low stress to φ m < φ RCP for frictional particles at high stress. Full flow curves η(σ) can be calculated using the phenomenological Wyart-Cates (WC) model [11], in which the stressdependent φ J (σ) interpolates between these two limits, reproducing both simulations and experimental results [12,13,14].
The importance of static friction and constraints in shear thickening [15] suggests an intimate connection with the physics of dry granular media. Similar features in the spectrum of vibrational modes (density of states) in suspensions and jammed grains have been shown in simplified simulations of frictionless suspensions [16,17,18], with low frequency 'soft modes' that vanish at a critical coordination number Z → Z J , driving the viscosity divergence. This result and the reduction of Z J in frictional packings [15,19] underpins the initial WC formulation [11], with φ m in suspensions corresponding to the limit of random loose packed frictional grains.
Motivated by a lack of obvious structural signatures of jamming in traditional real-space measures, there has been considerable effort devoted to characterizing the statistics and structure of the contact force networks of dry granular materials [20,21,22,23,24,25,26]. With limited exceptions [27], contact networks and force distributions in suspensions have received less attention. However, in a 'granular view' of shear thickening, force distributions are key [5,7], particularly at the thicken-ing transition where the form of the force distribution governs the fraction of frictional contacts at a given stress [12,28].
Here we probe this connection between suspensions and dry grains near jamming in detail, using DEM simulations to model suspensions of spherical particles at varying volume fractions φ and inter-particle friction coefficients µ. Computing the mean coordination number and force distributions in these simulations, we uncover numerous similarities between flowing suspensions near jamming and dry granular packs, supporting the 'granular view' of suspension rheology.

Methodology
We simulated 3D suspensions of neutrally-buoyant non-Brownian spheres in a Newtonian background fluid of viscosity η f . Hydrodynamic forces are calculated using the discrete element method (DEM) [29,5,30], including the Stokes drag and short-ranged lubrication interactions. The latter are regularized below a particle separation of ξ min where ξ = 2r/(d 1 + d 2 ) for two particles of diameter d 1 , d 2 at a center-to-center distance r [31,32]. We neglect long-ranged hydrodynamic interactions between particles because we only consider high φ.
The regularization of lubrication forces means that particles can come into contact. Such contacts are modeled by a linear Hookean spring F c n = k n h, with spring constant k n and the extent of particle overlap h. The tangential spring force is F c t = k t h t , where k t = (2/7)k n and the incremental tangential stretch h t is initialized to 0 at contact and updated following Silbert et al. [33]. The maximum F c t is set to satisfy the Coulomb criterion, |F c t | ≤ µ|F c n | , where µ is the friction coefficient. The critical load model (CLM), in which F c t = 0 when F c n ≤ F CLM , is used to simulate shear thickening [5,7]. F CLM mimics the repulsion that prevents facile particle contact [10,34], and sets a force scale for transition between frictionless to frictional flow.
Homogeneous simple shear at rateγ was imposed by affine deformation with Lees-Edwards periodic boundary condition in LAMMPS [35]. 1872 bi-disperse spheres at diameter ratio 1:1.4 (to prevent crystallization) mixed in equal volumes were simulated [5] at constant φ anḋ γ such that the Stokes' number ργd 2 /η f < 1, where ρ is the fluid or particle density.
The total stress was found by summing contributions from the contact, hydrodynamic forces between particles and isolated particle stresslets. The relative viscosity is taken to be η s = σ/(γη f ), where σ is the xy component of the stress tensor. All the results presented were obtained by averaging over at least 10 strain units  in steady state. To ensure hard-sphere like behavior, in all simulations we verified that the stiffness was set sufficiently high, k n σd, to avoid spurious shear thinning from particle overlaps [36].

Results
At fixed µ, the suspension's shear rheology is quasi Newtonian. Itsγ-independent viscosity increases with φ and diverges at a µ-dependent jamming volume fraction φ µ J , Fig. 1 (a), which we obtain by fitting where, following literature [14], we fix p = −2, leaving α µ and φ µ J as fitting parameters. As µ varies from 10 −4 to 10, φ µ J moved from 0.65 in the low friction limit to 0.57 in the high friction limit, agreeing with previous work [7,14]. The pre-factor varies weakly with µ and, despite previous suggestion of p varying somewhat from −2 in the frictionless limit [37], our results are consistent with p = −2 at all µ ( Fig. S2 [38]). Fits yielded nearly identical φ µ J even if the p is allowed to vary. As φ → φ µ J , the average per-particle contact (or coordination) number Z increases, Fig. 1 (b). We estimate Z at jamming, Z µ J , by linearly extrapolating data for φ µ J − φ ≤ 0.011 towards zero. Large fluctuations occurred for µ > 0.75 and φ µ J − φ 0.01 due to jamming and unjamming transitions in a finite-size system [39]. Thus, we do not extrapolate to find Z µ J above µ = 0.75. The plots of φ µ J and Z µ J against µ, Fig. 2, show a remarkably similarity in form compared to simulations of isotropically-compressed packings [19] and simple shear [40] of monodisperse granular spheres. Our φ µ J values are slightly higher than those in monodisperse systems, as expected [41]. The same is true of Z µ J even in the 0.56 This work [19] [40] Fig. 2 Dependence of (a) the critical jamming volume fraction φ µ J , and (b) the critical contact number Z µ J of bidisperse suspensions with friction coefficient µ compared with granular monodisperse spheres [19,40]. The inset shows the relation between Z µ J and φ µ J . The dashed line shows Eq. 2 with q = 1.676.
low-µ limit, where Z µ=0 iso = 6 is expected. It is not clear if the finite particle softness in DEM simulations or shear-induced anisotropic structures is responsible for this difference, which has also been observed for quasistatically sheared bidisperse granular spheres [42].
For µ 0.1, the low ∆φ data for different µ still collapse onto a power-law behavior with n ≈ 1. Now, however, we no longer find collapse at high ∆φ. Instead, data at increasing µ deviate from the n = 1 power law sooner and drop to a lower asymptote. Data at µ 2 is subject to significant uncertainty because of the impossibility of extrapolating to find Z µ J in Fig. 1(b). Instead, we solve for Z µ J using Eq. 2 by substituting values of φ µ J obtained previously from fitting Fig. 1(a). Including the points so obtained, we find that the data sets for µ 0.1 appear to converge toward a n = 0.3 power law, the same exponent as the low-µ data sets, but now with a lower amplitude. It is yet unknown if the power law differs at lower ∆φ at the highest µ.
The rigidity of granular packs is mediated through grain-grain contacts, which are distributed heterogeneously in space [47], with some grains bearing many times the mean force while others bear almost none. This is captured by the probability distribution P (θ), where θ ≡ F n / F n is the normal force between neighbors F n normalized by the mean. In dry granular packing, normal forces arise exclusively from direct contacts. In suspensions, inter-particle hydrodynamic interactions also contribute, although their exclusion does not change our results (Fig. S6 [38]). We calculate P (θ) from snapshots of the force network collected every 0.1 strain units averaged over 10 strain units, with error bars indicating the standard deviation.
In Fig. 4(a) we plot P (θ) for sheared suspensions close to jamming with ∆φ 5 × 10 −3 for a range of µ. In common with dry granular packings, we find exponential-like tails at high forces. At µ = 10 −4 , P (θ) is peaked at θ ≈ 0.5. As µ increases, the probability of low-force contacts increases, until at µ = 0.5 the peak in P (θ) is no longer evident. An increase in P (θ → 0) with increasing µ was also reported in simulations of dry granular packings [19]. We observe a similar trend for fixed µ = 0.2 at varying φ, with a slight increase in P (θ → 0) as ∆φ decreases (Fig. S3 [38]). Distributions of the tangential contact forces and the fraction of mobilized contacts are likewise similar to those obtained from dry-grain simulations (Fig. S4 [38]).
In compressed dry granular packings, P (θ) can be described by the empirical relation [21,23] We fit our data to this form, Table 1. The data and fit, Fig. 4(a) (full line), at µ = 10 −4 are more pronouncedly peaked than the P (θ) for compressed amorphous packings of smooth spheres, Fig. 4(a) (dot-dashed). At µ = 0.5, the peak in P (θ) is no longer evident: b = 0 and the distribution is purely exponential, Fig. 4(a) (dashed). Note that our results at high µ show, and our results at low µ are consistent with, a finite plateau rather than a power-law scaling as θ → 0 [17,20,45,48].
To examine the force network during shear thickening, we implement a critical load model (CLM), where frictional contacts are activated whenever F n exceeds a threshold F CLM . For a shear-thickening suspension with φ = 0.54, µ = 0.5, our data for P (θ), Fig. 4(b), are almost identical to those at fixed µ, Fig. 4(a), and can be fitted to Eq. 3 in the low-and high-stress limits using the previous low-and high-µ limit parameters (solid and dashed lines respectively). In the CLM, F CLM determines the onset stress σ * = F CLM /(1.5πd 2 ) for shear thickening [5,7,28]. The WC jamming volume fraction φ J (σ) shifts from φ 0 J to φ µ J as the fraction of frictional contacts f (σ) increases from f → 0 for σ σ * to f → 1 for σ σ * , giving a thickening flow curve η(σ). Typically, f (σ) = exp(−σ * /σ) in fitting WC-type models to experiments or simulations [12,14], wherẽ σ * ≈ σ * to O(1). In general, the fraction of frictional contacts is given by where θ CLM (σ) = F CLM / F n (σ) . Assuming F n ∝ σ, we can write θ CLM = ασ * /σ, where the pre-factor α 1.85 is found to be independent of φ (see Fig.  S7 and associated discussion [38]). For a simple exponential force distribution, corresponding to our high-µ form of Eq. 3 with b = 0 and β = 1, integrating Eq. 4 gives f (σ) = exp(−ασ * /σ), a result previously derived to motivate an exponential form for f (σ) [12,28]. This procedure can be repeated using the low-µ fitting parameters in Eq. 3, where b, c = 0. Interestingly, the resulting f (σ) obtained using either low-and high-µ fitted parameters in Eq. 4, corresponding to the solid and dashed lines in Fig. 4(c), differ little; each gives a reasonable account of the data. This reflects the dominant importance of the high-force exponential tail, which remains largely unchanged as either µ or φ is varied.

Discussion and Conclusions
We find that sheared dense suspensions exhibit a number of quantitative similarities to dry granular materials near jamming. The µ-dependence of the volume fraction and coordination number at jamming, φ µ J and Z µ J , are reminiscent of corresponding functions in isotropic sphere packings [19]. Thus, theoretical models relating φ µ J and Z µ J in the latter may be extendable to suspensions [44,45]. The force distribution in suspensions also recalls that in granular packings. An exponential tail dominates in the µ → 0 and ∞ limits, justifying the use of a simple exponential form for the fraction of frictional contacts in WC-type models of shear thickening.
Despite these similarities, we find that the relation between ∆Z and ∆φ is more complex in suspensions than in dry grains, undercutting a number of assumptions underpinning the WC model for shear thickening. As initially formulated [11], WC drew upon simulations of frictionless hard spheres [16,17], which found η ∼ ∆Z −ν with ν ∼ 2. The basic physics is that 'soft modes') are lost as the system approaches isostaticity (∆Z → 0). These soft modes are characterized by the vibrational density of states D(ω), which gives the number of modes per particle at a given frequency. In dry granular packings above jamming, D(ω) trasitions from classical Debye scaling (∝ ω 2 in three dimensions) to a constant low-frequency plateau as |Z − Z µ J | → 0 in both the frictionless and frictional case [15,49]. A similar change in the shape of D(ω) was seen in simulated frictionless suspensions as Z → Z µ=0 J [16].
WC assume that the same physics applies in wet frictional suspensions, so that the shape of D(ω) is likewise controlled by Z µ J and the η ∼ ∆Z −ν scaling applies in both the frictionless and frictional case. To make useful predictions, they then need to relate the 'natural variable' in dry grains, Z, to the 'natural variable' in suspensions, φ, for which they make the further assumption that ∆Z ∝ ∆φ in both the low-and high-µ limits. We find such proportionality only for µ 0.06 and very close to jamming (∆φ 0.03), Fig. 2. In these low-friction suspensions there is a rollover to a weaker power law at higher ∆φ, and for higher values of µ we never reach a regime where ∆Z ∝ ∆φ.
Despite these discrepancies, Eq. 1 fits our results over a relatively wide range of ∆φ for all values of µ, and the WC model formulated in terms of φ has proved successful in capturing experimental results [12,13]. This suggests that the empirical relation between η and ∆φ should in fact be viewed as more universal, while the 'soft mode' approach to understanding the viscosity divergence in suspensions may have limited applicability. Indeed, this is in line with work to develop empirical constitutive relations for frictional suspensions by analogy to empirical constitutive relations for flowing grains [50]. It has also recently been shown that viscosity divergence η ∝ ∆φ −2 can be derived for suspensions of frictionless spheres through a non-equilibrium kinetic theory approach [51], and this approach could perhaps be extended to frictional suspensions as well.
In our analysis of the force networks in these suspensions, we have neglected the contribution from contact anisotropy to the stress and possible microstructure evolution during shear thickening. However, given the similarities between the force networks in wet suspensions and dry grains, it is likely that more advanced methods used to characterize force networks in jammed packings [52,53,54] could be applied to open a new window into the rheology and dynamics of dense suspensions.