A constitutive model for simple shear of dense frictional suspensions

Dense non-Brownian suspensions of rigid particles exhibit diverse rheological behavior, including yielding, shear thinning or shear thickening, normal stress differences, particle migration and shear jamming [1–4]. Shear thickening of dense suspensions is a phenomenon in which, for a range of applied shear stress, the apparent viscosity increases, sometimes by an order of magnitude or more [5–7]. Strong continuous shear thickening (CST) is observed in concentrated suspensions, where the viscosity increases continuously with an increase in shear rate. At volume fractions $ϕ$ above a critical value $ϕC$⁠, an abrupt increase in viscosity may be observed and is termed discontinuous shear thickening (DST).

Recent experimental [8–13] and computational [8,14–18] work has demonstrated that shear thickening (both CST and DST) can arise due to frictional particle-particle contacts. In a suspension, a repulsive force is often present as a result of steric (e.g., due to adsorbed polymer) or electrostatic stabilization. When the shear forces acting to bring a pair of particles into contact exceed the repulsive force threshold, only fluid mechanical forces are available to keep the particle surfaces apart, and Melrose and Ball [19] have shown that the lubrication film can go to arbitrarily small values. We therefore assume the lubrication film can break, allowing for contact interactions of both normal and tangential (frictional) form, which most simply can be seen as representative of direct surface contacts due to roughness, but could also be a result of other phenomena, e.g., polymer brush interactions [8]. Thus the repulsive force F₀ gives rise to a stress scale $σ0=F0/6πa2$ for particles of radius a, which marks a crossover from lubricated (frictionless) contacts between particles to direct, frictional contacts. This concept can be idealized, by considering the stabilization forces to be of zero range so that they do not contribute a noncontact stress, but set a stress scale σ₀ for the onset of frictional interaction, as described in defining “the critical load model” [16]. The behavior is then divided into two regimes which are essentially rate independent, i.e., where the stresses are linear in the deformation rate. One is the case at low stress (⁠$σ≪σ0$⁠), where particle interactions are lubricated, with the viscosity diverging at the frictionless jamming point $ϕJ0$⁠, corresponding to random close packing of $ϕrcp≈0.64$ for monodisperse spheres. The second is at high stress (⁠$σ≫σ0$⁠) where almost all contacts are frictional, and the viscosity diverges at a volume fraction $ϕJμ<ϕJ0$⁠. Here, $ϕJμ=ϕJ(μ)$ is used as shorthand to denote the dependence of this jamming fraction on the interparticle friction coefficient μ [16,20–24]. With increase in σ, the crossover between these two states results in the shear thickening behavior. At large $ϕ$⁠, a finite-range stabilizing force becomes one of the major sources of stress at low stresses, which leads to a shear-thinning rheology. In this work, we focus on developing a constitutive model for the shear-thickening part of the flow curves. Our proposed model does not include the physics behind the shear-thinning at low stresses, but as this is due to an unrelated microscopic mechanism with its own set of parameters, our idealized model could in principle be augmented to include the relevant low stress physics and display both shear thinning and shear thickening.

Along this line of thought, aspects of which had been developed in [25], Wyart and Cates [26] (WC) have shown that under rather broad conditions a viscosity increasing with shear stress by interpolating between two rate-independent asymptotic rheologies can lead to three forms of shear stress curves as a function of shear rate (cf. Fig. 1), depending on the viscosity difference between the two states representative of unthickened and thickened suspensions. When this difference is small, i.e., for $ϕ≪ϕJμ$⁠, the shear thickening is continuous, with a monotonic relation between shear stress and rate. For a large enough but finite viscosity contrast, i.e., when $ϕC<ϕ<ϕJμ$⁠, the flow curve $σ(γ̇)$ becomes nonmonotonic, S-shaped, and the thickening becomes discontinuous. Finally, when the thickened branch is actually jammed, i.e., for $ϕ≥ϕJμ$⁠, the system can flow only for low stresses, while at high stresses frictional contacts cause the system to shear jam, and in this case the flow curve at large stress tends toward a zero shear rate state. Nonmonotonic flow curves have been reported under controlled stress conditions in several subsequent experimental and simulation studies [10,17,27,28].

