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 $\delta =1.4$ (nearly monodisperse) to $\delta =4$, and the large particle fraction of the total solid loading with values $\zeta =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 $\varphi J$ with shear stress $\sigma $ are reported. At a fixed volume fraction $\varphi $, bidispersity decreases the suspension relative viscosity $\eta r=\eta s/\eta 0$, where $\eta s$ is the suspension viscosity and $\eta 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 $\zeta \u22480.5$, but instead decreases with further increase of $\zeta $ 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 $\eta r$ of bidisperse rate-dependent suspensions can be predicted by a power law linking it to $\varphi J$, $\eta r=(1\u2212\varphi /\varphi J)\u22122$ 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 $\varphi /\varphi J$.

## I. INTRODUCTION

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 $\varphi $ but also on shear rate $\gamma \u02d9$ or shear stress $\sigma $ [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.

### A. Rate-dependency

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 $\sigma 0=FR/6\pi a2$, where $FR$ is a repulsive interparticle force and $a$ is the particle radius. At low shear stress ($\sigma \u226a\sigma 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 $\eta r=\eta s/\eta 0$, where $\eta s$ is the suspension viscosity and $\eta 0$ is the suspending fluid viscosity) diverges as $\varphi \u2192\varphi J$, with $\varphi J=\varphi 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 $\sigma 0$), the repulsive forces acting between the particles are overwhelmed and particle contacts proliferate. In the high shear stress regime ($\sigma \u226b\sigma 0$), the interactions between the particles are mostly frictional and the viscosity diverges at the *frictional* maximum packing fraction [5,7,8,11] $\varphi J\mu <\varphi J0$. We note that, in various literature studies, the maximum packing fraction, or jamming fraction, is denoted as $\varphi J$, $\varphi max$, or $\varphi m$. For consistency with the work of Mari *et al.* [5] and Singh *et al.* [12] we use $\varphi J$, but use the term maximum packing fraction. The relative viscosity divergence is taken in the classical Maron–Pierce [13] form $\eta r\u223c(1\u2212\varphi /\varphi J)\u22122$ 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 $\u2212$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] $\Pi $ ($=\u2212[\Sigma 11+\Sigma 22+\Sigma 33]/3$), as was also assumed by Bashkirtseva *et al.* [15] Since $\Pi $ and $\sigma $ for non-colloidal suspensions are related [16], $f$ may also be taken as a function of shear stress. According to Wyart and Cates, $\varphi J$ is a function of $f$ and decreases from $\varphi J0$ in the frictionless state ($f=0$) to $\varphi J\mu $ in the frictional state ($f\u21921$). Based on more recent experimental work [9], it is argued that the fraction of frictional contacts may depend on both $\varphi $ 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: $\varphi $, $\sigma $, and interparticle friction coefficient $\mu $. 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 $\delta =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 $\delta =4$ are clearly different from the conclusions of Singh *et al.* [12], who performed simulations with $\delta =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.

### B. Bidispersity

