Discrete-particle simulations of bidisperse shear thickening suspensions are reported. The work considers two packing parameters, the large-to-small particle radius ratio ranging from δ=1.4 (nearly monodisperse) to δ=4, and the large particle fraction of the total solid loading with values ζ=0.15, 0.5, and 0.85. Particle-scale simulations are performed over a broad range of shear stresses using a simulation model for spherical particles accounting for short-range lubrication forces, frictional interaction, and repulsion between particles. The variation of rheological properties and the maximum packing fraction ϕJ with shear stress σ are reported. At a fixed volume fraction ϕ, bidispersity decreases the suspension relative viscosity ηr=ηs/η0, where ηs is the suspension viscosity and η0 is the suspending fluid viscosity, over the entire range of shear stresses studied. However, under low shear stress conditions, the suspension exhibits an unusual rheological behavior: the minimum viscosity does not occur as expected at ζ0.5, but instead decreases with further increase of ζ to 0.85. The second normal stress difference N2 acts similarly. This behavior is caused by particles ordering into a layered structure, as is also reflected by the zero slope with respect to time of the mean-square displacement in the velocity gradient direction. The relative viscosity ηr of bidisperse rate-dependent suspensions can be predicted by a power law linking it to ϕJ, ηr=(1ϕ/ϕJ)2 in both low and high shear stress regimes. The agreement between the power law and experimental data from literature demonstrates that the model captures well the effect of particle size distribution, showing that viscosity roughly collapses onto a single master curve when plotted against the reduced volume fraction ϕ/ϕJ.

Dense suspensions are found in a number of industrial (cement, ceramics, concrete casting) and geophysical (muds, sediments, lava flows, slurries) settings. Controlling the flow of dense suspensions is of major importance. In the cement industry, for instance, it is desirable for concrete to be of high concentration while being able to flow easily [1].

Particles suspended in fluid often cause non-Newtonian flow behavior, where rheological properties become dependent not only on particle volume fraction ϕ but also on shear rate γ˙ or shear stress σ [2]. Non-Brownian dense suspensions, which we study here, are known to undergo a range of rheological behaviors, such as shear thinning, shear thickening, shear jamming, and normal stress differences [3,4]. Shear thickening, a phenomenon where viscosity of the suspension increases with shear rate, can vary from a mild form when viscosity increases continuously (continuous shear thickening, or CST), to severe thickening, when the increase in viscosity occurs abruptly (discontinuous shear thickening, or DST), sometimes even by orders of magnitude [5,6]. Accurately accounting for the conditions that give rise to shear thickening has been the focus of many studies. Yet, predicting the appearance of shear thickening behavior remains a fundamental challenge in this field [2], especially for polydisperse suspensions, even in the simplest case of bidispersity. Only limited studies have explored rate dependence for bidisperse systems, and this is particularly true with respect to recent understanding based on considerations of frictional interactions in suspension rheology. In the next section we provide background on the lubricated-to-frictional transition scenario, followed by a brief review of relevant work on bidisperse suspension rheology.

Recent efforts to explain the mechanism of shear thickening have credited the shear-induced transition from lubricated (frictionless) to frictional contact between particles [5,7–10]. In the simulation studies of Seto et al. [7] and Mari et al. [5], the transition between the lubricated and frictionally dominated states is defined by the characteristic stress σ0=FR/6πa2, where FR is a repulsive interparticle force and a is the particle radius. At low shear stress (σσ0), the interactions of particles are governed by lubricational forces. The repulsive forces, which maintain the separation of particle surfaces, lead to a larger apparent particle size, which results in shear thinning rheology. In the low stress regime, the relative viscosity (defined ηr=ηs/η0, where ηs is the suspension viscosity and η0 is the suspending fluid viscosity) diverges as ϕϕJ, with ϕJ=ϕJ0 being the frictionless maximum packing fraction, corresponding to a limiting value of volume fraction for purely lubricated interactions. As shear stress increases and reaches the critical value (comparable to σ0), the repulsive forces acting between the particles are overwhelmed and particle contacts proliferate. In the high shear stress regime (σσ0), the interactions between the particles are mostly frictional and the viscosity diverges at the frictional maximum packing fraction [5,7,8,11] ϕJμ<ϕJ0. We note that, in various literature studies, the maximum packing fraction, or jamming fraction, is denoted as ϕJ, ϕmax, or ϕm. For consistency with the work of Mari et al. [5] and Singh et al. [12] we use ϕJ, but use the term maximum packing fraction. The relative viscosity divergence is taken in the classical Maron–Pierce [13] form ηr(1ϕ/ϕJ)2 used by Wyart and Cates [8]. Singh et al. [12] found this form to fit well most of the data they generated from nearly monodisperse particle simulations using the method applied here. However, Ramaswamy et al. [9] recently demonstrated that this form systematically deviates, to an exponent smaller in magnitude at 3/2, from experimental data in the near-jamming frictional regime.

The rheological model proposed by Wyart and Cates [8] predicts shear thickening based on a microscopic parameter f, the fraction of frictional contacts, that is taken as a function of particle pressure [14] Π (=[Σ11+Σ22+Σ33]/3), as was also assumed by Bashkirtseva et al. [15] Since Π and σ for non-colloidal suspensions are related [16], f may also be taken as a function of shear stress. According to Wyart and Cates, ϕJ is a function of f and decreases from ϕJ0 in the frictionless state (f=0) to ϕJμ in the frictional state (f1). Based on more recent experimental work [9], it is argued that the fraction of frictional contacts may depend on both ϕ and stress. Singh et al. [12] extended the Wyart and Cates model to a constitutive law, which predicts well viscosity, the second normal stress difference N2 and particle pressure based on three scaling parameters: ϕ, σ, and interparticle friction coefficient μ. Testing the Wyart and Cates model on bidisperse suspensions, Guy et al. [17] showed that binary mixtures of spheres with large-to-small particles size ratio δ=4 fail to follow the model, which as they note is “blind” to distinguish between the types of contacts in bidisperse systems. The different types of particle contacts (large-large, large-small, and small-small) apparently contribute differently to the development of stress. One sees a dependence on the level of the size difference, as the findings of Guy et al. regarding the shortcomings of the model at δ=4 are clearly different from the conclusions of Singh et al. [12], who performed simulations with δ=1.4 and validated the Wyart and Cates model. This leaves open the question of the effect of particle size distribution on rate dependence of the rheology.

As seen from the brief review above and more extensive discussions elsewhere [18,19], the majority of studies available on rate-dependent suspensions have been focused on nearly monodisperse systems. In a number of materials, for example, in cement [20], crystal-bearing magma [21], muds [22], ceramic precursors [23], hydrocolloid gums [24], or chocolate during the refining stages [25], however, the particles are not identical and exhibit a certain size distribution among other complexities such as particle shape. Here, we limit ourselves to spherical particles and the simplest of polydisperse suspensions by considering the bidisperse case. A recent study has demonstrated that polydisperse suspensions can be represented well by an equivalent binary system [26,27].

In this work, we perform a study of rate-dependent bidisperse suspensions, for which the rheology picture is incomplete. We first introduce the parameters used to describe the bidisperse systems and provide background on known effects of bidispersity on the rheology of suspensions.

In addition to the solid volume fraction ϕ used to describe monodisperse systems, two other parameters are needed to characterize a binary mixture. We choose these to be a particle size ratio δ=al/as, where al and as are the radii of the large and small particles, respectively, and the fraction of large particles ζ=ϕl/ϕ, where ϕl is the volume fraction of large particles. The works in this section discuss the effect of bidispersity on the rheology of rate-independent suspensions, except for Bender and Wagner [28] who also consider shear stress dependence.

