A constraint-based approach to granular dispersion rheology

We present a phenomenological model for granular suspension rheology in which particle interactions enter as constraints to relative particle motion. By considering constraints that are formed and released by stress respectively, we derive a range of experimental flow curves in a single treatment and predict singularities in viscosity and yield stress consistent with literature data. Fundamentally, we offer a generic description of suspension flow that is independent of bespoke microphysics.

We present a phenomenological model for granular suspension rheology in which particle interactions enter as constraints to relative particle motion. By considering constraints that are formed and released by stress respectively, we derive a range of experimental flow curves in a single treatment and predict singularities in viscosity and yield stress consistent with literature data. Fundamentally, we offer a generic description of suspension flow that is independent of bespoke microphysics.
Such behavior likely stems from details of the particle interactions [2] set by, e.g., surface chemistry [19] or roughness [20]. Most models incorporate such interactions in a bespoke manner. Notably, a phenomenological model by Wyart and Cates (WC) [21] predicts thickening (class 2) due to a transition from frictionless (static friction coefficient µ ≈ 0) to frictional (µ > 0) particle contacts above a critical "onset stress". Atomic force microscopy confirms this picture for several systems [15,22] and the WC model fits a number of experimental flow curves [7,8,18]; although, quantitative discrepancies with microscopic simulations remain [23].
To recast the WC model within a more general framework, recall that frictional contacts constrain interparticle sliding. Crucially, the WC model is agnostic to the exact mechanism by which sliding is constrained, so that disparate microphysics, e.g., stress-induced interlocking of asperities [20,24], hydrogen bonding [25] or 'traditional' Coulomb friction can all give rise to the same macroscopic, shear-thickening phenomenology.
In this broader framework, the WC model deals with a single type of constraint: sliding. Rolling (rotations about axes perpendicular to the line of centres) and twisting (rotations about the line of centres) degrees of freedom remain unconstrained. By assuming that sliding constraints are formed at increasing stress, the WC model accounts for class 2 behavior, which, however, is rare in practice. Real systems are typically class 1 or 3, for which current explanations involve the ad hoc "bolting together" of different kinds of bespoke physics [10].
Here, we generalize the WC model to two constraint types, A and B ≠ A. For example, A = sliding and B = rolling (e.g., due to adhesive contact [26]). Constraint A is formed by stress, while constraint B is released by stress. Now, one single model predicts all observed classes of experimental flow curves. This includes class 3b, whose rheology we find is sensitive to the exact interplay of constraint formation and release, thus explaining the variability of such systems observed experimentally. Moreover, we make non-trivial predictions for the emergence of singularities that are consistent with literature data, e.g., re-entrant jamming in class 3a and a yield stress diverging below random close packing for class 1.
We begin by considering the case of a single constraint type, A, which is formed by increasing stress σ. In the original WC model, A = sliding. There are two possible contact states: one in which A is constrained, A, and one in which A is unconstrained, A. (Thus, e.g.,Ā='non-sliding' and A='sliding'.) These states are associated with different jamming points, φ A and φĀ; thus, e.g., a suspension with all contacts in state A jams (η → ∞) at φ A . Existing literature suggests that φĀ < φ A [27][28][29]. Thus, for monodisperse spheres with A = sliding, φ A ≈ φ rcp ≈ 0.64 (random close packing) and φĀ ≈ φ rlp ≈ 0.55 < φ A [27] (random loose packing [30]).
Under shear, A is constrained in a σ-dependent fashion, so that a fraction fĀ of the contacts are in theĀ state, with dfĀ dσ > 0. A suspension with fĀ = 0 (all contacts in state A) jams at φ A ; when fĀ = 1 (all contacts in stateĀ), it jams at φĀ. For 0 < fĀ < 1, the system jams at an intermediate volume fraction, φ J (fĀ), the functional form of which is not known for all constraint types. We use WC's form for constrained sliding [21]: (1) The viscosity η increases as the distance to jamming, ∆φ = φ J (fĀ) − φ, decreases, diverging as φ → φ J . Again, we use WC's form for sliding constraints [21]: where η 0 is the viscosity of the suspending medium. This single-constraint scenario leads to shear thickening: increasing σ increases fĀ, decreasing φ J (and ∆φ), arXiv:1807.11356v1 [cond-mat.soft] 30 Jul 2018 Eq. (1), driving up η, Eq. (2). The exact form of the flow curve η(σ) depends on fĀ(σ), which encapsulates the stress-dependent microphysics between particles. Now introduce a second constraint type B ≠ A that is released by stress, giving four contact states: AB, AB, AB andĀB (whereB/B means B is/is not constrained), each with an associated jamming point: φ AB , φĀ B , φ AB and φĀB. Random close packing at φ AB = φ rcp is the same for all identities of A and B, since both are unconstrained in the AB state. The other jamming points depend on the nature of A and B, and are unknown in general. In static packings of dry grains, combining multiple constraint types typically lowers the minimum packing fraction for mechanical stability [29,31,32]. We suppose that this also applies to the jamming point of sheared suspensions, so that φĀB < φĀ B , φ AB < φ AB .
The rheology for the two-constraint case is far richer than that for a single constraint. Even rather bland choices for fĀ(σ) and fB(σ) readily lead to flow curves of all classes, Fig. 1, including the little-understood class 3b. Within each class, particularly class 3, the exact phenomenology is sensitive to the detailed interplay between fĀ and fB, i.e., to the relative formation of type-A and release of type-B constraints with stress.
We present predictions using for fĀ(σ): where σ A is the characteristic stress for the formation of type-A constraints and α controls the rapidity of type-A constraint formation with σ, and for fB(σ): where σ B is the characteristic stress for the release of type-B constraints and β controls the rapidity of type-B constraint release. The rheological class of the system is now controlled by σ A σ B and α β. The exact scenario depends on the relative values of the jamming points; however, the qualitative phenomenology is the same provided φĀB < φĀ B , φ AB < φ AB . To make concrete predictions, Class 1 (shear thinning) arises whenever σ A σ B ≪ 1, independent of α β, Fig. 1(a-b) (see caption for parameters). Since σ A is readily exceeded, A is always constrained and fĀ = 1, Fig. 1(b). The shape of the flow curve η(σ), Fig. 1(a), reflects fB(σ): it shear thins as  type-B constraints are progressively released, reaching a Newtonian plateau when they are all removed.
Class 3b (thickening then thinning) occurs only when σ A σ B ∼ 1, Fig. 1(e-f). Now, the form of η(σ) cannot be easily deduced from the evolution of fĀ(σ) and fB(σ), and the viscosity "peak" is profoundly sensitive to changes in σ A σ B . Thus, for α β = 1.4 and φ φ AB = 0.69, halving σ A σ B drops the peak by a factor of ten. There is a similar sensitivity to α β at fixed σ A σ B , which we explore in the Supplementary Information (SI) [33]. In each case, however, the extended region of shear thinning after the peak is controlled by the progressive release of constraints on B, fB → 0, with constraints on A almost saturated, fĀ ≈ 1.
Such sensitivity is consistent with experiments, which find class 3b η(σ) sensitive to small changes in, e.g., particle size [11] and suspending medium [12]. Presumably, the resulting alterations to particle surface properties perturb the stress-dependent microphysics and hence fĀ and fB, leading to O(1) variations in η(σ).
Our model also predicts a plethora of class 3c flow curves (mapped out fully in the SI [33]), including a commonly-observed variation of class 3a in which shear thickening is followed by further thinning [6,10,14,34]. Such behavior arises if constraints on B are released more slowly than constraints on A are formed (α β > 1).
Thus, considering two constraint types gives all three classes of flow curve. In the SI [33], we discuss the transition between the classes as σ A σ B is varied. In our scheme, system-specific interactions affect the rheology only in how they constrain sliding, rolling and twisting, and how the resulting constraints depend on σ. Constraints control the jamming point, φ J , Eq. (3), which is therefore itself σ dependent, φ J (σ) ≡ σ J [fĀ(σ), fB(σ)]. At fixed φ, η(σ) is determined by the (σ-dependent) distance to jamming ∆φ(σ) = φ J (σ) − φ, Eq. (2). Some interactions do not impose constraints on interparticle sliding, rolling or twisting, e.g., conservative repulsive or attractive interactions (electrostatic, depletion etc.), but still lead to aγ-or σ-dependent rheology. How, then, does one determine whether the rheology is driven by constraints, or by other physics?
To answer this, note that at fixed φ the system can be in one of two rheological states, depending on σ. If φ J (σ) > φ, the system flows. If φ J (σ) ≤ φ, it jams, i.e., η → ∞. Although non-constraint physics can give rise to such bipartite behavior, our model makes specific predictions for the emergence of viscosity singularities that are not expected to arise generically otherwise. Whenever confirmed, these would rule in constraint physics.
We explore these predictions for class 1 and class 3a (see SI [33] for class 3b) and compare them to literature data. Figure 2(a) shows a typical σ-φ "phase diagram" for class 1 (see caption for parameters). Red states are jammed, φ J (σ) ≤ φ, and white states are flowing, φ J (σ) > φ. The boundary of jammed states, σ jam (φ), is defined by φ J (σ) = φ, which we solve for numerically using Eq. (3)(4)(5). Figure 2(b) (lines) shows flow curves at different φ generated using the same model parameters. The form of η(σ) (vis-á-vis singularities) at any given φ can be derived by tracing a vertical path with increasing stress in Fig. 2(a). Thus, at some φ < φĀB the system flows at all σ. At φ ≥ φĀB, the system is jammed until σ exceeds σ jam , whereupon it flows, i.e., the system has a yield stress σ y = σ jam . This σ y increases with φ, diverging at φĀ B < φ AB . Thus, we predict jamming at any σ if φĀ B ≤ φ < φ AB .
Here, σ y [ Fig. 2(a) (symbols)] emerges at φ ≈ 0.2, which we take as φĀB, and diverges at ≈ 0.55, which we take as φĀ B . Since φĀ B ≈ φ rlp , A = sliding in this system. The phase diagram for class 3a is more complex. We plot a representative example for φ AB < φĀ B in Fig. 3(a) (see SI [33] for the qualitatively similar cases of φ AB = φĀ B and φ AB > φĀ B ) and a corresponding set of flow curves at different φ in Fig. 3(b) (lines). Note first that there is a φ window between φĀ B and some φ max in which σ jam (φ), Fig. 3(a) (solid line), is double-valued. For fixed φ in this window, we predict re-entrant jamming: the system un-jams above a yield stress σ y [=lower part of σ jam (φ), green] and subsequently re-jams at a higher stress [=upper part of σ jam (φ), blue] [35]. Secondly, there is a critical φ max < φ AB above which the system is jammed at all σ. Note that φ max is not associated with any divergence; rather this is the φ at which the yield stress (green) and re-jamming stress (blue) become equal. Thus, as φ → φ − max , the stress window of flowing states shrinks, Fig. 3(b), and vanishes at φ max .
There is a region in Fig. 3(a) (gray) in which dσ dγ < 0, corresponding to d log η d log σ > 1 in Fig. 3(b) (gray lines). Here, shear flow is unstable [36]. The boundary of this region, σ uns (φ) (black dashes), corresponds to solutions of d log η d log σ = 1 and meets the yield (green) and re-jamming (blue) stresses at φ max . Under imposeḋ γ, σ uns is the maximum stress for stable flow.
Most systematic work on class 3a singular behavior has been done for aqueous cornstarch dispersions (R ≈ 7 µm) [16][17][18]37]. Under imposed σ, the only study of the entire phase diagram [16] found pronounced shear thinning followed by a completely rigid state, consistent with our prediction. Typical phenomenology under imposeḋ γ is shown in Fig. 3(b) (symbols) [17]. Cornstarch flows steadily above a yield stress [◻ in Fig. 3(a)] until some higherγ limit is reached, above which the suspension separates into shear bands, Fig. 3(b) (gray symbols). The stress at this upper limit plausibly corresponds to our σ uns (φ) line. Indeed, the experimental σ uns [○ in Fig. 3(a)] and σ y (◻) appear to converge at φ max ≈ 0.45 ≪ φ rcp . As we would predict, flow was found impossible for φ beyond this limit. The existence of a continuum of fully jammed states at φ max < φ < φ rcp , a surprise for the authors of Ref. [17], emerges naturally from our model.
Note that all classes of granular suspension flow curve can be obtained by appealing to a single constraint, but with f (σ) mimicking η(σ) in each case, e.g., a peaked f (σ) will give class 3b flow curves [38]. This contrasts with our appeal to two constraints with simple, generic forms for f , Eqs. (4) and (5), to give all classes of flow curves. A single constraint can lead to phase diagrams that are topologically similar to Fig. 2(a) and Fig. 3(a); however, singular behaviour can arise only in a relatively narrow φ window, φĀ ≤ φ < φ A = φ rcp .
We should emphasize that our framework has been constructed for granular suspensions. Recently, class 3a flow curves, Fig. 3(b), have been observed in simulations of Brownian particles with inter-particle friction (= constrained sliding) and conservative potential attraction [39], giving a yield stress that depends much more weakly on φ (roughly ∝ φ 2 ) than is observed in granular systems, Fig. 3(a) (◻). If the potential attraction is strong enough, the yield stress "masks" shear thickening and the system becomes class 1, Fig. 2(b). However, the physics of yielding in a system where constraints act in the presence of Brownian motion is little understood and may be quite distinct from that explored here.
Before concluding, we turn briefly to microphysics. While much is known about the tribology of dry grains, little is known about the effect of a solvent on sliding, rolling and twisting resistance. If, like others [21,40], we assume that 'dry tribology' is possible in a suspension, then one realization of our scheme is the finite-area, adhesive contact with asymmetric pressure distribution between frictional particles [26,41]. Such a contact is pinned by surface roughness, leading to a critical torque M * for "peeling" particles apart and σ B ∼ M * R 3 . Interestingly, simulations of adhesive, frictional dry grains [29] are consistent with the observation of φĀ B ≈ φ AB for cornstarch in Fig. 3(a). On the other hand, our AB state (sliding unconstrained, rolling constrained) has no obvious dry-granular analogue, suggesting that the contact physics of dispersed particles holds surprises.
To conclude, we have presented a constraint-based phenomenological model for granular dispersion flow that predicts all known classes of experimental flow curve. Several non-trivial predictions for the emergence of singular behavior follow, which are borne out by literature data. Our notion that system-specific microphysics enters the rheology only on the level of constraints is a powerful one. It allows disparate systems to be captured in a single treatment agnostic of microphysics, much like the way multifarious details are subsumed into a "tube constraint" in polymer rheology [42]. Many challenges remain, including making a formal link between constraints, jamming and dissipation [43,44] and probing the microphysics of constraints in the presence of solvent.