Wyart and Cates [26] have also proposed a rheological model for shear thickening exhibiting the features noted above. This model is based on an interpolation between two diverging stress-independent rheologies, where the interpolation depends on a unique microscopic state parameter identified as the “fraction of frictional contacts,” f. Wyart and Cates described f as a function of Π, the particle pressure, but in standard rheometric experiments, where $ϕ$ is fixed, Π and σ are directly related [29]. Since σ is more readily controlled, we find it more convenient to consider $f(σ̃)$⁠, where $σ̃=σ/σ0$⁠. In rate-controlled simulations, f has been found to be a function of $σ̃$ [16], while in pressure-controlled simulations, f is not found to be a unique function of Π [30]. The quantity f serves, in essence, as an order parameter for the shear thickening transition, assuming low values in the low viscosity state under small stress, and asymptoting to $f≈1$ in the large viscosity state at large stress.

In this article, we explore the WC model by comparison to extensive numerical simulations we have performed by varying the volume fraction and friction coefficient of the microscopic model and exploring the resulting stress (or rate) dependence. Informed by the numerical results, we then introduce empirical expressions for the relations between the parameters μ, $ϕ$⁠, the order parameter f, and the steady-state stress. We find that upon increasing the friction coefficient, the frictional jamming point $ϕJμ$ decreases, but the power-law exponent of divergence for both shear and normal stresses, expressed as $(ϕJ−ϕ)−β$ remains well-described by β = 2. We show that the WC model is highly effective at predicting the shear stress for a wide range of μ and $ϕ$⁠. An extension following the same model structure provides predictions for the normal stress response to allow for a complete description of the viscometric functions in steady simple-shear flow of shear thickening suspensions.

We simulate an assembly of inertialess frictional spheres immersed in a Newtonian fluid under an imposed shear stress σ, giving rise to an imposed velocity field $v=γ̇(t)v̂(x)=γ̇(t)(x2,0,0)$⁠. We use Lees-Edwards periodic boundary conditions with N = 500 (or to test size dependence N = 2000 in a few cases) particles in a unit cell. To avoid ordering, we use bidisperse particles, with radii a and $1.4a$ mixed at equal volume fractions. The particles interact through short-range hydrodynamic forces (lubrication), a short-ranged repulsive force and frictional contacts; this simulation model has been shown to reproduce accurately many features of the experimentally measured rheology for dense shear-thickening suspensions [14,16], although discrepancies have been noted [31] in the small first normal stress difference relative to some experimental observations [32].

The equation of motion for N spheres is the 6N-dimensional force/torque balance between hydrodynamic (⁠$FH$⁠), repulsive (⁠$FR$⁠), and contact (⁠$FC$⁠) interactions

where the particle positions are denoted by X and their velocities/angular velocities by U. $FR$ is a conservative force and can be determined based on the positions X of the particles, while the calculation of the tangential component of the contact force $FC$ is more involved as it depends on the deformation history of the contact.

We make the translational velocities dimensionless with $γ̇a$ and the shear rate and rotation rates by $γ̇$⁠. Decomposing the dimensionless velocity as $v̂(r)=ω̂×r+ê·r$ in rotational $ω̂=(0,0,−1/2)$ and extensional $ê12=ê21=1/2$ parts, the hydrodynamic force and torque vector takes the form

