Suspensions of solid particles in liquids, under highly concentrated or ‘dense’ conditions in which the volume fraction ϕ approaches the jamming fraction, can exhibit strong shear thickening. This is seen as an abrupt increase in apparent viscosity η with increasing shear stress σ or shear rate γ ̇ [1–5]. Strong shear thickening has been linked to a shear-induced transition from lubricated or ‘frictionless’ particle interactions at low shear stress, where the interaction between the particles is hydrodynamic, to predominately contact interactions at high shear stress [6–9], with the critical stress determined by interparticle repulsive forces that tend to maintain the lubricated state at low stress [10–12].

Based on the lubrication-to-frictional transition scenario [12], Wyart and Cates (WC) [10] developed a model where shear thickening viscosity η is controlled by divergences at two jamming volume fractions: ϕ J 0 for the frictionless regime at low shear stress (in monodisperse systems, corresponding to random close packing ϕ J 0 = ϕ RCP ∼ 0.64 [13]) and ϕ J μ < ϕ J 0 for the frictional regime at high shear stress. The suspension viscosity was expressed as η ∼ ( ϕ J − ϕ ) − 2 , with the transition between ϕ J 0 and ϕ J μ implying a stress-dependent jamming fraction ϕ _J (σ). This was modeled using a parameter f, the fraction of close-pair interactions that are frictional, as ϕ J ( σ ) = ϕ J 0 ( 1 − f ( σ ) ) + ϕ J μ f ( σ ) . Here, 0 ≤ f ≤ 1, with f = 0 when neighboring particles have lubricated (frictionless) interactions at low stress, and f → 1 when the neighboring particles are predominately in the frictional state at high stress. In their work, f was described as a function of the particle pressure Π; since Π and σ for non-colloidal suspensions are related by an O (1) factor, μ _B = σ/Π [14], f can be written in terms of shear stress f(σ). Mari et al. [11] have validated elements of the WC theory by carrying out simulations of nearly monodisperse suspensions with particle size ratio δ = 1.4, where they have described f as f(σ) = exp (σ*/σ). Singh et al. [15] adopted this expression and extended the model to develop a constitutive description for simple shear flow. This provides a description of relative viscosity η/η ₀ (η ₀ is the suspending fluid viscosity), the second normal stress difference N ₂, and the particle pressure Π as a function of parameters ϕ, σ, and interparticle friction coefficient μ. As seen in other work [16], the first normal stress difference N ₁ = Σ₁₁ − Σ₂₂ was found to change sign from negative to positive with increasing ϕ and σ and was not well-described. We confine our consideration of the normal stress response to N ₂ = Σ₂₂ − Σ₃₃ and Π = −(Σ₁₁ + Σ₂₂ + Σ₃₃)/3, where 1, 2, and 3 (or x, y, and z) are, respectively, the flow, gradient, and vorticity directions of a viscometric flow. We study only simple shear here, of form u x = γ ̇ y .

Rheological experiments were carried out by Guy et al. [17] on monodisperse suspensions of poly-methylmethacrylate particles sterically stabilized by poly-12-hydroxystearic acid, where the data are shown to be fitted by the WC model. This model’s predictions were, however, found [18] to exhibit significant discrepancies from the experimental data for binary mixtures of spheres with large-to-small particle size ratio δ = 4. This was especially pronounced for cases with a predominance of particle volume from the large particles. Unlike monodisperse systems, where there is only one type of frictional contact, the bidisperse systems have three different types of frictional contact coming from the large-large, the large-small, and the small-small particles contact that contribute to stress development differently. This motivates our consideration of suspensions with significant bidispersity.

