Open Access Article
This Open Access Article is licensed under a
Creative Commons Attribution 3.0 Unported Licence

Testing the Wyart–Cates model for non-Brownian shear thickening using bidisperse suspensions

Ben M. Guy*a, Christopher Ness*bc, Michiel Hermesad, Laura J. Sawiaka, Jin Sunb and Wilson C. K. Poona
aSchool of Physics and Astronomy, University of Edinburgh, King's Buildings, Peter Guthrie Tait Road, Edinburgh EH9 3FD, UK. E-mail:
bSchool of Engineering, University of Edinburgh, King's Buildings, Peter Guthrie Tait Road, Edinburgh EH9 3JL, UK. E-mail:
cDepartment of Chemical Engineering and Biotechnology, University of Cambridge, Cambridge CB3 0AS, UK
dSoft Condensed Matter, Debye Institute for Nanomaterials Science, Utrecht University, Princetonplein 5, 3584 CC Utrecht, The Netherlands

Received 7th January 2019 , Accepted 15th November 2019

First published on 18th November 2019

There is a growing consensus that shear thickening of concentrated dispersions is driven by the formation of stress-induced frictional contacts. The Wyart–Cates (WC) model of this phenomenon, in which the microphysics of the contacts enters solely via the fraction f of contacts that are frictional, can successfully fit flow curves for suspensions of weakly polydisperse spheres. However, its validity for “real-life”, polydisperse suspensions has yet to be seriously tested. By performing systematic simulations on bidisperse mixtures of spheres, we show that the WC model applies only in the monodisperse limit and fails when substantial bidispersity is introduced. We trace the failure of the model to its inability to distinguish large–large, large–small and small–small frictional contacts. By fitting our data using a polydisperse analogue of f that depends separately on the fraction of each of these contact types, we show that the WC picture of shear thickening is incomplete. Systematic experiments on model shear-thickening suspensions corroborate our findings, but highlight important challenges in rigorously testing the WC model with real systems. Our results prompt new questions about the microphysics of thickening for both monodisperse and polydisperse systems.

1 Introduction

Shear thickening, the increase in viscosity η with shear stress σ or rate [small gamma, Greek, dot above], is ubiquitous in concentrated suspensions.1 Its microscopic origin has been hotly debated.2 Recent experiments,3–7 simulations8,9 and theoretical modelling10 point to a σ-dependent transition from frictionless (sliding) to frictional (rolling) inter-particle contacts. A phenomenological model by Wyart and Cates10 (WC) predicts thickening based on a single microphysical parameter, the fraction of frictional contacts, f. It fits well the rheology of model systems;3,5 however, its validity for complex industrial suspensions remains untested. We systematically explore the conditions under which the WC model breaks down for one kind of complexity: size polydispersity, and reveal important shortcomings in our current understanding of shear thickening.

The phenomenology is generic.3,5,11 Fig. 1(a) shows literature flow curves,3 η(σ), for buoyancy-matched suspensions of polymethylmethacrylate (PMMA) spheres with diameter d ≈ 4 μm at different volume fractions ϕ. At any fixed σ (vertical lines), the viscosity increases with ϕ, Fig. 1(b) (symbols). The viscosity “branches” at different σ are well described by

η/ηf = (1 − ϕ/ϕJ)−2, (1)
which diverges at a σ-dependent jamming volume fraction ϕJ(σ); ηf is the solvent viscosity. Fig. 1(b) shows example fits of eqn (1) (lines) with ϕJ as a free parameter. The fitted ϕJ(σ), Fig. 1(c), is a decreasing function of σ; so, increasing σ at fixed ϕ, i.e., traversing a vertical path in Fig. 1(b) (arrow), decreases ϕJ and causes η to increase: the suspension shear thickens, Fig. 1(a). The limiting low- and high-σ viscosity plateaux, η0 and ηm [blue and red in Fig. 1(b)], diverge at ϕ0 and ϕm < ϕ0, respectively.

image file: c9sm00041k-f1.tif
Fig. 1 Experimental shear thickening phenomenology. (a) Relative viscosity η/ηf as a function of shear stress σ at different volume fractions ϕ (as labelled) for d = 3.78 μm, PHSA-stabilised, PMMA spheres in a cyclohexylbromide–decalin mixture of viscosity ηf = 2.83 × 10−3 Pa s (taken from Guy et al.3). The grey region is inaccessible due to inertial edge fracture. (b) Symbols, viscosity “branches” for different (fixed) values of σ in (a). Lines correspond to fits by eye to eqn (1). The jamming volume fraction at which each viscosity branch diverges, ϕJ, depends on σ. Blue and red lines and symbols correspond respectively to the limiting low-σ and high-σ viscosities, η0 and ηm. (c) ϕJ(σ), obtained from the fits shown in (b). ϕJ(σ) decreases smoothly from ϕ0, the ϕ at which η0 diverges, to ϕm, the ϕ at which ηm diverges. In all parts: shear thickening arises at any ϕ, e.g., ϕ = 0.51, because increasing σ decreases ϕJ [black arrow in (c)], shifting the viscosity branch in (b) to the left and so increasing η [black arrows in (b) and (a)].

There is evidence of this scenario in a range of experimental systems.3,5,11 The precise values of ϕ0 and ϕm, and the form of ϕJ(σ) [and hence η(σ)], depend on details of particle shape,12 size polydispersity11 and surface roughness.13,14 In all systems, shear-induced jamming,15 inhomogeneous flow (shear banding)16 or unsteady flow17 are observed for ϕmϕ < ϕ0, where conditions exist for which ϕ > ϕJ, Fig. 1(c), i.e., the system can exhibit solid-like behaviour.

In the WC model,10 inter-particle contacts are either lubricated, with static friction coefficient μ = 0, or frictional, with μ > 0. The fraction of the latter, f, increases with σ, Fig. 2(a). WC's jamming volume fraction is a function of f only:

ϕWCJ = m + (1 − f)ϕ0, (2)
changing linearly from random close packing, ϕ0, at f = 0 (all lubricated contacts) to frictional jamming, ϕm, at f = 1 (all frictional contacts), Fig. 2(b). Thus, ϕWCJ(f(σ)) decreases with σ, Fig. 2(c), and determines η via some empirical form, e.g., eqn (1),3,5 leading to shear-thickening flow curves, Fig. 2(d) (line).

image file: c9sm00041k-f2.tif
Fig. 2 Logic of the WC shear thickening model. (a) The fraction of frictional contacts, f, takes a sigmoidal form. (b) ϕWCJ(f) is linearly interpolated between ϕ0 at f = 0 to ϕm at f = 1. (c) The previous two plots directly give ϕWCJ(f(σ)), which is inverse sigmoidal. (d) Using ϕWCJ(f(σ)) in eqn (1) gives η(σ), which shows shear thickening (line). Testing the WC model using simulations. The plot in (a) is calculated directly using contact forces from simulations of pure small-sphere suspension at ϕ = 0.53. See the text for how we obtain values for ϕ0 and fm to calculate the ϕWCJ(f) plotted in (b) from eqn (2). These two plots directly give the ϕWCJ in (c), which, when used in eqn (1) gives the flow curve in (d), ηWC(σ) (line). The symbols in (d) are the computed viscosity from simulations.

The WC model, eqn (1) and (2), predicts the σ- and ϕ-dependent viscosity, ηWC(σ,ϕ), from three inputs: the limiting frictionless and frictional jamming points, ϕ0 and ϕm, and the σ-dependent fraction of frictional contacts, f. ϕm and ϕ0, can be obtained by fitting viscosity branches at different σ, as done in Fig. 1(b). They are not related to shear thickening per se. On the other hand, f, which determines the shape of the flow curve, is currently inaccessible in experiments. Thus, various ansatzs are used to fit the WC model to experiments. For sterically-stabilised PMMA spheres, Guy et al.3 used a ϕ-independent sigmoidal form:

f(σ) = exp[−(σ*/σ)β], (3)
with β = 0.85. The single stress scale, σ*, scales as the “engineering” onset stress at which η(σ) visibly begins to increase.

Importantly for this study, we note that the particle size does not appear in WC model. On the other hand, and perhaps significantly in light of the findings we will report, Guy et al. found that the onset stress decreases with particle size, with σ* ∝ d−2, suggesting σ* ∝ F*/d2 for their PMMA particles, where F* ∼ kBT/nm is a constant force.3,10 Royer et al.5 used a similar form to fit data for dispersions of charge-stabilised silica.

In discrete-element (DEM) simulations, f can be calculated directly from particle coordinates. A popular choice is to use the “critical-load model” (CLM), in which μ jumps from 0 to >0 when the normal contact force between particles, F, exceeds a threshold value, F* (the critical load). This model reproduces9 the phenomenology of Fig. 1 and unstable flow at high ϕ. For a bidisperse mixture of spheres with diameter d1 and d2 = d1/1.4, Mari et al. found9 a ϕ-independent f(σ) of the form eqn (3), with β = 1.1 and σ* ≈ F*/[6π(d2/2)2], and later used it to fit flow curves at a range of ϕ.18 Thus, in this one case, the WC model is fully validated: using the measured fraction of frictional contacts in eqn (2) correctly predicts the viscosity. The similarity between the forms of f(σ) used to fit experiments and measured in simulations suggests that f in mildly polydisperse experimental systems can indeed be well described by eqn (3) or some similar form.

In these studies, a weak polydispersity (= standard deviation normalised by the mean of the particle size distribution) of s ≲ 20% was used to inhibit shear-induced crystallisation.9 Industrial dispersions typically have broad, often multimodal, size distributions with s ≳ 100%. Nevertheless, the low-s phenomenology in Fig. 1 continues to hold.11,19–21 However, the validity of the WC model in such higher-s suspensions has not been tested.

Indeed, it is a surprise for the WC model to work, and work well, even for low-s systems. In suspension rheology, details of the microstructure, fabric of the contact network and distribution of forces matter. Cates pointed out long ago that the relatively small number of nearest neighbours, ∼[scr O, script letter O](10), usually precludes any mean-field description.22 The success of his shear-thickening model with Wyart contradicts this norm. In the WC model, η(σ) is controlled primarily by a single scalar parameter f that is agnostic to exact microstructural details. For this reason alone, it is important to test the limitations of the WC model.

Here, we do so in suspensions of strongly bidisperse spheres. As before,9,18 we use DEM simulations to extract f(σ) for different mixtures and compare the predictions of the WC model to bulk flow curves. We then use the same f(σ) to test the model against experimental data for bidisperse PMMA spheres. We find that WC works for nearly-monodisperse suspensions [i.e., the simulated f(σ) correctly predicts η(σ)], but fails in general for bidisperse suspensions. We show that, nevertheless, the model can be extended to fit our data if the fractions of each contact type (large–large, large–small and small–small) are taken into account separately. Our results indicate that, in its original form, the WC model is at least incomplete, and highlight a number of unresolved issues in the current understanding of shear thickening and in the use of the WC framework to make inferences about microphysics. We propose directions for future research to address these issues.

2 Methods

A binary mixture of spheres is characterised by four parameters: d1, d2 < d1, the fraction of small particles ξ = V2/(V1 + V2) (where V1 and V2 are the total volumes of large and small particles, respectively) and the total volume fraction ϕ = (V1 + V2)/V (where V is the total volume). We fix ϕ and the size ratio αd2/d1 ≈ 0.25, and vary ξ. At our α, the small spheres are slightly too large to fit inside four touching large spheres (for which αmax = 0.225).