It is well known that bidispersity significantly affects the suspension rheology. For instance, a bidisperse suspension of spheres has lower relative viscosity as compared to the monodisperse suspension at the same solid volume fraction. This can be attributed to the ability to achieve a larger maximum packing fraction ϕJ in bidisperse suspensions, as shown by numerous experimental [29,30] and numerical [27,31] studies. Specifically, the relative viscosity is well-described as a function of the reduced volume fraction ϕ/ϕJ, with ϕJ being an empirically determined function of the packing parameters [32] and contact friction. [5,12] Using the simulation model employed here, Pednekar et al. [27] showed a viscosity collapse for the high-stress limit of polydisperse suspensions of frictionally-interacting particles as well as several datasets from the literature. They found that the second normal stress difference of the simulations was similarly well-described for nearly monodisperse and polydisperse suspensions with different packing parameters when plotted as a function of ϕ/ϕJ. The collapse of all data regardless of the packing parameters suggests that ϕ/ϕJ is one of the most important parameters in suspension flow with ϕJ encoding information on the particle size distribution. Note that the work of Pednekar et al. [27] did not consider a stress-dependent jamming fraction, and this will be a central consideration in the current work.

Early work on multimodal suspensions of rigid particles was done by Farris [33], who proposed a model for the calculation of the viscosity of a bidisperse suspension based on monodisperse viscosity data. Assuming no interaction between small and large particles, the relative viscosity of a bidisperse suspension can be represented as the product of relative viscosities of each particle size. Chong et al. [29], working with glass beads of various δ and ζ suspended in polyisobutylene, showed that the Farris model works for δ of up to 10, but found further increase of δ to decrease the viscosity in an inconsistent way. Investigating the effect of bidispersity on the rheology of suspensions of various ϕ, these authors observed that small changes in ζ caused the viscosity to vary significantly with the minimum viscosity being achieved at ζ=0.650.75. The most pronounced drop in viscosity compared to a monodisperse system was observed for the largest solid volume fraction (ϕ=0.65) with the largest particle size ratio studied (δ=22.22). Their relative viscosity data plotted against the reduced volume fraction ϕ/ϕJ collapsed onto a single curve, leading the authors to suggest that the relative viscosity of suspensions is only a function of ϕ/ϕJ. This concept was corroborated by Chang and Powell [34], who worked with the results of several experimental and numerical studies.

Studying glass spheres in glycerine, Shapiro and Probstein [35] demonstrated that the viscosity of bidisperse suspensions is a decreasing function of the maximum packing fraction, which, in turn, is dependent on the particle size distribution. The most pronounced increase of maximum packing, and, associated decrease of the viscosity, was observed at ζ=0.75 for the largest δ studied (δ=4). These trends were further supported by Chang and Powell [34], as well as by Gondret and Petit [36] who interpreted the increase of dynamic viscosity through the decrease of ϕJ. The above findings show that the increase of ϕJ in bidisperse suspensions relative to monodisperse suspensions and thus, at a given ϕ, increase of the “distance from jamming” ϕJϕ, is the basis for the lower viscosity of the bidisperse material. Since defining ϕJ for any system provides a valuable insight into how the system will behave, we will use the maximum packing fraction to frame our results.

Bender and Wagner [28], working with monodisperse and bidisperse sheared suspensions of colloidal silica spheres (ranging from 160 to 400 nm in diameter), showed the dependence of viscosity on ζ and σ. The viscosity noticeably decreased when just a small fraction of smaller particles was introduced to the suspension of larger particles. For instance, the viscosity of a ϕ=0.64 bidisperse suspension with ζ=0.90 and δ=2.1 is one-fifth that of a monodisperse suspension of the same concentration. By decreasing ζ further to 0.50, the reduction in viscosity became even stronger: compared to ζ=0.90, it decreased by more than an order of magnitude at high shear stress and a factor of 4 when moving to lower shear stress. Also, the onset of shear thickening shifts to higher values of shear stress for bidisperse systems. The simulation work of Pednekar et al. [27] explored the effect of δ and ζ of concentrated bidisperse suspensions in the high shear limit, and showed that the viscosity reaches its minimum at ζ0.65.

In our study on rate-dependent bidisperse suspensions, we investigate how particle size ratio δ and large particle fraction ζ affect the rheology over a wide range of shear stresses, using a simulation tool that we describe in Sec. II A. The simulation results are discussed in Sec. III. We demonstrate the relative viscosity dependence on ϕ, δ, ζ, and σ and interpret the set of results for various bidispersity parameters through the values of maximum packing. We present a surprising gradual decrease of viscosity as large particle fraction increases (within the range of 00.85 considered in the current work) under low shear rate conditions and attribute this to the suspended particles undergoing ordering. Our conclusions follow in Sec. IV.

We simulate 3D suspensions using the lubrication flow—discrete element model (LF-DEM) [7], which has been found to successfully describe the flow of dense suspensions [5] and closely reproduce the essential features of DST and shear jamming seen experimentally [37]. We simulate neutrally buoyant particles of different radii, as for the smaller and al for the larger particles, immersed in a viscous incompressible Newtonian fluid of density ρ and viscosity η0 under an imposed shear stress σ. We consider motion in the limit of vanishing particle Reynolds number, Re=ργ˙a2/η00, such that inertia can be neglected. Hence, the equation of motion reduces to the balance between hydrodynamic FH, repulsive FR, and contact FC forces (and torques) on each particle, 0=FH+FR+FC.

The hydrodynamic interactions are written as

FH=RFU(UU)+RFE:E,

where U is the particle (translational and angular) velocity and U=γ˙ye^x is the imposed shear flow. RFU and RFE are resistance matrices corresponding to short-range lubrication forces [38] on the particles due to their motion and to the imposed shear flow, respectively; long-range hydrodynamics are neglected. Here, E=(e^xe^z+e^ze^x)γ˙/2 is the rate-of-strain tensor. The dominant squeeze mode resistance to motion along the line of centers for the pair, RFU, scales as (h+ϵ)1, while the resistance to tangential motions scales as log[1/(h+ϵ)], where h=2(raiaj)/(ai+aj) is the dimensionless interparticle gap; ϵ=103 is included such that the lubrication resistance remains finite at contact (h=0), which could be interpreted as representative of the influence of particle surface roughness. Contact may occur when the applied stress overcomes the repulsive forces. The tangential and normal contact forces satisfy the Coulomb law FC,tanμFC,nor with components of the contact force computed using linear spring models as FC,tan=ktξ and FC,nor=knhn, where ξ is the tangential spring stretch vector, which is initialized to zero when a pair first makes contact (i.e., when h=0), and n is the normal vector along the pair line of centers. A damping is applied to the normal contact motion, as described in the full elaboration of the model [5]. In this study, the friction coefficient is taken as μ=1; when the Coulomb threshold is exceeded, i.e., when FC,tan>μFC,nor=FC,nor, the contact point can slide with the same friction coefficient. A repulsive electrical double layer (EDL) force maintains the particle surface separation at low stress. This force decays exponentially with the interparticle gap h as FR=FR0eκh, where κ1 is the Debye length, taken in this work as κ1=0.05as. The repulsive force introduces a stress scale σ0=FR/6πas2. This form of the simulation is called the electrostatic repulsion model (ERM) [5]. When the applied stress σσ0, the interactions between the particles are lubricated (frictionless). When σσ0, the interactions are predominately frictional contacts.

At each time step, we evaluate FR and FC and solve the equation of motion for the particle velocities, UU=RFU1(RFE:E+FR+FC). The model employs Lees–Edwards boundary conditions and is periodic in all directions. The flow, velocity gradient, and vorticity directions are x, z, and y, respectively. We omit Brownian motion, i.e., the Péclet number Pe=η0γ˙a3/kT.

