Particles suspended in a Newtonian fluid raise the viscosity and also generally give rise to a shear-rate dependent rheology. In particular, pronounced shear thickening may be observed at large solid volume fractions. In a recent article [R. Seto et al., Phys. Rev. Lett. 111, 218301 (2013)], we have considered the minimum set of components to reproduce the experimentally observed shear thickening behavior, including discontinuous shear thickening. We have found frictional contact forces to be essential and were able to reproduce the experimental behavior by a simulation including this physical ingredient along with viscous lubrication. In the present article, we thoroughly investigate the effect of friction and express it in the framework of the jamming transition. The viscosity divergence at the jamming transition has been a well known phenomenon in suspension rheology, as reflected in many empirical laws for the viscosity. Friction can affect this divergence, and in particular the jamming packing fraction is reduced if particles are frictional. Within the physical description proposed here, shear thickening is a direct consequence of this effect: As the shear rate increases, friction is increasingly incorporated as more contacts form, leading to a transition from a mostly frictionless to a mostly frictional rheology. This result is significant because it shifts the emphasis from lubrication hydrodynamics and detailed microscopic interactions to geometry and steric constraints close to the jamming transition.

Suspensions (solid particles immersed in a fluid) exhibit a wide range of rheological behaviors, including shear thinning, shear thickening, and finite normal stress differences. Shear thickening [Barnes (1989); Mewis and Wagner (2011); Brown and Jaeger (2014)], where the viscosity increases with shear rate even if the suspending fluid is Newtonian, is a particularly intriguing phenomenon. In the extreme case of discontinuous shear thickening [DST; for early descriptions see the work of Williamson (1930); Williamson and Hecker (1931); Freundlich and Roder (1938)], which is observed at high volume fractions of solid material, the viscosity can increase by several orders of magnitude at a critical shear rate. Note that we do not include here the irreversible shear thickening occurring due to particle aggregation {as seen for example in fumed silica [Crawford et al. (2013)]} under the name DST.

The variety of systems showing a DST suggests that this may be a universal behavior that is obtained with a minimum set of physical ingredients. Even though most of the data available are for Brownian suspensions (that is, submicrometer particles) [Metzner and Whitlock (1958); Hoffman (1972); Bender and Wagner (1995, 1996); Frith et al. (1996); Fagan and Zukoski (1997); Boersma et al. (1990); D'Haene et al. (1993); O'Brien and Mackay (2000); Maranzano and Wagner (2001b,a, 2002)], thermal motion does not seem necessary to observe DST, as experiments with micrometer scale [Boersma et al. (1990); Lootens et al. (2004); Lootens et al. (2005); Larsen et al. (2010)] or larger particles [Freundlich and Roder (1938); Bertrand et al. (2002); Brown and Jaeger (2009, 2012); Fall et al. (2010, 2012)] show. Inertia has also been associated with the existence of shear thickening in suspensions [Lemaître et al. (2009); Fall et al. (2010); Trulsson et al. (2012); Fernandez et al. (2013); Kawasaki et al. (e-print)], based on the initial ideas of Bagnold (1954), but the Bagnoldian scaling of a viscosity proportional to the shear rate is clearly milder than the abrupt shear thickening observed in experiments. More importantly, the thickening due to inertia arises at Stokes number S t ρ a 2 γ ̇ / η 0 (with γ ̇ being the shear rate, η0 being the viscosity of the fluid phase, and ρ and a are the mass density and size of the solid particles, respectively) of order S t 1 in the simulations [Trulsson et al. (2012); Fernandez et al. (2013); Kawasaki et al. (e-print)], which contrasts with the values of at most S t 10 3 found in experiments {see for example [Maranzano and Wagner (2001b); Fall et al. (2010); Brown and Jaeger (2012)]}. Thus inertia is not necessary for DST in the conditions probed by rheometric flows. (Flows involving significantly larger particles, say in the mm range, or significantly larger shear rates γ ̇ 1 s 1 , can probe the regime where inertia matters).

Indeed, the apparently simple experimental system of nearly rigid non-Brownian neutrally buoyant particles immersed in a Newtonian fluid exhibits DST in the Stokes regime [Brown and Jaeger (2009, 2012)], indicating that the phenomenon stems from a rather restricted set of simple ingredients. The puzzle has very few pieces. In the Stokes regime, however, non-Brownian neutrally buoyant hard spheres in a Newtonian fluid will create a suspension whose rheology is independent of shear rate, as there is only one force scale, the hydrodynamic one (we will further explain this fact in Sec. II B 1). So there cannot be too few pieces.

The flow of suspensions has historically been studied from a fluid mechanics perspective, and the emphasis has been on a description based on hydrodynamic interactions. Suspensions are usually described as hard particles immersed in a Newtonian fluid, interacting through hydrodynamics (including Brownian forces) and sometimes an additional soft repulsive potential (approximating an electric double-layer or mimicking polymer coating, for instance). They are typically studied in the Stokes flow regime, at vanishing particle Reynolds number R e ρ 0 a 2 γ ̇ / η 0 (with ρ0 being the mass density of the fluid phase) and Stokes number. The key point is that the hard cores of the particles are treated as boundary conditions for the Stokes equations that describe the fluid phase; the particles never directly generate forces through contacts. This treatment is self-consistently justified by the fact that the Stokes flow between two rigid surfaces leads to a lubrication force whose resistance coefficient diverges at contact, effectively preventing two particles from colliding. Within this framework, shear thickening is explained by the creation at large shear rates of locally denser clusters of particles (hydro-clusters), which are highly dissipative due to the singular lubrication flows between the particles [Brady and Bossis (1985); Bender and Wagner (1996); Melrose and Ball (2004a); Wagner and Brady (2009)]. While this purely fluid mechanical point of view is able to describe the flow and rheology of moderately concentrated suspensions, it seems unable to explain the abrupt shear thickening observed in highly concentrated suspensions. Simulations by Stokesian Dynamics give a weak logarithmic shear thickening [Brady and Bossis (1988); Bossis and Brady (1989); Phung et al. (1996); Foss and Brady (2000); Melrose and Ball (2004b)], thereby raising the issue of the validity of this approach at high volume fractions. Moreover, DST is also observed in athermal suspensions, where the shear rate dependence introduced by the Brownian force scale disappears: In this regime, a purely hydrodynamic perspective would predict a shear-rate independent rheology.

In the past few years, new ideas have emerged from the granular rheology perspective. Boyer et al. (2011) developed an analogy between the rheology of suspensions and the rheology of granular materials. In dry granular flow, the rheology depends on the ratio between inertia and particle pressure [Deboeuf et al. (2009)]: If the inertia dominates, particles bounce around and the stresses are essentially due to momentum exchange based on collisions (a “granular gas” regime), whereas if the pressure dominates, particles are forced to stay at contact and the stresses are dominated by contact force chains (a “granular packing” regime). This perspective can be incorporated in a dimensionless inertial number I γ ̇ d ρ / Π [da Cruz et al. (2005)], where d is the diameter of the solid particles, ρ their density, and Π is the particle pressure. Similarly, for suspensions in the Stokes regime, there is a viscous number I v η 0 γ ̇ / Π , where η0 is the viscosity of the suspending fluid, which compares viscous dissipation to the particle pressure [Boyer et al. (2011)]. The stresses can be dominated either by viscous dissipation if I v 1 (when particles are far apart) or by a contact network if I v 1 . This results in constitutive laws for the viscosity and other stress components and for the volume fraction that are unique functions of Iv. One can anticipate that the regimes I 1 and I v 1 , both of which are dominated by the approach to the jamming transition, will share some similarities. Indeed, they do share a power-law scaling (typical of the jamming transition) between particle pressure and distance to the jamming point Π ( ϕ J ϕ ) λ [Forterre and Pouliquen (2008); Boyer et al. (2011)]. The scaling is the same for the other stress components, including the important case of shear stress σ ( ϕ J ϕ ) λ . These power laws are now argued to be a generic behavior close to jamming [Lerner et al. (2012); Andreotti et al. (2012)], though the value of the exponent λ might be system dependent (it might depend on shape or friction for instance). With these ideas, the emphasis has shifted from hydrodynamics and detailed microscopic interactions to geometry and steric constraints close to the jamming transition.

While those new ideas have proven to be successful in explaining some of the rheology at high volume fractions [Boyer et al. (2011); Trulsson et al. (2012); Lerner et al. (2012)], they cannot account for non-Newtonian behavior in the Stokes regime. The case of DST is particularly puzzling: Whereas intuition suggests that DST is a manifestation of jamming [Cates et al. (1998); Bertrand et al. (2002); Lootens et al. (2003); Hébraud and Lootens (2005); Brown and Jaeger (2009)], DST is not captured by the Iv-based rheology, which predicts no shear rate dependence in the Stokes regime [Lerner et al. (2012); Trulsson et al. (2012)], and reduces to a ϕ -dependent rheology [Boyer et al. (2011)]. Simply stated, this rheology predicts the correct volume fraction dependence of the viscosity at fixed shear rate but completely fails at predicting the shear rate dependence at fixed volume fraction. Again, it is worth noticing that this absence of shear rate dependence hints at a missing force (or time) scale in the description of dense suspensions.

In this article, which extends and complements a previous publication [Seto et al. (2013)], we show how jamming is related to shear thickening and how to escape the apparent contradiction described above. The starting point is to note that the jamming transition, while it exhibits some generic features that are independent of the microscopic interactions like the scalings described above, retains a clear signature of the microscopic details in the volume fraction at which it occurs. For instance, it is well known that jamming can occur anywhere between ϕ J μ = 0.55 and ϕ J μ = 0 0.64 depending on the friction coefficient μ between the grains [Silbert et al. (2002); Jerkins et al. (2008); Song et al. (2008); Silbert (2010)]. This means that the relation ϕ ( I v ) actually depends on μ, which in turn implies that the viscosity η ( ϕ ) depends on μ. Therefore, if there exists a mechanism such that μ (or an effective friction coefficient) varies with the shear rate γ ̇ , one can obtain a non-Newtonian rheology within the framework of Boyer et al. (2011). This mechanism has to be associated with an extra force scale besides the hydrodynamic one in order to provide the shear-rate dependence.

We introduced such a mechanism in the recent article [Seto et al. (2013)], and here we thoroughly explore its consequences through numerical simulations of dense (i.e., highly concentrated) non-Brownian suspensions under simple shear flow. We note that two other recent papers have independently introduced similar but somewhat different mechanisms to achieve a rate-dependent friction coefficient μ ( γ ̇ ) in simple shear flow simulations [Fernandez et al. (2013); Heussinger (2013)]. Both of them prove that interparticle friction is essential to shear thickening of dense suspensions. In contrast to those works, the present work however shows that inertia is inessential to this phenomenon, as we obtain it in an overdamped simulation (see Sec. II). This is important because results coming from simulations including inertia give shear thickening at shear rates that are several orders of magnitude larger than the experimental values (see previous Sec. I A). We show here that our inertialess model gives a quantitative description of shear thickening (see Sec. V).

The present work shows the importance of frictional contacts in dense suspension rheology. Of course, particles must form contacts during flow in order for the rheology to depend on the friction coefficient, and here we come up against one consequence of the fluid mechanics perspective. In fact, it is known that particles do come in contact under flow, even in the dilute limit, perhaps due to surface roughness [Davis (1992); Zhao and Davis (2002); Lootens et al. (2005); Blanc et al. (2011)]. One simply has to accept that those contacts, though not essential to the rheology of dilute or semidilute suspensions (as the successes of Stokesian Dynamics demonstrate), become important at large volume fractions. But contacts are a necessity if one thinks of the large volume-fraction limit of jamming ( I v 0 ) : There is no more room in this limit for lubrication films between particles, and contacts must proliferate and dominate the rheology. Besides enforcing the geometric constraint of no overlapping particles, those contacts bring a new restriction to reorganization at the microscopic level: They carry significant tangential forces due to friction. It is worth noting that even though the close-range hydrodynamics (i.e., lubrication forces) also generate some tangential forces, they are of a fundamentally different nature. For a relative motion between two nearby particles, the normal lubrication force diverges as the inverse of the interparticle gap, whereas the tangential lubrication force diverges only logarithmically with vanishing gap, which means that the effective friction coefficient vanishes in the limit of zero gap for smooth bodies. Thus lubrication, as in its literal meaning, provides very little resistance to relative tangential motion compared to relative normal motion. This stands in contrast to the frictional contact forces, for which the friction coefficient is finite; i.e., tangential and normal contact forces are of comparable scale. It follows that the constraint introduced by contact friction is much more efficient at increasing the viscosity than is lubrication.

Describing the physics of contact and motion when the gap between two particles in a suspension is smaller than, say, 1 nm, is not an easy task. Moreover, one cannot expect to find generic behavior, as interactions will vary depending on the nature of the solid particles and of the suspending fluid. However, this level of detail may not be essential to understanding the physics that emerges at a macroscopic level. At high volume fractions, it is reasonable to expect that the singularity associated with the jamming transition will dominate the rheology of suspensions. Our approach in this work is thus to study the rheology of minimal model systems that include hydrodynamics, a repulsive interaction, and frictional contacts between particles. These interactions, while being realistic, will be simplified to their essence.