We sheared N = 2000 bidisperse, repulsive spheres at fixed ϕ in a periodic cell with Lees–Edwards boundary conditions. Short-range lubrication and repulsive contact forces described by linear springs of stiffness k were resolved using a classical DEM code that allows marginal overlaps δ between the surfaces of pairs of particles.23 We employed a contact model9 in which Coulomb friction with static friction coefficient μ = 1 appears beyond a critical overlap δ*, corresponding to a critical normal load F* = *. For simplicity, and consistency with experiments for nearly-monodisperse systems,3 we use a constant F* (and hence δ*) that is independent of d1 and d2.

Our unit of stress is σ0 = F*/(3πd22/2), at which purely small particles (ξ = 1) shear thicken;9 pure large spheres shear thicken at ∼α2σ0. For bidisperse mixtures, we averaged σ over the strain interval γ ⊂ [1.5,3] or ⊂ [1.5,2], in which the system had reached steady state for all ξ. For monodisperse suspensions, we averaged over γ ⊂ [0.7,1] to avoid the onset of large-scale crystallization.24

We performed experiments on binary suspensions of PMMA spheres stabilised with poly-12-hydroxystearic acid (PHSA) in a near-density-matching mixture of cyclohexylbromide and decalin (density ≈ 1.18 g cm−3, viscosity ηf = 2.4 mPa s). We varied ξ at a fixed ϕ = 0.51 by mixing together monomodal suspensions with mean particle diameters d2 = 0.712 μm and d1 = 2.76 μm and s ≈ 10% (from static light scattering). Flow curves were measured using an Anton Paar MCR 301 rheometer with sandblasted steel cone (angle 1°; diameter 50 mm; truncation 100 μm; roughness ∼10 μm) and roughened aluminium base plate (roughness ∼10 μm) at 20 °C. A solvent trap minimised evaporation. Details of experiments are given in the ESI.

3 Results

3.1 Bidisperse shear thickening phenomenology

We first present the simulated relative viscosity η/ηf as a function of shear stress σ/σ0 for fixed ϕ = 0.53 and α = 0.25 at various fractions of small particles, Fig. 3(a), ξ = 0 (pentagon), 0.2 (□), 0.5 (▽), 0.65 (△), 0.8 (○) and 1 (⋄).
image file: c9sm00041k-f3.tif
Fig. 3 Bidisperse shear thickening phenomenology. (a) η/ηf as a function of σ/σ0 from simulations at α = 0.25 and ϕ = 0.53, with ξ = 0 (pentagon), 0.2 (□), 0.5 (▽), 0.65 (△), 0.8 (○) and 1 (⋄). (b) Frictionless relative viscosity η0/ηf from simulations (image file: c9sm00041k-u7.tif) and experiments (image file: c9sm00041k-u8.tif), and frictional relative viscosity ηm/ηf from simulations (image file: c9sm00041k-u9.tif). (c) Limiting jamming volume fractions, ϕ0 (blue) and ϕm < ϕ0 (red), versus ξ from simulations (image file: c9sm00041k-u10.tif, image file: c9sm00041k-u11.tif) and experiments (image file: c9sm00041k-u12.tif, image file: c9sm00041k-u13.tif). (d) Experimental η/ηf versus σ for PMMA spheres at α = 0.26 and ϕ = 0.51, with ξ = 0 (pentagon), 0.2 (□), 0.5 (▽), 0.65 (△), 0.8 (○) and 1 (⋄). Inertial fracture at [small gamma, Greek, dot above] ≈ 8 × 103 s−1 renders the grey-shaded region inaccessible.3

Bidisperse and monodisperse flow curves are qualitatively similar, showing shear thickening between two Newtonian plateaux. Fig. 3(b) shows the ξ-dependent plateau viscosities, η0(ξ) (image file: c9sm00041k-u1.tif) and ηm(ξ) (image file: c9sm00041k-u2.tif), estimated by eye from Fig. 3(a). Mixing particles reduces both limiting viscosities relative to the values for single-sized spheres. Such a “Farris effect”25 has been widely observed in fixed-friction (i.e., Newtonian) suspensions.25,26 The limiting volume fractions, ϕ0(ξ) and ϕm(ξ), Fig. 3(c) (image file: c9sm00041k-u3.tif and image file: c9sm00041k-u4.tif respectively), are calculated using the simulated η0 and ηm values in eqn (1) with ϕ = 0.53. The non-monotonic behaviour directly mirrors that of η0(ξ) and ηm(ξ).

Experimental flow curves for binary mixtures of PMMA with α = 0.26 and ϕ = 0.51, Fig. 3(d), show similar phenomenology, except that the limiting high-σ behaviour is preempted by edge fracture due either to an inertial instability (grey region), or a different fracture mechanism when σ exceeds ≈103 Pa (ESI). Thus, we cannot access ηm(ξ) directly for all ξ. Shear thickening is preceded by shear thinning, presumably due to residual Brownian motion;§3 so, we estimate η0 by the viscosity minimum before the onset of thickening, Fig. 3(b) (image file: c9sm00041k-u5.tif). The experimental η0(ξ) show the same non-monotonicity as the simulated values, but are always too high, by up to a factor of ≲2 for the two end members (ξ = 0 or 1). Using the experimental values of η0(ξ) in eqn (1) with ϕ = 0.53 gives us an experimental estimate of ϕ0(ξ), Fig. 3(c) (image file: c9sm00041k-u6.tif). Consistent with the experimental viscosities η0(ξ) being higher than simulated values, the experimentally deduced ϕ0(ξ) are somewhat lower than the simulated values, i.e. the experimental system at ϕ = 0.53 is closer to jamming than the corresponding simulated system.

3.2 Comparing simulations to the WC model

We test the WC model by comparing simulated flow curves, η(σ), with those calculated using the fraction of frictional contacts, f(σ), measured from the simulations, ηWC(σ). To explain our procedure, consider data for monodisperse small particles (ξ = 1). First, we calculate f(σ) directly from inter-particle forces by counting, at each σ, the fraction of contacts with F > F*. The f(σ) so obtained, Fig. 2(a), is sigmoidal, similar to the f(σ) in bidisperse mixtures9 with α = 0.71 and eqn (3).