A recent study [19] showed that the WC model can be framed in the language of crossover scaling. Specifically, the shear thickening transition was described in terms of a crossover between two critical points, at each of which the viscosity diverges, the frictionless ϕ J 0 and the frictional ϕ J μ maximum packing fractions. In this scaling analysis, the viscosity was expressed as a function of the fraction of frictional contacts f(σ) and the distance to the frictionless maximum packing fraction, η r ϕ J 0 − ϕ 2 ∼ F W C f σ ϕ J 0 − ϕ , F W C ∼ 1 ϕ J 0 − ϕ J μ − f σ ϕ J 0 − ϕ − 2 , (1) where the scaling variable x W C = f ( σ ) / ( ϕ J 0 − ϕ ) is the part of the scaling function associated with the frictionless jamming critical point, and x c = 1 / ( ϕ J 0 − ϕ J μ ) is its critical value at which there is a divergence of the scaling function F W C . More generally, the scaling may be written η r ϕ J 0 − ϕ 2 ∼ F g σ , ϕ ϕ J 0 − ϕ , (2) with g allowing dependence on variables beyond the stress dependence of f(σ), and F W C a particular form of F ( x ) . With these generalizations, an excellent collapse to a single universal curve was achieved for experimental viscosity data for both cornstarch and silica particle suspensions measured over a wide range of volume fractions and shear stresses. In our work, we will not make allowance for ϕ-dependence in g, but will consider the scaling based on the WC scaling variable, x _WC , with f(σ) describing the microscopic interaction state, and identify any limitations of the generalized scaling function, F ( x W C ) , with regard to describing the bidisperse suspension rheology.

The scaling of the viscosity presented in Eq. 2 resembles a well-known form called crossover scaling. An alternative way of representing this scaling, which is better suited to identifying any change in the nature of the divergence with g (σ, ϕ) is the Cardy scaling form [20]. Cardy scaling is typically applied in study of critical points and phase transitions in thermal equilibrium systems, among which a well-known example is that of magnetic systems described by Heisenberg and Ising models. These differ in their symmetry properties [21], as the Heisenberg model is isotropic, whereas the Ising model has uniaxial symmetry as it considers a single component of the vector spin. Based on these different symmetry properties, the two models for magnetic systems show different universality classes characterized by different critical exponents and scaling functions [22], and a crossover in the dominant singularity may be observed. Analogously, the dense suspension that undergoes shear thickening by the lubricated-to-frictional mechanism also consists of different states, namely frictionless and frictional. Note that the function g (σ, ϕ) plays the same role as the uniaxial anisotropy plays in the Heisenberg-Ising crossover, but the suspension crossover does not involve a clear change of symmetry. Ramaswamy et al. [19] expressed the scaling function for the viscosity in the Cardy scaling form η r g 2 σ , ϕ ∼ H | 1 / x c − 1 / x | , (3) where H is a scaling function. The authors found that the viscosity divergence weakened to a -3/2 exponent, i.e., H ∼ | 1 / x c − 1 / x | − 3 / 2 , in the frictional regime, distinctly different from H ∼ | 1 / x c − 1 / x | − 2 in the frictionless limit. The two divergences of F suggest critical points of different universality classes.

In this work, we apply the scaling approach presented in Ramaswamy et al. [19] to dense suspension rheological data obtained from numerical simulation. We consider nearly-monodisperse systems, as well as more strongly bidisperse systems of varying large particle fraction. The analyses of the simulation data are carried out by two different approaches. We first consider the critical exponents and scaling functions. Following this, we implement Cardy scaling to study the crossover behavior.

Shear thickening has been observed in molecular dynamics simulations [23], including approaches which also introduce friction to the particle interactions [24]. Here, we prefer to consider the viscous fluid effects, but in a simplified fashion relative to the more rigorous Stokesian Dynamics [25, 26]. To this end, we use the lubrication flow-discrete element model (LF-DEM) [6] to simulate shear flow of dense suspensions of non-Brownian frictional spheres immersed in a viscous Newtonian fluid under an imposed shear stress σ. With a repulsive interaction and contact forces including friction, this method has been found to successfully describe the flow of dense suspensions [11], and to closely reproduce essential features of DST and shear jamming seen experimentally [5]. Recent work [27] has extended this approach to include a rolling friction [28], but here we consider only frictional resistance to slipping of the contact.