In this section, we will describe our numerical model of a dense non-Brownian suspension under shear flow [Denn and Morris (2014)]. We intend to have a minimal model containing as few physical elements as possible while exhibiting a DST. Even if the model is intended to be as generic as possible, a few parameter choices are necessary. For those choices, we try to set values that correspond to a realistic interpretation of a suspension of hard (say, silica) spheres with radius in the 1–10 μm range, density matched (or such that the sedimentation occurs over times much larger than the inverse shear rate) and charge stabilized. The only parameter that takes an unrealistic value is the particles' stiffness, which is smaller in the simulations than in a typical experimental system. A stiffness comparable to experimental values for silica (or even Polymethyl methacrylate (PMMA)), would require the use of very small time steps in the simulation in order to resolve the dynamics of the contacts which correspond to an overlap in the nm range. We however circumvent this problem by tuning the particle stiffness with applied shear rate such that we get rid of particle deformation effects in the rheology (see Sec. II A 3 and Appendix  B for details).

1. Equations of motion

In a suspension, the flow of the fluid is described by the Navier–Stokes equation, while the motion of the particles is given by Newton's equations of motion

m · d d t ( U Ω ) = α ( F α ( U , Ω , r ) T α ( U , Ω , r ) ) ,
(1)

where m is the mass/moment-of-inertia matrix, U and Ω are the translational and rotational velocity vectors, and F α and T α are the force and torque vectors, respectively. The right-hand side consists of two different types of forces and torques: Some depend only on the configurations r N of the particles (e.g., forces derived from a potential), but some also depend on the velocities U and Ω , including inelastic contact forces (e.g., contacts with dashpot terms) and fluid-particle interactions (hydrodynamic interactions).

We will study the rheology of this suspension under an imposed simple shear flow field U ( r ) expressed using the vorticity Ω and the rate-of-strain tensor E as

U ( r ) = Ω × r + E · r .
(2)

At a shear rate γ ̇ , a simple shear flow corresponds to the following nonzero elements: Ω 3 = γ ̇ / 2 and E 12 = E 21 = γ ̇ / 2 .

In many experimental flows, including the ones for which shear thickening is typically observed, the particle Reynolds number and the Stokes number are very small due to the small size of the suspended particles. This means that the equation of motion for the particles can be studied in its overdamped limit, which is simply a quasistatic force balance equation. As will be discussed in Secs. II A 2, II A 3, and II B below, in this regime the velocity dependence coming from the hydrodynamic interactions and the contact forces has a linear form, and the force balance equation can be written as a linear algebraic relation having the form

0 = R · ( U U Ω Ω ) + ( F ( r N ) T ( r N ) ) .
(3)

Therefore, our simulation of a suspension under a simple shear flow in the Stokes regime requires solving (3) for the velocities, and obtaining the positions r N of the particles at any time t through time integration of these velocities.

2. Hydrodynamic interactions

In principle, hydrodynamic interactions can be determined by solving the Stokes equations, but this is extremely expensive from a computational view point. For dense suspensions, however, the dominant hydrodynamic interactions come from the fluid flow in the narrow gaps between nearby solid particles [Ball and Melrose (1997)], because the resistance to relative motion becomes singular when the gaps are vanishingly small. They give rise to pair-wise short-range lubrication forces (as opposed to the many-body nature of the full, long-range hydrodynamic interactions). As a consequence, the linear relations between velocities and hydrodynamic forces simply contain a contribution from the Stokes drag and a contribution from the lubrication

( F H T H ) = ( R Stokes + R Lub ) · ( U U Ω Ω ) + R Lub : E ,
(4)

where R Stokes is a diagonal matrix giving Stokes drag forces and torques, and R Lub and R Lub are sparse matrices [Ball and Melrose (1997)].

Consistent with our choice of keeping only physically relevant near-distance interactions, we use only the leading terms for the resistance matrices. Physically, the terms included in our model correspond to the squeeze, shear, and pump modes of Ball and Melrose (1997) (we do not consider their twist mode, as it is not associated with a divergence of the resistance at contact and is thus subdominant). Defining the nondimensional gap h(i,j) between particles i and j having radii ai and aj as h ( i , j ) 2 ( r ( i , j ) a i a j ) / ( a i + a j ) and the center-to-center unit vector n i j ( r ( j ) r ( i ) ) / r ( i , j ) with r ( i , j ) | r ( j ) r ( i ) | , the modes associated with relative displacements, respectively, along and tangential to n i j have singular behaviors with leading terms diverging as 1/h(i,j) and log ( 1 / h ( i , j ) ) [Jeffrey and Onishi (1984); Jeffrey (1992)]. The consequence of the singularity in the squeeze mode (i.e., along n i j ) is that there should not be any contact between particles [Ball and Melrose (1995)]. However, as mentioned in the introduction, this statement is only true for perfectly idealized situations and not realistic even for model experimental systems consisting of spherical particles due to, for example, a finite surface roughness [Davis (1992); Zhao and Davis (2002); Lootens et al. (2005); Blanc et al. (2011)] or a breakdown of the continuity assumption for the fluid at the molecular mean free path scale [Ho and Tai (1998)]. In order to mimic this reality, we regularize the singularities arising in the lubrication resistances by inserting a small cutoff length scale δ [Trulsson et al. (2012)], which can be thought of as the length scale of the particle surface roughness; the leading terms we use for normal and tangential displacements then behave as 1 / ( h ( i , j ) + δ ) and log ( 1 / ( h ( i , j ) + δ ) ) . In the simulations, we set δ = 10−3. This value is typical for non-Brownian particles [Davis (1992); Blanc et al. (2011)]. The detailed expressions for the resistance matrices R Lub and R Lub are given in Appendix  A.

3. Contact model

There are several models that one can use to describe frictional contacts between particles. We use a stick/slide friction model employing springs and dashpots that is commonly used in granular physics [Cundall and Strack (1979); Luding (2008)]. The normal and tangential components of the force and the torque for particles having radii ai and aj are obtained as

F C , nor ( i , j ) = k n h ( i , j ) n i j + γ n U n ( i , j ) , F C , tan ( i , j ) = k t ξ ( i , j ) , T C ( i , j ) = a i n i j × F C , tan ( i , j ) ,
(5)

and fulfill Coulomb's friction law | F C , tan ( i , j ) | μ | F C , nor ( i , j ) | . In the above expressions, kn and kt are the normal and tangential spring constants, respectively, and γn is the damping constant. The normal and tangential velocities are U n ( i , j ) n i j n i j · ( U ( j ) U ( i ) ) and U t ( i , j ) ( I n i j n i j ) · [ U ( j ) U ( i ) ( a i Ω ( i ) + a j Ω ( j ) ) × n i j ] . Finally, the quantity ξ ( i , j ) is the tangential spring stretch. This contact model could be made more general, but at the price of numerical difficulties; see Appendix  B.

The computation of the tangential spring stretch ξ ( i , j ) , described in the following, requires some care, as we have to impose Coulomb's law. We apply an algorithm described in Luding (2008). At the time t0 at which the contact (i,j) is created, we set an unstretched tangential spring ξ ( i , j ) ( t 0 ) = 0 . At any further time step t in the simulation, the tangential stretch ξ ( i , j ) ( t ) is incremented according to the value of a “test” force F C , tan ( i , j ) ( t + d t ) = k t ξ ( i , j ) ( t + d t ) with ξ ( i , j ) ( t + d t ) = ξ ( i , j ) ( t ) + U t ( i , j ) ( t ) d t . If | F C , tan ( i , j ) ( t + d t ) | μ | F C , nor ( i , j ) ( t + d t ) | , the contact is in a static friction state and we update the spring stretch and force as

ξ ( i , j ) ( t + d t ) = ξ ( i , j ) ( t + d t ) , F C , tan ( i , j ) ( t + d t ) = F C , tan ( i , j ) ( t + d t ) .
(6)

However, if | F C , tan ( i , j ) ( t + d t ) | > μ | F C , nor ( i , j ) ( t + d t ) | , the contact is in a sliding state and the spring and force are updated as

ξ ( i , j ) ( t + d t ) = μ k t | F C , nor ( i , j ) ( t + d t ) | t i j , F C , tan ( i , j ) ( t + d t ) = k t ξ ( i , j ) ( t + d t ) ,
(7)

where the direction t i j is the same as the one for the test force, i.e., t i j F C , tan ( i , j ) ( t + d t ) / | F C , tan ( i , j ) ( t + d t ) | . In this case, Coulomb's law is not violated, as | F C , tan ( i , j ) ( t + d t ) | = μ | F C , nor ( i , j ) ( t + d t ) | .

1. Shear rate dependence: Why Stokes hydrodynamics plus hard spheres is not enough

The shear-rate dependence of a suspension requires the existence of a time scale distinct from the inverse shear rate. A well known case is the one of Brownian suspensions, where the inverse shear rate competes with the Brownian diffusion time to give a shear-rate dependent rheology. This competition is adequately measured by a nondimensionalized shear rate, the Péclet number, which is the ratio of the two time scales: P e 6 π η 0 a 3 γ ̇ / k B T . This quantity can also be thought of as a ratio of two force scales: The typical hydrodynamic force and the typical Brownian force. Thus, a shear-rate dependent rheology generically requires at least one other force scale besides the hydrodynamic one.

Such a time scale does not exist for non-Brownian suspensions of hard spheres in the Stokes regime. The reason is straightforward: The hard sphere force has no typical value, and a hard-sphere contact between two particles can withstand any applied load. Therefore, there is nothing with which to compare the hydrodynamic force. This means that the solutions of the force/torque balance equations (3) (i.e., the particle trajectories) are independent of the shear rate, leading to a shear-rate independent rheology. Notice that this conclusion holds whether the hard spheres are frictionless or frictional: Coulomb's law does not introduce a force scale. There must be another force scale in the system to provide the shear-rate dependence.

An important point in our model is that, even if we try to mimic a hard sphere suspension, our contact forces actually are elastic forces, thus bringing a natural force scale kna. We preserve the shear-rate independent hard sphere behavior by selecting the particle stiffness kn such that the nondimensional number based on the stiffness 6 π η 0 a γ ̇ / k n remains much smaller than unity; i.e., we stay in the asymptotic regime of nearly hard particles. We describe the techniques we use to achieve this in our simulation in Appendix  B.

2. A minimal model: Critical-load friction

The first model we consider is an intentionally simplistic one, arguably the simplest model with an additional force scale besides the hydrodynamic force. Simplicity is the main motivation for this model, as it enables the clear physical discussion one can expect from a minimal model.

In this model (which we call critical load model, CLM), we introduce an extra force scale in the friction law itself, namely a threshold normal load FCL to activate friction. When the normal load is smaller than this threshold, particles interact as frictionless hard spheres. When the load is beyond the threshold, friction is activated. Overall, the friction law reads