To calculate ϕWCJ(f) from eqn (2), we need ϕ0 and ϕm, which could be obtained by simulating and fitting η(σ,ϕ) at a range of ϕ, as done in Fig. 1. Instead, we use our simulated values of the low- and high-σ viscosities, η0 and ηm, at ϕ = 0.53 in eqn (1) to obtain ϕ0 and ϕm, giving the ϕWCJ(f) in Fig. 2(b).

From f(σ), Fig. 2(a), and ϕWCJ(f), Fig. 2(b), we now calculate ϕWCJ(f(σ)), Fig. 2(c), which decreases smoothly from ≈ϕ0 at σ/σ0 ≪ 1 to ϕm at σ/σ0 ≫ 1. Finally, we calculate the viscosity by substituting ϕWCJ(f(σ)), Fig. 2(c), into eqn (1). The flow curve, ηWC(σ), Fig. 2(d) (solid line), increases smoothly from η0 to ηm.

We compare this flow curve calculated using the measured f(σ) with the simulated η(σ), Fig. 2(d) (symbols). The two calculated plateaux agree with the simulated values by construction. The WC model is judged instead by how well it captures the shear thickening process. It does this well for monodisperse spheres. Both model-predicted and simulation viscosities start to increase at σ/σ0 ≈ 1, reaching saturation at σ/σ0 ≳ 30.

We repeat this procedure for our bidisperse suspensions with ξ = 0.2, 0.5, 0.65 and 0.8. Again, the measured f(σ), Fig. 4(a), and linearly interpolated ϕWCJ(f), Fig. 4(b), are used to calculate ϕWCJ(σ), Fig. 4(c), from which we obtain flow curves, Fig. 4(d) (lines). We compare these with the simulated viscosities, Fig. 4(d) (symbols), recalling that the limiting viscosities are constrained to fit, and noting that data for different ξ have been shifted vertically to aid visibility (see caption for shift factors).

image file: c9sm00041k-f4.tif
Fig. 4 Failure of the WC model for bidisperse simulations. (a) Fraction of frictional contacts f as a function of σ/σ0, extracted from simulations at different ξ, as labelled. (b) WC jamming volume fraction, ϕWCJ, as a function of f for different ξ [colours as in part (a)]. (c) ϕWCJ as a function of σ calculated using (b) and f(σ) from (a). (d) Symbols: shifted flow curves for different ξ, as labelled. Shift factors are: ξ = 0, 0.025; ξ = 0.2, 0.35; ξ = 0.5, 1.6; ξ = 0.65, 5; ξ = 0.8, 10 and ξ = 1, 15. Lines: predictions of the WC model shifted by the same factors.

Note, first, that f(σ) for the monodisperse end members, ξ = 0 and 1, are identical in shape, Fig. 4(a), but with the former shifted to the left by a factor of (d2/d1)2 = α2 = 0.0625 due to a decrease in σ* by the same factor for the larger particles.3 Addition of just 20% of small spheres to a suspension of large spheres (ξ = 0.2) produces a dramatic effect, Fig. 4(a). While frictional contacts still start to form at σ1* ≈ 0.06σ0, the onset is now much more gradual, until σσ2* ≈ σ0, whereupon df/dσ abruptly becomes as large as the monodisperse case (either ξ = 0 or 1), before f saturates at a σ that is comparable to (but slightly later than) that of monodisperse small spheres, even though only 20% of these are present. By ξ = 0.5, f(σ) becomes very similar to that of that of monodisperse small spheres (ξ = 1); although, the onset is still clearly somewhat earlier and the saturation somewhat lower. These features become progressively less obvious at ξ = 0.65 and 0.80. The effect of bidispersity is therefore asymmetrical: the presence of 20% of large spheres in 80% of small spheres has far smaller an effect on f(σ) than the reverse situation.

Turning to the shear-thickening flow curves, Fig. 4(d), we see that, as expected, the WC model reproduces the simulated data for the two monodisperse end members. It gives a tolerable representation of the data at ξ = 0.8, i.e. when there are 20% of large spheres present in a predominantly small-sphere system; but, it fails badly in the reverse situation, when there are 20% of small spheres in a mainly large-sphere system (ξ = 0.2). The disagreement between the WC prediction and simulation data is still substantial at ξ = 0.5, and remains perceptible at ξ = 0.65.

3.3 Comparing experiments to the WC model

Testing the WC model against experimental data is more involved. Fig. 3(a) shows clearly that introducing bidispersity alters shear thickening; however, simulations9 and recent experiments6,27 have shown that, even for nearly monodisperse suspensions, thickening is also sensitive to the relationship between the particle static friction coefficient, μ, and the normal contact force, F. The function μ(F) is fully prescribed in our simulations: μ = 0 below a threshold force F*, and μ = μm > 0 above F*, allowing us to isolate the effect of bidispersity on the shear-thickening phenomenology. However, for our sterically stabilised PMMA particles the load-dependent friction μ(F) has not been measured; hence, we do not know a priori the role of the specific tribology of our particles.

For simplicity, we assume that the experimental μ(F) obeys the same contact model (CLM) as in simulations and treat the critical load, F*, as an unknown parameter; as in simulations, we take F* to be independent of the size of the contacting particles. F* defines a contact stress scale σ0F*/d22, the unit of stress in our simulations. For each σ and ξ in experiments we calculate σ/σ0 and read off the corresponding fraction of frictional contacts, f(σ/σ0), from simulations, Fig. 4(a).

Using σ0 as a global fitting parameter, and ϕm and ϕ0 as local fitting parameters, we use the f(σ/σ0) so obtained to calculate ϕJ(σ) via eqn (2), from which we compute ηWC(σ) with eqn (1). In Fig. 5, we plot measured flow curves (symbols) and WC flow curves (lines) for σ0 = 250 Pa and ϕm(ξ) = Λϕ0(ξ) with Λ = 0.89. Choosing a ξ-dependent Λ = ϕm/ϕ0 does not affect any of our conclusions. Data and fits have been shifted vertically for clarity (see caption).