with $Û∞=(v̂(y1),…,v̂(yN),ω̂(y1),…,ω̂(yN))$ and $Ê=(ê(y1),…,ê(yN))$⁠. The position dependent resistance tensors $RFU$ and $RFE$ include the “squeeze,” “shear,” and “pump” modes of pairwise lubrication [33], as well as one-body Stokes drag. The occurrence of contacts between particles due, for example, to surface roughness is mimicked by a regularization of the resistance divergence at vanishing interparticle gap $hij=2(rij−ai−aj)/(ai+aj)$⁠: The squeeze mode resistance is proportional to $1/(h+δ)$⁠, while the shear and pump mode resistances are proportional to $log(h+δ)$ [16], where h is the dimensionless gap scaled by particle radius a. Here, we take $δ=10−3$⁠. The lubrication force is upper limited, and negative h, i.e., particle overlaps (contacts) are allowed in our simulations.

The electrostatic repulsion force is taken to represent a simple electrostatic double layer interaction between particles, with the force decaying exponentially with the interparticle surface separation h as $|FR|=F0 exp(−h/λ)$⁠, with a Debye length λ.

Contacts are modeled by linear springs and dashpots. Tangential and normal components of the contact force $FC(ij)$ between two particles satisfy Coulomb's friction law $|FC,t(ij)|≤μ|FC,n(ij)|$⁠, where μ is the interparticle friction coefficient. Some softness is allowed in the contact, but the spring stiffnesses are tuned for each $(ϕ,μ,σ)$ (that is, stiffer springs at large $(ϕ,μ,σ)$⁠) such that the largest particle overlaps do not exceed 3% of the particle radius during the simulation, so stay near the rigid limit [16,34].

The equation of motion (1) is solved under the constraint of flow at constant shear stress σ. At any time, the shear stress in the suspension is given by

where η₀ is the suspending fluid viscosity, $ηHγ̇=γ̇V−1{(RSE−RSU·RFU−1·RFE):Ê∞}12$ is the contribution of hydrodynamic interactions to the stress, and $σR,C=V−1{XFR,C−RSU·RFU−1·FR,C}12$⁠, where $RSU$ and $RSE$ are position-dependent resistance matrices giving the lubrication stresses from the particle velocities and resistance to deformation, respectively [16,35], and V is the volume of the simulation box. Note that the resistance tensors are proportional to the suspending fluid viscosity. At fixed shear stress σ, the shear rate $γ̇$ is the dependent variable that is to be determined at each time step by [17]

The full solution of the equation of motion (1) under the constraint of fixed stress (3) is thus the velocity [17]

From these velocities, the positions are updated at each time step. Lastly, the unit scales are $γ̇0≡F0/6πη0a2$ for the strain rate and $σ0≡η0γ̇0=F0/6πa2$ for the stress. In the rest of the paper, we use scaled stress defined as $σ̃=σ/σ0$⁠.

In this section, we present simulation results for values of friction coefficient $μ=0.1,0.2,0.5,1,5$⁠, and 10. As the friction coefficient determines $ϕJμ$⁠, the results allow an exploration of the effect of the separation between $ϕJμ$ and $ϕJ0$ on the rheology. As previously shown [20,21,23], $0.1≤μ≤1.0$ corresponds to the moderate friction limit, where the friction coefficient affects the jamming point $ϕJμ$⁠; $μ>1$ corresponds to large friction, and we find that $ϕJμ$ saturates rapidly with $μ>1$⁠.

In the following, the shear stress σ, particle pressure Π, and normal stress differences N₁ and N₂ are defined as $σ≡Σ12$⁠, $Π≡−(Σ11+Σ22+Σ33)/3, N1≡Σ11−Σ22$⁠, and $N2≡Σ22−Σ33$⁠, respectively.

The basic assumption of the model [26] is that $ηr, Π/η0γ̇0$⁠, and $N2/η0γ̇0$ are in distinct stress-independent states both at low $(σ̃≪1)$ and high $(σ̃≫1)$ stress. Here, we model the $ϕ$ dependence of viscosity and normal stresses as $(ϕJ−ϕ)−2$ and $ϕ2(ϕJ−ϕ)−2$⁠, respectively. The forms are consistent with the proposed correlations for constant volume [29,36] as well as constant pressure conditions [37].