Our typical simulation is performed using 1,000 particles in a cubic box with periodic boundary conditions, to mimic an infinite system. Shear flow is imposed using Lees-Edwards boundary conditions. The suspensions studied in this work are bidisperse with a _s and a _l as the radii of the smaller and larger particles, respectively, with the size ratio defined as δ ≡ a _l /a _s . The motion of particles is considered inertialess, i.e. to be at particle Reynolds number R e = ρ γ ̇ a 2 / η 0 → 0 . The equation of motion is hence simply the force/torque balance 0 = F _H + F _R + F _C, where F _H is the finite-range hydrodynamic force/torque, F _R is the repulsive force, and F _C is the contact force/torque.

The hydrodynamic forces are of the form F _H = − R _FU ⋅ ( U − U ^∞ ) + R _FE: E ^∞ , where U is the particle translational and angular velocity, U ∞ = γ ̇ y e ̂ x is the flow due to imposed shear, and E ∞ = ( e ̂ x e ̂ z + e ̂ z e ̂ x ) γ ̇ / 2 is the rate-of-strain tensor. The resistance matrices R _FU and R _FE contain single-particle Stokes drag and the leading terms of the pairwise hydrodynamic interaction corresponding to short-range lubrication forces [29]. The divergence of the resistance matrix at vanishing dimensionless interparticle gap h = 2 (r − a _i − a _j )/(a _i + a _j ) is regularised to allow contact: the squeeze mode resistance is proportional to (h + ϵ)⁻¹ and shear mode resistance is proportional to log [(h + ϵ)⁻¹], where ϵ = 10⁻³ a _s . Including ϵ allows the lubrication resistance to be finite at contact (h = 0), which could be interpreted as representing the influence of particle surface roughness.

A repulsive electrical double layer (EDL) force maintains the particle surface separation at low stress. The force decays exponentially with the interparticle gap h over the length scale defined by Debye length κ ⁻¹ as F R = F R 0 e − κ h . For this work, based on prior experience [11], the Debye length is taken as κ ⁻¹ = 0.05a _s . The repulsive force introduces a stress scale of σ 0 = F R / 6 π a s 2 . This form of the simulation is called the electrostatic repulsion model (ERM) [11]. The contact between the particles occurs when the applied stress overcomes the repulsive forces. The contact force is modeled by linear springs in both normal and tangential directions; the friction prohibits slipping at the contact if the force components satisfy the Coulomb friction law, F _C, tan ≤ μF _C,nor , where μ is the interparticle friction coefficient. When the particles make contact, the friction is activated. The friction coefficient is fixed at μ = 1 for this work.

The repulsive force sets the level of the critical stress, such that when the applied stress σ ≪ σ ₀, the interactions between the particles are lubricated (frictionless). When σ ≫ σ ₀, the interactions are predominately frictional contacts. At each time step, we evaluate F _R and F _C and solve the equation of motion for the particle velocities, U − U ∞ = R FU − 1 ⋅ ( R FE : E ∞ + F R + F C ) . The flow, velocity gradient, and vorticity directions are x, z, and y, respectively. We do not consider Brownian motion, i.e. the Péclet number P e = η 0 γ ̇ a 3 / k T → ∞ .

The stress-controlled simulations were performed over a range of dimensionless shear stresses of 0.3 ≤ σ/σ ₀ ≤ 562. Each simulation was run to γ ̇ t = 30 strain units. The results of the transient period were not included, as we discard the initial five strain units. We report the suspension viscosity in the form of the relative viscosity η _r = η/η ₀, and for scaling purposes we will use the particle interaction viscosity η′ = η _r − (1 + 2.5ϕ), where (1 + 2.5ϕ) is the correction for the Einstein viscosity. The reported shear stresses and shear rates will be nondimensionalized by σ 0 = F R / 6 π a s 2 and γ 0 ̇ = F R / 6 π η 0 a s 2 , respectively. The bidisperse suspensions in this study are in the range of total solid volume fraction of 0.54 ≤ ϕ ≤ 0.60, with particle size ratio of δ = a _l /a _s = 1.4 and 3. For δ = 3, the large particle fraction of the total volume fraction is varied by taking ζ = ϕ _l /ϕ = 0.15, 0.50, and 0.85.