We simulate bidisperse suspensions of non-Brownian spherical frictional particles with a total of N=1000 in a unit cell. The total solid volume fraction range is 0.54ϕ0.60. The particle size ratio is 1.4δ4 and the fraction of large particles occupying the solid volume is ζ=0.15, 0.50, or 0.85. To help visualize the particle configurations, Fig. 1 shows a ϕ=0.57, δ=3 suspension at ζ=0.15, 0.50, and 0.85. The scale for shear stress is σ0=FR/6πas2, and hence the scale for shear rate is defined as γ˙0=FR/6πη0as2. Reported shear stresses and shear rates will be nondimensionalized by σ0 and γ˙0, respectively. We report the suspension viscosity ηs in the form of relative viscosity ηr=ηs/η0, while normal stresses are scaled by η0γ˙.

FIG. 1.

Packings of simulated bidisperse suspensions for various ζ=ϕl/ϕ=0.15, 0.50, and 0.85 with δ=al/as=3. The total volume fraction is ϕ=0.57. (a) ζ=0.15, (b) ζ=0.50, (c) ζ=0.85.

FIG. 1.

Packings of simulated bidisperse suspensions for various ζ=ϕl/ϕ=0.15, 0.50, and 0.85 with δ=al/as=3. The total volume fraction is ϕ=0.57. (a) ζ=0.15, (b) ζ=0.50, (c) ζ=0.85.

Close modal

The stress-controlled simulations were performed over a range of dimensionless shear stress of 0.01σ/σ0562. Each simulation in this study was run over a period of γ˙t=30 strain units. For reporting the rheological properties, we systematically discard the initial five strain units in each data set to account for the transient period as local structure develops, and report means based on a uniform sampling at strain steps of 0.01 of the remaining data. This accounts fully for the transient in almost all cases. However, the transient period extends longer for ζ=0.85 suspensions at low stress (σ/σ0<0.3) due to a large-scale ordering to be described in detail below, and while the bulk of the rheological property change occurs in the first five strain units, this procedure results in a modest (<10%) overestimate of the viscosity for these conditions.

In this section, we present the simulated viscosity data for all systems studied. Data for δ=1.4, ζ=0.50 were adopted from Mari et al. [5].

In Fig. 2, we plot the simulated relative viscosity ηr as a function of dimensionless shear stress σ/σ0 and dimensionless shear rate γ˙/γ˙0 for suspensions with δ=3. Large particle fractions ζ=0.15, 0.50, and 0.85 at several solid volume fractions ϕ0.54 were studied. In all cases, the flow curves share the pattern of shear thinning at low shear stress followed by shear thickening beginning at 0.1σ/σ01. As expected, the viscosity increases with increasing volume fraction. The onset of DST appears to occur at ϕ=0.57 for the systems with ζ=0.15 and ζ=0.85 and at ϕ>0.60 for ζ=0.50. This corresponds to the start of nonmonotonic behavior of the flow curves in Fig. 3 where γ˙/γ˙0 is plotted against σ/σ0. The critical shear stress, which defines the onset of shear thickening, is at σ/σ0=0.3 for ζ=0.15 and ζ=0.50 and σ/σ0=0.1 for ζ=0.85 and roughly constant for given ζ across all volume fractions studied [Fig. 2(a), 2(c), and 2(e)]. The second normal stress difference N2=Σ22Σ33=ΣzzΣyy is always found to be negative here (Fig. 4), which agrees with prior work [39]. It qualitatively follows the same dependence as the viscosity in suspensions, showing a decrease in magnitude when shear thinning occurs, followed by an increase upon shear thickening. In this work, we focus our attention in consideration of the normal stresses on N2. While reliable experimental data on N1 for bidisperse suspensions have been obtained [40], the first normal stress difference N1=Σ11Σ22=ΣxxΣzz is known to be noisy and difficult to reliably simulate under bidisperse conditions [27], a point that may be related to its tendency in shear thickening dense suspensions to change sign from negative to positive [41] with increasing ϕ and σ. We do, however, present N1 data in Fig. 5(b) to provide an indication of the scale and variability of this quantity based on the simulations.

FIG. 2.

Relative viscosity as a function of (a), (c), and (e) shear stress and (b), (d), and (f) shear rate for particle size ratio of δ=3 at varying solid volume fraction and large particle fraction of (upper row) ζ=0.15, (middle row) ζ=0.50, (bottom row) ζ=0.85. The legend in panel c is used in all plots. (a) ζ=0.15. (b) ζ=0.15. (c) ζ=0.50. (d) ζ=0.50. (e) ζ=0.85. (f) ζ=0.85.

FIG. 2.

Relative viscosity as a function of (a), (c), and (e) shear stress and (b), (d), and (f) shear rate for particle size ratio of δ=3 at varying solid volume fraction and large particle fraction of (upper row) ζ=0.15, (middle row) ζ=0.50, (bottom row) ζ=0.85. The legend in panel c is used in all plots. (a) ζ=0.15. (b) ζ=0.15. (c) ζ=0.50. (d) ζ=0.50. (e) ζ=0.85. (f) ζ=0.85.

Close modal
FIG. 3.

Shear rate as a function of shear stress for particle size ratio of δ=3 at varying solid volume fraction and large particle fraction of (a) ζ=0.15, (b) ζ=0.50, and (c) ζ=0.85. The legend in panel b is used in all plots. (a) ζ=0.15. (b) ζ=0.50. (c) ζ=0.85. https://doi.org/10.1122/8.0000495.1

FIG. 3.

Shear rate as a function of shear stress for particle size ratio of δ=3 at varying solid volume fraction and large particle fraction of (a) ζ=0.15, (b) ζ=0.50, and (c) ζ=0.85. The legend in panel b is used in all plots. (a) ζ=0.15. (b) ζ=0.50. (c) ζ=0.85. https://doi.org/10.1122/8.0000495.1

Close modal
FIG. 4.

Second normal stress difference (scaled by η0γ˙) as a function of shear rate for particle size ratio of δ=3 at varying solid volume fraction and large particle fraction of (a) ζ=0.15, (b) ζ=0.50, (c) ζ=0.85. The legend in panel b is used in all plots. (a) ζ=0.15. (b) ζ=0.50. (c) ζ=0.85. https://doi.org/10.1122/8.0000495.2

FIG. 4.

Second normal stress difference (scaled by η0γ˙) as a function of shear rate for particle size ratio of δ=3 at varying solid volume fraction and large particle fraction of (a) ζ=0.15, (b) ζ=0.50, (c) ζ=0.85. The legend in panel b is used in all plots. (a) ζ=0.15. (b) ζ=0.50. (c) ζ=0.85. https://doi.org/10.1122/8.0000495.2

Close modal
FIG. 5.

(a)-(c) Rheological functions as a function of shear stress for varying large particle fraction ζ at fixed solid volume fraction ϕ=0.57 and particle size ratio δ=3. (a) Relative viscosity, ηr (b) first normal stress difference, N1, and (c) second normal stress difference, N2; both N1 and N2 are scaled by η0γ˙. (d) Relative viscosity as a function of large particle fraction for δ=3 at ϕ=0.57 for different shear stresses. https://doi.org/10.1122/8.0000495.3

FIG. 5.

(a)-(c) Rheological functions as a function of shear stress for varying large particle fraction ζ at fixed solid volume fraction ϕ=0.57 and particle size ratio δ=3. (a) Relative viscosity, ηr (b) first normal stress difference, N1, and (c) second normal stress difference, N2; both N1 and N2 are scaled by η0γ˙. (d) Relative viscosity as a function of large particle fraction for δ=3 at ϕ=0.57 for different shear stresses. https://doi.org/10.1122/8.0000495.3

Close modal