The viscosity, particle pressure, and second normal stress difference as functions of $ϕ$ in low and high stress states can be expressed as

where $α0,μ, β0,μ, and K20,μ$ are friction dependent constant coefficients, independent of $ϕ$ (and σ). Recall that $ϕJ0$ and $ϕJμ$ denote the jamming volume fraction for μ = 0 (frictionless state) and nonzero values of μ (frictional states), respectively. Different functional forms for shear stress and particle pressure [note the leading $ϕ2$ term in Eq. (6c)] leads to a stress ratio $μbulk=σ/Π$ (or $q=Π/σ$ in suspension flow modeling [29]) being volume fraction dependent, consistent with the experimental results of Boyer et al. [37].

The friction-dependent parameters are found empirically from our simulations to be expressible as functions of μ,

as shown by the fits in Fig. 3. The fitting parameters $μϕ$ and $μα$ are reported in Table I; note that this table contains all friction-independent parameters of the model.

Next, we specify the flow curves utilizing a stress-dependent jamming volume fraction, using an expression similar to that proposed by Wyart and Cates [26],

where $f(σ̃)$ denotes the fraction of close particle interactions in which shear forces have overcome the stabilizing repulsive force F₀ to achieve contact. We propose stress-dependent coefficients

Finally, using Eqs. (6)–(8) we propose the dependences of the rheological functions on $σ̃,ϕ$⁠, and μ,

The divergences of the rheological functions described above—viscosity, N₂, and Π—have the same algebraic sign at low and high stress. By contrast, N₁ presents a special case, in that it appears to have different signs under conditions dominated by lubrication and friction [38–41]. We model N₁ as

Now the stress-and volume fraction-dependent N₁ can be written as

where $K1m(σ̃)$ is given by

The transition between the lubricated and frictionally dominated stress states is captured by the fraction of frictional interactions, $f(σ̃)$⁠, which we model as $f(σ̃)=exp[−σ̃∗/σ̃]$⁠, with $σ∗=1.45σ0$ based on simulations here and previously published results [11,16,18,28,42]. We assume that $f(σ̃)$ does not depend on μ.

Before turning to stress dependence, we show in Fig. 2 our simulation results for $ηr$⁠, Π and N₂ for rate-independent states. These agree well with Eq. (6) for all values of μ. The viscosity at large $ϕ$ is well represented by $(ϕJ−ϕ)−2$ independent of μ, as shown in Fig. 2(b); $ϕJ$ is obtained by a least-squares fit of Eqs. (6a) and (6b) to the volume fraction dependence of the viscosity at low $(0.1<σ̃<0.3)$ and high $(σ̃>10)$ stresses, respectively.

The friction-dependent constants are found to fit well to exponential functions of the form proposed in Eq. (7) as shown in Fig. 3; values of these parameters are presented in Table I. Because N₁ proves more difficult to reliably simulate (or experimentally measure [32,43,44]), we defer its consideration to Sec. IV C.

To develop a sense of the entire flow behavior, we present in Fig. 4 the viscosity data with interparticle friction coefficient μ = 1. This value of μ is comparable to the experimentally measured values of Fernandez et al. [8], where μ is in the range 0.6–1.1 for polymer brush-coated quartz particles of diameter $2a∼10$ $μm$⁠, but is higher than the value of 0.5 reported by Comtet et al. [13]. The data are presented in the forms $γ̇(σ̃)/γ̇0, ηr(γ̇/γ̇0)$⁠, and $ηr(σ̃)$ for a range of values of $ϕ≥0.45$⁠. Figure 4(a) shows that for volume fractions $ϕ<0.56$ the $γ̇(σ̃)/γ̇0$ curves are monotonic and show CST for $0.3≤σ̃≤10$⁠. At $ϕ=0.56$⁠, the curve exhibits the first sign of nonmonotonicity: We define $ϕC≐0.56$ to signify the DST onset volume fraction. For $ϕ=0.57$ and $ϕ=0.58$⁠, the slope is negative (i.e., $dγ̇/dσ̃<0$⁠) for intermediate stress but crosses over to a positive slope for $σ̃>10$⁠, corresponding to the S-shaped $ηr(γ̇/γ̇0)$ curves shown in Fig. 4(b). The viscosity as modeled in Eq. (9a) is shown by solid lines and agrees with the simulation data except at high stresses at $ϕ=0.58$⁠. This discrepancy is due to the closeness to $ϕJ(μ)$ given in Eq. (7a), which slightly underestimates the jamming volume fraction for μ = 1.