We present the results of the scaling analysis for all suspensions studied, first applying it to prior work on nearly monodisperse suspensions [15] with δ = 1.4 and ζ = 0.50 for volume fractions 0.52 ≤ ϕ ≤ 0.63 and shear stresses 0.1 ≤ σ/σ ₀ ≤ 100. For brevity, we call these suspensions monodisperse, as the slight bidispersity avoids the layering or string ordering of particles subjected to shear flow [30] but the shear thickening at δ = 1.4 is found to differ little from smaller δ = 1.2 [11]. We then present the scaling analysis applied to our rheological property data for bidisperse suspensions with δ = 3 and varying ζ = 0.15, 0.50, and 0.85 for 0.54 ≤ ϕ ≤ 0.60 and 0.3 ≤ σ/σ ₀ ≤ 562. Unless otherwise stated, we model the fraction of frictional contacts as f ( σ ) = e − σ * / σ , with σ*/σ ₀ = 1.45 based on previous work [11, 15, 31]. The rheology data that is considered in the scaling analysis, for both the monodisperse and bidisperse suspensions, is shown in Figure 1. Note that we do not consider data for cases where the viscosity is shear thinning, and thus for the monodisperse case we use data σ/σ ₀ ≥ 0.2 and for the bidisperse suspensions σ/σ ₀ ≥ 0.32, data which is to the right of the dashed lines shown in Figures 1A,G. This figure provides the relative viscosity η _r , normal stress difference as N ̂ 2 = N 2 / η 0 γ ̇ , and particle pressure as Π ̂ = Π / η 0 γ ̇ ; plotting the normal stresses normalized by η 0 γ ̇ produces a property analogous to the viscosity and proportional to the well-known ‘normal stress viscosity’ [14]. All quantities are shown as functions of the dimensionless shear stress σ/σ ₀. As noted, we do not consider the first normal stress difference in the scaling analysis.

As the two scaling methods we apply in this work have not previously been applied to simulated suspension rheology data, we provide some guidance before presenting our results. The two methods (which we will call crossover and Cardy scaling) have the same content, but the insights gained from the two are different. In the crossover scaling analysis (presented in Figure 2, and 4–6), the scaling approach can be seen based on the WC [10] form of the viscosity η r ∼ ( ϕ J ( σ ) − ϕ ) − 2 with ϕ J ( σ ) = f ( σ ) ϕ J μ + ( 1 − f ( σ ) ) ϕ J 0 using the fraction of frictional contacts f. This yields η r ∼ f ϕ J μ + 1 − f ϕ J 0 − ϕ − 2 = ϕ J 0 − ϕ − 2 1 − x / x c − 2 where, as defined previously, x = f / ϕ J 0 − ϕ a n d x c = 1 / ϕ J 0 − ϕ J μ .

Thus we see that if x → x _c , there are two possible divergences. From the crossover scaling analysis we find the data presented as η ′ ( ϕ J 0 − ϕ ) 2 as a function of x collapse to the generalized scaling form, F of Eq. 2 (not necessarily the precise form of WC), and confirm the presence of a second divergence, i.e. a divergence in the scaled data. Note that η′ is the singular part of the viscosity, which will be defined in the following paragraph. The form of the second divergence, at the frictional jamming fraction, revealed by crossover scaling and assumed to have exponent -2 in the WC form, is more effectively probed by Cardy scaling. For this purpose, we manipulate the expressions of the viscosity and other rheological functions differently, again using the WC form for concreteness, η r ∼ f ϕ J μ + 1 − f ϕ J 0 − ϕ − 2 = f − 2 ϕ J 0 − ϕ / f − ϕ J 0 − ϕ J μ − 2 = f − 2 x − 1 − x c − 1 − 2 .