image file: c9sm00041k-f5.tif
Fig. 5 Failure of the WC model for PMMA spheres. Symbols, shifted flow curves from experiments for different ξ, as labelled. Lines, shifted WC model predictions. Shift factors are: ξ = 0, 0.035; ξ = 0.2, 0.2; ξ = 0.5, 1.8; ξ = 0.65, 5; ξ = 0.8, 10 and ξ = 1, 15.

The all-large (ξ = 0) flow curve is well fit by the WC model, in agreement with our simulations. This justifies a posteriori our use of the CLM for μ(F) for this sample. The model should equally well describe the all-small (ξ = 1) flow curve; however, this is not the case. Although the onset of shear thickening is correctly predicted, the rise in η(σ) is overestimated by the model, implying a different μ(F) for the small particles, e.g., CLM with a lower μm than for the large particles.9 As a consequence, the present map between σ and f is not reliable for our bidisperse suspensions; to calculate f(σ) properly, one would have to independently characterise μ(F) experimentally for the different contact types (large–large, small–small and large–small) and compare with simulations employing a representative DEM contact model. We do not do so here; however, we point out that all the same trends noted when we compared simulations with the WC model are clearly reproduced in our bidisperse experiments ([0 < ξ < 1 in Fig. 5]). In particular, there is a striking disagreement between model and experiment at ξ = 0.2, for which, as in Fig. 4(d), the onset of shear thickening is grossly underestimated. Based on this similarity, we can infer already that the same shortcomings of the WC model as applied to bidisperse simulations should also apply to bidisperse experiments.

4 How the WC model fails

Previous experiments and simulations find, and we confirm, that the WC model works well in the quasi-monodisperse limit (s ≲ 0.2).18 This phenomenological model is designed to reveal the consequences of a simplified set of assumptions in the most perspicuous way, and (according to its authors10) not meant for the fitting of data. Thus, that it works quantitatively in the small-s limit is itself non-trivial, especially given its mean-field nature.22

It is perhaps unsurprising that we find the WC model fails to account for a binary mixture with α = 0.25 (size ratio 1[thin space (1/6-em)]:[thin space (1/6-em)]4), which translates, using a previously-proposed definition,28 to an effective polydispersity seff = (d1d2)/(d1 + d2) ≈ 60%.|| The pertinent question is: precisely where is the WC model failing in this case?

In a monodisperse system, there is a single kind of frictional contact. In a bidisperse system, such contacts are of three kinds: large–large (‘11’), large–small (‘12’), and small–small (‘22’). Fig. 6(a) shows how the three types of frictional contact develop with stress, f11(σ), f12(σ) and f22(σ), in our simulated ξ = 0.2 system, where we observe maximal discrepancy with the WC model. Not surprisingly, frictional contacts first form amongst the large species, at σ1* ≈ 0.06σ0; however, this contribution rapidly saturates to f11 ≲ 0.1. The latest frictional contacts to form are the small–small ones: f22 does not start to increase until σ2* ≈ σ0; but, these saturate to about f22 ≈ 0.4 ≳ 4f11. Ultimately, the largest contribution is from ‘mixed’ contacts, f12 ≈ 0.5 ≳ 5f11, which start to form at σ12* ≈ 0.2σ0 (perhaps fortuitously close to image file: c9sm00041k-t1.tif).

image file: c9sm00041k-f6.tif
Fig. 6 WC model fails due to equal weighting of contact types. (a) Simulated σ-dependent fraction of frictional contacts, f(σ) assumed implicitly by WC in their original model (black line), for ξ = 0.2 decomposed into contributions from large–large, f11; large–small, f12 and small–small, f22, contacts, as labelled. (b) Weighted total fraction of frictional contacts fpoly needed to fit the WC model to our data, eqn (4) (black line), and individual weighted fractions, fpoly11 (red), fpoly12 (green) and fpoly22 (blue), as defined in eqn (5). (c) Symbols, simulation flow curve. Dashed line, prediction of the WC model using f. Solid line, prediction of the WC model using fpoly.

In its original form, the WC model is agnostic to particle size. Consistency with this feature requires that, when applied to our biphasic system, we take f = f11 + f12 + f22, so that any new frictional contact formed as stress builds up contributes equally to the lowering of ϕWCJ, and therefore to the viscosity increment via eqn (1) and (2). Thus, because f12f22f11 at ξ = 0.2, the WC flow curve at this composition is much more similar in shape to that for the all-small (ξ = 1) system than that for the ξ = 0 system. In reality, the simulated flow curves start to shear thicken at σ1*, which is where f11 starts to increase; i.e., large–large contacts dominate despite the smallness of f11, and many small–large and small–small contacts seem to contribute little to the shift in ϕJ.

This suggests that we should write ϕJ a function of f11, f12 and f22, separately. A simple ansatz is to retain the functional form of ϕJ in eqn (2) and replace f with a polydisperse crossover function, fpoly = fpoly11 + fpoly12 + fpoly22, where the weighted fraction of frictional contacts for contacts of type (ij) is fpolyij = κij/fij. The coefficient κij corresponds to the large-σ limit of fpolyij; fij denotes the corresponding large-σ limit of fij in Fig. 6(a). We choose κ11 + κ12 + κ22 = 1 to ensure fpoly(σ/σ0 ≪ 1) = 0 and fpoly(σ/σ0 ≫ 1) = 1. So, our extended WC model reads

ϕJ = fpolyϕm + (1 − fpoly)ϕ0, (4)
image file: c9sm00041k-t2.tif(5)
We continue to use eqn (1) for the relative viscosity.