Although plotted as a function of shear rate, it is important to note that these simulations were performed at fixed shear stress: DST would be observed at fixed rate for this range of volume fraction (⁠$0.56≤ϕ≤0.58$⁠) as S-shaped curves are not accessible in a rate-controlled scenario [17,26]. For this range of $ϕ$⁠, where DST is observed between two flowing states, we term the discontinuous shear thickening as “pure DST.” However for $ϕ≥0.59$⁠, the upper, i.e., high stress, branch of the S-shaped $γ̇(σ̃)$ curves [Fig. 4(a)] is not accessible, signifying that the suspension enters a shear jammed (SJ) state above a critical stress depending on $ϕ$⁠, which in concept is similar to the one observed in dry granular systems [45] (a difference is that steady flow cannot occur at low stress for granular systems). The curve defining these critical stresses is explained in Sec. IV B. We term such DST, in which the thickening continues until reaching a shear-jammed state, as “DST-SJ.” When the traditional flow curve $ηr(γ̇/γ̇0)$ is plotted for the same data in Fig. 4(b), we observe nonmonotonicity for $ϕ≥0.56$⁠. These data, when presented in the form $ηr(σ̃)$⁠, show that the onset stress for shear thickening $σ̃ST≈0.3$ is roughly independent of volume fraction, as observed in previous studies [7,16,46,47], except close to jamming where it steadily decreases. To characterize shear thickening as CST or DST, we fit $ηr(σ̃)$ in the thickening regime to $ηr∼σ̃ζ$⁠: $ζ<1$ implies CST, while ζ = 1 implies the onset of DST, and larger ζ are in the DST region. We also observe some shear thinning at low stress as a result of the short-range electrostatic repulsion [16].

The viscosity obtained from Eq. (9) is shown in Fig. 5 along with simulation data for different values of μ. The model is in excellent agreement with the data. The phenomenology is the same for any value of μ expect for μ = 0. With increasing μ, the DST onset volume fraction $ϕC$ decreases. It is also important to mention that the simulation results presented here do not show any noticeable system size dependence, as shown in supplementary Fig. S3 [57]. Increasing the number of particles by a factor of four only decreases the fluctuation in the data, with no strong effect on the mean viscosity. In the proximity of jamming (⁠$ϕJ−ϕ<0.01$⁠), we do observe some size dependence, as is usually observed with systems showing continuous phase transitions [48,49].

The shear rheology described above is controlled by three dimensionless parameters, namely, the solid volume fraction $ϕ$⁠, the dimensionless shear stress $σ̃$⁠, and the interparticle friction coefficient μ.