Plotting η′f ² as a function of | x c − 1 − x − 1 | (as indicated by Eq. 3 with the more general g in place of f), a change in slope of the curve provides insight to whether the two divergences have the same exponent. In fact, we find evidence they are different.

We have carried out a scaling analysis on the rheology determined from simulations of shear-thickening suspensions. This is the first application of the crossover and Cardy scaling to simulated suspension rheology data. Our analysis includes nearly monodisperse suspensions studied in prior work [15] and bidisperse suspensions with particle size ratio δ = 3 for fractions of the particle volume occupied by large particles given by ζ = 0.15, 0.50, and 0.85. We considered first a critical exponent approach and then Cardy scaling. We showed that the singular part of the suspension relative viscosity given by η′ = η _r − (1 + 5ϕ/2) exhibits critical scaling, with the data collapsing to a universal curve when plotted as η ′ ( ϕ J 0 − ϕ ) 2 as a function of f / ( ϕ J 0 − ϕ ) for both monodisperse and bidisperse suspensions studied. Here, ϕ J 0 is the frictionless (or low stress) jamming fraction, and the fraction of frictional contacts f(σ) plays a key role. Note that while the data collapses with that of the monodisperse case for the smaller ζ bidisperse suspensions and the flow state diagrams [5, 10, 15] would be qualitatively similar, the forms would differ quantitatively because jamming fractions and the separation between them differ with the packing parameters δ and ζ.

The normal stress response was also found to follow the scaling, albeit with more scatter. Adapting Eq. 1 for N ̂ 2 = N 2 / η 0 γ ̇ , − N ̂ 2 ( ϕ J 0 − ϕ ) 2 showed reasonably good collapse. The particle pressure, however, exhibited a different exponent in order for the data to achieve a collapse, namely Π ̂ ( ϕ J 0 − ϕ ) 2.75 collapsed the data, again as a function of f / ( ϕ J 0 − ϕ ) ; we cannot at present offer a mechanistic basis for this difference in form for the particle pressure. All functions plotted in this fashion had a divergence at x c = 1 / ( ϕ J μ − ϕ J 0 ) , with ϕ J μ the frictional or high stress jamming fraction. As noted, ϕ J μ − ϕ J 0 differs depending upon the packing parameters.

The Cardy scaling analysis revealed a change in the exponent for the divergences of η′, N ̂ 2 , and Π ̂ . We find the frictionless limit to have a −2 exponent for η′ and N ̂ 2 and −2.75 for Π ̂ , and the magnitude decreases to an exponent of ≈ 1.5 as frictional jamming is approached. The change in slope occurred in the range |1/x _c − 1/x| of 10⁻¹ to 10⁻², similar to the location of the slope change observed in the experimental work of Ramaswamy et al. for two different systems. In essence, we demonstrated that all viscosity data, as well as N ₂ and Π data could be related to each other through the scaling parameter x = f / ( ϕ J 0 − ϕ ) .

We have also demonstrated that the fraction of frictional contacts f necessary to obtain collapse differs in form depending on the suspension packing parameters. Specifically, for monodisperse suspensions and the bidisperse ζ = 0.15-suspensions, f follows the exponential form f ( σ ) = e ( σ * / σ ) which has been previously described [15]. In order to achieve a collapse for the bidisperse ζ = 0.50 and 0.85 suspensions, however, we found that f had to be modified and this was done by fitting the function manually for each set of packing parameters (δ and ζ), but otherwise retaining only stress dependence, such that a scaling collapse could be obtained. This was due in part to the fact that the different size particles are driven into contact across a range of stresses, such that the single critical stress is incomplete. However, consistent with the arguments of Guy et al. [18], the need for different forms of f depending on ζ appears to also arise from the different contributions to the stress of frictional contacts between different types of particles, i.e. small-small, small-large, or large-large particle interaction, and this topic warrants further investigation.