The weighted fraction of frictional contacts of type (ij), fpolyij(σ) needed to fit our data, shown in Fig. 6(b) for ξ = 0.2, has the same shape as fij(σ), Fig. 6(a), but scaled up by a factor of κij, which sets the limiting value of fpolyij as σ/σ0 → ∞. [Note that the slight non-monotonicity of f11(σ) in Fig. 6(a) (red) means that fpoly exceeds 1 using our normalisation scheme.] By varying the free parameters κ11 and κ12 (with κ22 = 1 − κ11κ12), we can readily fit all of our bidisperse simulation flow curves; Fig. 6(c) (solid line) shows the best-fit curve, obtained by eye, for ξ = 0.2. In Fig. 7(a), we plot the best-fit coefficients for ξ = 0.2, along with the coefficients for other ξ (for the full fits, see the ESI).

image file: c9sm00041k-f7.tif
Fig. 7 (a) Coefficients, κ11 (red), κ12 (green) and κ22 (blue), obtained from fitting bidisperse simulation flow curves using fpoly, eqn (5), as a function of ξ. (b) Fractions of frictional contacts in the large-stress limit, f11 (11) (red), f12 (green) and f22 (blue), as a function of ξ.

For comparison, we consider the special case in which κ11 = f11, κ12 = f12 and κ22 = f22, so that fpoly reduces to f, the unweighted fraction of frictional contacts, Fig. 4(a), and the original WC model is recovered, Fig. 6(c) (dashed line). We plot the coefficients for this case, f11(ξ), f12(ξ) and f22(ξ), in Fig. 7(b).

For ξ = 0.2, our set of fitted {κij}, Fig. 7(a), reaffirms quantitatively what we proposed qualitatively earlier. The largest contribution to changes in fpoly (and hence ϕJ) is from large–large contacts, for which κ11 = 0.7, while there is a negligible contribution from small–small contacts, κ22 ≈ 0. In contrast, in the original WC model κ11 = f11 ≈ 0.1 and κ22 ≈ 0.4, Fig. 7(b). Increasing ξ to 0.5 sees the increasing importance of large–small contacts and decreasing importance of large–large contacts, while small–small contacts remain irrelevant. At ξ = 0.65, fpoly is determined almost entirely by large–small contact formation. Only for ξ = 0.8 do small particles have a measurable contribution, where they enter on equal footing with large–small contacts; large–large contacts are, by this point, irrelevant.

5 Discussion

The WC coefficients, fij, in Fig. 7(b) correspond to the relative numbers of each kind of contact (‘11’, ‘12’ or ‘22’) in the high-σ/σ0 limit. Thus, the ratio Δijκij/fij measures the relative contribution to ϕWCJ (and hence to η) due to the formation of a single frictional contact of type ij. In the WC model, all contact types give rise to the same increment in ϕJ and Δ11 = Δ12 = Δ22 = 1. In Fig. 8, we plot the Δs as a function of ξ (colours) and overlay the WC prediction (horizontal dashed line). Strikingly, Δ11Δ12Δ22 for all our bidisperse mixtures. Thus, at ξ = 0.2, for example, a large–large contact contributes over an order of magnitude more than a large–small one; while, the effect of forming a small–small contact is negligible (Δ22 ≈ 0); i.e., small particles are largely redundant for stress transmission. Only at ξ = 0.8, where the fraction of small particles is largest and there are no large–large contacts (Δ11 is undefined here; so, we do not plot it), do small–small contacts have an appreciable contribution (Δ22 becomes non-zero). Even then, a single large–small contact contributes the same as ∼[scr O, script letter O](10) small–small contacts.
image file: c9sm00041k-f8.tif
Fig. 8 ϕJ-Contribution per contact Δ = κij/fij for ‘11’ (red), ‘12’ (green) and ‘22’ (blue) contacts, obtained by taking the ratio of the data in Fig. 7(a and b). The horizontal dashed line is the WC prediction, Δ11 = Δ12 = Δ22 = 1.

Importantly, none of the bidisperse Δs follow the WC prediction (dashed line). Since different contact types do not contribute equally to changes in ϕWCJ, one needs to know not only the total number of frictional contacts, but the sizes of the particles participating in all those contacts to predict η. Thus, f, which assumes that all frictional contacts contribute the same, is inherently unsuitable for modelling bidisperse systems.

The particle-size-dependence in Fig. 8 is reminiscent of sheared polydisperse dry granular packings, in which stress transmission is strongly spatially heterogeneous with contacts between larger particles carrying a higher load on average than those between smaller particles.28 In bidisperse dry granular systems, the detailed partition of stress is sensitive to both ξ and the size ratio α; e.g., as the size disparity grows (α decreases), the contribution of contacts involving small particles progressively decreases.29 While this problem has been studied at length in dry systems under imposed particle pressure (⇒ varying ϕ),30–33 it has received relatively little attention for fixed-ϕ, granular suspensions, in which particles interact not only through contact-, but also hydrodynamic forces. Presumably, the trend with α is similar to the dry-grain case, so that, as α → 1, the disparity between different contact types vanishes, i.e., Δ11Δ12Δ22 = 1, and hence f eventually becomes a reasonable approximation for fpoly, which would explain the success of the WC model for weakly bidisperse mixtures.18 Clearly, focussed work is needed to understand the details and origins of stress partitioning, and its relation to shear thickening, before a realistic model can be constructed. Minimally, the relative weights of different contact types should be allowed to vary, e.g., like our eqn (5).

Alternative existing models of frictional shear thickening could prove more successful in capturing bidisperse flowcurves. A recent idea is that thickening is not driven by the formation of frictional contacts per se, but by the changes in anisotropy of stress transmission that this induces.**[thin space (1/6-em)]34,35 In particular, Thomas et al. proposed an ab initio model for two-dimensional systems based on the ratio of the shear stress and particle pressure, σxz/P, which they relate to the anisotropy of contact forces. Interestingly, for simulations of dry grains with a uniform particle size distribution, σxz/P is found to be independent of polydispersity (controlled via the span of the distribution),28,33 which suggests that their approach may account for polydispersity effects naturally in a way that the WC model, which is agnostic to the spatial distribution of contact forces, does not. This merits a thorough study of the role of stress anisotropy in bidisperse systems, and the extension of Thomas et al.'s theory to 3-d and polydisperse systems.