Viscosity data of the ϕ=0.57 suspension at different ζ are plotted in Fig. 5(a). This shows that at ζ=0.85 the suspension has a distinctive behavior: its flow curve ηr(σ/σ0) does not repeat the other curves’ features, but instead crosses other curves and shows a very low viscosity at σ/σ0=0.1. The same behavior extends to the second normal stress difference [Fig. 5(c)].

The unusual low-shear-stress behavior of the ζ=0.85 suspension is also reflected in Fig. 5(d), where relative viscosity is plotted as a function of large particle fraction at different shear stresses. At σ/σ0=0.1, the viscosity continuously decreases with increasing ζ, reaching its minimum among the simulated conditions at ζ=0.85. As will be discussed in Sec. III C, this low viscosity value at ζ=0.85 is attributed to the ordering of large particles into layers with normal along the velocity gradient direction.

In the high shear stress regime, σ/σ03, the viscosity reaches its lowest value at ζ=0.50, with its intermediate value at ζ=0.85, and the highest value at ζ=0.15. These values follow the trends reported in the simulation study of Pednekar et al. [27], which considered the high shear limit, as well as a number of experimental studies described in Sec. I B. However, our values are shown to be one order of magnitude higher than those reported by Pednekar et al. [27], due to our use of μ=1 rather than μ=0.2 in that work. The minimum viscosity values at ζ=0.50 in comparison with the points ζ=0.15 and ζ=0.85 (or the monodisperse suspensions, not studied here) can be related to the ζ=0.50 suspension having a larger ϕJ in the high shear stress regime (see Sec. III D). This is in agreement with the work of Shapiro and Probstein [35], who obtained the largest ϕJ for suspension[s] with the fraction of large particles of about ζ=0.6 to 0.75.

It is established that in the high shear limit, the viscosity of bidisperse suspensions decreases significantly compared to the monodisperse case as the particle size ratio increases [27,35]. In Fig. 6, our ϕ=0.57 suspension data shown for two particle size ratios, δ=3 and 4, are consistent with the trend: the viscosity at δ=4 is slightly lower than that for δ=3. In addition, our rate-dependent data show that this is true across the entire range of shear stresses studied. Note that in our stress-controlled simulations, shear stress is scaled by σ0=FR/6πas2, which is expected to be appropriate if the force to push two particles together is related to the small particle radius. In suspensions where the fraction of large particles is high, the meaningful length scale is expected to be the large particle radius. The alternative scaling of stress as σ0=FR/6πal2 is shown in the inset in Fig. 6(c): the curves related to large ζ are shifted to the right by 32 and 42 units for δ=3 and δ=4, respectively, and the onset of shear thickening occurs at the same value of σ/σ0 in this alternative scaling. Further consideration of this rescaling is presented in the supplementary material [42].

FIG. 6.

Relative viscosity as a function of shear stress for particle size ratio of δ=1.4, δ=3, and δ=4 at solid volume fraction ϕ=0.57 (and ϕ=0.55 for δ=1.4) and large particle fraction of (a) ζ=0.15, (b) ζ=0.50, (c) ζ=0.85. The inset in (c) shows rescaled data, when switching to large particle radius length scale σ0=FR/(6πal2), where al=3as and al=4as for δ=3 and δ=4 suspensions, respectively. Data for δ=1.4 are adapted from Mari et al. [5]. The legend in panel b is used in all plots. (a) ζ=0.15. (b) ζ=0.50. (c) ζ=0.85.

FIG. 6.

Relative viscosity as a function of shear stress for particle size ratio of δ=1.4, δ=3, and δ=4 at solid volume fraction ϕ=0.57 (and ϕ=0.55 for δ=1.4) and large particle fraction of (a) ζ=0.15, (b) ζ=0.50, (c) ζ=0.85. The inset in (c) shows rescaled data, when switching to large particle radius length scale σ0=FR/(6πal2), where al=3as and al=4as for δ=3 and δ=4 suspensions, respectively. Data for δ=1.4 are adapted from Mari et al. [5]. The legend in panel b is used in all plots. (a) ζ=0.15. (b) ζ=0.50. (c) ζ=0.85.

Close modal

For the purpose of comparison, we plot our results for δ=3 and 4 at ζ=0.50 against those of Mari et al. [5] for δ=1.4 at ζ=0.50 [Fig. 6(b)]. We observe that the most significant drop in viscosity takes place when moving from δ=1.4 to 3 as the suspension response to shear stress changes from sharp DST to CST. Figure 6(b) illustrates the severity of shear thickening at δ=1.4. For instance, at ϕ=0.55 and δ=1.4 we still observe higher viscosity than the one at ϕ=0.57 and δ=3. When moving from δ=3 to 4 the decrease in viscosity is weaker but the suspension exhibits a smoother CST for δ=4. The decrease of viscosity and, hence, severity of shear thickening with δ is a result of increasing ϕJ, and is in agreement with experimental work of Shapiro and Probstein [35].

It is known that suspensions can exhibit a shear-induced ordering. In their study on sheared non-Brownian suspensions, Sierou and Brady [43] found an unexpected decrease in viscosity with increasing ϕ for 0.50<ϕ<0.55 before the viscosity again increased at higher ϕ. It was suggested that the decrease in viscosity was related to the suspension undergoing flow-induced ordering. Kulkarni and Morris [44], studying the ordering transition by Stokesian Dynamics simulation, observed that monodisperse Brownian suspensions with ϕ>0.50 have a tendency for layered or string-like ordering at elevated shear rates, Pe>1. Similar ordering has also been observed in shear flow of very dense soft-particle suspensions [45].

In order to understand the nature of the unusually low viscosity point (ζ=0.85, σ/σ0=0.1) in Fig. 5, we examined the structure of the simulated suspension with ϕ=0.57 and δ=3. The snapshots in Fig. 7 show the initial (γ˙t=0) and the final (γ˙t=30) states. At the initial state, the small and large particles are well-mixed. Under shear, however, the large particles order. In the ordered state the large particles slip past each other, which results in the decrease of viscosity. The progress of ordering phenomenon can be easily seen in a video presented in the supplementary material [42]. We demonstrate the time dependence this introduces to the suspension properties in Fig. 8 by plotting the viscosity variation across the entire simulation, showing that the viscosity decreases with time (strain) as the large particles order. The characteristic strain scale is found to be γ˙t7, as shown by the fitting ηr=AeBγ˙t+C to the viscosity data, with approximate parameter values A=12,B=0.15, and C=13. The behavior is reminiscent of the absorbing state dynamics seen in oscillatory shear flow of a simplified suspension model [46], where the fraction of contacts relaxes as (eγ˙t/τ)/(γ˙t)λ; a fit of this form in which we require monotonic decrease of the viscosity consistent with the trend observed results in λ=0 and the fit is the same as that shown.

FIG. 7.

Microstructure of a ϕ=0.57 suspension with δ=3 and ζ=0.85 at the (a) initial γ˙t=0 and (b) final γ˙t=30 states at dimensionless shear stress of σ/σ0=0.1. See the supplementary material [42] for the animation, and Fig. 8 for the associated viscosity evolution. (a) initial state γ˙t=0. (b) final state γ˙t=30.

FIG. 7.

Microstructure of a ϕ=0.57 suspension with δ=3 and ζ=0.85 at the (a) initial γ˙t=0 and (b) final γ˙t=30 states at dimensionless shear stress of σ/σ0=0.1. See the supplementary material [42] for the animation, and Fig. 8 for the associated viscosity evolution. (a) initial state γ˙t=0. (b) final state γ˙t=30.

Close modal
FIG. 8.

Relative viscosity, ηr, as a function of strain, γ˙t (between the beginning and ending conditions of Fig. 7) for ϕ=0.57 with the large particle fraction of ζ=0.85 and the particle size ratio of δ=3 at shear stress of σ/σ0=0.1. Solid black line: exponential fit ηr=AeBγ˙t+C, with approximate parameters A=12, B=0.15, and C=13.