In addition to the solid volume fraction $\varphi $ used to describe monodisperse systems, two other parameters are needed to characterize a binary mixture. We choose these to be a particle size ratio $\delta =al/as$, where $al$ and $as$ are the radii of the large and small particles, respectively, and the fraction of large particles $\zeta =\varphi l/\varphi $, where $\varphi 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 $\varphi 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 $\varphi /\varphi J$, with $\varphi 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 $\varphi /\varphi J$. The collapse of all data regardless of the packing parameters suggests that $\varphi /\varphi J$ is one of the most important parameters in suspension flow with $\varphi 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 $\delta $ and $\zeta $ suspended in polyisobutylene, showed that the Farris model works for $\delta $ of up to 10, but found further increase of $\delta $ to decrease the viscosity in an inconsistent way. Investigating the effect of bidispersity on the rheology of suspensions of various $\varphi $, these authors observed that small changes in $\zeta $ caused the viscosity to vary significantly with the minimum viscosity being achieved at $\zeta =0.65$–$0.75$. The most pronounced drop in viscosity compared to a monodisperse system was observed for the largest solid volume fraction ($\varphi =0.65$) with the largest particle size ratio studied ($\delta =22.22$). Their relative viscosity data plotted against the reduced volume fraction $\varphi /\varphi J$ collapsed onto a single curve, leading the authors to suggest that the relative viscosity of suspensions is only a function of $\varphi /\varphi 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 $\zeta =0.75$ for the largest $\delta $ studied ($\delta =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 $\varphi J$. The above findings show that the increase of $\varphi J$ in bidisperse suspensions relative to monodisperse suspensions and thus, at a given $\varphi $, increase of the “distance from jamming” $\varphi J\u2212\varphi $, is the basis for the lower viscosity of the bidisperse material. Since defining $\varphi 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 $\zeta $ and $\sigma $. 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 $\varphi =0.64$ bidisperse suspension with $\zeta =0.90$ and $\delta =2.1$ is one-fifth that of a monodisperse suspension of the same concentration. By decreasing $\zeta $ further to $0.50$, the reduction in viscosity became even stronger: compared to $\zeta =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 $\delta $ and $\zeta $ of concentrated bidisperse suspensions in the high shear limit, and showed that the viscosity reaches its minimum at $\zeta \u22480.65$.

In our study on rate-dependent bidisperse suspensions, we investigate how particle size ratio $\delta $ and large particle fraction $\zeta $ 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 $\varphi $, $\delta $, $\zeta $, and $\sigma $ 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 $0$–$0.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.

## II. MODEL AND METHODOLOGY

### A. Simulation tool

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 $\rho $ and viscosity $\eta 0$ under an imposed shear stress $\sigma $. We consider motion in the limit of vanishing particle Reynolds number, $Re=\rho \gamma \u02d9a2/\eta 0\u21920$, 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

where $U$ is the particle (translational and angular) velocity and $U\u221e=\gamma \u02d9ye^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\u221e=(e^xe^z+e^ze^x)\gamma \u02d9/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+\u03f5)\u22121$, while the resistance to tangential motions scales as $log\u2061[1/(h+\u03f5)]$, where $h=2(r\u2212ai\u2212aj)/(ai+aj)$ is the dimensionless interparticle gap; $\u03f5=10\u22123$ 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\u2264\mu FC,nor$ with components of the contact force computed using linear spring models as $FC,tan=kt\xi $ and $FC,nor=knhn$, where $\xi $ 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 $\mu =1$; when the Coulomb threshold is exceeded, i.e., when $FC,tan>\mu 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\u2212\kappa h$, where $\kappa \u22121$ is the Debye length, taken in this work as $\kappa \u22121=0.05as$. The repulsive force introduces a stress scale $\sigma 0=FR/6\pi as2$. This form of the simulation is called the electrostatic repulsion model (ERM) [5]. When the applied stress $\sigma \u226a\sigma 0$, the interactions between the particles are lubricated (frictionless). When $\sigma \u226b\sigma 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, $U\u2212U\u221e=RFU\u22121\u22c5(RFE:E\u221e+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=\eta 0\gamma \u02d9a3/kT\u2192\u221e$.

### B. Suspension parameters

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\u2264\varphi \u22640.60$. The particle size ratio is $1.4\u2264\delta \u22644$ and the fraction of large particles occupying the solid volume is $\zeta =0.15$, $0.50$, or $0.85$. To help visualize the particle configurations, Fig. 1 shows a $\varphi =0.57$, $\delta =3$ suspension at $\zeta =0.15$, $0.50$, and $0.85$. The scale for shear stress is $\sigma 0=FR/6\pi as2$, and hence the scale for shear rate is defined as $\gamma \u02d90=FR/6\pi \eta 0as2$. Reported shear stresses and shear rates will be nondimensionalized by $\sigma 0$ and $\gamma \u02d90$, respectively. We report the suspension viscosity $\eta s$ in the form of relative viscosity $\eta r=\eta s/\eta 0$, while normal stresses are scaled by $\eta 0\gamma \u02d9$.

The stress-controlled simulations were performed over a range of dimensionless shear stress of $0.01\u2264\sigma /\sigma 0\u2264562$. Each simulation in this study was run over a period of $\gamma \u02d9t=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 $\zeta =0.85$ suspensions at low stress ($\sigma /\sigma 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.

## III. RESULTS

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

### A. Rate-dependent viscosity

In Fig. 2, we plot the simulated relative viscosity $\eta r$ as a function of dimensionless shear stress $\sigma /\sigma 0$ and dimensionless shear rate $\gamma \u02d9/\gamma \u02d90$ for suspensions with $\delta =3$. Large particle fractions $\zeta =0.15$, $0.50$, and $0.85$ at several solid volume fractions $\varphi \u22650.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\u2264\sigma /\sigma 0\u22641$. As expected, the viscosity increases with increasing volume fraction. The onset of DST appears to occur at $\varphi =0.57$ for the systems with $\zeta =0.15$ and $\zeta =0.85$ and at $\varphi >0.60$ for $\zeta =0.50$. This corresponds to the start of nonmonotonic behavior of the flow curves in Fig. 3 where $\gamma \u02d9/\gamma \u02d90$ is plotted against $\sigma /\sigma 0$. The critical shear stress, which defines the onset of shear thickening, is at $\sigma /\sigma 0=0.3$ for $\zeta =0.15$ and $\zeta =0.50$ and $\sigma /\sigma 0=0.1$ for $\zeta =0.85$ and roughly constant for given $\zeta $ across all volume fractions studied [Fig. 2(a), 2(c), and 2(e)]. The second normal stress difference $N2=\Sigma 22\u2212\Sigma 33=\Sigma zz\u2212\Sigma 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=\Sigma 11\u2212\Sigma 22=\Sigma xx\u2212\Sigma 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 $\varphi $ and $\sigma $. 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.

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

The unusual low-shear-stress behavior of the $\zeta =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 $\sigma /\sigma 0=0.1$, the viscosity continuously decreases with increasing $\zeta $, reaching its minimum among the simulated conditions at $\zeta =0.85$. As will be discussed in Sec. III C, this low viscosity value at $\zeta =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, $\sigma /\sigma 0\u22653$, the viscosity reaches its lowest value at $\zeta =0.50$, with its intermediate value at $\zeta =0.85$, and the highest value at $\zeta =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 $\mu =1$ rather than $\mu =0.2$ in that work. The minimum viscosity values at $\zeta =0.50$ in comparison with the points $\zeta =0.15$ and $\zeta =0.85$ (or the monodisperse suspensions, not studied here) can be related to the $\zeta =0.50$ suspension having a larger $\varphi 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 $\varphi J$ for suspension[s] with the fraction of large particles of about $\zeta =0.6$ to $0.75$.

### B. Particle size ratio

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 $\varphi =0.57$ suspension data shown for two particle size ratios, $\delta =3$ and $4$, are consistent with the trend: the viscosity at $\delta =4$ is slightly lower than that for $\delta =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 $\sigma 0=FR/6\pi 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 $\sigma 0=FR/6\pi al2$ is shown in the inset in Fig. 6(c): the curves related to large $\zeta $ are shifted to the right by $32$ and $42$ units for $\delta =3$ and $\delta =4$, respectively, and the onset of shear thickening occurs at the same value of $\sigma /\sigma 0$ in this alternative scaling. Further consideration of this rescaling is presented in the supplementary material [42].

For the purpose of comparison, we plot our results for $\delta =3$ and $4$ at $\zeta =0.50$ against those of Mari *et al.* [5] for $\delta =1.4$ at $\zeta =0.50$ [Fig. 6(b)]. We observe that the most significant drop in viscosity takes place when moving from $\delta =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 $\delta =1.4$. For instance, at $\varphi =0.55$ and $\delta =1.4$ we still observe higher viscosity than the one at $\varphi =0.57$ and $\delta =3$. When moving from $\delta =3$ to $4$ the decrease in viscosity is weaker but the suspension exhibits a smoother CST for $\delta =4$. The decrease of viscosity and, hence, severity of shear thickening with $\delta $ is a result of increasing $\varphi J$, and is in agreement with experimental work of Shapiro and Probstein [35].

### C. Ordering phenomena

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 $\varphi $ for $0.50<\varphi <0.55$ before the viscosity again increased at higher $\varphi $. 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 $\varphi >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 ($\zeta =0.85$, $\sigma /\sigma 0=0.1$) in Fig. 5, we examined the structure of the simulated suspension with $\varphi =0.57$ and $\delta =3$. The snapshots in Fig. 7 show the initial ($\gamma \u02d9t=0$) and the final ($\gamma \u02d9t=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 $\gamma \u02d9t\u22487$, as shown by the fitting $\eta r=Ae\u2212B\gamma \u02d9t+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\u2212\gamma \u02d9t/\tau )/(\gamma \u02d9t)\lambda $; a fit of this form in which we require monotonic decrease of the viscosity consistent with the trend observed results in $\lambda =0$ and the fit is the same as that shown.

We have also found this type of ordering at conditions near $\zeta =0.85$ and $\delta =3$, under low stress. As shown in Fig. 9, the ordering and time-dependent viscosity appear in suspensions at different $\varphi $; the ordering also occurs at $\delta =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, $\zeta =0.15$ and $\zeta =0.50$, show disordered structure and time-independent rheology with properties fluctuating around their mean values, as shown by Fig. 10.

#### 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)\u2211i=1N[ri(t)\u2212ri(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, $\u27e8\Delta y2\u27e9$, and velocity gradient, $\u27e8\Delta z2\u27e9$, directions, and we do not consider the dispersion in the mean flow direction ($\u27e8\Delta x2\u27e9$). The MSD is scaled by $as2$.

Figure 11 shows the time variation of $\u27e8\Delta y2\u27e9$ and $\u27e8\Delta z2\u27e9$ for small and large particles for a suspension with $\varphi =0.57$, $\delta =3$, and for different $\zeta $ at low ($0.1\u2264\sigma /\sigma 0\u22640.32$) and high ($\sigma /\sigma 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 $\zeta =0.85$. At low shear stress, the $\u27e8\Delta z2\u27e9$ curve for large particles displays a near-zero slope after $\gamma \u02d9t\u22487$ with $\u27e8\Delta z2\u27e9/as2\u22485$. At $\gamma \u02d9t\u22487$, the viscosity is seen to approach a statistically steady state for $\varphi =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 $\u27e8\Delta y2\u27e9$ and $\u27e8\Delta z2\u27e9$, aside from short-term decreases due to fluctuations, for both small and large particles indicates that ordering is no longer present.

In the $\zeta =0.15$ suspension at low shear stress, the average displacement of small particles systematically increased, while that of large particles remained small, $\u27e8\Delta y2\u27e9\u22480.3$ and $\u27e8\Delta z2\u27e9\u22482$. 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 $\u27e8\Delta y2\u27e9$ and $\u27e8\Delta z2\u27e9$ 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 $\zeta =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 $\zeta =0.50$ suspensions with varying $\delta $. The data for a suspension where $\delta =1.4$, $\varphi =0.56$, $\zeta =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 $\u27e8\Delta y2\u27e9$ and $\u27e8\Delta z2\u27e9$ at $\delta =4$, where the values at high shear stress are slightly lower. In the long-time limit, the $\u27e8\Delta z2\u27e9$ 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 $\gamma \u02d9t\u224815$ at low shear stress and $\gamma \u02d9t\u22489$ at high shear stress. Reconsidering the viscosity data in Fig. 6(b), we find that the $\delta =4$ values at high stress are lower than what the high-stress data of Pednekar *et al.* [27] would suggest. For instance, their $\delta =4$ suspension viscosity values are $25%$ lower than those of $\delta =3$. In our case this gap is almost $50%$, a difference which it appears can be explained by the large particle ordering.

### D. Maximum packing fraction

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, $\varphi J$. To determine $\varphi J$, we used the least-square error method to fit the power law

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

δ
. | ζ
. | $\varphi J0$ . | $\varphi J\mu $ . |
---|---|---|---|

1 | 0 or 1 | 0.62 | 0.61 |

1.4 | 0.50 | 0.667 | 0.572 |

3 | 0.15 | 0.675 | 0.599 |

0.50 | 0.723 | 0.626 | |

0.85 | 0.785 | 0.605 | |

4 | 0.15 | 0.684 | 0.602 |

0.50 | 0.741 | 0.645 | |

0.85 | 0.771 | 0.613 |

δ
. | ζ
. | $\varphi J0$ . | $\varphi J\mu $ . |
---|---|---|---|

1 | 0 or 1 | 0.62 | 0.61 |

1.4 | 0.50 | 0.667 | 0.572 |

3 | 0.15 | 0.675 | 0.599 |

0.50 | 0.723 | 0.626 | |

0.85 | 0.785 | 0.605 | |

4 | 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 $\zeta =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 ($\varphi =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 $\varphi =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 $\zeta =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 $\varphi J0=0.785$, although the order may break down before this solid fraction.

In Fig. 14(a) we plot the relative viscosity $\eta r$ as a function of the reduced volume fraction $\varphi /\varphi J$. We observe that the viscosity data of all suspensions studied roughly collapse onto a single master curve but with poorer fit for $\delta =1.4$ generally, and for $\zeta =0.15$ at $\delta =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.

In Fig. 14(b) the same viscosity data are plotted in the form of $\eta r$ as a function of $1\u2212\varphi /\varphi J$. This shows excellent agreement between our data and the power law up until the two leftmost points on the plot. These correspond to $\varphi =0.59$ suspensions with $\zeta =0.15$ and $\zeta =0.85$ at high shear stress, $\sigma /\sigma 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 $\varphi 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 $\u2212$2 to $\u2212$3/2 when transitioning from frictionless to frictional shear jamming regimes.

Pednekar et al. [27] | |||||

ϕ | δ | ζ | η | ϕ_{J}, our method | ϕ_{J} in original work |

0.60 | 70 | ||||

0.62 | 120 | ||||

0.64 | 3 | 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.66^{a} |

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 | 4 | 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 | 3 | 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.66^{a} |

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 | 4 | 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 $\zeta =0$ and $1$, where the packing is monodisperse, we also included $\varphi 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 $\varphi J0$ is larger than at high shear stress $\varphi J\mu $. The values of $\varphi J0$ peak at $\zeta =0.85$ for $\delta =3$ and $4$. This is where we observed the lowest viscosity and the ordering. By contrast, $\varphi J\mu $ is largest among our simulated conditions at $\zeta =0.50$, with this peak higher as the size ratio $\delta $ increases. An increase of $\varphi J$ with $\delta $ is observed, with the exception of the $\zeta =0.85$ suspension at low shear stress where the $\delta =4$ value of $\varphi J0=0.771$ is lower than for $\delta =3$, $\varphi J0=0.785$; both are ordered states. Note, the elevated values of $\varphi J$ at $\zeta =0$ and $1$ are explained by the lower values of viscosities reported in Pednekar *et al.* [27], due to their use of $\mu =0.2$ (compared to our use of $\mu =1$). Those results agree with experimental [35] and simulation work [27] with regard to the influence of $\zeta $ and $\delta $ at high shear stress.

## IV. CONCLUSIONS

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 $\zeta $ studied. The thickening is due to frictional contacts proliferating for $\sigma /\sigma 0\u226b1$, where $\sigma 0=FR/6\pi as2$ is the characteristic stress associated with the maximum of the EDL force. The $\delta =3$ suspensions with $\zeta =0.15$ and $0.85$ exhibited both CST ($\varphi <0.57$) and DST ($\varphi \u22650.57$), while those with $\zeta =0.50$ showed no DST up to $\varphi =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 $\zeta =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=\gamma \u02d9z$. Somewhat surprisingly, we also found some large-particle ordering at high shear stress for $\delta =4$ and $\zeta =0.50$, which requires further investigation.

Moreover, we found that the power law $\eta r=(1\u2212\varphi /\varphi J)\u22122$ 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 $\varphi /\varphi J$, when we include the effect of stress on $\varphi J$ (the low stress maximum packing fraction $\varphi J0$ and high stress maximum packing fraction $\varphi J\mu $) as shown in Table II; the exception is the case $\zeta =0.15$ and low stress, for which the suspension underwent a partial ordering of the small particles to full disorder with increase of $\varphi $. 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 $\varphi $ 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 $\delta $ and $\zeta $, 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 $\varphi 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.

## ACKNOWLEDGMENTS

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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## REFERENCES

*Proceedings, Shotcrete for Underground Support XI*, ECI Symposium Series, edited by F. Amberg and K. F. Garshol (ECI, ECI Digital Archives, Zurich, 2009).