Perhaps as important as the contact-type-dependent contributions to ϕJ is our observation that the original WC model fits our monodisperse simulation data. This result is non-trivial: it implies that the microphysics of shear thickening can be captured by a single scalar parameter (f) that is agnostic to the spatial distribution of contacts. Remarkably, there is evidence in the literature that this scenario may be true, at least in the quasi-monodisperse limit. By simulating a weakly bidisperse mixture of particles interacting via the critical-load model and the same mixture containing particles with different but load-independent μ, Dong and Trulsson37 showed that η is uniquely defined by f, even though the microstructure for both setups is very different. For strongly bidisperse suspensions, the roles of microstructure and stress partitioning in shear thickening remain to be disentangled. If they bear similarity to polydisperse dry grains, then the two should be strongly correlated. Specifically, we would expect big particles, which carry the largest loads, to align with the compressive axis; whereas small particles, which carry a negligible load, would form an almost isotropic background of “spectator” particles.28,38 Studying the spatial distribution of contacts systematically in these systems, e.g., in the vein of Dong and Trulsson,37 would shed light on this issue and help to establish whether the independence to microstructure in the monodisperse case has deeper physical meaning, or if it is entirely fortuitous.

Before concluding, we comment on the utility of experimental data in testing the WC model. Although our experiments qualitatively support the notion that the WC model fails for bidisperse suspensions, the inability of the model to describe the all-small η(σ) based on simulation f(σ) highlights an important obstacle to rigorous testing: the F-dependent tribology of interparticle contacts is a priori unknown. Our work with PHSA-stabilised PMMA dispersions suggests that, even for particles with ostensibly well-controlled surface properties, μ(F) may vary considerably from batch to batch. For example, in the ESI we show flow curves for quasi-monodisperse PMMA spheres showing “two-stage” shear-thickening, with two distinct onset stresses. Such behaviour clearly cannot arise from the single-stress-scale CLM used here. Thus, microtribology experiments6,39,40 must play a central role in future model testing. Indeed, the scarce tribology measurements that already exist for sterically-stabilised particles indicate a rich behaviour, particularly at large normal loads.27 Finally, we note that even in experimental systems where μ(F) is described by the CLM, “fitting” the WC model to experimental data with a presumed form for f(σ) will result in a f(σ) that is not the fraction of frictional contacts except in the monodisperse limit.3,5

6 Summary and conclusions

We have simulated and experimentally measured the rheology of a bidisperse suspension of repulsive spheres to test the validity of the WC model of shear thickening. By using the fraction of frictional contacts f extracted directly from simulations, we showed that the WC model works in the special case of monodisperse particles, but is incomplete when applied to bidisperse mixtures. While our study focusses on continuous shear thickening, we expect all the same conclusions to apply at higher volume fractions, where discontinuous shear thickening is observed.

In practical terms, our results suggest caution when using the WC model as anything other than an empirical fitting tool. Specifically, little, if any, meaning can be ascribed to f extracted from fits to flow curves. On a fundamental level, our work highlights the need for a focussed effort to understand the link between σ-dependent frictional contact formation and dissipation. Existing studies of shear thickening consider either bulk rheology3,5 or ex situ, two-particle properties,6 with little or no concerted effort to bridge the two regimes. We believe that any serious effort to make this link should consider polydispersity from the outset in its own right, rather than merely as a means of mitigating crystallisation. Indeed, our work hints that the monodisperse limit is a singular one, and so cannot be used as a guide to developing models for the flow of polydisperse systems.

Data from this article are available at

Conflicts of interest

There are no conflicts to declare.


BMG and MH were funded by EPSRC EP/J007404/1. CJN was funded by EPSRC EP/N025318/1 and the Maudslay-Butler Research Fellowship at Pembroke College, Cambridge. LJS was funded by EPSRC SOFI CDT (EP/L015536/1). JS was funded by EPSRC EP/N025318/1 and The Royal Academy of Engineering/The Leverhulme Trust Senior Research Fellowship LTSRF1617/13/2. WCKP was funded by EPSRC EP/J007404/1 and EP/N025318/1. We thank Andrew Schofield for synthesising the particles, and John Royer, Dan Hodgson and an anonymous referee for helpful discussions. The simulation makes use of the LF-DEM code published in Mari et al.9 and available at as well as LAMMPS.41