The results discussed are presented in a flow state diagram. As this depends on three variables, we present two views: In the $ϕ$ − $σ̃$ plane for μ = 1, and in the μ − $ϕ$ plane for stress $σ̃C$⁠. Figure 6(a) displays the observed flow state diagram in the $ϕ$ − $σ̃$ plane, and here we identify three volume fractions: $ϕC, ϕJ(μ)$⁠, and $ϕJ0$⁠. Vertical lines represent frictional $ϕJ(μ)$ and frictionless $ϕJ0$ jamming points. In the lower part of the diagram, where the stress is too low to overcome the interparticle repulsive force, friction between close particles is not activated and hence the rheology diverges at $ϕJ0$⁠. However, in the upper part of the flow state diagram where the stress is large, most of the close interactions (or “contacts”) are frictional, which leads to divergence of the viscosity and other rheological functions at $ϕJ(μ)<ϕJ0$⁠. In the two extremes, the viscosity in the model is rate independent. However, in the simulations at low stress, the finite range of repulsion leads to a larger apparent particle size, and the competition between short-range repulsion and external applied stress creates a shear thinning behavior at the conditions indicated by the + symbols. For intermediate stress, CST is observed in the range of $ϕ<ϕC$⁠. For $ϕC≤ϕ<ϕJ(μ)$⁠, “pure DST” is observed (shown by triangles). In this range of $ϕ$⁠, the dashed line is the envelope of the pure DST states, with $(ϕC,σ̃C)$ being the point with the minimum $ϕ$ value along this line. This line is determined as the locus of points for which $dγ̇/dσ̃=0$ in a flow curve $γ̇(σ̃)$ as shown in Fig. 4(a): There are two such points on a curve for any $ϕ>ϕC$⁠, and coalescence of these two points occurs at a critical point (⁠$ϕC,σ̃C$⁠).

For $ϕ>ϕJ(μ)$⁠, the upper boundary of DST states is the stress-dependent jamming line $ϕm(σ̃)$⁠. The jamming line separates the DST regime (diamonds) from conditions yielding a shear-jammed state (squares). The shear rates in these states are vanishingly small, i.e., the suspension no longer flows. The distinction between two types of DST regimes is based on differences in the high stress state, which is flowable in the pure DST regime (for $ϕC≤ϕ<ϕJ(μ)$⁠) and jammed for DST-SJ (for $ϕJ(μ)≤ϕ<ϕJ0$⁠). The minimum stress required to observe DST and SJ states decreases with increasing $ϕ$⁠, and eventually these curves converge and the minimum stress for jamming tends to zero as the frictionless jamming point $ϕJ0$ is approached. Our findings are in qualitative agreement with the recent experimental study on shear jamming by Peters et al. [4], which shows that a minimum volume fraction is required to observe SJ states, and that the shear stress required to shear jam the system decreases with an increase in volume fraction and vanishes as $ϕ$ approaches the isotropic jamming point. The general structure of the flow-state diagram presented here is similar to the one presented by Peters et al. [4], with the details like the shape of $ϕm(σ̃)$ being different in the two cases.

Figure 6(b) displays the flow state diagram in the μ − $ϕ$ plane for a constant stress $σ̃C$⁠, the onset stress for DST for minimum volume fraction $ϕC$⁠; we note that $σ̃C$ is roughly independent of μ. Here, we see a demonstration of how the volume fractions $ϕC, ϕJ(μ)$⁠, and $ϕm(σ̃C)$ decrease as functions of μ. The region enclosed between $ϕC$ and $ϕJ(μ)$ broadens at larger μ, illustrating that the range of $ϕ$ over which pure DST is observed broadens with increasing interparticle friction. For the range of volume fractions $ϕJ(μ)<ϕ<ϕm(σ̃C)$⁠, the suspension is in the DST-SJ region, and above $ϕm(σ̃C)$⁠, the system is in the shear-jammed state at the imposed stress.

The simulation data along with the model predictions for the particle pressure $Π/η0γ̇$ and second normal stress difference $N2/η0γ̇$ are presented in Fig. 7. The proposed model is in good agreement with the simulations. We observe that $N2/η0γ̇$ is always negative, and is comparable to but smaller than $ηr$⁠. For volume fraction $ϕ≤0.45, Π/η0γ̇$ is smaller than $ηr$⁠. With increasing $ϕ$ the particle pressure increases faster than the shear stress, and for $ϕ$ approaching $ϕJ(μ)$⁠, $Π/η0γ̇$ becomes larger than the relative viscosity $ηr$⁠, as deduced in modeling based in part on particle migration data by Morris and Boulay [29]. The experimental data by Boyer et al. [37] also show similar decrease in bulk friction coefficient as the jamming volume fraction is approached.