FIG. 8.

Relative viscosity, ηr, as a function of strain, γ˙t (between the beginning and ending conditions of Fig. 7) for ϕ=0.57 with the large particle fraction of ζ=0.85 and the particle size ratio of δ=3 at shear stress of σ/σ0=0.1. Solid black line: exponential fit ηr=AeBγ˙t+C, with approximate parameters A=12, B=0.15, and C=13.

Close modal

We have also found this type of ordering at conditions near ζ=0.85 and δ=3, under low stress. As shown in Fig. 9, the ordering and time-dependent viscosity appear in suspensions at different ϕ; the ordering also occurs at δ=4 (data not shown). The magnitude of the drop in viscosity with time becomes more pronounced as volume fraction increases. Suspensions with other large particle fractions studied, ζ=0.15 and ζ=0.50, show disordered structure and time-independent rheology with properties fluctuating around their mean values, as shown by Fig. 10.

FIG. 9.

Relative viscosity ηr as a function of strain γ˙t for large particle fraction of ζ=0.85 and particle size ratio of δ=3 for varying volume fraction in the range of ϕ=0.540.59 at dimensionless shear stress of σ/σ0=0.1. Viscosity results for ϕ=0.57 are averaged over 3 simulation runs.

FIG. 9.

Relative viscosity ηr as a function of strain γ˙t for large particle fraction of ζ=0.85 and particle size ratio of δ=3 for varying volume fraction in the range of ϕ=0.540.59 at dimensionless shear stress of σ/σ0=0.1. Viscosity results for ϕ=0.57 are averaged over 3 simulation runs.

Close modal
FIG. 10.

Relative viscosity ηr as a function of strain γ˙t for the ϕ=0.57δ=3 suspensions with varying ζ at low shear stress 0.1σ/σ00.32.

FIG. 10.

Relative viscosity ηr as a function of strain γ˙t for the ϕ=0.57δ=3 suspensions with varying ζ at low shear stress 0.1σ/σ00.32.

Close modal

1. Mean-square displacement:

We studied the influence of the suspension structure on transport properties of suspensions through the mean-square particle displacement, MSD = (1/N)i=1N[ri(t)ri(0)]2, where N is the number of particles, t is the simulation time, ri(t) and ri(0) are the positions of particle i at time t and the initial time, respectively. As we are interested in the MSD as it relates to the observed ordering, the MSD is calculated in the vorticity, Δy2, and velocity gradient, Δz2, directions, and we do not consider the dispersion in the mean flow direction (Δx2). The MSD is scaled by as2.

Figure 11 shows the time variation of Δy2 and Δz2 for small and large particles for a suspension with ϕ=0.57, δ=3, and for different ζ at low (0.1σ/σ00.32) and high (σ/σ0=100) shear stresses. It is evident that the long-time high stress MSD values are larger than those at low stress. The most noticeable difference is seen for ζ=0.85. At low shear stress, the Δz2 curve for large particles displays a near-zero slope after γ˙t7 with Δz2/as25. At γ˙t7, the viscosity is seen to approach a statistically steady state for ϕ=0.57 in Fig. 8, although there is a slight further relaxation. The flattening of both viscosity and mean-square z-displacement curves is related to the large particles ordering into layers in the flow direction (with normal in the z direction). At high stress, the small and large particle MSD curves have comparable behavior. The systematic increase of Δy2 and Δz2, aside from short-term decreases due to fluctuations, for both small and large particles indicates that ordering is no longer present.

FIG. 11.

Mean-square displacement in the vorticity Δy2 and velocity gradient Δz2 directions as a function of time (presented as a strain unit γ˙t) for small and large particles for a suspension with ϕ=0.57, δ=3 for varying large particle fraction ζ=0.15,0.50,0.85 at low (LSS, 0σ/σ00.32) and high (HSS, σ/σ0=100) shear stresses. The legend in panel a is used in all plots. (a) ζ=0.15. (b) ζ=0.15. (c) ζ=0.50. (d) ζ=0.50. (e) ζ=0.85. (f) ζ=0.85.

FIG. 11.

Mean-square displacement in the vorticity Δy2 and velocity gradient Δz2 directions as a function of time (presented as a strain unit γ˙t) for small and large particles for a suspension with ϕ=0.57, δ=3 for varying large particle fraction ζ=0.15,0.50,0.85 at low (LSS, 0σ/σ00.32) and high (HSS, σ/σ0=100) shear stresses. The legend in panel a is used in all plots. (a) ζ=0.15. (b) ζ=0.15. (c) ζ=0.50. (d) ζ=0.50. (e) ζ=0.85. (f) ζ=0.85.

Close modal

In the ζ=0.15 suspension at low shear stress, the average displacement of small particles systematically increased, while that of large particles remained small, Δy20.3 and Δz22. In this suspension with a few large particles, the large particles seem to behave as single bodies immersed in a bath of “fluid” of small particles subjected to shear flow. At high shear stress, the Δy2 and Δz2 curves display larger values, with the small particles curves growing approximately linearly with time. The average y and z displacements for small and large particles in the ζ=0.50 suspension increase significantly with time, indicating the absence of ordering, for both low and high stress.

The ordering described to this point was for large particles at low shear stress. However, when analyzing the MSD results, we also found ordering occurring at high shear stress. Figure 12 shows MSD for small and large particles for ζ=0.50 suspensions with varying δ. The data for a suspension where δ=1.4, ϕ=0.56, ζ=0.50 were adopted from Singh et al. [12] The high stress data clearly have larger MSD values than the low stress ones, in both directions, except for the large particle Δy2 and Δz2 at δ=4, where the values at high shear stress are slightly lower. In the long-time limit, the Δz2 curves fluctuate around their mean values for both stress regimes, although this fluctuation starts earlier at high shear stress. The examination of the microstructure of the simulated suspension (see the supplementary material [42]) revealed the large particle ordering starting at γ˙t15 at low shear stress and γ˙t9 at high shear stress. Reconsidering the viscosity data in Fig. 6(b), we find that the δ=4 values at high stress are lower than what the high-stress data of Pednekar et al. [27] would suggest. For instance, their δ=4 suspension viscosity values are 25% lower than those of δ=3. In our case this gap is almost 50%, a difference which it appears can be explained by the large particle ordering.

FIG. 12.

Mean-square displacement in the vorticity Δy2 and velocity gradient Δz2 directions as a function of time (presented as a strain unit γ˙t) for small and large particles for a suspension with fixed ζ=0.50: for ϕ=0.56 and δ=1.4; for ϕ=0.57 and δ=3; for ϕ=0.57 and δ=4. The plots presented are for (a) and (c) low 0σ/σ00.32 and (b) and (d) high σ/σ050 shear stresses. The data for δ=1.4 were adopted from Singh et al. [12] The legend in panel b is used in all plots.

FIG. 12.

Mean-square displacement in the vorticity Δy2 and velocity gradient Δz2 directions as a function of time (presented as a strain unit γ˙t) for small and large particles for a suspension with fixed ζ=0.50: for ϕ=0.56 and δ=1.4; for ϕ=0.57 and δ=3; for ϕ=0.57 and δ=4. The plots presented are for (a) and (c) low 0σ/σ00.32 and (b) and (d) high σ/σ050 shear stresses. The data for δ=1.4 were adopted from Singh et al. [12] The legend in panel b is used in all plots.

Close modal

In the present section, we show that the effect of the particle size distribution can be successfully described by means of the maximum packing fraction, ϕJ. To determine ϕJ, we used the least-square error method to fit the power law

ηr=(1ϕ/ϕJ)2
(1)