Notes and references

  1. H. A. Barnes, J. Rheol., 1989, 33, 329–366 CrossRef CAS.
  2. E. Brown and H. M. Jaeger, Rep. Prog. Phys., 2014, 77, 046602 CrossRef PubMed.
  3. B. M. Guy, M. Hermes and W. C. K. Poon, Phys. Rev. Lett., 2015, 115, 088304 CrossRef CAS PubMed.
  4. N. Y. C. Lin, B. M. Guy, M. Hermes, C. Ness, J. Sun, W. C. K. Poon and I. Cohen, Phys. Rev. Lett., 2015, 115, 228304 CrossRef PubMed.
  5. J. R. Royer, D. L. Blair and S. D. Hudson, Phys. Rev. Lett., 2016, 116, 188301 CrossRef PubMed.
  6. J. Comtet, G. Chatté, A. Niguès, L. Bocquet, A. Siria and A. Colin, Nat. Commun., 2017, 8, 15633 CrossRef CAS PubMed.
  7. C. Clavaud, A. Bérut, B. Metzger and Y. Forterre, Proc. Natl. Acad. Sci. U. S. A., 2017, 114, 5147–5152 CrossRef CAS PubMed.
  8. R. Seto, R. Mari, J. F. Morris and M. M. Denn, Phys. Rev. Lett., 2013, 111, 218301 CrossRef PubMed.
  9. R. Mari, R. Seto, J. F. Morris and M. M. Denn, J. Rheol., 2014, 58, 1693–1724 CrossRef CAS.
  10. M. Wyart and M. E. Cates, Phys. Rev. Lett., 2014, 112, 098302 CrossRef CAS PubMed.
  11. D. J. Hodgson, M. Hermes and W. C. K. Poon, 2015, arXiv preprint, arXiv:1507.08098.
  12. E. Brown, H. Zhang, N. A. Forman, B. W. Maynor, D. E. Betts, J. M. DeSimone and H. M. Jaeger, J. Rheol., 2010, 54, 1023–1046 CrossRef CAS.
  13. L. C. Hsiao, S. Jamali, E. Glynos, P. F. Green, R. G. Larson and M. J. Solomon, Phys. Rev. Lett., 2017, 119, 158001 CrossRef PubMed.
  14. C.-P. Hsu, S. N. Ramakrishna, M. Zanini, N. D. Spencer and L. Isa, Proc. Natl. Acad. Sci. U. S. A., 2018, 115, 5117–5122 CrossRef CAS PubMed.
  15. I. Peters, S. Majumdar and H. Jaeger, Nature, 2016, 532, 214 CrossRef CAS PubMed.
  16. A. Fall, F. Bertrand, D. Hautemayou, C. Mézière, P. Moucheront, A. Lematre and G. Ovarlez, Phys. Rev. Lett., 2015, 114, 098301 CrossRef CAS PubMed.
  17. M. Hermes, B. M. Guy, W. C. K. Poon, G. Poy, M. E. Cates and M. Wyart, J. Rheol., 2016, 60, 905–916 CrossRef CAS.
  18. A. Singh, R. Mari, M. M. Denn and J. F. Morris, J. Rheol., 2018, 62, 457–468 CrossRef CAS.
  19. P. D'Haene and J. Mewis, Rheol. Acta, 1994, 33, 165–174 CrossRef.
  20. J. Bender and N. J. Wagner, J. Rheol., 1996, 40, 899–916 CrossRef CAS.
  21. B. M. Guy, PhD thesis, The University of Edinburgh, 2017.
  22. M. E. Cates, Annal. Henri Poincaré Suppl. 2, 2003, 4, S647–S661 CrossRef.
  23. C. Ness and J. Sun, Soft Matter, 2016, 12, 914–924 RSC.
  24. A. Sierou and J. Brady, J. Rheol., 2002, 46, 1031–1056 CrossRef CAS.
  25. R. Farris, Trans. Soc. Rheol., 1968, 12, 281–301 Search PubMed.
  26. S. Pednekar, J. Chun and J. F. Morris, J. Rheol., 2018, 62, 513–526 CrossRef CAS.
  27. G. Chatté, J. Comtet, A. Niguès, L. Bocquet, A. Siria, G. Ducouret, F. Lequeux, N. Lenoir, G. Ovarlez and A. Colin, Soft Matter, 2018, 14, 879–893 RSC.
  28. C. Voivret, F. Radjai, J.-Y. Delenne and M. S. El Youssoufi, Phys. Rev. Lett., 2009, 102, 178001 CrossRef CAS PubMed.
  29. T. Shire, C. O'sullivan and K. Hanley, Granular Matter, 2016, 18, 52 CrossRef.
  30. J. Gray and A. Thornton, Proc. R. Soc. A, 2005, 461, 1447–1473 CrossRef CAS.
  31. J. Gray and V. Chugunov, J. Fluid Mech., 2006, 569, 365–398 CrossRef.
  32. T. Weinhart, S. Luding and A. R. Thornton, AIP Conf. Proc., 2013, 1202–1205 CrossRef.
  33. D. Cantor, E. Azéma, P. Sornay and F. Radjai, Phys. Rev. E, 2018, 98, 052910 CrossRef CAS.
  34. J. E. Thomas, K. Ramola, A. Singh, R. Mari, J. F. Morris and B. Chakraborty, Phys. Rev. Lett., 2018, 121, 128002 CrossRef CAS PubMed.
  35. M. Otsuki and H. Hayakawa, 2018, arXiv preprint arXiv:1810.03846.
  36. J. F. Brady and G. Bossis, J. Fluid Mech., 1985, 155, 105–129 CrossRef.
  37. J. Dong and M. Trulsson, Phys. Rev. Fluids, 2017, 2, 081301 CrossRef.
  38. F. Radjai, D. E. Wolf, M. Jean and J.-J. Moreau, Phys. Rev. Lett., 1998, 80, 61 CrossRef CAS.
  39. L. O. Gálvez, S. de Beer, D. van der Meer and A. Pons, Phys. Rev. E, 2017, 95, 030602 CrossRef PubMed.
  40. N. M. James, C.-P. Hsu, N. D. Spencer, H. M. Jaeger and L. Isa, J. Phys. Chem. Lett., 2019, 10, 1663–1668 CrossRef CAS PubMed.
  41. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19 CrossRef CAS.


Electronic supplementary information (ESI) available: For details on experimental methods and artefacts; an example of complex shear thickening of monodisperse spheres; and fits of the extended WC model for all our simulated binary mixtures. See DOI: 10.1039/c9sm00041k
By number, small particles dominate the large particles for all the ξ we study. For the data in Fig. 3(a), the number fractions of small particles, x = 1/[1 + α3(1/ξ − 1)], are: x = 0 (pentagon), 0.941 (□), 0.985 (▽), 0.992 (△), 0.996 (○) and 1 (⋄).
§ Thus, the viscosity of the small spheres is greater than that of the large spheres below the onset of thickening, e.g., at σ = 1 Pa.
We expect this to be a reasonable approximation, since eqn (1) has previously been used to fit η(ϕ) for various frictional, bidisperse sphere mixtures.26
|| A more natural definition would normalise to the average size to give s ≈ 120%.
** This is reminiscent of the “hydrocluster” driven thickening of lubricated spheres observed in Stokesian dynamics simulations,36 which, however, is distinct from the contact-driven thickening we observe.4

This journal is © The Royal Society of Chemistry 2019