Finally, $N1/η0γ̇$ from Eq. (11) together with the simulation data are shown in Fig. 8. Figure 8(a) displays the divergences of $N1/η0γ̇$ in stress-independent states, where we choose the divergent volume fraction to be the same as that of $ηr$ (as well as particle pressure $Π/η0γ̇$ and second normal stress difference $N2/η0γ̇$⁠) for μ = 1. We use $K10=0.055$ and $K1(μ)=0.045$ to match the model with the simulation data for $ϕ>0.54$⁠. The predictions of the model for the stress-dependent $N1/η0γ̇$ are plotted along with the simulation data in Figs. 8(b) and 8(c). The modeled $N1/η0γ̇$ exhibits several features: For all volume fractions $ϕ, N1/η0γ̇$ is small and negative at small stress, and becomes increasingly negative with increasing stress, reaching a minimum for $σ̃≈1$⁠. The magnitude of this negative $N1/η0γ̇$ becomes larger with increase in $ϕ$⁠. At stress values $σ̃>1, N1/η0γ̇$ tends toward positive values, crossing zero at $σ̃p$⁠, which is independent of $ϕ$⁠. The proposed model is in good agreement with the simulation data for volume fraction $ϕ≥0.54$⁠, where the simulation data show positive N₁ at high stress. On the other hand, for volume fraction $ϕ<0.54$⁠, the model does not agree with the data at high stress, where simulations show negative N₁, while the model predicts N₁ to be positive. The simulation data and model predictions for N₁ at larger volume fractions $ϕ≥0.54$ agree in being negative at small stress (where lubrication films between most particles remain) and becoming positive at large stress (where most of the contacts are frictional). This also agrees in part with observations of Lootens et al. [44], Dbouk et al. [50], and Royer et al. [42], but not with the data of Cwalina and Wagner [32]. However, at lower volume fractions $ϕ<0.54$⁠, experiments which have shown $N1<0$ for all stresses [32,42,44] are in agreement with our simulations, but are not captured by the model. This suggests that there is a difference in microstructure between the lower and higher particle fractions such that behavior consistent with the lubricated regime is observed in the lower-$ϕ$ suspension even when contacts are frictional.

The rheology of dense frictional suspensions determined by the extensive numerical simulations presented here displays continuous and discontinuous shear thickening (CST and DST, respectively) and, at sufficiently large $ϕ$⁠, shear-induced jamming. All of these behaviors are predicted by the model structure of Wyart and Cates [26]. We provide a thorough examination of the influence of the interparticle friction coefficient on the jamming fraction in the near-hard-sphere limit. We find that the approaches to the jamming point of both the shear and normal stresses scale well with volume fraction as $(ϕJ−ϕ)−2$ for either the lubricated or frictional case. The comprehensive and coherent database from these simulations allows us to make detailed comparisons against a constitutive model incorporating stress-dependent frictional effects in dense suspensions.

We find that a model defined by three parameters—solid volume fraction $ϕ$⁠, dimensionless stress $σ̃=σ/σ0$⁠, and interparticle friction coefficient μ—captures well the extreme rate-dependence of the rheology of these materials. Here, $σ0=F0/6πa2$ is a stress scale determined by a stabilizing repulsive force of magnitude F₀ at contact for particles of radius a. The central concept is that this stress scale divides the material response into low stress and high stress regions: When the stress is small, $σ̃≪1$⁠, particle surfaces remain separated by lubrication films and the viscosity and normal stresses are relatively small. When the stress overwhelms the repulsive force, i.e., when $σ̃≫1$⁠, frictional contacts dominate and the rheological functions are much larger. This leads to two limiting jamming fractions: $ϕJ0$ in the frictionless states at low stress, and $ϕJ(μ)<ϕJ0$ at high stress. A stress-dependent jamming fraction $ϕm(σ)$ can be defined by interpolating between these two jamming fractions as a function of the applied stress, as shown in Eq. (8a), in the manner proposed by Wyart and Cates [26] using the fraction of frictional contacts. The divergence of the stresses approaching the interpolated jamming fraction $ϕm$ follows the form of the two limits, with the stresses growing as $(ϕm−ϕ)−2$⁠.