to the simulated viscosity as a function of volume fraction η(ϕ). This was done for low shear stress 0.1<σ/σ0<0.3 to obtain a frictionless maximum packing fraction, ϕJ0, and high shear stress σ/σ0>70 to obtain a frictional maximum packing fraction, ϕJμ. The optimal values of ϕJ for each suspension with δ=1.4, δ=3, and δ=4 at both low and high stress are reported in Table I.

TABLE I.

Maximum packing fractions at low (ϕJ0) and high (ϕJμ) shear stresses for suspensions with various δ and ζ. The monodisperse ϕJ (δ = 1, ζ = 0 or 1) were obtained after treating data of Chong et al. [29] for low shear stress and of Pednekar et al. [27] for high shear stress.

δζϕJ0ϕJμ
0 or 1 0.62 0.61 
1.4 0.50 0.667 0.572 
0.15 0.675 0.599 
 0.50 0.723 0.626 
 0.85 0.785 0.605 
0.15 0.684 0.602 
 0.50 0.741 0.645 
 0.85 0.771 0.613 
δζϕJ0ϕJμ
0 or 1 0.62 0.61 
1.4 0.50 0.667 0.572 
0.15 0.675 0.599 
 0.50 0.723 0.626 
 0.85 0.785 0.605 
0.15 0.684 0.602 
 0.50 0.741 0.645 
 0.85 0.771 0.613 

Figure 13 illustrates that Eq. (1) represents relatively well the viscosity of simulated suspension in all cases but one: for ζ=0.15 at low shear stress. In this case, where the small particles make up most of the suspension, we observe a change in microstructural behavior from an ordered to a disordered state as the volume fraction increases. At low solid volume fractions (ϕ=0.54), the small particles tend to order in layers, with only sparse disorder regions around the few large particles, which results in low viscosity. As we increase the volume fraction, this ordering becomes less pronounced, while viscosity gradually increases. Finally, at ϕ=0.59, we observe that all ordering has disappeared, and the viscosity becomes noticeably higher. This transition from ordered to disordered state makes it challenging to describe the viscosity with this power law, which describes well the disordered states. It is worth noting that the low-stress ζ=0.85 suspension exhibits ordering for all volume fractions studied and, hence, is well-described by Eq. (1); the simulated conditions are predicted to be well below the maximum packing fraction of ϕJ0=0.785, although the order may break down before this solid fraction.

FIG. 13.

Relative viscosity plotted as a function of solid volume fraction, ηr(ϕ), at (a) low shear stress 0.1σ/σ00.32 and (b) high shear stress σ/σ070. The solid lines are the fit to power law Eq. (1) with maximum packing fractions as mentioned. The legend in panel a is used in all plots.

FIG. 13.

Relative viscosity plotted as a function of solid volume fraction, ηr(ϕ), at (a) low shear stress 0.1σ/σ00.32 and (b) high shear stress σ/σ070. The solid lines are the fit to power law Eq. (1) with maximum packing fractions as mentioned. The legend in panel a is used in all plots.

Close modal

In Fig. 14(a) we plot the relative viscosity ηr as a function of the reduced volume fraction ϕ/ϕJ. We observe that the viscosity data of all suspensions studied roughly collapse onto a single master curve but with poorer fit for δ=1.4 generally, and for ζ=0.15 at δ=3. The power law [Eq. (1)] satisfactorily predicts the effect of particle size ratio and large particle fraction on the value of maximum packing fraction at both low and high shear stresses.

FIG. 14.

Relative viscosity as a function of reduced volume fraction (a) ϕ/ϕJ and (b) (1ϕ/ϕJ) for different bidisperse suspensions at low shear stress (LSS) 0.1σ/σ00.32, and high shear stress (HSS) σ/σ070. Viscosity data of Chong et al. [29] were obtained at near zero shear rate. The solid lines are Eq. (1). The legend in panel a is used in both plots.

FIG. 14.

Relative viscosity as a function of reduced volume fraction (a) ϕ/ϕJ and (b) (1ϕ/ϕJ) for different bidisperse suspensions at low shear stress (LSS) 0.1σ/σ00.32, and high shear stress (HSS) σ/σ070. Viscosity data of Chong et al. [29] were obtained at near zero shear rate. The solid lines are Eq. (1). The legend in panel a is used in both plots.

Close modal

In Fig. 14(b) the same viscosity data are plotted in the form of ηr as a function of 1ϕ/ϕJ. This shows excellent agreement between our data and the power law up until the two leftmost points on the plot. These correspond to ϕ=0.59 suspensions with ζ=0.15 and ζ=0.85 at high shear stress, σ/σ0=100. The discrepancy for these conditions between the simulated data and Eq. (1) is most likely due to the contact rigidity threshold when the suspension concentration is close to ϕJ. Specifically, when the particles are constantly in contact, the interparticle interaction plays a major role and this might slightly underestimate the maximum packing fraction.

Experimental data (Chong et al. [29], Shapiro and Probstein [35]) and other simulations (Pednekar et al. [27]) also collapse onto the master curve (Fig. 14). Applying to these data the methodology previously described in this paper, we obtain maximum packing fractions that differ slightly from the value reported by Chong et al. and more significantly from the value reported by Shapiro and Probstein, with the values we determine being larger. When we plot their viscosity data against the reduced volume fraction using the maximum packing fraction we obtained, the fit appears consistent. Both the maximum packing fractions originally reported and the ones obtained here are included in Table II. The agreement between our data, the experimental data and Eq. (1) suggests that the power law captures the essential features of the systems studied. However, a recent experimental work of Ramaswamy et al. [9] found that the power law exponent changes from 2 to 3/2 when transitioning from frictionless to frictional shear jamming regimes.

TABLE II.

Data obtained from Shapiro and Probstein [35], Chong et al. [29], Pednekar et al. [27] along with ϕJ obtained with the method associated with Eq. (1).