| F C , tan ( i , j ) | { μ ( | F C , nor ( i , j ) | F CL ) for | F C , nor ( i , j ) | F CL , 0 otherwise .
(8)

We may consider the CLM as the short Debye length limit of the electrostatic repulsion model (ERM) presented below (Sec. II B 3).

3. Electrostatic repulsion model

An electrostatic repulsion is a simple and plausible interaction. Many experimental systems have such a force, as it stabilizes the suspension. Thus, we will consider an ERM including this ingredient in the simplest form of a repulsive double-layer electrostatic force F R ( i , j ) [Israelachvili (2011)]. If the Debye length κ−1 is small compared to the radius of the particles, the approximate form

F R ( i , j ) ( h ) { 2 F ER a i a j / a ( a i + a j ) e κ ( r ( i , j ) a i a j ) n i j if h ( i , j ) 0 , 2 F ER a i a j / a ( a i + a j ) n i j if h ( i , j ) < 0 ,
(9)

can be used. In the simulations, we set κ−1 = 0.05a, meaning a Debye length in the 0.1 μm range for particles of a few μm; see Mewis and Wagner (2011).

The mechanical stress applied to the suspension arises from the several interactions included in the model: Particles disturb the flow, creating hydrodynamic stresses, and they develop force chains via contacts (and/or electrostatic repulsion for the ERM).

The hydrodynamic stresslets acting on the particles are given by [Batchelor (1970); Brady and Bossis (1988)]

S H = ( R Stokes S + R Lub S ) · ( U U Ω Ω ) + R Lub S : E ,
(10)

where S H ( S H ( 1 ) , , S H ( n ) ) and the matrices R Lub S and R Lub S contain leading terms of the lubrication resistances [Jeffrey and Onishi (1984); Jeffrey (1992)] (see Appendix  A) in a manner consistent with the hydrodynamic forces considered in Sec. II A 2.

The stress due to the contact or repulsive force between particles i and j is simply written as

S C ( i , j ) = ( r ( j ) r ( i ) ) F C ( i , j ) or S R ( i , j ) = ( r ( j ) r ( i ) ) F R ( i , j ) ,
(11)

respectively.

The bulk stress, in which the isotropic part of the fluid pressure is omitted, is the sum of the above contributions

Σ 2 η 0 E + 1 V ( i S H ( i ) + i > j S C ( i , j ) + i > j S R ( i , j ) ) ,
(12)

where V is the volume of the simulation box (the electrostatic stress S R appears only for the ERM). Shear stress σ, normal stress differences N1 and N2, and particle pressure Π [Yurkovetsky and Morris (2008)] are defined as σ Σ 12 , N 1 Σ 11 Σ 22 , N 2 Σ 22 Σ 33 , and Π ( Σ 11 + Σ 22 + Σ 33 ) / 3 , respectively. The relative viscosity is given by η r σ / η 0 γ ̇ .

  • Bidispersity: We consider a suspension of bidisperse frictional hard spheres immersed in a Newtonian fluid with viscosity η0. The bidispersity is introduced to hinder the formation of the ordered phase observed for dense monodisperse suspensions under shear. An effective choice for the size ratio is a2/a1 = 1.4 (where a1 = a), with the two populations occupying equal volumes; i.e., ϕ 1 = ϕ 2 . Other effective size ratios are possible, as is polydispersity, with the same qualitative effects. The choice of size ratio 1:1.4 is however the most common one in the granular matter literature, where bidispersity is used to avoid crystallization [see for example O'Hern et al. (2003)]. Indeed, with a smaller size ratio 1:1.2, slow and weak strain thinning due to ordering is seen in the low viscosity state.

  • Dimensionless shear rate: We discussed in Sec. II B 1 that another force scale besides the hydrodynamic one is required to yield shear-rate dependence. In CLM, the threshold value gives the force scale F = F CL , and in ERM, the force at contact gives F = F ER . Therefore, the shear rate dependence is given by the ratio γ ̇ / γ ̇ 0 , with γ ̇ 0 F / 6 π η 0 a 2 .

  • Periodic boundary conditions and fixed-volume simulation: Rheology is a bulk property. If solid walls are used for the boundaries of the system, we need a large system size to reduce the influence of walls. We can avoid the use of solid walls in a simulation by using periodic boundary conditions. For the strain-controlled simple shear, we use the Lees–Edwards [Lees and Edwards (1972)] boundary condition.

Particle migration never develops under these periodic boundary conditions. Although some experimental observations suggest that global migration may be a cause of shear thickening [Fall et al. (2010)], our simulation is free of this effect.

There is no way for the system to dilate with these periodic boundary conditions. Such a fixed volume condition is expected in most rheology measurements. For shear thickening fluids, however, the open edges may have some influence [Brown and Jaeger (2012)]. In granular physics, systems are often sheared under a given normal stress, hence the volume is not fixed [Boyer et al. (2011)].

This section presents the results from the simulation of the two models, the “minimal” CLM and the ERM. In Secs. III A–III E, whenever possible, the two models are treated at the same time in the text, and we will show data plots for the CLM. However, when the two models give slightly different results, the results of the CLM will be described first, as it allows a simple and clear understanding of the underlying physics, and the more realistic model ERM will be described afterward. The differences are essentially in the low shear-rate limit for the ERM (existence of a shear-thinning regime), due to the repulsive potential.

Although friction is the key factor in this work, few experimental estimates are available for the friction coefficient. We select a friction coefficient μ = 1 [which is comparable to the measurements of Fernandez et al. (2013) for ∼10 μm quartz particles coated by polymer brushes] for most of the simulations, except when the dependence on μ is specifically investigated.

In the data plots shown in this section, the error bars represent the standard deviation, which we define for an observable A(t) as A 2 A 2 , where · is the time (strain) average T 1 0 T · d t over the simulated units of strain. (We perform simulations over 50 strain units, hence T = 50 except when the time averages are performed over subsets of the whole simulation in the case of intermittent data around DST.)

Rheological data are plotted versus shear rates and shear stresses nondimensionalized by the characteristic shear rate γ ̇ 0 = F / 6 π η 0 a 2 and stress η 0 γ ̇ 0 , respectively, where η0 is the viscosity of the suspending fluid.

In the CLM, due to the threshold force, the friction between grains is absent at low shear rates and activated at high shear rates. Because of this, we expect that the low shear-rate limit for concentrated suspensions will have a rheology η ( ϕ , γ ̇ 0 ) typical of a system close to the jamming transition for frictionless particles, while the high shear-rate limit shows a rheology η ( ϕ , γ ̇ ) typical of a system close to the jamming transition for particles with a friction coefficient μ = 1.

These two limiting viscosities are shown in Fig. 1, where we also show the high shear-rate behavior for the infinite friction case (μ = ) for reference. Each diverges at a different volume fraction, thus friction shifts the jamming point [Otsuki and Hayakawa (2011)]. We fit our data with power-law divergences η C ( 1 ϕ / ϕ J ) λ , with parameters ( ϕ J , λ , C ) as detailed in the caption of Fig. 1.

FIG. 1.

Relative viscosity ηr as a function of the volume fraction ϕ in the two limits γ ̇ 0 and γ ̇ (left). The γ ̇ 0 viscosity (circles) is independent of the friction coefficient μ as the friction is not activated at low stresses, which leads to a relatively lower viscosity diverging at a higher volume fraction ϕ J 0 (which is the jamming point for frictionless systems). The γ ̇ viscosity however directly depends on μ, as is seen from the difference between μ = 1 (squares) and μ =  (diamonds) plots. In particular, the jamming volume fraction decreases with increasing μ. We fit our data with power laws η r = C ( 1 ϕ / ϕ J ) λ (right). The best fitting parameters are ( ϕ J 0 , λ 0 , C 0 ) ( 0.66 , 1.6 , 1.40 ) , ( ϕ J μ = 1 , λ μ = 1 , C μ = 1 ) ( 0.58 , 2.3 , 0.71 ) , and ( ϕ J μ = , λ μ = , C μ = ) ( 0.56 , 2.4 , 0.63 ) .

FIG. 1.

Relative viscosity ηr as a function of the volume fraction ϕ in the two limits γ ̇ 0 and γ ̇ (left). The γ ̇ 0 viscosity (circles) is independent of the friction coefficient μ as the friction is not activated at low stresses, which leads to a relatively lower viscosity diverging at a higher volume fraction ϕ J 0 (which is the jamming point for frictionless systems). The γ ̇ viscosity however directly depends on μ, as is seen from the difference between μ = 1 (squares) and μ =  (diamonds) plots. In particular, the jamming volume fraction decreases with increasing μ. We fit our data with power laws η r = C ( 1 ϕ / ϕ J ) λ (right). The best fitting parameters are ( ϕ J 0 , λ 0 , C 0 ) ( 0.66 , 1.6 , 1.40 ) , ( ϕ J μ = 1 , λ μ = 1 , C μ = 1 ) ( 0.58 , 2.3 , 0.71 ) , and ( ϕ J μ = , λ μ = , C μ = ) ( 0.56 , 2.4 , 0.63 ) .

Close modal

For the ERM, the situation is the same at high shear rates but differs at low shear rates. While friction is not felt for γ ̇ 0 , as particles do not contact, the finite range of the repulsive potential creates a shear-thinning behavior from which we could not obtain the low shear-rate limit η ( ϕ , γ ̇ 0 ) . The shear thinning comes from the fact that the system behaves essentially as soft particles at low shear stresses, with an apparent size that includes the hard sphere and a part of the surrounding soft repulsive potential. Indeed, a simple force balance argument at the scale of a particle says that if a particle is subject to a driving shear force σa2, the minimum gap hmin between this particle and its neighbors will be such that F R ( h min ) σ a 2 [with the limitation that it must respect the geometrical constraint, i.e., h min < ( 1 ϕ / ϕ J 0 ) 1 / 3 ]. Then, at this shear stress the minimum center-to-center distance between two particles is 2 a ( 1 + h min ( σ ) / 2 ) , which means that the particles have an apparent radius larger than the one of their hard core. The jamming transition of these effectively bigger frictionless particles is at an apparent packing fraction (that is, based on the apparent diameter) ϕ = ϕ J 0 0.66 , which corresponds to an actual packing fraction ϕ significantly lower than ϕ J 0 . Because of this low volume fraction jamming transition at γ ̇ = 0 , the viscosity is larger than at finite (but small) γ ̇ , leading to the shear thinning we observe. Such thinning is actually known to occur due to repulsive forces in charge stabilized suspensions [Krieger (1972); Maranzano and Wagner (2001a)].

We can switch from one rheology to the other by varying the shear rate. Physically, the transmitted stress increases as the shear rate increases, which triggers the formation of frictional contacts between particles. Thus, by increasing the shear rate, the viscosity interpolates between the frictionless and frictional rheology curves, which means we can observe shear thickening. All this should be a natural consequence of the existence of two distinct rheologies at γ ̇ = 0 and γ ̇ = . What we cannot anticipate a priori is the way in which the system switches from the low viscosity state to the high viscosity one: Do we observe a continuous shear thickening (CST) or a DST?

The shear rate dependence of the viscosity, shown in Fig. 2 for the CLM, demonstrates the existence of both CST and DST in our system, depending on the volume fraction. As in experiments, when ϕ < ϕ c the shear thickening is continuous, getting steeper and steeper as we approach ϕ c , at which point it becomes discontinuous and keeps this behavior for ϕ > ϕ c and up to ϕ J 0 . Note that there appears to be a real discontinuity in our data for these volume fractions: The time series of the viscosity show an intermittent behavior switching between two states, which we address in Sec. III C. In Fig. 2, the intermittent data are split between low and high viscosity states: Two points that correspond to a separate time average for each of the two states appear at the same shear rate.

FIG. 2.

Shear rate γ ̇ (left) and shear stress σ (right) dependences of the relative viscosity ηr for the CLM, with friction coefficient μ = 1, for volume fractions 0.45 ϕ 0.56 . The system size is N = 1000.

FIG. 2.

Shear rate γ ̇ (left) and shear stress σ (right) dependences of the relative viscosity ηr for the CLM, with friction coefficient μ = 1, for volume fractions 0.45 ϕ 0.56 . The system size is N = 1000.

Close modal

When plotted against stress in Fig. 2, the viscosity curves show another interesting feature, namely that the onset of shear thickening occurs at a stress σ ST / ( η 0 γ ̇ 0 ) 0.3 that is roughly independent of the volume fraction, as is observed in many experiments [Frith et al. (1996); Bender and Wagner (1996); Maranzano and Wagner (2001b,a); Lootens et al. (2005); Fall et al. (2010); Larsen et al. (2010); Brown and Jaeger (2012, 2014)]. Similarly, the transition toward the shear-rate independent plateau at high viscosity also occurs at a stress independent of ϕ . The stress range over which thickening occurs remains constant from a mild shear thickening to a marked DST, as expected from a simple balance argument between the driving stress σ = η γ ̇ and the stress required to create a frictional contact F CL / a 2 (for the CLM) with the area a2 meaning that we require almost every particle to have a contact. It implies the thickening should roughly take place around σ F CL / a 2 = 6 π η 0 γ ̇ 0 , which gives σ / ( η 0 γ ̇ 0 ) 6 π . This value falls at the upper end of the thickening regime as shown in Fig. 2, which is consistent with our assumption that at this stress scale every particle should have contacts with its neighbors, hence the system is in the frictional rheology state. The inverse argument applied to the onset stress σST finds that the area AST over which one unit of critical load applies is σST = FCL/AST, with A ST ( 7 a ) 2 . That means that at the onset of shear thickening, the contact chains must be sparsely distributed.

Friction is an essential ingredient for the shear thickening to occur, as is illustrated in Fig. 3, which shows the effect of a reduction of the friction coefficient on the η r ( γ ̇ ) curve in the CST regime. If the friction is removed (μ = 0), the effect is completely suppressed (which is expected, as the threshold force scale FCL disappears from the equations of motion for μ = 0, and the simulation is then rigorously the same for all γ ̇ ). If the friction is only reduced (μ = 0.2), the shear thickening is milder than with μ = 1 friction at the same volume fraction. In that case, the CST however becomes more pronounced and turns to DST when one increases the volume fraction, just as for the μ = 1 case.

FIG. 3.

Effect of the friction coefficient μ on the shear rate dependence of the relative viscosity η r ( γ ̇ ) in the CLM. The friction is essential to the shear thickening, as is illustrated by the reduction of the effect for μ = 0.2 and its complete suppression for μ = 0.

FIG. 3.

Effect of the friction coefficient μ on the shear rate dependence of the relative viscosity η r ( γ ̇ ) in the CLM. The friction is essential to the shear thickening, as is illustrated by the reduction of the effect for μ = 0.2 and its complete suppression for μ = 0.

Close modal

The ERM shown in Fig. 4 behaves in a similar manner, only adding a shear-thinning regime at low shear rates as discussed in the last paragraph of Sec. III A.

FIG. 4.

Shear rate γ ̇ (left) and shear stress σ (right) dependences of the relative viscosity ηr for the ERM, with friction coefficient μ = 1, for volume fractions 0.48 ϕ 0.57 . The system size is N = 512.

FIG. 4.

Shear rate γ ̇ (left) and shear stress σ (right) dependences of the relative viscosity ηr for the ERM, with friction coefficient μ = 1, for volume fractions 0.48 ϕ 0.57 . The system size is N = 512.

Close modal

The existence of two distinct flowing states suggests that one might expect to observe a hysteresis associated with DST. It is thus remarkable that, aside from suspensions that are hysteretic at the microscopic level (with attractive interactions, for example, since once particles are aggregated under flow they stay aggregated upon cessation of the flow), there are very few experimental reports of hysteresis in systems showing DST [one of the best examples is provided by Bender and Wagner (1996)].

What is commonly observed, however, is a switching behavior, where the time series of the stress shows that the system shares its time between the two states, erratically switching in what looks like activated events [Boersma et al. (1991); D'Haene et al. (1993); Bender and Wagner (1996); Lootens et al. (2003); Lootens et al. (2004); Lootens et al. (2005); Hébraud and Lootens (2005)]. We observe this switching behavior in both models that we studied. In Fig. 5, we show time series of the stress and the corresponding histograms around the DST. These data show a typical “coexistence” behavior: Significantly below and above shear thickening, the system, respectively, stays in the low and high viscosity state, but close to DST, the two states coexist in the same time series, the proportion of time spent in the high viscosity state increasing with γ ̇ . While it is tempting to use the analogy of DST with an equilibrium first-order transition to conclude that the switchings are finite-size “activated” events, we need to be careful in doing so, because intermittency seems to survive for unusually large numbers of particles. The intermittency is still present in simulations with N = 1000 particles, and the switching frequency does not seem to decrease relative to simulations with N = 512 particles, whereas in a first-order thermodynamic transition, one would expect the switching rate to be drastically reduced by a doubling of the system size. In fact, this effect might persist for much larger particle numbers, as is suggested by the experimentally observed intermittency [Boersma et al. (1991); D'Haene et al. (1993); Bender and Wagner (1996); Lootens et al. (2003); Lootens et al. (2004); Lootens et al. (2005); Hébraud and Lootens (2005)] in systems where N is considerably larger. The simulations of Heussinger (2013) also show intermittency for systems of several thousands of particles. Moreover, should this phenomenon disappear in the thermodynamic limit, the fact that this limit is not observed until there are very large numbers of particles might be valuable information about the nature of this out-of-equilibrium phase transition, just as strong finite-size effects have a profound physical significance near a critical point.

FIG. 5.

Strain (or dimensionless time) series of the shear stress close to the DST (left), and the corresponding histograms (right). The simulations are performed at ϕ = 0.57 with the ERM, with shear rates γ ̇ / γ ̇ 0 = 0.0094 (just below the DST region), γ ̇ / γ ̇ 0 = 0.0096 (in the DST region), and γ ̇ / γ ̇ 0 = 0.0098 (just above the DST region). Data at γ ̇ / γ ̇ 0 = 0.0096 show a superposition of the two states. This is very similar to what is seen in experiments [Boersma et al. (1991)].

FIG. 5.

Strain (or dimensionless time) series of the shear stress close to the DST (left), and the corresponding histograms (right). The simulations are performed at ϕ = 0.57 with the ERM, with shear rates γ ̇ / γ ̇ 0 = 0.0094 (just below the DST region), γ ̇ / γ ̇ 0 = 0.0096 (in the DST region), and γ ̇ / γ ̇ 0 = 0.0098 (just above the DST region). Data at γ ̇ / γ ̇ 0 = 0.0096 show a superposition of the two states. This is very similar to what is seen in experiments [Boersma et al. (1991)].

Close modal

In order to observe hysteresis, one must have a system where the relaxation time τrelax (the length of the transient) is much smaller than the activation time τact (the typical time between two switches). One then does a proper measurement with a shear rate sweep on a time scale τsw such that τ relax τ sw τ act : The first inequality ensures that the system is always in a steady state during the sweep and that the measurement is done with sufficient time averaging, while the second inequality ensures that we are not averaging data over two distinct states. In our simulations, as can be seen by looking at the time series in Fig. 5, a clear separation of time scales is not achieved (we have τ relax τ act , which does not permit finding a proper τsw) and hysteresis cannot be observed unambiguously. Many experimental data actually suffer from the same limitations, and this might be the reason for the very small number of hysteretic flow curves actually reported. One exception is the noted work of Bender and Wagner (1996).

One common rheological characteristic of dense suspensions is the appearance of finite normal stress differences. Our measurements of the two normal stress differences N1 and N2 are shown in Figs. 7 and 6 for the CLM. The ERM behaves very similarly and is not shown.

FIG. 7.

Shear rate dependence of the normalized first normal stress difference N 1 / η 0 γ ̇ (left) and shear stress dependence of N1σ (right) for the CLM.

FIG. 7.

Shear rate dependence of the normalized first normal stress difference N 1 / η 0 γ ̇ (left) and shear stress dependence of N1σ (right) for the CLM.

Close modal
FIG. 6.

Shear rate dependence of the normalized second normal stress difference N 2 / η 0 γ ̇ (left) and shear stress dependence of N2σ (right) for the CLM. The stress σST at which thickening starts is marked by an arrow.

FIG. 6.

Shear rate dependence of the normalized second normal stress difference N 2 / η 0 γ ̇ (left) and shear stress dependence of N2σ (right) for the CLM. The stress σST at which thickening starts is marked by an arrow.

Close modal

We obtain for N2 a behavior consistent with most experimental data available: It is negative, comparable to the shear stress (and much larger than N1 that we discuss in the next paragraph), and its behavior is reminiscent of that of the shear stress. This comes from the fact that most of the stress is transmitted by forces in the shear plane (flow-gradient), not along the vorticity direction. Indeed the ratio N2σ is rather insensitive to the stress magnitude, increasing at most by a factor of two across shear thickening while the stresses increase by almost two orders of magnitude. What is remarkable is the independence of N2σ on the volume fraction: All volume fractions investigated here collapse on a master curve for both models.

There has been debate in the community regarding the value (and even the algebraic sign) of N1. The first notable feature of our data concerning N1 is that it is dominated by the fluctuations. In one time series, a cursory look indicates that its value fluctuates around zero. Long time averaging reveals more structure, however, as shown in Fig. 7 for the CLM. Prior to shear thickening, the average value of N1 is nearly zero (slightly negative for CLM, slightly positive for ERM). The shear thickening transition is marked by two different behaviors, depending on the volume fraction: At the lowest volume fractions studied, N1 decreases across shear thickening, while at larger ϕ there is a clear upturn toward positive values at the shear thickening transition. For the CLM, the behavior can be systematized by plotting N1 as a function of the stress, as in Fig. 7: N1σ is an increasing function of σ, negative at low stresses, with a qualitative behavior roughly independent of ϕ . The same plot for the ERM is similar but less systematic. In any case, even above shear thickening, the ratio N1σ is always small, never exceeding 0.1.

Overall, the behavior we observe for N1 is in reasonable agreement with most experiments [Lootens et al. (2005); Lee et al. (2006); Larsen et al. (2010); Couturier et al. (2011); Dbouk et al. (2013)]. Zarraga et al. (2000); Singh and Nott (2003); Dai et al. (2013) report a negative value for N1. Lee et al. (2006) agree with our trend at low shear rates but observe a subsequent change toward negative N1 at very high shear rates, while Dbouk et al. (2013) find a positive but significantly larger N1. In any case, the amplitude of the fluctuations seems to be the relevant physical information about N1, as the average, positive, or negative, is buried in the fluctuations in the time series, such that at any time N1 can be either positive or negative, even at large ϕ and γ ̇ .

Non-Newtonian behavior of Brownian suspensions is often discussed with shear induced microstructure; the particle configuration is nearly at equilibrium at low Péclet number, while an anisotropic microstructure is induced by shear at high Péclet number [Morris (2009)]. Ideas associated with the order-disorder transition scenario [Hoffman (1972, 1974, 1998)] also predicted a clear structural change with shear thickening. Experimental observations, however, revealed that the low viscosity state was not always ordered and that the shear thickening can occur with only a subtle signature in the microstructure [Bender and Wagner (1996); Watanabe et al. (1997); Watanabe et al. (1998); Newstein et al. (1999); Maranzano and Wagner (2002)]. One suggestion comes from the hydro-clustering scenario, which assumes the existence of hydrodynamically created extended density fluctuations. Some experimental data are indeed in qualitative agreement with this scenario [Maranzano and Wagner (2002); Cheng et al. (2011)].

In the scenario we suggest in this article, the main structural modifications across the shear thickening transition should be sought in the contacts between particles. At low shear rates, frictional contacts are avoided because the particle pressure is too small to overcome the threshold (for the CLM) or the repulsion between particles (for the ERM). Those changes should thus be detected by measures specifically sensitive to the contacts. The pair correlation function should show few dramatic modifications across shear thickening. This is what we show in this section.

We define the pair correlation function

g all ( r ) V N 2 i , j δ ( r r i j ) .
(13)

The system being bidisperse, we can also define partial pair correlation functions restricted to, e.g., pairs of small-small particles g SS ( r ) ( V / N S 2 ) i , j S S δ ( r r i j ) where S S is the subset of NS small particles. In the same way, we define the structure factor (which is related to the pair correlation function via a Fourier Transform) as

S all ( k ) 1 N i , j S , i j e i k · r i j .
(14)

In Figs. 8 and 9, we show the pair correlation function restricted to the shear plane (velocity/gradient) at ϕ = 0.57 for four values of the shear rate: γ ̇ / γ ̇ 0 = 0.005 , well below DST; γ ̇ / γ ̇ 0 = 0.015 , just below DST; γ ̇ / γ ̇ 0 = 0.02 , just above DST and γ ̇ / γ ̇ 0 = 0.05 , well above DST. As expected, the evolution of those plots seems weak around DST. The main feature is a loss of contrast, unveiling a more isotropic structure above shear thickening.

FIG. 8.

Pair correlation functions restricted in the shear plane g all ( r , θ ) for the CLM at ϕ = 0.56 . The four plots correspond to shear rates γ ̇ = 0.005 (well below DST), γ ̇ = 0.015 (just below DST), γ ̇ = 0.02 (just above DST), and γ ̇ = 0.05 (well above DST).

FIG. 8.

Pair correlation functions restricted in the shear plane g all ( r , θ ) for the CLM at ϕ = 0.56 . The four plots correspond to shear rates γ ̇ = 0.005 (well below DST), γ ̇ = 0.015 (just below DST), γ ̇ = 0.02 (just above DST), and γ ̇ = 0.05 (well above DST).

Close modal
FIG. 9.

Pair correlation functions restricted to pairs of small particles in the shear plane g SS ( r , θ ) for the CLM at ϕ = 0.56 . The two plots correspond to shear rates γ ̇ = 0.015 (just below DST) and γ ̇ = 0.02 (just above DST).

FIG. 9.

Pair correlation functions restricted to pairs of small particles in the shear plane g SS ( r , θ ) for the CLM at ϕ = 0.56 . The two plots correspond to shear rates γ ̇ = 0.015 (just below DST) and γ ̇ = 0.02 (just above DST).

Close modal

The structure factor shown in Fig. 10 reveals another interesting feature: There is a clear anisotropy in the shear plane in the low viscosity state, with a peak in the gradient direction that is strongly reduced (but not absent, see the right of Fig. 10) in the high viscosity state, above DST. This anisotropy is absent in the flow-vorticity plane; see Fig. 11. Looking back at the pair correlation functions in Fig. 8, we see that the peak of S all ( k ) along the gradient direction is the signature of a slight ordering of the low viscosity phase: There are small but visible downstream ripples aligned with the flow direction at small shear rate.

FIG. 10.

Structure factor in the shear plane S all ( k 1 , k 2 ) for the CLM at ϕ = 0.56 . Left: γ ̇ = 0.015 (just below DST). Right: γ ̇ = 0.02 (just above DST). The low viscosity state shows a peak in the gradient direction (around k 1 2 π / d ) associated with a small amount of layering, which is strongly attenuated in the high viscosity phase.

FIG. 10.

Structure factor in the shear plane S all ( k 1 , k 2 ) for the CLM at ϕ = 0.56 . Left: γ ̇ = 0.015 (just below DST). Right: γ ̇ = 0.02 (just above DST). The low viscosity state shows a peak in the gradient direction (around k 1 2 π / d ) associated with a small amount of layering, which is strongly attenuated in the high viscosity phase.

Close modal
FIG. 11.

Structure factor in the flow-vorticity plane S all ( k 1 , k 3 ) for the CLM at ϕ = 0.56 . Left: γ ̇ = 0.015 (just below DST). Right: γ ̇ = 0.02 (just above DST). In this plane, the structure is isotropic in both states.

FIG. 11.

Structure factor in the flow-vorticity plane S all ( k 1 , k 3 ) for the CLM at ϕ = 0.56 . Left: γ ̇ = 0.015 (just below DST). Right: γ ̇ = 0.02 (just above DST). In this plane, the structure is isotropic in both states.

Close modal

The slight layering we observe is however quite far from the string order usually associated with the order-disorder scenario. In supplementary material, we show movies of the system projected on the gradient-vorticity plane (i.e., looking down the stream direction) for shear rates just below and just above shear thickening, where we reduced the size of the particles to a third of their actual radius to allow a better visualization. Those movies show warped layers along the flow-vorticity plane in some parts of the system, while in other parts ordered layers are completely absent. In any case, there is no further string ordering within the layers. This is consistent with our structure factor data in the gradient-vorticity plane (not shown), which exhibits only some structure in the gradient direction and no six-fold symmetry patterns typically observed when string ordering takes place [Laun et al. (1992); Kulkarni and Morris (2009)].

In our simulation, the shear thickening transition is due to the appearance of a growing number of frictional contacts as the shear rate (or more precisely the shear stress) increases. The contact network that is built across shear thickening is shown in Fig. 12 for the CLM for the two cases of CST and DST, respectively. (Two movies showing γ ̇ / γ ̇ 0 = 0.015 and 0.018 of DST are available in supplementary material. The time in these movies is scaled with 1 / γ ̇ .)

FIG. 12.

Snapshots of the contact network for the CLM at ϕ = 0.5 (top) and 0.56 (bottom). Frictionless contacts (with normal force below the threshold FCL) are drawn in gray segments joining the centers of the two involved particles, while frictional contacts are drawn in red. Each row corresponds to a single shear rate, ranging from the low viscosity state to the high viscosity state and across CST (top) and DST (bottom). For each shear rate, we show four snapshots at different times (or equivalently strains) along the same simulation.

FIG. 12.

Snapshots of the contact network for the CLM at ϕ = 0.5 (top) and 0.56 (bottom). Frictionless contacts (with normal force below the threshold FCL) are drawn in gray segments joining the centers of the two involved particles, while frictional contacts are drawn in red. Each row corresponds to a single shear rate, ranging from the low viscosity state to the high viscosity state and across CST (top) and DST (bottom). For each shear rate, we show four snapshots at different times (or equivalently strains) along the same simulation.

Close modal

At low shear rate, in the low viscosity state, frictional contacts only seldomly appear, being concentrated in small force chains along the compressional axis. At high shear rate, frictional contacts are the norm, creating a frictional contact network that is very close to being jammed; i.e., having only a few (collective) degrees of freedom left to reorganize under the applied stress.

Between those two extreme situations, a whole continuum of gradually denser frictional contact networks is seen across the CST transition in Fig. 12 (top). But in the case of DST in Fig. 12 (bottom), the contact network discontinuously changes from sparse to dense at the transition, never showing configurations of intermediate densities. Another interesting point is the occurrence of intermittency in the contact network, which is the immediate structural origin of the intermittent stress behavior shown in Fig. 5: Close to DST, the contact network suddenly switches from mostly frictionless to mostly frictional, showing a large sensitivity to fluctuations. When this happens, the viscosity immediately follows by switching to high/low viscosity when the contact network, respectively, densifies/loosens.

We can make the link between frictional contacts and shear stress explicit by looking at the fraction of frictional contacts f [Wyart and Cates (2014)]. This quantity is unambiguously defined in the CLM model as the ratio of the number of contacts N C f that are in the frictional state (those for which the normal force exceeds the critical load FCL) divided by the total number of contacts NC. This ratio has a direct relation to the shear stress, as the function f(σ) plotted in Fig. 13 shows. This single relation demonstrates that, more than structural aspects, the shear thickening is above all a manifestation of the proliferation of frictional contacts in the system. We note one remarkable aspect of the relation f(σ): It is independent of the volume fraction, at least in the range of ϕ that we have studied (which is rather large, owing to the fact that it should be understood as a range of distance to jamming ϕ ϕ J μ rather than a bare range of ϕ ). This peculiar aspect will be studied in a separate article.

FIG. 13.

The fraction of frictional contacts f as a function of the shear stress σ, for several volume fractions. There is a direct relation between the stress and the proportion of frictional contacts, which emphasizes the fundamental role of friction in shear thickening.

FIG. 13.

The fraction of frictional contacts f as a function of the shear stress σ, for several volume fractions. There is a direct relation between the stress and the proportion of frictional contacts, which emphasizes the fundamental role of friction in shear thickening.

Close modal

While the jamming transition is the critical phenomenon underlying DST, it should be noted that the jamming transition and DST are distinct, and in particular ϕ c ϕ J μ . According to a recent theoretical argument by Wyart and Cates (2014), a scenario based on two diverging rheologies like the one we present in this work implies under reasonable assumptions that ϕ c < ϕ J μ ; i.e., DST could then occur between two flowable (unjammed) states for volume fractions ϕ c < ϕ < ϕ J μ .

In our simulations, we estimate ϕ J μ = 1 0.58 (see Fig. 1), while we observe DST for ϕ = 0.56 (see Fig. 2), which indeed implies ϕ c < ϕ J μ = 1 . For a volume fraction ϕ > ϕ J μ = 1 only the low viscosity state would be visible under shear, because the high viscosity frictional system is jammed at this volume fraction. This sets a maximum shear rate for the shear of the system at this volume fraction.

Physically, the two qualitatively different behaviors CST and DST, and the fact that ϕ c < ϕ J μ = 1 can be explained in our model as stemming from the level of contrast that exists between frictionless and frictional rheologies, as suggested by Wyart and Cates (2014) [part of this explanation, when the contrast is diverging, was also suggested by Nakanishi et al. (2012)]. Recalling that the fraction of frictional contacts f(σ) essentially depends only on the applied shear stress (see the previous Sec. IV A), and noting that the viscosity interpolates between the two rheologies according to the number of frictional contacts, we see that we can write a direct relation η(σ) between the viscosity and the applied stress. If the difference between the two rheologies is small—i.e., η ( σ ) η ( σ 0 ) is small—the curve η(σ) is a sufficiently mildly increasing function such that γ ̇ = σ / η ( σ ) is also an increasing function of σ; this corresponds to a single-valued curve σ ( γ ̇ ) , and thus to CST. But if the difference η ( σ ) η ( σ 0 ) becomes too large, η(σ) increases faster than σ in some interval, which means that γ ̇ ( σ ) is a decreasing function in this same interval. This corresponds to a multivalued curve σ ( γ ̇ ) , which is unstable, and shows up experimentally as a discontinuous curve, i.e., DST.

An increasing difference η ( σ ) η ( σ 0 ) is naturally provided in our model by the fact that the frictional rheology diverges at a jamming volume fraction ϕ J μ that is smaller than the frictionless rheology jamming point ϕ J 0 . Then, at low volume fraction, the difference is small, but there must be a point ϕ c below ϕ J μ where the difference becomes large enough to observe a DST.

These ideas can be summed up in a phase diagram. The rheology is essentially frictionless in the lower part of the diagram Fig. 14, as the stress is too small to activate friction between grains. This rheology diverges at the frictionless jamming point ϕ J 0 . The rheology is frictional in the upper part of the diagram, as friction is activated under the applied stress. This rheology diverges at the frictional jamming point ϕ J μ and thus shows a larger viscosity than the frictionless rheology. Thus two rheologies coexist on the shear rate-volume fraction plane. Those rheologies are separated by a shear thickening that is continuous for ϕ < ϕ c and discontinuous for ϕ c < ϕ < ϕ J μ .

FIG. 14.

Schematic phase diagram in the shear rate-volume fraction plane. The viscosity is color coded, from small (dark) to large values (white). Along the ϕ axis, ϕ c is the point above which we can observe DST, ϕ J μ is the jamming point for frictional spheres with friction coefficient μ, and ϕ J 0 is the jamming point for frictionless spheres. Intermittency is observed in the region delimited by a dashed line.

FIG. 14.

Schematic phase diagram in the shear rate-volume fraction plane. The viscosity is color coded, from small (dark) to large values (white). Along the ϕ axis, ϕ c is the point above which we can observe DST, ϕ J μ is the jamming point for frictional spheres with friction coefficient μ, and ϕ J 0 is the jamming point for frictionless spheres. Intermittency is observed in the region delimited by a dashed line.

Close modal

The discontinuity is actually related to the coexistence of the two rheologies in the triangle delimited by a dashed line. In this region, we observe one consequence of the coexistence in the intermittency of the flow: The system switches from one state to the other through activation events, as is shown by the time series of the stress (see Fig. 5). For ϕ J μ < ϕ < ϕ J 0 , DST is strictly speaking no longer observed, as the high viscosity flowable state does not exist any more. In this region, DST is actually replaced by shear jamming [Bi et al. (2011)]: If one applies a shear stress larger than σST, the system goes to a solid state and cannot flow. This implies the existence of a forbidden region, in gray, where no flow is possible for hard spheres. In the lower shear-rate domain above ϕ J μ , the low viscosity state is stable, but the activated events responsible for the intermittency at lower ϕ in small systems might lead to complete jamming in a finite time. Lastly, the dot in the upper left corner of the DST region is a critical point separating CST from DST. It shares some features with a critical point of a second order phase transition, for instance, a diverging susceptibility d σ / d γ ̇ .

Of course, as we numerically work with a shear-rate controlled scheme, the system always flows at any shear rate and any volume fraction, even in the forbidden region of the phase diagram. But in the latter region flow is only achieved by creating large overlaps between the particles (i.e., compressing them), hence violating the criteria we set to mimic hard sphere suspensions. At high shear rate and for ϕ > ϕ J μ , a real hard sphere system would not flow, but it does in the simulation because of our inability to enforce the hard sphere condition in a high (possibly infinite) stress state. The spheres we simulate in that case cannot be considered hard any more, and their stiffness sets a cutoff to the stress scales under shear. This situation is somehow similar to the one observed in experiments, where it is observed that the measured stress in the high viscosity state is set by the weakest stress scale in the sample (which might be a large stress scale for the rheometer usage range), whether it is the particles' stiffness or the surface tension on the edge of the rheometer [Brown and Jaeger (2012)]. Therefore, in an experiment, even if the idealized system would be jammed under the investigated conditions, the real system still flows, perhaps thanks to dilation, whereas in our simulation, under the same conditions, it still flows thanks to the deformability of particles.

The difference between ϕ c and ϕ J μ may have been observed in some stress drop experiments [O'Brien and Mackay (2000); Larsen et al. (2010)]. For ϕ c < ϕ < ϕ J μ , the phase diagram of Fig. 14 predicts that the stress in the thickened state would relax quickly upon flow cessation, as the contact network is unjammed. By contrast, for ϕ J μ < ϕ < ϕ J 0 the stress would not relax entirely, as part of it would be stored elastically in the jammed contact network. This scenario is actually consistent with the data from [O'Brien and Mackay (2000); Larsen et al. (2010)].

We have introduced a frictional-viscous model of dense suspensions under shear that is built on the framework used by previous standard models of suspensions: Simulations are performed in the zero particle-inertia and zero fluid-inertia limits ( S t 0 and R e 0 ) and include relevant hydrodynamic interactions and a short-range repulsive potential. The critical innovation is that, besides this purely hydrodynamic basis, it includes frictional contacts between solid particles. This model can be used to simulate sheared suspensions over a range of volume fractions, from mildly dense, where the short-range hydrodynamic interactions dominate, to the jamming point, where contact forces dominate.

The effect of adding frictional contacts is striking at large volume fractions: Tangential forces due to friction restrict microscopic rearrangements in the system, resulting in a large shear stress. The number of contacts created during the flow is the result of a competition between applied shear forces and the short-range repulsive forces. Therefore, alongside the contactless (hence frictionless) rheology at low shear rates a new contact dominated frictional rheology appears at high shear rates. Shear thickening is then simply the transition from the contactless rheology to the contact dominated rheology with increasing shear stress (and shear rate). The strongest evidence that this stress-induced friction scenario is correct is the ability of our model to reproduce shear thickening, including its discontinuous form, and its volume-fraction dependence. Therefore, we conclude that friction is a key element for shear thickening.

The shear thickening mechanism we describe is qualitatively independent of our parameter choices. Quantitatively, we can summarize the effects of the various parameters as follows. The friction coefficient μ only affects the volume fraction range over which the phenomenon happens: If μ increases, the critical volume fraction ϕ c decreases. Any friction coefficient μ > 0 leads to a DST, with ϕ c ( μ ) ϕ J 0 in the limit μ 0 . The lubrication cutoff δ mainly sets the dissipation level in the systems and thus sets an overall prefactor in the entire curve η(σ): The viscosity is larger when δ decreases. For instance, we performed simulations with a cutoff δ = 10 2 and compared the results to the ones shown in the present work obtained with δ = 10 3 : The curve η δ = 10 2 ( σ ) can be superimposed to η δ = 10 3 ( σ ) by a mere multiplicative factor η δ = 10 2 ( σ ) 0.8 η δ = 10 3 ( σ ) . This has the effect of pushing the shear thickening to lower shear rates, as the onset stress σST is independent of δ. The particle stiffness constants kn and kt are chosen such that the simulations are run in the asymptotic regime where the dimensionless numbers 6 π η 0 a γ ̇ / k n , t are very small. In this asymptotic regime, the rheology is independent of these stiffness constants. The Debye length κ–1 governs the shear-thinning behavior at low shear rates. Our CLM can be thought as having a vanishing Debye length. Comparing with the ERM with κ– 1 = 0.05a (with a the radius of the small particles), we observe that besides the shear-thinning regime, a finite κ–1 has virtually no effect on the shear thickening and the high viscosity state.

We now turn to a quantitative comparison with experimental data. The main tuning parameter in our description is the force scale FCL for CLM (or FER for the ERM). This parameter sets the onset stress at which shear thickening occurs. We can estimate this force for a system of 1 μm silica spheres in water. We find an electrostatic repulsion at contact of order F ER 0.2 nN [as evaluated in Behrens and Grier (2001) with a charge density 5 × 10 4 C m 2 and a Debye length κ 1 0.1 μ m ]. This gives for our model the stress scale η 0 γ ̇ 0 = ( 6 π a 2 ) 1 F ER , which together with the nondimensionalized onset stress σ ST / ( η 0 γ ̇ 0 ) 0.3 obtained the ERM model (see Fig. 4) gives a predicted onset stress of σ ST 8 Pa . For the very same system, Lootens et al. (2005) experimentally find an onset stress of σ ST 5 Pa , with which our simulation is thus in good agreement. Thus, our mechanism seems to give a quantitative description of shear thickening in dense suspensions, even if it would require comparison to more experimental data to be definitely conclusive. (Obtaining a reliable value of the microscopic force scale is rarely feasible in the published data on DST, and this restricts the number of comparisons we can make at the moment.)

Although our model assumes non-Brownian suspensions, we may expect that the same mechanism explains shear thickening of Brownian colloidal suspensions. Indeed, the Brownian forces qualitatively play a role similar to the repulsive potential that we use in our ERM: At low shear stress, the Brownian forces are large enough to randomly open gaps between contacting particles, while at high shear stress this becomes much less likely, as relative trajectories of pairs of particles become smooth at high Peclet numbers [Nazockdast and Morris (2013)]. The fact that a random noise can relax friction between grains is a known phenomenon in granular matter [Gao et al. (2009); Karim and Corwin (2014)]. Friction would then develop in a Brownian suspension in a manner very similar to the present work.

The authors would like to thank Mike Cates and Matthieu Wyart for the many discussions concerning our respective work and Philippe Claudin for discussions about the μ(Iv) rheology. The authors' code makes use of the CHOLMOD library by Tim Davis (https://www.cise.ufl.edu/research/sparse/cholmod/) for direct Cholesky factorization of the sparse resistance matrix. This research was supported in part by a grant of computer time from the City University of New York High Performance Computing Center under NSF Grant Nos. CNS-0855217, CNS-0958379, and ACI-1126113. J.F.M. was supported in part by NSF PREM (DMR 0934206).

The hydrodynamic interactions arising from an imposed background flow and relative motions between nearby particles are given by

( F H T H ) = ( R Stokes + R Lub ) · ( U U Ω Ω ) + R Lub : E ,
(A1)

where a diagonal matrix R Stokes comes from the (one-body) Stokes drag and sparse matrices R Lub and R Lub come from the (two-body) lubrication.

Using the basic units L 0 a for lengths, U 0 L 0 γ ̇ for velocities, and F 0 6 π η 0 L 0 U 0 for forces, the elements of R Stokes give the Stokes drag through

( F Stokes ( i ) T Stokes ( i ) ) = R Stokes ( i , i ) · ( U ( i ) U ( r ( i ) ) Ω ( i ) Ω ) = ( U ( i ) U ( r ( i ) ) ( 4 / 3 ) ( Ω ( i ) Ω ) ) ,
(A2)

while R Lub and R Lub consist of off-diagonal blocks giving the lubrication forces and torques for a pair (i, j) through

( F H ( i , j ) F H ( j , i ) T H ( i , j ) T H ( j , i ) ) = R Lub ( i , j ) · ( U ( i ) U ( r ( i ) ) U ( j ) U ( r ( j ) ) Ω ( i ) Ω Ω ( j ) Ω ) + R H ( i , j ) : ( E E ) .
(A3)

Following the notation of Jeffrey and Onishi (1984) and Jeffrey (1992), the matrices R Lub ( i , j ) and R Lub ( i , j ) are obtained from the particle positions as (we give only the upper triangular part of the symmetric R Lub ( i , j ) )

R Lub ( i , j ) ( X i i A P n i j + Y i i A P n i j X i j A P n i j + Y i j A P n i j Y i i B P n i j r Y j i B P n i j r . X j j A P n i j + Y j j A P n i j Y i j B P n i j r Y j j B P n i j r . . Y i i C P n i j Y i j C P n i j . . . Y j j C P n i j ) ,
(A4)
R Lub ( i , j ) ( X i i G Q n i j + Y i i G Q n i j X j i G Q n i j + Y j i G Q n i j X i j G Q n i j + Y i j G Q n i j X j j G Q n i j + Y j j G Q n i j Y i i H Q n i j r Y j i H Q n i j r Y i j H Q n i j r Y j j H Q n i j r ) .
(A5)

In these expressions, we have introduced the normal projection operator P n i j n i j n i j , the tangential projection operator P n i j I n i j n i j , and the “cross product” operator P n i j r defined as P n i j r q n i j × q . We also used the operators Q n i j , Q n i j , and Q n i j r defined for an arbitrary matrix M as

Q n i j M ( M : n i j n i j 1 3 Tr M ) n i j , Q n i j M ( M + M T ) · n i j 2 ( M : n i j n i j ) n i j , Q n i j r M 2 n i j × [ ( M + M T ) · n i j ] .
(A6)

The scalar coefficients X and Y have an explicit dependence on the nondimensional interparticle gap h ( i , j ) 2 ( r ( i , j ) a i a j ) / ( a i + a j ) , and we use only the terms of leading order

X g X 1 h ( i , j ) + δ , Y g Y log 1 h ( i , j ) + δ .
(A7)

With λ a j / a i , the coefficients gX and gY appearing in R Lub ( i , j ) are written

g X i i A ( λ ) = 2 a i λ 2 ( 1 + λ ) 3 , g X j j A ( λ ) = λ g X i i A ( λ 1 ) , g X i j A ( λ ) = 2 ( a i + a j ) λ 2 ( 1 + λ ) 4 , g X j i A ( λ ) = g X i j A ( λ 1 ) = g X i j A ( λ ) , g Y i i A ( λ ) = 4 a i 15 λ ( 2 + λ + 2 λ 2 ) ( 1 + λ ) 3 , g Y j j A ( λ ) = λ g Y i i A ( λ 1 ) , g Y i j A ( λ ) = 4 ( a i + a j ) 15 λ ( 2 + λ + 2 λ 2 ) ( 1 + λ ) 4 , g Y j i A ( λ ) = g Y i j A ( λ 1 ) = g Y i j A ( λ ) , g Y i i B ( λ ) = 2 a i 2 15 λ ( 4 + λ ) ( 1 + λ ) 2 , g Y j j B ( λ ) = λ 2 g Y i i B ( λ 1 ) , g Y i j B ( λ ) = 2 ( a i + a j ) 2 15 λ ( 4 + λ ) ( 1 + λ ) 4 , g Y j i B ( λ ) = g Y i j B ( λ 1 ) , g Y i i C ( λ ) = 8 a i 3 15 λ 1 + λ , g Y j j C ( λ ) = λ 3 g Y i i C ( λ 1 ) = λ 2 g Y i i C ( λ ) , g Y i j C ( λ ) = 2 ( a i + a j ) 3 15 λ 2 ( 1 + λ ) 4 , g Y j i C ( λ ) = g Y i j C ( λ 1 ) = g Y i j C ( λ ) .
(A8)

Similarly, the terms appearing in R Lub ( i , j ) are

g X i i G ( λ ) = 2 a i 2 λ 2 ( 1 + λ ) 3 , g X j j G ( λ ) = λ 2 g X i i G ( λ 1 ) , g X i j G ( λ ) = 2 ( a i + a j ) 2 λ 2 ( 1 + λ ) 5 , g X j i G ( λ ) = g X i j G ( λ 1 ) , g Y i i G ( λ ) = a i 2 15 4 λ λ 2 + 7 λ 3 ( 1 + λ ) 3 , g Y j j G ( λ ) = λ 2 g Y i i G ( λ 1 ) , g Y i j G ( λ ) = ( a i + a j ) 2 15 4 λ λ 2 + 7 λ 3 ( 1 + λ ) 5 , g Y j i G ( λ ) = g Y i j G ( λ 1 ) , g Y i i H ( λ ) = 2 a i 3 15 2 λ λ 2 ( 1 + λ ) 2 , g Y j j H ( λ ) = λ 3 g Y i i G ( λ 1 ) , g Y i j H ( λ ) = ( a i + a j ) 3 15 λ 2 + 7 λ 3 ( 1 + λ ) 5 , g Y j i H ( λ ) = g Y i j H ( λ 1 ) .
(A9)
1. Tuning the spring constants

Our objective is to study the rheology of hard-sphere suspensions. We have no way to mimic rigorously hard spheres by using a contact model with springs. To provide the best approximation to hard spheres, the elastic constants appearing in the contact model must satisfy a simple constraint: They should be large enough to generate as little geometric overlap as possible between the particles. (An equivalent way to state this is through the requirement that the contact based nondimensionalized shear rate 6 π η 0 a γ ̇ / k n 1 introduced in Sec. II B 1 is sufficiently small.)

We therefore set a criterion for the overlap: The largest overlap between any two particles during the simulation should not be larger than hmax (5% of the particle radius in this simulation). As the overlap | h | depends on the shear stress, an estimate being | h | ( 1 / k n ) σ ( ϕ , γ ̇ ) a 2 , the spring constant has to be tuned for every ϕ and γ ̇ so that max ( | h | ) < h max . We introduce a similar criterion for the tangential spring stretch mimicking static friction: We pick k t ( ϕ , γ ̇ ) so that ξ is smaller than 5% of the particle radius. We detail below how we fulfill these two conditions and retain hard-sphere like behavior in our simulation.

For a given ϕ , the largest stress is obtained in the high shear-rate limit. Hence, we start by determining high shear-rate spring constants ( k n ( ϕ ) , k t ( ϕ ) ) by running preliminary simulations at γ ̇ , where these values are regularly updated with a certain interval until the criteria on the overlap and the tangential spring stretch are fulfilled.

Second, a trivial shear-rate dependence comes in if the contact model introduces a force scale in addition to the hydrodynamic one. Avoiding this requires picking the shear rate dependence of k n ( ϕ , γ ̇ ) and k n ( ϕ , γ ̇ ) by scaling these parameters with the shear rate; i.e., k n ( ϕ , γ ̇ ) = γ ̇ k n ( ϕ ) and k t ( ϕ , γ ̇ ) = γ ̇ k t ( ϕ ) . With this scaling, there is no competition between hydrodynamic and contact forces, as they are proportional, so an additional explicit force scale (as the one in Sec. II B) is required to have a shear-rate dependence in our simulation. Also note that as the high shear-rate limit is always the largest viscosity at a given ϕ , with this choice of scaling k n ( ϕ , γ ̇ ) and k t ( ϕ , γ ̇ ) always fulfill the criteria we set for the overlap and tangential spring stretch.

2. Combining the overdamped dynamics and the contact model

In this appendix, we present a slightly more general contact model than the one presented in the main text. While the contact model may not deserve a lengthy description in itself, immersing an orthodox contact model in overdamped dynamics leads to nontrivial difficulties, which justifies our use of a simpler version of the model.

A more general stick/slide friction model using springs and dashpots is the following [Cundall and Strack (1979); Luding (2008)]:

F C , nor ( i , j ) = k n h ( i , j ) n i j + γ n U n ( i , j ) , F C , tan ( i , j ) = k t ξ ( i , j ) + γ t U t ( i , j ) , T C , tan ( i , j ) = a i n i j × F C , tan ( i , j ) .
(B1)

These forces must fulfill Coulomb's law of friction | F C , tan ( i , j ) | μ | F C , nor ( i , j ) | . In the above expression, kn and kt are, respectively, the normal and tangential spring constants, and γn and γ t are the damping constants. The normal and tangential relative velocities between two particles i and j are U n ( i , j ) P n i j ( U ( j ) U ( i ) ) and U t ( i , j ) P n i j [ U ( j ) U ( i ) ( a i Ω ( i ) + a j Ω ( j ) ) × n i j ] . Finally, the quantity ξ ( i , j ) is the tangential spring stretch.

The computation of the tangential spring stretch ξ ( i , j ) , described in the following, requires some care as we have to impose Coulomb's law, which is made difficult by the overdamped dynamics. At the time t0 at which the contact (i, j) is created, we set an unstretched tangential spring ξ ( i , j ) ( t 0 ) = 0 . At any further time step t in the simulation, the tangential stretch ξ ( i , j ) ( t ) is incremented according to the value of a test force, F C , tan ( i , j ) ( t + d t ) = k t ξ ( i , j ) ( t + d t ) + γ t U t ( i , j ) ( t + d t ) with the tentative update of stretch ξ ( i , j ) ( t + d t ) = ξ ( i , j ) ( t ) + U t ( i , j ) ( t ) d t . If | F C , tan ( i , j ) ( t + d t ) | μ | F C , nor ( i , j ) ( t + d t ) | , the contact is in a static friction state and we update the spring stretch as ξ ( i , j ) ( t + d t ) = ξ ( i , j ) ( t + d t ) and the corresponding tangential force is F C , tan ( i , j ) ( t + d t ) = F C , tan ( i , j ) ( t + d t ) . However, if | F C , tan ( i , j ) ( t + d t ) | > μ | F C , nor ( i , j ) ( t + d t ) | , the contact is in sliding state and the spring is updated as ξ ( i , j ) ( t + d t ) = ( 1 / k t ) ( μ | F C , nor ( i , j ) ( t + d t ) | t i j γ t U t ( i , j ) ( t ) ) , where the direction is the same as the test force, i.e., t i j F C , tan ( i , j ) ( t + d t ) / | F C , tan ( i , j ) ( t + d t ) | . Assuming that the velocities are continuous in time, this ensures that the contact force will at most weakly violate the Coulomb friction law for the next time step, as the force effectively used in the equation of motion will be F C , tan ( i , j ) ( t + d t ) = k t ξ ( i , j ) ( t + d t ) + γ t U t ( i , j ) ( t + d t ) = μ | F C , nor ( i , j ) ( t + d t ) | t i j + γ t U . t ( i , j ) d t , so that | F C , tan ( i , j ) ( t + d t ) | = μ | F C , nor ( i , j ) ( t + d t ) | + O ( d t ) .

Of course, if velocities happen to be discontinuous, the Coulomb's law might be violated significantly, and this is the source of one problem when merging this contact model with the hydrodynamic interaction model. Indeed, when contacts are created and destroyed during the flow (which happens very frequently), the overall resistance matrix changes discontinuously, as dashpots are switched on and off. But, as other position-dependent forces are continuous in time, solving the force balance equation will lead to a discontinuity in velocities when a contact forms or breaks. This, in turn, leads to large violations of Coulomb's law, and hence to numerical instabilities where contacts keep switching between static and dynamic cases at every time step.

Another problem occurs when a contact forms (i.e., ξ ( i , j ) ( t ) = 0 ): As by definition the normal load on the new contact | F C , nor ( i , j ) ( t + d t ) | vanishes (or is very small in the simulation owing to the finite time step), the finite test force F C , tan ( i , j ) ( t + d t ) γ t U t ( i , j ) ( t ) will lead to an immediate rescaling of the tangential spring as ξ ( i , j ) ( t + d t ) ( γ t / k t ) U t ( i , j ) ( t ) ; i.e., an entirely new finite force appears right at contact time.

These two problems actually also exist in simulations of dry granular matter [Walton (1993b)], but they are more acute with an overdamped dynamics because of the direct relation between forces and velocities. In this case, any discontinuity in forces or velocities rapidly causes numerical instabilities. A solution is to use a continuously varying damping γ t ( h i j ) such that γ t ( h i j = 0 ) = 0 , which ensures the continuity of forces and velocities [Walton (1993a)] at the cost of increased complexity of the model. The solution that we prefer using here is to eliminate of the tangential dashpot by setting γ t = 0 . In any case, this dashpot has no physical significance, and it is only used in granular matter simulations as an efficient numerical stabilizer. In our already overdamped context, as long as the hydrodynamic resistance associated with tangential motion is not too small, we do not need such an extra resistive stabilizer, and we do not face any major difficulty by simply dropping this dashpot.

We can quantify this assertion by looking at the relaxation times associated with a contact. Both normal and tangential contacts have relaxation times. For the normal part, it is the one of a spring and dashpot system, τ n = γ n / k n . For the tangential part, the damping is provided by the hydrodynamic resistance, which we here simply denote γ t H , so that τ t = γ t H / k t . On the one hand, in order to have a stable numerical scheme, one should choose these relaxation times large enough compared to the time step τ n , τ t d t . On the other hand, contacts of hard spheres should react instantaneously to any external load, so physics imposes relaxation times smaller than other typical times τ n , τ t 1 / γ ̇ . For the normal part of the contact, we are free to choose the damping γn to achieve this. For the tangential part however, we check that those scale separations are verified.

1.
Andreotti
,
B.
,
J.-L.
Barrat
, and
C.
Heussinger
, “
Shear flow of non-Brownian suspensions close to jamming
,”
Phys. Rev. Lett.
109
,
105901
(
2012
).
2.
Bagnold
,
R. A.
, “
Experiments on a gravity-free dispersion of large solid spheres in a Newtonian fluid under shear
,”
Proc. R. Soc. London, Ser. A
225
,
49
63
(
1954
).
3.
Ball
,
R. C.
, and
J. R.
Melrose
, “
Lubrication breakdown in hydrodynamic simulations of concentrated colloids
,”
Adv. Colloid Interface Sci.
59
,
19
30
(
1995
).
4.
Ball
,
R. C.
, and
J. R.
Melrose
, “
A simulation technique for many spheres in quasi-static motion under frame-invariant pair drag and Brownian forces
,”
Physica A
247
,
444
472
(
1997
).
5.
Barnes
,
H. A.
, “
Shear-thickening (“dilatancy”) in suspensions of nonaggregating solid particles dispersed in Newtonian liquids
,”
J. Rheol.
33
,
329
366
(
1989
).
6.
Batchelor
,
G. K.
, “
The stress system in a suspension of force-free particles
,”
J. Fluid Mech.
41
,
545
570
(
1970
).
7.
Behrens
,
S. H.
, and
D. G.
Grier
, “
The charge of glass and silica surfaces
,”
J. Chem. Phys.
115
,
6716
6721
(
2001
).
8.
Bender
,
J.
, and
N. J.
Wagner
, “
Reversible shear thickening in monodisperse and bidisperse colloidal dispersions
,”
J. Rheol.
40
,
899
916
(
1996
).
9.
Bender
,
J. W.
, and
N. J.
Wagner
, “
Optical measurement of the contributions of colloidal forces to the rheology of concentrated suspensions
,”
J. Colloid Interface Sci.
172
,
171
184
(
1995
).
10.
Bertrand
,
E.
,
J.
Bibette
, and
V.
Schmitt
, “
From shear thickening to shear-induced jamming
,”
Phys. Rev. E
66
,
060401
(
2002
).
11.
Bi
,
D.
,
J.
Zhang
,
B.
Chakraborty
, and
R. P.
Behringer
, “
Jamming by shear
,”
Nature
480
,
355
358
(
2011
).
12.
Blanc
,
F.
,
F.
Peters
, and
E.
Lemaire
, “
Experimental signature of the pair trajectories of rough spheres in the shear-induced microstructure in noncolloidal suspensions
,”
Phys. Rev. Lett.
107
,
208302
(
2011
).
13.
Boersma
,
W. H.
,
J.
Laven
, and
H. N.
Stein
, “
Shear thickening (dilatancy) in concentrated dispersions
,”
AIChE J.
36
,
321
332
(
1990
).
14.
Boersma
,
W. H.
,
P. J. M.
Baets
,
J.
Laven
, and
H. N.
Stein
, “
Time-dependent behavior and wall slip in concentrated shear thickening dispersions
,”
J. Rheol.
35
,
1093
1120
(
1991
).
15.
Bossis
,
G.
, and
J. F.
Brady
, “
The rheology of Brownian suspensions
,”
J. Chem. Phys.
91
,
1866
1874
(
1989
).
16.
Boyer
,
F.
,
É.
Guazzelli
, and
O.
Pouliquen
, “
Unifying suspension and granular rheology
,”
Phys. Rev. Lett.
107
,
188301
(
2011
).
17.
Brady
,
J. F.
, and
G.
Bossis
, “
The rheology of concentrated suspensions of spheres in simple shear flow by numerical simulation
,”
J. Fluid Mech.
155
,
105
129
(
1985
).
18.
Brady
,
J. F.
, and
G.
Bossis
, “
Stokesian dynamics
,”
Annu. Rev. Fluid Mech.
20
,
111
157
(
1988
).
19.
Brown
,
E.
, and
H. M.
Jaeger
, “
Dynamic jamming point for shear thickening suspensions
,”
Phys. Rev. Lett.
103
,
086001
(
2009
).
20.
Brown
,
E.
, and
H. M.
Jaeger
, “
The role of dilation and confining stresses in shear thickening of dense suspensions
,”
J. Rheol.
56
,
875
923
(
2012
).
21.
Brown
,
E.
, and
H. M.
Jaeger
, “
Shear thickening in concentrated suspensions: Phenomenology, mechanisms and relations to jamming
,”
Rep. Prog. Phys.
77
,
046602
(
2014
).
22.
Cates
,
M. E.
,
J. P.
Wittmer
,
J.-P.
Bouchaud
, and
P.
Claudin
, “
Jamming, force chains, and fragile matter
,”
Phys. Rev. Lett.
81
,
1841
1844
(
1998
).
23.
Cheng
,
X.
,
J. H.
McCoy
,
J. N.
Israelachvili
, and
I.
Cohen
, “
Imaging the microscopic structure of shear thinning and thickening colloidal suspensions
,”
Science
333
,
1276
1279
(
2011
).
24.
Couturier
,
É.
,
F.
Boyer
,
O.
Pouliquen
, and
É.
Guazzelli
, “
Suspensions in a tilted trough: Second normal stress difference
,”
J. Fluid Mech.
686
,
26
39
(
2011
).
25.
Crawford
,
N. C.
,
S. K. R.
Williams
,
D.
Boldridge
, and
M. W.
Liberatore
, “
Shear-induced structures and thickening in fumed silica slurries
,”
Langmuir
29
,
12915
12923
(
2013
).
26.
Cundall
,
P. A.
, and
O. D. L.
Strack
, “
A discrete numerical model for granular assemblies
,”
Geotechnique
29
,
47
65
(
1979
).
27.
da Cruz
,
F.
,
S.
Emam
,
M.
Prochnow
,
J.-N.
Roux
, and
F.
Chevoir
, “
Rheophysics of dense granular materials: Discrete simulation of plane shear flows
,”
Phys. Rev. E
72
,
021309
(
2005
).
28.
Dai
,
S.-C.
,
E.
Bertevas
,
F.
Qi
, and
R. I.
Tanner
, “
Viscometric functions for noncolloidal sphere suspensions with Newtonian matrices
,”
J. Rheol.
57
,
493
510
(
2013
).
29.
Davis
,
R. H.
, “
Effects of surface roughness on a sphere sedimenting through a dilute suspension of neutrally buoyant spheres
,”
Phys. Fluids A
4
,
2607
2619
(
1992
).
30.
Dbouk
,
T.
,
L.
Lobry
, and
E.
Lemaire
, “
Normal stresses in concentrated non-Brownian suspensions
,”
J. Fluid Mech.
715
,
239
272
(
2013
).
31.
Deboeuf
,
A.
,
G.
Gauthier
,
J.
Martin
,
Y.
Yurkovetsky
, and
J. F.
Morris
, “
Particle pressure in a sheared suspension: A bridge from osmosis to granular dilatancy
,”
Phys. Rev. Lett.
102
,
108301
(
2009
).
32.
Denn
,
M. M.
, and
J. F.
Morris
, “
Rheology of non-Brownian suspensions
,”
Annu. Rev. Chem. Biomol. Eng.
5
,
203–228
(
2014
).
33.
D'Haene
,
P.
,
J.
Mewis
, and
G. G.
Fuller
, “
Scattering dichroism measurements of flow-induced structure of a shear thickening suspension
,”
J. Colloid Interface Sci.
156
,
350
358
(
1993
).
34.
Fagan
,
M. E.
, and
C. F.
Zukoski
, “
The rheology of charge stabilized silica suspensions
,”
J. Rheol.
41
,
373
397
(
1997
).
35.
Fall
,
A.
,
A.
Lemaître
,
F.
Bertrand
,
D.
Bonn
, and
G.
Ovarlez
, “
Shear thickening and migration in granular suspensions
,”
Phys. Rev. Lett.
105
,
268303
(
2010
).
36.
Fall
,
A.
,
F.
Bertrand
,
G.
Ovarlez
, and
D.
Bonn
, “
Shear thickening of cornstarch suspensions
,”
J. Rheol.
56
,
575
591
(
2012
).
37.
Fernandez
,
N.
,
R.
Mani
,
D.
Rinaldi
,
D.
Kadau
,
M.
Mosquet
,
H.
Lombois-Burger
,
J.
Cayer-Barrioz
,
H. J.
Herrmann
,
N. D.
Spencer
, and
L.
Isa
, “
Microscopic mechanism for shear thickening of non-Brownian suspensions
,”
Phys. Rev. Lett.
111
,
108301
(
2013
).
38.
Forterre
,
Y.
, and
O.
Pouliquen
, “
Flows of dense granular media
,”
Annu. Rev. Fluid Mech.
40
,
1
24
(
2008
).
39.
Foss
,
D. R.
, and
J. F.
Brady
, “
Structure, diffusion and rheology of Brownian suspensions by stokesian dynamics simulation
,”
J. Fluid Mech.
407
,
167
200
(
2000
).
40.
Freundlich
,
H.
, and
H. L.
Roder
, “
Dilatancy and its relation to thixotropy
,”
Trans. Faraday Soc.
34
,
308
316
(
1938
).
41.
Frith
,
W. J.
,
P.
d'Haene
,
R.
Buscall
, and
J.
Mewis
, “
Shear thickening in model suspensions of sterically stabilized particles
,”
J. Rheol.
40
,
531
548
(
1996
).
42.
Gao
,
G.-J.
,
J.
Blawzdziewicz
,
C. S.
O'Hern
, and
M.
Shattuck
, “
Experimental demonstration of nonuniform frequency distributions of granular packings
,”
Phys. Rev. E
80
,
061304
(
2009
).
43.
Hébraud
,
P.
, and
D.
Lootens
, “
Concentrated suspensions under flow: Shear-thickening and jamming
,”
Mod. Phys. Lett. B
19
,
613
624
(
2005
).
44.
Heussinger
,
C.
, “
Shear thickening in granular suspensions: Inter-particle friction and dynamically correlated clusters
,”
Phys. Rev. E
88
,
050201(R)
(
2013
).
45.
Ho
,
C.-M.
, and
Y.-C.
Tai
, “
Micro-electro-mechanical-systems (mems) and fluid flows
,”
Annu. Rev. Fluid Mech.
30
,
579
612
(
1998
).
46.
Hoffman
,
R. L.
, “
Discontinuous and dilatant viscosity behavior in concentrated suspensions. I. Observation of a flow instability
,”
Trans. Soc. Rheol.
16
,
155
173
(
1972
).
47.
Hoffman
,
R. L.
, “
Discontinuous and dilatant viscosity behavior in concentrated suspensions. II. Theory and experimental tests
,”
J. Colloid Interface Sci.
46
,
491
506
(
1974
).
48.
Hoffman
,
R. L.
, “
Explanations for the cause of shear thickening in concentrated colloidal suspensions
,”
J. Rheol.
42
,
111
123
(
1998
).
49.
Israelachvili
,
J. N.
,
Intermolecular and Surface Forces
(
Academic Press
, New York,
2011
).
50.
Jeffrey
,
D. J.
, “
The calculation of the low Reynolds number resistance functions for two unequal spheres
,”
Phys. Fluids A
4
,
16
29
(
1992
).
51.
Jeffrey
,
D. J.
, and
Y.
Onishi
, “
Calculation of the resistance and mobility functions for two unequal rigid spheres in low-Reynolds-number flow
,”
J. Fluid Mech.
139
,
261
290
(
1984
).
52.
Jerkins
,
M.
,
M.
Schröter
,
H. L.
Swinney
,
T. J.
Senden
,
M.
Saadatfar
, and
T.
Aste
, “
Onset of mechanical stability in random packings of frictional spheres
,”
Phys. Rev. Lett.
101
,
018301
(
2008
).
53.
Karim
,
M. Y.
, and
E. I.
Corwin
, “
Eliminating friction with friction: 2d Janssen effect in a friction-driven system
,”
Phys. Rev. Lett.
112
,
188001
(
2014
).
54.
Kawasaki
,
T.
,
A.
Ikeda
, and
L.
Berthier
, “
Thinning or thickening? multiple rheological regimes in dense suspensions of soft particles
,” Europhysics Letters (to be published); e-print arXiv:1404.4778[cond-mat].
55.
Krieger
,
I. M.
, “
Rheology of monodisperse lattices
,”
Adv. Colloid Interface Sci.
3
,
111
136
(
1972
).
56.
Kulkarni
,
S. D.
, and
J. F.
Morris
, “
Ordering transition and structural evolution under shear in Brownian suspensions
,”
J. Rheol.
53
,
417
439
(
2009
).
57.
Larsen
,
R. J.
,
J.-W.
Kim
,
C. F.
Zukoski
, and
D. A.
Weitz
, “
Elasticity of dilatant particle suspensions during flow
,”
Phys. Rev. E
81
,
011502
(
2010
).
58.
Laun
,
H. M.
,
R.
Bung
,
S.
Hess
,
W.
Loose
,
O.
Hess
,
K.
Hahn
,
HäE.
dicke
,
R.
Hingmann
,
F.
Schmidt
, and
P.
Lindnerm
, “
Rheological and small angle neutron scattering investigation of shear-induced particle structures of concentrated polymer dispersions submitted to plane Poiseuille and Couette flow
,”
J. Rheol.
36
,
743
787
(
1992
).
59.
Lee
,
M.
,
M.
Alcoutlabi
,
J. J.
Magda
,
C.
Dibble
,
M. J.
Solomon
,
X.
Shi
, and
G. B.
McKenna
, “
The effect of the shear-thickening transition of model colloidal spheres on the sign of n1 and on the radial pressure profile in torsional shear flows
,”
J. Rheol.
50
,
293
311
(
2006
).
60.
Lees
,
A. W.
, and
S. F.
Edwards
, “
The computer study of transport processes under extreme conditions
,”
J. Phys. C.
5
,
1921
1930
(
1972
).
61.
Lemaître
,
A.
,
J.-N.
Roux
, and
F.
Chevoir
,
“What do dry granular flows tell us about dense non-Brownian suspension rheology?”
Rheol. Acta
48
,
925
942
(
2009
).
62.
Lerner
,
E.
,
G.
Düring
, and
M.
Wyart
, “
A unified framework for non-Brownian suspension flows and soft amorphous solids
,”
Proc. Natl. Acad. Sci. U.S.A.
109
,
4798
4803
(
2012
).
63.
Lootens
,
D.
,
H.
Van Damme
, and
P.
Hébraud
, “
Giant stress fluctuations at the jamming transition
,”
Phys. Rev. Lett.
90
,
178301
(
2003
).
64.
Lootens
,
D.
,
H.
van Damme
,
Y.
Hémar
, and
P.
Hébraud
, “
Dilatant flow of concentrated suspensions of rough particles
,”
Phys. Rev. Lett.
95
,
268302
(
2005
).
65.
Lootens
,
D.
,
P.
Hébraud
,
E.
Lécolier
, and
H.
Van Damme
, “
Gelation, shear-thinning and shear-thickening in cement slurries
,”
Oil Gas Sci. Technol.- Rev. IFP.
59
,
31
40
(
2004
).
66.
Luding
,
S.
, “
Cohesive, frictional powders: Contact models for tension
,”
Granular Matter
10
,
235
246
(
2008
).
67.
Maranzano
,
B. J.
, and
N. J.
Wagner
, “
The effects of interparticle interactions and particle size on reversible shear thickening: Hard-sphere colloidal dispersions
,”
J. Rheol.
45
,
1205
1222
(
2001a
).
68.
Maranzano
,
B. J.
, and
N. J.
Wagner
, “
The effects of particle size on reversible shear thickening of concentrated colloidal dispersions
,”
J. Chem. Phys.
114
,
10514
10527
(
2001b
).
69.
Maranzano
,
B. J.
, and
N. J.
Wagner
, “
Flow-small angle neutron scattering measurements of colloidal dispersion microstructure evolution through the shear thickening transition
,”
J. Chem. Phys.
117
,
10291
10302
(
2002
).
70.
Melrose
,
J. R.
, and
R. C.
Ball
, “
“contact networks” in continuously shear thickening colloids
,”
J. Rheol.
48
,
961
978
(
2004a
).
71.
Melrose
,
J. R.
, and
R. C.
Ball
, “
Continuous shear thickening transitions in model concentrated colloids—The role of interparticle forces
,”
J. Rheol.
48
,
937
960
(
2004b
).
72.
Metzner
,
A. B.
, and
M.
Whitlock
, “
Flow behavior of concentrated (dilatant) suspensions
,”
Trans. Soc. Rheol.
2
,
239
253
(
1958
).
73.
Mewis
,
J.
, and
N. J.
Wagner
,
Colloidal Suspension Rheology
(
Cambridge University Press
, New York,
2011
).
74.
Morris
,
J. F.
, “
A review of microstructure in concentrated suspensions and its implications for rheology and bulk flow
,”
Rheol. Acta
48
,
909
923
(
2009
).
75.
Nakanishi
,
H.
,
S.
Nagahiro
, and
N.
Mitarai
, “
Fluid dynamics of dilatant fluids
,”
Phys. Rev. E
85
,
011401
(
2012
).
76.
Nazockdast
,
E.
, and
J. F.
Morris
, “
Pair-particle dynamics and microstructure in sheared colloidal suspensions: Simulation and Smoluchowski theory
,”
Phys. Fluids
25
,
070601
(
2013
).
77.
Newstein
,
M. C.
,
H.
Wang
,
N. P.
Balsara
,
A. A.
Lefebvre
,
Y.
Shnidman
,
H.
Watanabe
,
K.
Osaki
,
T.
Shikata
,
H.
Niwa
, and
Y.
Morishima
, “
Microstructural changes in a colloidal liquid in the shear thinning and shear thickening regimes
,”
J. Chem. Phys.
111
,
4827
4838
(
1999
).
78.
O'Brien
,
V. T.
, and
M. E.
Mackay
, “
Stress components and shear thickening of concentrated hard sphere suspensions
,”
Langmuir
16
,
7931
7938
(
2000
).
79.
O'Hern
,
C. S.
,
L. E.
Silbert
,
A. J.
Liu
, and
S. R.
Nagel
, “
Jamming at zero temperature and zero applied stress: The epitome of disorder
,”
Phys. Rev. E
68
,
011306
(
2003
).
80.
Otsuki
,
M.
, and
H.
Hayakawa
, “
Critical scaling near jamming transition for frictional granular particles
,”
Phys. Rev. E
83
,
051301
(
2011
).
81.
Phung
,
T. N.
,
J. F.
Brady
, and
G.
Bossis
, “
Stokesian dynamics simulation of Brownian suspensions
,”
J. Fluid Mech.
313
,
181
207
(
1996
).
82.
Seto
,
R.
,
R.
Mari
,
J. F.
Morris
, and
M. M.
Denn
, “
Discontinuous shear thickening of frictional hard-sphere suspensions
,”
Phys. Rev. Lett.
111
,
218301
(
2013
).
83.
Silbert
,
L. E.
, “
Jamming of frictional spheres and random loose packing
,”
Soft Matter
6
,
2918
2924
(
2010
).
84.
Silbert
,
L. E.
,
D.
Ertaş
,
G.
Grest
,
T. C.
Halsey
, and
D.
Levine
, “
Geometry of frictionless and frictional sphere packings
,”
Phys. Rev. E
65
,
1
6
(
2002
).
85.
Singh
,
A.
, and
P. R.
Nott
, “
Experimental measurements of the normal stresses in sheared stokesian suspensions
,”
J. Fluid Mech.
490
,
293
320
(
2003
).
86.
Song
,
C.
,
P.
Wang
, and
H. A.
Makse
, “
A phase diagram for jammed matter
,”
Nature
453
,
629
632
(
2008
).
87.
Trulsson
,
M.
,
B.
Andreotti
, and
P.
Claudin
, “
Transition from the viscous to inertial regime in dense suspensions
,”
Phys. Rev. Lett.
109
,
118305
(
2012
).
88.
Wagner
,
N. J.
, and
J. F.
Brady
, “
Shear thickening in colloidal dispersions
,”
Phys. Today
62
(
10
),
27
32
(
2009
).
89.
Walton
,
O. R.
, “
Numerical simulation of inclined chute flows of monodisperse, inelastic, frictional spheres
,”
Mech. Mater.
16
,
239
247
(
1993a
).
90.
Walton
,
O. R.
, “
Numerical simulation of inelastic, frictional particle-particle interactions
,”
in Particulate Two-Phaseflow (Butterworth-Heinemann
Boston,
1993b
), Chap. 25
, pp.
884
911
.
91.
Watanabe
,
H.
,
M.-L.
Yao
,
K.
Osaki
,
T.
Shikata
,
H.
Niwa
, and
Y.
Morishima
, “
Nonlinear rheology of a concentrated spherical silica suspension
,”
Rheol. Acta
36
,
524
533
(
1997
).
92.
Watanabe
,
H.
,
M.-L.
Yao
,
K.
Osaki
,
T.
Shikata
,
H.
Niwa
,
Y.
Morishima
,
N. P.
Balsara
, and
H.
Wang
, “
Nonlinear rheology and flow-induced structure in a concentrated spherical silica suspension
,”
Rheol. Acta
37
,
1
6
(
1998
).
93.
Williamson
,
R. V.
, “
Some unusual properties of colloidal dispersions
,”
J. Phys. Chem.
35
,
354
359
(
1930
).
94.
Williamson
,
R. V.
, and
W. W.
Hecker
, “
Some properties of dispersions of the quicksand type
,”
Ind. Eng. Chem.
23
,
667
670
(
1931
).
95.
Wyart
,
M.
, and
M. E.
Cates
, “
Discontinuous shear thickening without inertia in dense non-Brownian suspensions
,”
Phys. Rev. Lett.
112
,
098302
(
2014
).
96.
Yurkovetsky
,
Y.
, and
J. F.
Morris
, “
Particle pressure in sheared Brownian suspensions
,”
J. Rheol.
52
,
141
164
(
2008
).
97.
Zarraga
,
I. E.
,
D. A.
Hill
, and
D. T.
Leighton
, “
The characterization of the total stress of concentrated suspensions of noncolloidal spheres in Newtonian fluids
,”
J. Rheol.
44
,
185
220
(
2000
).
98.
Zhao
,
Y.
, and
R. H.
Davis
, “
Interaction of two touching spheres in a viscous fluid
,”
Chem. Eng. Sci.
57
,
1997
2006
(
2002
).
99.
See supplementary material (movies) at http://dx.doi.org/10.1122/1.4890747 for (a) the contact network at volume fraction 0.56 around the DST transition and (b) the configuration under flow seen from upstream in the flow direction.

Supplementary Material