Based on these concepts, a constitutive model for dense frictional suspensions in steady, simple shear flow has been proposed. Comparison of the model predictions, e.g., relative viscosity $ηr(ϕ,σ̃,μ)$⁠, agree well over the full range of parameters with the simulations reported here. The overall behavior is described by a flow-state diagram given in Fig. 6. This diagram, particularly in the $ϕ$ − $σ̃$ plane, displays the various regions of material behavior obtainable at a fixed value of μ. At smaller $ϕ$⁠, the material shear thickens continuously (CST), while above a critical solid fraction $ϕC$ the shear thickening becomes discontinuous. We find two regimes of DST: (i) A pure DST regime between two flowing states for $ϕC<ϕ<ϕJ(μ)$⁠, and (ii) a DST-SJ regime where with the increase of $σ̃$⁠, DST gives way to shear-jamming for $ϕJ(μ)<ϕ<ϕJ0$⁠. In both the scenarios (i) and (ii), upon increase in $σ̃$ the suspension under shear goes through CST over a range of stress before entering DST. With increasing μ, $ϕJ(μ)$ decreases, and as a consequence the range of $ϕ$ over which shear jamming is observed, i.e., $ϕJ(μ)<ϕ<ϕJ0$⁠, increases. The range of volume fraction $ϕC<ϕ<ϕJ(μ)$ for which DST is observed also broadens with increase of μ. The onset stress to the shear-jammed state decreases with increase in $ϕ$⁠, tending to zero as $ϕ→ϕJ0$⁠. This can be seen by considering $ϕm(σ̃)$ curve in the flow state diagram.

Once the key parameters in the model are fitted the entire flow-state diagram can be constructed. To achieve the fitting, only a measure of the two jamming volume fractions (⁠$ϕJ0$ and $ϕJ(μ)$⁠) and a stress ramp $ηr(σ)$ at one volume fraction $ϕ$ are required. The two jamming fractions can also be extracted from stress ramps at several $ϕ$ provided these are sufficiently concentrated. It is important to mention that the model is expected to be applicable only for volume fractions where the frictional contribution to the viscosity at high stress is greater than the lubrication contribution.

The three material functions are sufficient to define the rate-dependent coefficients in the CEF (Criminale, Ericksen, Filbey) formalism [51,52], which is exact for steady viscometric flows of a nonlinear material whose stress depends only on the history of the rate of deformation. The CEF equation could be used for continuum simulations of nonviscometric flows in conjunction with a model for particle migration like that of Morris and Boulay [29]. The current formulation should be applicable to concentrated Brownian suspensions, since the Brownian force term can be replaced by an additive term in the interparticle force [17,53].

Finally, we expect that the formulation of the model itself should be robust to changes of particle properties such as polydispersity, particle shape, and other surface properties. However, the values of the parameters such as $ϕJ$ are known to be sensitive to these details [54–56].

The author's code makes use of the CHOLMOD library by Tim Davis (http://faculty.cse.tamu.edu/davis/suitesparse.html) for direct Cholesky factorization of the sparse resistance matrix. This work was supported, in part, under the National Science Foundation Grant Nos. CNS-0958379, CNS- 0855217, and ACI-1126113 and the City University of New York High Performance Computing Center at the College of Staten Island. J.F.M. was supported by NSF 1605283.