Pednekar et al. [27
ϕ δ ζ η ϕJ, our method ϕJ in original work 
0.60   70   
0.62   120   
0.64 0.65 220 0.69 0.69 
0.65   320   
0.66   560   
Chong et al. [29
ϕ δ ζ η ϕJ, our method ϕJ in original work 
0.550   13   
0.575   19   
0.585   20   
0.605 3.2 0.75 60 0.68 0.66a 
0.615   70   
0.630   170   
0.650   700   
Shapiro and Probstein [35
ϕ δ ζ η ϕJ, our method ϕJ in original work 
0.30   2.8   
0.35   4.0   
0.40 0.50 6.0 0.64 0.60 
0.45   11.0   
0.50   21.0   
0.55   0.50   
Pednekar et al. [27
ϕ δ ζ η ϕJ, our method ϕJ in original work 
0.60   70   
0.62   120   
0.64 0.65 220 0.69 0.69 
0.65   320   
0.66   560   
Chong et al. [29
ϕ δ ζ η ϕJ, our method ϕJ in original work 
0.550   13   
0.575   19   
0.585   20   
0.605 3.2 0.75 60 0.68 0.66a 
0.615   70   
0.630   170   
0.650   700   
Shapiro and Probstein [35
ϕ δ ζ η ϕJ, our method ϕJ in original work 
0.30   2.8   
0.35   4.0   
0.40 0.50 6.0 0.64 0.60 
0.45   11.0   
0.50   21.0   
0.55   0.50   
a

Based on the method of Chong et al. associated with their Fig. 10.

To extend the data in Table I to ζ=0 and 1, where the packing is monodisperse, we also included ϕJ for monodisperse data. Those were obtained after treating data found in works of Chong et al. [29] (experimental, for low shear stress) and Pednekar et al. [27] (simulations in high shear limit). As expected from previous works [8,12,17], the maximum packing fraction at low shear stress ϕJ0 is larger than at high shear stress ϕJμ. The values of ϕJ0 peak at ζ=0.85 for δ=3 and 4. This is where we observed the lowest viscosity and the ordering. By contrast, ϕJμ is largest among our simulated conditions at ζ=0.50, with this peak higher as the size ratio δ increases. An increase of ϕJ with δ is observed, with the exception of the ζ=0.85 suspension at low shear stress where the δ=4 value of ϕJ0=0.771 is lower than for δ=3, ϕJ0=0.785; both are ordered states. Note, the elevated values of ϕJ at ζ=0 and 1 are explained by the lower values of viscosities reported in Pednekar et al. [27], due to their use of μ=0.2 (compared to our use of μ=1). Those results agree with experimental [35] and simulation work [27] with regard to the influence of ζ and δ at high shear stress.

In summary, we numerically studied the rheological properties of bidisperse rate-dependent dense suspensions with a focus on the role of particle size distribution on shear thickening. We have shown that the bidisperse suspensions we simulate, with electrical double layer (EDL) repulsion and contact friction, undergo shear thinning followed by shear thickening with increasing shear rate, for all values of large particle fractions ζ studied. The thickening is due to frictional contacts proliferating for σ/σ01, where σ0=FR/6πas2 is the characteristic stress associated with the maximum of the EDL force. The δ=3 suspensions with ζ=0.15 and 0.85 exhibited both CST (ϕ<0.57) and DST (ϕ0.57), while those with ζ=0.50 showed no DST up to ϕ=0.6, the largest volume fraction studied here. The second normal stress difference N2 was found to be negative (common for non-Brownian suspensions), and qualitatively similar to the viscosity in showing a decrease in magnitude during shear thinning and an increase during shear thickening. A significant reduction in viscosity as particle size ratio increases agrees well with known results, and is shown to extend to the entire range of shear stresses studied.

We have also demonstrated that under low shear rate conditions the ζ=0.85 suspension viscosity decreases over a long period. As it reaches a steady state, the large-particle mean-square displacement in the velocity gradient direction shows a zero slope, or near zero shear-induced self-diffusivity. The reduction in viscosity with time and a zero slope of the mean-square displacement are related to the large particles forming layers with normal in the velocity gradient direction, i.e., z for ux=γ˙z. Somewhat surprisingly, we also found some large-particle ordering at high shear stress for δ=4 and ζ=0.50, which requires further investigation.

Moreover, we found that the power law ηr=(1ϕ/ϕJ)2 captures well the influence on viscosity of changes in particle size ratio and the fraction of large particles. We showed that viscosities of all suspensions with one exception collapsed onto a single master curve when plotted against the reduced volume fraction ϕ/ϕJ, when we include the effect of stress on ϕJ (the low stress maximum packing fraction ϕJ0 and high stress maximum packing fraction ϕJμ) as shown in Table II; the exception is the case ζ=0.15 and low stress, for which the suspension underwent a partial ordering of the small particles to full disorder with increase of ϕ. We find that the cases with large particle ordering are captured by this power law because they were ordered at all conditions studied, and of course the fit may break down at higher ϕ if the ordering is disturbed. It is a quite simplifying feature that just these two stress-dependent jamming fractions, with both dependent on packing parameters δ and ζ, allow collapse of the rate-dependent viscosity data (Fig. 14). Our ability to collapse the viscosities found in literature (when treating their data using our approach in determining ϕJ) onto the same master curve shows that this simple power law can be used as a reliable tool in predicting the viscosity of bidisperse rate-dependent suspensions in both low and high shear stress regimes.

The authors gratefully acknowledge the discussions with Sidhant Pednekar, Omer Sedes, and Art Diky. This work was supported by a collaborative grant, NSF CBET-1916877 (B.C.) and NSF CBET-1916879 (J.F.M.). Computations in this work were supported, in part, under the NSF Grant Nos. CNS-0958379, CNS-0855217, and ACI-1126113 to the City University of New York High Performance Computing Center at the College of Staten Island.

The authors have no conflicts to disclose.

1.
Jolin
,
M.
,
D.
Burns
,
B.
Bissonnette
,
F.
Gagnon
, and
L.
Bolduc
, “Understanding the pumpability of concrete,” in Proceedings, Shotcrete for Underground Support XI, ECI Symposium Series, edited by F. Amberg and K. F. Garshol (ECI, ECI Digital Archives, Zurich, 2009).
2.
Brown
,
E.
, and
H. M.
Jaeger
, “
Shear thickening in concentrated suspensions: Phenomenology, mechanisms and relations to jamming
,”
Rep. Prog. Phys.
77
(
4
),
046602
(
2014
).
3.
Wagner
,
N. J.
, and
J. F.
Brady
, “
Shear thickening in colloidal dispersions
,”
Phys. Today
62
(
10
),
27
32
(
2009
).
4.
Guazzelli
,
É.
, and
J. F.
Morris
,
A Physical Introduction to Suspension Dynamics
(
Cambridge University Press
,
Cambridge
,
2011
).
5.
Mari
,
R.
,
R.
Seto
,
J. F.
Morris
, and
M. M.
Denn
, “
Shear thickening, frictionless and frictional rheologies in non-Brownian suspensions
,”
J. Rheol.
58
(
6
),
1693
1724
(
2014
).
6.
Egres
,
R. G.
, and
N. J.
Wagner
, “
The rheology and microstructure of acicular precipitated calcium carbonate colloidal suspensions through the shear thickening transition
,”
J. Rheol.
49
(
3
),
719
746
(
2005
).
7.
Seto
,
R.
,
R.
Mari
,
J. F.
Morris
, and
M. M.
Denn
, “
Discontinuous shear thickening of frictional hard-sphere suspensions
,”
Phys. Rev. Lett.
111
(
21
),
218301
(
2013
).
8.
Wyart
,
M.
, and
M. E.
Cates
, “
Discontinuous shear thickening without inertia in dense non-Brownian suspensions
,”
Phys. Rev. Lett.
112
(
9
),
098302
(
2014
).
9.
Ramaswamy
,
M.
,
I.
Griniasty
,
D. B.
Liarte
,
A.
Shetty
,
E.
Katifori
,
E.
Del Gado
,
J. P.
Sethna
,
B.
Chakraborty
, and
I.
Cohen
, “Universal scaling of shear thickening transitions” arXiv:2107.13338 (2021).
10.
Lin
,
N. Y. C.
,
B. M.
Guy
,
M.
Hermes
,
C.
Ness
,
J.
Sun
,
W. C. K.
Poon
, and
I.
Cohen
, “
Hydrodynamic and contact contributions to continuous shear thickening in colloidal suspensions
,”
Phys. Rev. Lett.
115
(
22
),
228304
(
2015
).
11.
Guy
,
B. M.
,
M.
Hermes
, and
W. C. K.
Poon
, “
Towards a unified description of the rheology of hard-particle suspensions
,”
Phys. Rev. Lett.
115
(
8
),
088304
(
2015
).
12.
Singh
,
A.
,
R.
Mari
,
M. M.
Denn
, and
J. F.
Morris
, “
A constitutive model for simple shear of dense frictional suspensions
,”
J. Rheol.
62
(
2
),
457
468
(
2018
).
13.
Maron
,
S. H.
, and
P. E.
Pierce
, “
Application of Ree–Eyring generalized flow theory to suspensions of spherical particles
,”
J. Colloid Sci.
11
(
1
),
80
95
(
1956
).
14.
Ad Hoc Committee on Official Nomenclature and Symbols
, “
Official symbols and nomenclature of the Society of Rheology
,”
J. Rheol.
57
(
4
),
1047
1055
(
2013
).
15.
Bashkirtseva
,
I. A.
,
A. Y.
Zubarev
,
L. Y.
Iskakova
, and
L. B.
Ryashko
, “
On rheophysics of high-concentrated suspensions
,”
Colloid J.
71
(
4
),
446
454
(
2009
).
16.
Morris
,
J. F.
, and
F.
Boulay
, “
Curvilinear flows of noncolloidal suspensions: The role of normal stresses
,”
J. Rheol.
43
(
5
),
1213
1237
(
1999
).
17.
Guy
,
B. M.
,
C.
Ness
,
M.
Hermes
,
L. J.
Sawiak
,
J.
Sun
, and
W. C. K.
Poon
, “
Testing the Wyart–Cates model for non-Brownian shear thickening using bidisperse suspensions
,”
Soft Matter
16
(
1
),
229
237
(
2020
).
18.
Mari
,
R.
,
R.
Seto
,
J. F.
Morris
, and
M. M.
Denn
, “
Nonmonotonic flow curves of shear thickening suspensions
,”
Phys. Rev. E
91
(
5
),
052302
(
2015
).
19.
Pan
,
Z.
,
H.
de Cagny
,
B.
Weber
, and
D.
Bonn
, “
S-shaped flow curves of shear thickening suspensions: Direct observation of frictional rheology
,”
Phys. Rev. E
92
(
3
),
032202
(
2015
).
20.
Neville
,
A. M.
,
Properties of Concrete
, 4th ed. (
Prentice-Hall
,
London
,
1995
).
21.
Klein
,
J.
,
S. P.
Mueller
, and
J. M.
Castro
, “
The influence of crystal size distributions on the rheology of magmas: New insights from analog experiments
,”
Geochem. Geophys. Geosyst.
18
(
11
),
4055
4073
, (
2017
).
22.
Nascimento
,
D. R.
,
B. R.
Oliveira
,
V. G. P.
Saide
,
S. C.
Magalhães
,
C. M.
Scheid
, and
L. A.
Calçada
, “
Effects of particle-size distribution and solid additives in the apparent viscosity of drilling fluids
,”
J. Pet. Sci. Eng.
182
(
106275
),
106275
(
2019
).
23.
Rasteiro
,
M. G.
, and
I.
Salgueiros
, “
Rheology of particulate suspensions in ceramic industry
,”
Part. Sci. Technol.
23
(
2
),
145
157
(
2005
).
24.
Huang
,
X.
,
Y.
Kakuda
, and
W.
Cui
, “
Hydrocolloids in emulsions: Particle size distribution and interfacial activity
,”
Food Hydrocoll.
15
(
4–6
),
533
542
(
2001
).
25.
Beckett
,
S. T.
,
The Science of Chocolate
, 3rd ed. (
Royal Society of Chemistry
,
Cambridge
,
2018
).
26.
Ogarko
,
V.
, and
S.
Luding
, “
Prediction of polydisperse hard-sphere mixture behavior using tridisperse systems
,”
Soft Matter
9
(
40
),
9530
9534
(
2013
).
27.
Pednekar
,
S.
,
J.
Chun
, and
J. F.
Morris
, “
Bidisperse and polydisperse suspension rheology at large solid fraction
,”
J. Rheol.
62
(
2
),
513
526
(
2018
).
28.
Bender
,
J.
, and
N. J.
Wagner
, “
Reversible shear thickening in monodisperse and bidisperse colloidal dispersions
,”
J. Rheol.
40
(
5
),
899
916
(
1996
).
29.
Chong
,
J. S.
,
E. B.
Christiansen
, and
A. D.
Baer
, “
Rheology of concentrated suspensions
,”
J. Appl. Polym. Sci.
15
(
8
),
2007
2021
(
1971
).
30.
Poslinski
,
A. J.
,
M. E.
Ryan
,
R. K.
Gupta
,
S. G.
Seshadri
, and
F. J.
Frechette
, “
Rheological behavior of filled polymeric systems II. the effect of a bimodal size distribution of particulates
,”
J. Rheol.
32
(
8
),
751
771
(
1988
).
31.
Chang
,
C.
, and
R. L.
Powell
, “
Dynamic simulation of bimodal suspensions of hydrodynamically interacting spherical particles
,”
J. Fluid Mech.
253
(
1
),
1
25
(
1993
).
32.
Hodgson
,
D. J. M.
,
M.
Hermes
,
E.
Blanco
, and
W. C. K.
Poon
, “Granulation and suspension rheology: A unified treatment,”
J. Rheol.
66
(
5
),
853
858
(
2022
).
33.
Farris
,
R. J.
, “
Prediction of the viscosity of multimodal suspensions from unimodal viscosity data
,”
Trans. Soc. Rheol.
12
(
2
),
281
301
(
1968
).
34.
Chang
,
C.
, and
R. L.
Powell
, “
Effect of particle size distributions on the rheology of concentrated bimodal suspensions
,”
J. Rheol.
38
(
1
),
85
98
(
1994
).
35.
Shapiro
,
A. P.
, and
R. F.
Probstein
, “
Random packings of spheres and fluidity limits of monodisperse and bidisperse suspensions
,”
Phys. Rev. Lett.
68
(
9
),
1422
1425
(
1992
).
36.
Gondret
,
P.
, and
L.
Petit
, “
Dynamic viscosity of macroscopic suspensions of bimodal sized solid spheres
,”
J. Rheol.
41
(
6
),
1261
1274
(
1997
).
37.
Mari
,
R.
,
R.
Seto
,
J. F.
Morris
, and
M. M.
Denn
, “
Discontinuous shear thickening in Brownian suspensions by dynamic simulation
,”
Proc. Natl. Acad. Sci. U.S.A.
112
(
50
),
15326
15330
(
2015
).
38.
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
(
1–4
),
444
472
(
1997
).
39.
Denn
,
M. M.
, and
J. F.
Morris
, “
Rheology of non-Brownian suspensions
,”
Annu. Rev. Chem. Biomol. Eng.
5
(
1
),
203
228
(
2014
).
40.
Gamonpilas
,
C.
,
J. F.
Morris
, and
M. M.
Denn
, “
Shear and normal stress measurements in non-Brownian monodisperse and bidisperse suspensions
,”
J. Rheol.
60
(
2
),
289
296
(
2016
).
41.
Royer
,
J. R.
,
D. L.
Blair
, and
S. D.
Hudson
, “
Rheological signature of frictional interactions in shear thickening suspensions
,”
Phys. Rev. Lett.
116
(
18
),
188301
(
2016
).
42.
See the supplementary material at https://www.scitation.org/doi/suppl/10.1122/8.0000495 for the ζ=0.50 suspension data after switching to the large particle size length scale σ0=FR/6πal2, the animation videos of the particle motions of our simulated suspensions, relaxation data for the normal stresses in the case of ordering, and further N1 data.
43.
Sierou
,
A.
, and
J. F.
Brady
, “
Rheology and microstructure in concentrated noncolloidal suspensions
,”
J. Rheol.
46
(
5
),
1031
1056
(
2002
).
44.
Kulkarni
,
S. D.
, and
J. F.
Morris
, “
Ordering transition and structural evolution under shear in Brownian suspensions
,”
J. Rheol.
53
(
2
),
417
439
(
2009
).
45.
Khabaz
,
F.
,
T.
Liu
,
M.
Cloitre
, and
R. T.
Bonnecaze
, “
Shear-induced ordering and crystallization of jammed suspensions of soft particles glasses
,”
Phys. Rev. Fluids
2
(
9
),
093301
(
2017
).
46.
Corte
,
L.
,
P. M.
Chaikin
,
J. P.
Gollub
, and
D. J.
Pine
, “
Random organization in periodically driven systems
,”
Nat. Phys.
4
(
5
),
420
424
(
2008
).