1. Introduction
Suspensions containing mixed Brownian colloids and larger, non-Brownian grains are prevalent across industry, with important examples being cement (Aitcin Reference Aitcin2000) and slurries (Chen et al. Reference Chen, Duan, Zhao and Yang2009). Obtaining a complete description of their constitutive behaviour is a fundamental challenge that embodies the intersection between fluid mechanics and rheology. When prepared with moderate to high solids volume fraction
$\phi$
, which is often desirable for function, these materials exhibit complex rheological behaviour including shear thinning, thickening and, as
$\phi$
approaches its maximal value
$\phi _m$
, flow arrest (de Kruif et al. Reference de Kruif, van Iersel, Vrij and Russel1985; Foss & Brady Reference Foss and Brady2000; Brown & Jaeger Reference Brown and Jaeger2014; Ness, Seto & Mari Reference Ness, Seto and Mari2022). In fully non-Brownian systems with bidisperse grains, the suspension viscosity diverges as a power law at large
$\phi$
, with (Ikeda, Berthier & Sollich Reference Ikeda, Berthier and Sollich2012; Pednekar, Chun & Morris Reference Pednekar, Chun and Morris2018)
\begin{align} \eta(\phi) \sim (1-(\phi/\phi_m))^{-\beta}. \end{align}
The jamming point
$\phi _m$
depends on the particle friction coefficient
$\mu$
and is also a non-monotonic function of the nominal size ratio
$\varDelta$
and volumetric mixing ratio
$\alpha$
of the two species (see e.g. Yu & Standish Reference Yu and Standish1991; Shapiro & Probstein Reference Shapiro and Probstein1992; Anzivino et al. Reference Anzivino, Casiulis, Zhang, Moussa, Martiniani and Zaccone2023). The exponent is usually reported as
$\beta \approx 2$
and, though weakly dependent on
$\varDelta$
and
$\alpha$
, it is likely sensitive to
$\mu$
and to particle shape. While geometric packing models are sufficient to predict the changes in constitutive curves observed when small particles are added to a system of larger ones (Farris Reference Farris1968; Singh et al. Reference Singh, Ness, Sharma, de, Juan and Jaeger2024), they neglect the additional complexity present when the smaller species is colloidal. Here we focus on non-attractive colloids with short-ranged repulsion, isolating the important role of the Brownian nature of the colloids. Understanding the physics of such systems is likely crucial for characterising phenomenology such as superplasticisation reported in concrete formulation and elsewhere (Flatt Reference Flatt2004; Roussel et al. Reference Roussel, Lemaître, Flatt and Coussot2010; Van Damme Reference Van Damme2018; Ito, Matsunaga & Imai Reference Ito, Matsunaga and Imai2019; Thiévenaz et al. Reference Thiévenaz, Rajesh and Sauret2021). Moreover, experimental measurements in bidisperse colloidal and granular suspensions (Cwalina & Wagner Reference Cwalina and Wagner2016; Madraki et al. Reference Madraki, Hormozi, Ovarlez, Guazzelli and Pouliquen2017, Reference Madraki, Ovarlez and Hormozi2018) report trends similar to those we explore here, highlighting the practical importance of understanding the interplay of Brownian motion, particle size and rheology.
In mixed colloidal and granular systems, taking
$\alpha$
as the fraction of the solid volume occupied by grains, the time scale introduced by thermal motion
${6 \pi \eta _0 a^3}/({k_B T})$
leads to rate-dependence in
$\phi _m$
and influences the functional form of
$\eta (\phi )$
. Here
$\eta _0$
is the solvent viscosity,
$a$
is a characteristic particle size,
$k_B$
is the Boltzmann constant and
$T$
is the temperature. Crucially, changing
$\alpha$
even at fixed shear rate
$\dot {\gamma }$
affects not only the geometry as in non-Brownian bidisperse suspensions but also the relative importance of diffusion, through the
$\alpha$
-dependence of the characteristic particle size used when calculating the Péclet number
$\textit{Pe} = {6 \pi \eta \dot \gamma a^3}/({k_B T})$
. For (nearly) monodisperse spheres at very large
$\textit{Pe}$
, the particles are effectively non-Brownian grains, so that the relevant
$\phi _m$
is the frictional random loose packing limit (e.g.
$\phi _m^{\mu =1}\approx 0.57$
). At
$\textit{Pe}\approx 1$
, recent work shows that thermal motion releases frictional constraints so that the relevant
$\phi _m$
becomes frictionless random close packing
$\phi _m^{\mu =0}\approx 0.64$
(Li, Royer & Ness Reference Li, Royer and Ness2024). Meanwhile in the limit of vanishing
$\textit{Pe}$
, the viscosity diverges at the glass transition
$\phi _G\approx 0.58$
with differing functional form
$\eta (\phi )\sim \exp {(k/(\phi _G-\phi ))}$
(Hunter & Weeks Reference Hunter and Weeks2012; Ikeda et al. Reference Ikeda, Berthier and Sollich2012). The nature of the crossover between this divergence and the power law is not well established, and our data reported later suggest an intermediate form present over a large range of
$\textit{Pe}$
. In characterising flowing colloidal-granular suspensions, it is thus crucial to understand the interplay between the frictionless and frictional character as functions of the flow rate and composition (Guy et al. Reference Guy, Hermes and Poon2015, Reference Guy, Richards, Hodgson, Blanco and Poon2018).
Real-world suspensions have, in practice, continuous mono- or multimodal particle-size distributions often spanning orders of magnitude in
$a$
from nano- to millimetric (Tari, Ferreira & Fonseca Reference Tari, Ferreira and Fonseca1999; Vlasak & Chara Reference Vlasak and Chara2011; Bauer & Nötzel Reference Bauer and Nötzel2014; Mehdipour & Khayat Reference Mehdipour and Khayat2017; Hlobil, Kumpová & Hlobilová Reference Hlobil, Kumpová and Hlobilová2022), obscuring the basic physics that link composition to constitutive behaviour. To make this problem tractable, we study by particle-based simulation a model bidisperse colloidal(
$c$
)-granular(
$g$
) system of frictional spheres with size ratio
$\varDelta \equiv a_g/a_c=5$
. This choice ensures that the separation in diffusive time scales
$\sim \! \varDelta ^3$
is large and that the variation in jamming point with
$\alpha$
qualitatively matches that for larger size ratios, so that the key features described in the following should be robust for larger
$\varDelta$
. The particles experience the same set of governing forces, so the system now has two separate dimensionless shear rates (Péclet numbers), with
$\textit{Pe}_c={6 \pi \eta \dot \gamma a_c^3}/({k_B T})$
for the colloids and likewise
$\textit{Pe}_g={6 \pi \eta \dot \gamma a_g^3}/({k_B T})$
for the grains. Varying
$\alpha$
from
$0$
to
$1$
over a wide range of dimensionless shear rates, our model allows us to map the full rheology of this bidisperse system at
$\phi$
close to the varying
$\phi _m$
.
In this article, we do three things. (i) We present the first systematic simulation study of the rheology of a colloid-grain mixture, showing separately the impact of adding colloids (grains) to a granular (colloidal) suspension. The rheology shows characteristic shear thinning followed by thickening with increasing
$\textit{Pe}$
, with the onset and extent of thickening being sensitive to
$\alpha$
. (ii) We provide the first simulation demonstration that adding colloids to a granular suspension can enhance flowability even as
$\phi$
increases. This can manifest as a modest viscosity reduction, or even as unjamming of a granular packing by adding a small volume of colloids that disrupt grain contacts. (iii) We explain the phenomenology by mapping
$\phi _m$
as a function of both
$\alpha$
and a suitably defined dimensionless shear rate. Collectively, these results offer a rationale for predicting the constitutive behaviour of suspensions of mixed colloids and grains, and provide a starting point for addressing the rheology of more complex multi-species suspensions, a major challenge in geophysics, formulation science and other application areas.
2. Simulation method overview
We simulate the trajectories of spheres with motion governed by Langevin equations comprising three pairwise force and torque contributions: direct contacts, short-ranged hydrodynamics and Brownian noise. Contact forces are modelled as linear springs with a tangential component regulated by a friction coefficient
$\mu =1$
. Following standard convention for dense suspensions, hydrodynamics enters through single-body drag and pairwise lubrication (Ball & Melrose Reference Ball and Melrose1997; Seto et al. Reference Seto, Mari, Morris and Denn2013; Cheal & Ness Reference Cheal and Ness2018), with long-range hydrodynamics expected to be largely screened and thus omitted from the model. Brownian forces and torques are generated on each time step to satisfy fluctuation–dissipation theorem as described in Li et al. (Reference Li, Royer and Ness2024). Our model thus incorporates sufficient microscopic physics to capture the suspension rheology around the colloidal-to-granular crossover as a function of
$\textit{Pe}$
for systems of non-attractive particles with short-ranged repulsion. A summary of the governing equations and parameter values is given in Appendix A, and a comprehensive model rationale and description is given in Li et al. (Reference Li, Royer and Ness2024).
The model is implemented in the open source molecular dynamics package LAMMPS, with shear applied by deforming the periodic simulation box at fixed volume
$V$
, updating the particle displacements following the velocity-Verlet algorithm. Simulations are initialised in randomised non-overlapping configurations of up to
$10^4$
particles, having verified that this is sufficient to mitigate finite-size effects by testing the model up to
$10^5$
particles. To prevent crystallisation but otherwise likely having a minimal effect on our results, our nominally bidisperse system comprises particles with size ratio
$1 : 1.4 : 5 : 7$
, with colloid and grain volumes given by
$V_c=V_1 + V_{1.4}$
and
$V_g=V_5 + V_{7}$
. The total solids volume fraction is
$\phi =(V_g+V_c)/V$
, with the volumetric mixing ratio
$\alpha =V_g/(V_c + V_g)$
and
$V_l=V-V_c-V_g$
the free ‘liquid’ volume. Snapshots in figure 1(a) show the simulated system at a range of
$\alpha$
. All particles experience the same forces, including the Brownian noise, but we refer to the smaller particles as colloids since their diffusive time scale is
$125\times$
shorter than that of the larger grains. We subject the system to constant
$\dot \gamma$
for sufficiently long run time
$t$
so that the shear stress
$\varSigma _{\textit{xy}}$
reaches a steady state at large strains
$\dot {\gamma }t$
, before taking the suspension viscosity as
$\eta = \Sigma _{xy}/ (\eta _0 \dot {\gamma})$
. Time-averaged velocity profiles follow closely the imposed shear rate and there is no evidence of heterogeneity or banding within the domain. To obtain constitutive curves, we adjust both
$\dot {\gamma }$
and
$k_BT$
to cover a wide range of
$\textit{Pe}_c=10^{-2}{-}10^6$
while maintaining both near inertia-free and hard-sphere conditions (see Appendix A).
3. Results: mixing colloids and grains
We first present two sets of constitutive curves in which particles of one species are systematically added to a suspension initially comprising only the other, figures 1(b) and 1(c). Here we plot the suspension viscosity
$\eta$
against
$\textit{Pe}_c$
, with corresponding values of
$\textit{Pe}_g$
along the top axis, so that by reading vertically upward in figures 1(b) and 1(c) one isolates the effect of adding colloids or grains at a fixed processing condition. The data represent averages of up to 50 realisations, with the errors bars representing the realisation-to-realisation variation of the time-averaged viscosities being smaller than the plot markers and thus not shown. The fluctuations in instantaneous viscosities are comparable to those reported earlier for near-monodisperse suspensions (Li et al. Reference Li, Royer and Ness2024).
Rheological impact of adding grains and colloids. (a) Snapshots of the simulated system at a range of compositions
$\alpha$
, showing particle sizes
$a$
(dark blue),
$1.4a$
(light blue),
$5a$
(grey) and
$7a$
(white). (b,c) Constitutive flow curves showing the effect of adding grains (b) or colloids (c) to a colloidal (b) or granular (c) suspension initially at
$\phi =0.55$
. Particle addition is done holding both the volume of the initial species,
$V_c$
(b) or
$V_g$
(c), and the liquid volume
$V_l$
fixed, so that both
$\phi$
and
$\alpha$
vary as detailed later. Vertical lines at selected
$\textit{Pe}_c$
highlight the variation of
$\eta (\phi )$
at selected shear rates replotted in (d) when adding grains (blue) and colloids (orange). Symbols for volume ratios
$\alpha$
and volume fractions
$\phi$
:
$\alpha =0,\phi = 0.55$
;
$\alpha =0.1, \phi = 0.576$
;
$\alpha =0.2,\phi = 0.604$
;
$\alpha =0.3,\phi = 0.636$
;
$\alpha =0.4,\phi = 0.671$
;
$\alpha =0.5,\phi = 0.710$
;
$\alpha =0.6,\phi = 0.671$
;
$\alpha =0.7,\phi = 0.636$
;
$\alpha =0.8,\phi = 0.604$
;
$\alpha =0.9,\phi = 0.576$
;
$\alpha =1,\phi = 0.55$
.
Starting from a pure colloidal system at
$\phi = 0.55$
and
$\alpha =0$
(green circles, figure 1
b), the rheology
$\eta (\textit{Pe}_c)$
is consistent with earlier works (Mari et al. Reference Mari, Seto, Morris and Denn2015; Li et al. Reference Li, Royer and Ness2024), with a shear-thinning regime at small
$\textit{Pe}_c\lesssim 1$
, followed by shear thickening to a frictional, rate-independent regime at large
$\textit{Pe}_c\gtrsim 10^3$
. Systematically introducing grains by increasing
$V_g$
while keeping both
$V_c$
and
$V_l$
constant, we find an increase in
$\eta$
with
$\alpha$
across all
$\textit{Pe}_c$
. This is as might be expected, since the addition of grains in this manner also increases the total volume fraction
$\phi$
. Additionally, the shear-thickening onset decreases monotonically to lower
$\textit{Pe}_c$
with increasing
$\alpha$
, accompanied by a steeper increase of
$\eta$
above the thickening onset. This mirrors experimental results, demonstrating enhanced thickening when adding larger non-Brownian grains to suspensions of smaller colloids (Cwalina & Wagner Reference Cwalina and Wagner2016; Madraki et al. Reference Madraki, Hormozi, Ovarlez, Guazzelli and Pouliquen2017, Reference Madraki, Ovarlez and Hormozi2018).
Repeating this process, now adding colloids to an initially granular system at
$\phi =0.55$
and
$\alpha =1$
, produces a markedly different result. The rheology of the pure granular system (red circles, figure 1
c) is similar to the pure colloidal system, but controlled by the granular Péclet number
$\textit{Pe}_g$
. However, the effect of adding colloidal particles by increasing
$V_c$
and now keeping both
$V_g$
and
$V_l$
constant, is
$\textit{Pe}_c$
-dependent. At low shear rates (
$\textit{Pe}_c \lesssim 0.1$
), the viscosity
$\eta$
increases monotonically with
$\phi$
, similar to our results for adding grains to colloidal suspensions. However, at higher
$\textit{Pe}_c$
, adding colloids can actually decrease the suspension viscosity, despite the increase in
$\phi$
. This viscosity reduction is most pronounced around
$\textit{Pe}_c \sim 1$
(
$\textit{Pe}_g\sim 10^2)$
, where the viscosity decreases by up to a factor of
${\sim}7$
at
$\alpha =0.7$
and
$\phi =0.636$
before increasing and eventually jamming at higher
$\phi$
. At this shear rate, the pure granular suspension is in the high-viscosity frictional flow regime while the pure colloids are close to their viscosity minimum. This non-monotonic behaviour persists at higher flow rates, where both the granular and colloidal suspensions have reached the frictional-flow regime, though the magnitude of the viscosity reduction is notably reduced. While there are numerous experimental reports of bidisperse suspensions having lower viscosity than their monodisperse counterpart at the same
$\phi$
(Chong, Christiansen & Baer Reference Chong, Christiansen and Baer1971; Gondret & Petit Reference Gondret and Petit1997; Zaman & Dutcher Reference Zaman and Dutcher2006; Guy et al. Reference Guy, Ness, Hermes, Sawiak, Sun and Poon2020), our simulation results are the first to show a significant viscosity reduction at higher
$\phi$
with the addition of colloidal solids. We further find that the range of shear rates over which shear thickening occurs is non-monotonic in
$\alpha$
, being narrower at
$\alpha =1$
and
$0.5$
compared with at intermediate values.
Re-plotting the data along the vertical lines of figures 1(b) and 1(c), showing
$\eta (\phi )$
at fixed
$\textit{Pe}_c$
, figure 1(d), highlights the differing response to the addition of either colloids or grains. (Note that
$\alpha$
varies along each of these lines, as indicated by the symbols which match those in figures 1
b and 1
c.) The addition of grains to colloids (blue lines) leads to a monotonically increasing
$\eta (\phi )$
, with the most rapid increase at intermediate
$\textit{Pe}_c$
. In contrast, when adding colloids to grains (orange lines)
$\eta (\phi )$
is non-monotonic for
$\textit{Pe}_c\gtrsim 0.1$
, with the viscosity first decreasing as
$\phi$
increases before eventually increasing again as the suspension approaches the jamming point. This viscosity reduction persists over a wide range of volume fraction, and is most notable for
$\textit{Pe}_c=1$
but also present for
$\textit{Pe}_c=10^3$
, indicating that thermal motion enhances this effect but is not crucial.
To rationalise the above mentioned results in terms of changes in the critical jamming point
$\phi _m$
, which should now depend on both
$\alpha$
and the shear rate, we must first reconsider how to define an appropriate Péclet number for these colloid–grain mixtures. Defining the Péclet number in terms of a single particle size (either
$\textit{Pe}_c$
or
$\textit{Pe}_g$
) allows for a systematic comparison of the constitutive curves at a fixed absolute shear rate
$\dot \gamma$
, but it becomes increasingly poor as a measure of the competition between diffusion and convection in mixtures at intermediate
$\alpha$
. This difficulty can be seen in constitutive flow curves
$\eta (\textit{Pe}_c)$
measured at fixed
$\phi$
and varying
$\alpha$
, where the transition from shear thinning to shear thickening shifts by over two orders of magnitude in
$\textit{Pe}_c$
, figure 2(a). Here the trend in high shear-rate viscosity is qualitatively consistent with (Pednekar et al. Reference Pednekar, Chun and Morris2018), with its range spanning an order of magnitude and having a minimum near
$\alpha =0.7$
.
Mapping between Péclet number definitions. Constitutive flow curves at various
$\alpha$
with fixed
$\phi = 0.565$
, reported as functions of (a) the colloidal
$\textit{Pe}_c$
and (b) the averaged
$\textit{Pe}_{\textit{cg}}$
. Coloured sections in (a) and (b) represent intermediate (blue), frictionless (purple) and frictional (green) states. (c) Sketch of the putative dependence of
$\phi _m$
on
$\alpha$
at the viscosity minimum where
$\textit{Pe}_{\textit{cg}}=5$
(purple) and in the granular limit where
$\textit{Pe}_{\textit{cg}}=10^5$
(green). Blue line interpolates between these as colloids are added at fixed
$\textit{Pe}_c$
, so that the effective jamming point
$\phi _m$
moves between green and purple lines as
$\alpha$
decreases from 1 to 0. Red lines sketch the change in
$\phi$
when adding either grains or colloids to a suspension at initial volume fraction
$\phi _0$
, highlighting the asymmetry in the change in the distance to jamming. Dotted blue line in (c) shows the volume fraction used in (a) and (b).
We find empirically that plotting viscosity against an averaged Péclet number
$\textit{Pe}_{\textit{cg}}={6 \pi \eta \dot \gamma a_{\textit{cg}}^3}/(k_B T)$
, defined using an
$\alpha$
-weighted mean particle size
$a_{\textit{cg}}=(1-\alpha )a_c + \alpha a_g$
, roughly aligns these flow curves, figure 2(b). This suggests that colloids contribute to the Brownian character in proportion to their volume rather than their number. The
$\eta ({Pe}_{\textit{cg}})$
flow curves now all have minima over a narrow range of
$\textit{Pe}_{\textit{cg}}\sim 10$
, and reasonable overlap in the shear thinning regime at lower
$\textit{Pe}_{\textit{cg}}$
. At higher
$\textit{Pe}_{\textit{cg}}$
, the shear-thickening dynamics show a stronger
$\alpha$
-dependence. Suspensions closer to the monodisperse limit (either
$\alpha \to 0$
or
$\alpha \to 1$
) have both an increased frictional (high-shear) viscosity and a more abrupt shear-thickening transition (
$10\lesssim \textit{Pe}_{\textit{cg}}\lesssim 10^3$
), while more evenly mixed suspensions (
$\alpha \approx 0.5$
) exhibit a broader thickening transition (
$10\lesssim \textit{Pe}_{\textit{cg}}\lesssim 10^4$
) to a lower high-shear viscosity despite all having the same
$\phi$
.
The aligned shear-thickening flow curves in figure 2(b) indicate that we can regard the intermediate- (pink) and high-
$\textit{Pe}_{\textit{cg}}$
(green) states as frictionless and frictional, respectively. Drawing from prior studies of binary sphere packings (Srivastava et al. Reference Srivastava, Roberts, Clemmer, Silbert, Lechman and Grest2021; Singh et al. Reference Singh, Ness, Sharma, de, Juan and Jaeger2024), we sketch
$\phi _m(\alpha )$
in these two limits in figure 2(c). We compute these curves using an analytic packing model (Yu & Standish Reference Yu and Standish1991), but note that their key features are insensitive to the specific choice of model (Anzivino et al. Reference Anzivino, Casiulis, Zhang, Moussa, Martiniani and Zaccone2023; Roquier Reference Roquier2024). The limits for the monodisperse cases at
$\alpha =(0,1)$
are set to
$\phi _m^{\mu =0}=0.64$
and
$\phi _m^{\mu =1}=0.57$
. In binary mixtures,
$\phi _m(\alpha )$
is non-monotonic and asymmetric, and for the limiting case of
$\mu =0$
and
$\varDelta \to \infty$
one finds the maximal
$\phi ^{\textit{max}}_m=\phi _m^{\mu =0} + (1-\phi _m^{\mu =0}) \phi _m^{\mu =0}=0.870$
occurring at
$\alpha = \phi _m^{\mu =0} \! /\phi ^{\textit{max}}_m=0.735$
. For finite
$\varDelta$
, the pure species jamming points are not reached in the mixture, so that
$\phi ^{\textit{max}}_m$
and its corresponding
$\alpha$
are reduced. At
$\varDelta =5$
, the maximum remains near
$\alpha \approx 0.7$
in both the frictional and frictionless cases. Due to this asymmetry, both the frictional and frictionless
$\phi _m$
curves increase more rapidly when
$\alpha$
decreases from an initially granular system (
$\alpha =1$
) compared with increasing
$\alpha$
from an initially colloidal system (
$\alpha =0$
).
Varying
$\alpha$
at fixed
$\phi =\phi _0$
in either the intermediate- or high-
$\textit{Pe}_{\textit{cg}}$
regimes alters the ratio
$\phi /\phi _m$
between our fixed packing fraction (horizontal dotted line in figure 2
c) and the relevant
$\alpha$
-dependent
$\phi _m$
. Close to jamming, we expect that reducing this ratio reduces the suspension viscosity (1.1). This is consistent with our viscosity measurements in these two highlighted flow regimes, with the viscosity reduction at intermediate
$\alpha$
most pronounced in the frictional high-
$\textit{Pe}_{\textit{cg}}$
regime due to the closer proximity to jamming.
If we vary
$\alpha$
at fixed
$\phi$
now at fixed intermediate
$\textit{Pe}_c=\mathcal{O}(1)$
within the shear-thickening transition, the relevant jamming point shifts between the two branches as we vary
$\alpha$
, with
$\phi _m\to \phi _m^{\mu =0}$
in the colloidal (
$\alpha \to 0$
) limit and
$\phi _m\to \phi _m^{\mu =1}$
in the granular (
$\alpha \to 1$
) limit. This is sketched schematically by the solid blue line in figure 2(c). As in the fixed
$\textit{Pe}_{\textit{cg}}$
examples, the viscosity again varies non-monotonically with
$\alpha$
but is now notably higher in the granular limit, reflecting the asymmetry in
$\phi _m(\alpha )$
.
Adding particles of a different species to an initially monodisperse suspension (fixing both the liquid volume and volume of the initial species, as in figure 1) increases
$\phi$
(red arrows in figure 2
c). Comparing the sketched
$\phi (\alpha )$
with the respective
$\phi _m$
lines, adding grains moves
$\phi$
closer to
$\phi _m$
for any choice of
$\textit{Pe}$
, whereas adding colloids initially moves
$\phi$
further from
$\phi _m$
. This holds even in the large-
$\textit{Pe}_{\textit{cg}}$
limit, implying an initial viscosity reduction when adding smaller particles to a suspension of larger ones even when both species are non-Brownian (and fully frictional). At intermediate
$\textit{Pe}_c=\mathcal{O}(1)$
, the shift in
$\phi _m$
between frictional and frictionless jamming points marked by the blue line in figure 2(c) suggests an enhanced viscosity reduction when adding colloids to a granular suspension near the colloidal-to-granular transition.
Jamming volume fraction at the colloidal-to-granular transition. (a) Viscosity divergence with
$\phi$
for
$\alpha = 0.4$
, at a range of
$\textit{Pe}_{\textit{cg}}$
. Data at
$\textit{Pe}_{\textit{cg}} = 10^6$
are fit to
$\eta(\phi) \sim (1-(\phi/\phi_m))^{-\beta}$
, shown by the black dashed line. Three representative
$\phi _m$
for
$\textit{Pe}_{\textit{cg}}=100,10,0.1$
are marked left-to-right by black vertical lines, and correspond to the black points in (b). Panel (b) shows
$\phi _m$
measured across a range of
$\textit{Pe}_{\textit{cg}}$
and
$\alpha$
, with the former indicated by the same colour legend as (a), and the latter with the same marker shapes as in figure 1. The solid and dashed black lines represent model predictions (Yu & Standish Reference Yu and Standish1991) and provide geometric
$\phi _m$
estimates taking the values at
$\alpha =(0,1)$
as inputs. (c) Example of
$\phi _m(\alpha )$
at fixed
$\textit{Pe}_c=1$
, obtained by interpolating through the fixed
$\textit{Pe}_{\textit{cg}}$
data in (b). Shown for comparison is
$\phi _m(\alpha )$
measured at fixed
$\textit{Pe}_{\textit{cg}}=1$
.
4. Results: mapping the jamming point
We now test this schematic picture using our simulation model. To do so, we first map the jamming point
$\phi _m(\alpha , \textit{Pe}_{\textit{cg}})$
over a wide range of our mixed Péclet number
$10^{-2} \leq \textit{Pe}_{\textit{cg}} \leq 10^6$
. This is done from steady-state viscosity measurements on systems spanning a range of
$\alpha$
with narrow increments of
$\phi$
to clearly resolve the viscosity divergence. Example results for
$\alpha = 0.4$
are shown in figure 3(a). At large
$\textit{Pe}_{\textit{cg}}$
, the viscosity follows a power-law divergence as
$\phi \to \phi _m$
(1.1, black dashed line), as expected for an athermal suspension. Reducing the Péclet number
$\textit{Pe}_{\textit{cg}} \leq 10^3$
changes the form of
$\eta (\phi )$
, with the smooth power-law divergence replaced by an abrupt viscosity jump. We identify
$\phi _m$
at these higher
$\textit{Pe}_{\textit{cg}}$
from the maximum of the gradient
$\partial \eta /\partial \phi$
(black vertical lines in figure 3
a and corresponding black points in figure 3
b), as the behaviour of
$\eta (\phi )$
below this point gives no indication of imminent jamming. This change in behaviour can be seen clearly in the data for
$\textit{Pe}_{\textit{cg}}=10^2$
in figure 3(a), where the viscosity is dramatically reduced relative to that at
$\textit{Pe}_{\textit{cg}}=10^6$
though there is no discernible change in the jamming point at
$\phi _m=\phi _m^{\mu =1}$
. This suggests that the weak Brownian motion is sufficient to decrease the contact stress but not to significantly change the jamming point.
Further reducing
$10^{-2} \lt {Pe}_{\textit{cg}} \leq 10^{2}$
, the jamming point
$\phi _m$
shifts upwards approaching the frictionless limit
$\phi _m^{\mu =0}$
, suggesting that within this range of Péclet number Brownian motion is now sufficient to mobilise and release frictional contacts and thus shift the jamming point. At small
$\textit{Pe}_{\textit{cg}}\lesssim 10$
, the viscosity curves shift vertically upward as the shear rate is reduced, reflecting an increasing contribution from the Brownian stress.
The full
$\phi _m(\alpha )$
curves at fixed values of
$\textit{Pe}_{\textit{cg}}$
shown in figure 3(b) validate the schematic picture sketched in figure 2(c), with the high-shear (
$\textit{Pe}_{\textit{cg}}\geq 10^2$
) and low-shear (
$\textit{Pe}_{\textit{cg}}\sim 10^{-2}$
) branches matching our estimates based on the monodisperse frictional and frictionless jamming points. At intermediate
$\textit{Pe}_{\textit{cg}}$
,
$\phi _m(\alpha )$
falls between these two limiting branches, with symmetric monodisperse limits
$\phi _m(\alpha =0, \textit{Pe}_{\textit{cg}})=\phi _m (\alpha =1, \textit{Pe}_{\textit{cg}})$
. Here
$\textit{Pe}_{\textit{cg}}$
is constant along each plotted line, so that the effective roles of Brownian versus convective transport are constant despite changes in composition. The transition between these two limiting branches gives the
$\textit{Pe}_{\textit{cg}}$
-controlled shear thickening shown in figure 2(b). This increase in
$\phi _m$
with decreasing
$\textit{Pe}_{\textit{cg}}$
mirrors the shift from frictional to frictionless
$\phi _m$
in bidisperse non-Brownian suspensions with a similar size ratio (Singh et al. Reference Singh, Ness, Sharma, de, Juan and Jaeger2024), indicating that thermal motion here plays a key role freeing frictional constraints (but see Hofmann et al. Reference Hofmann, Dormann, Liebchen and Schmid2025).
We next map these jamming curves back into
$\textit{Pe}_c$
space to describe flows where the absolute shear rate is fixed, independent of composition. To obtain
$\phi _m$
for fixed
$\textit{Pe}_c$
and
$\alpha$
, we first compute the corresponding
$\textit{Pe}_{\textit{cg}}$
and then use nearest-neighbour linear interpolation with the measured
$\phi _m(\alpha ,\textit{Pe}_{\textit{cg}})$
results shown in figure 3(b) to estimate
$\phi _m(\alpha ,\textit{Pe}_c)$
. As initially sketched in figure 2(c), these
$\phi _m(\alpha )$
profiles at fixed
$\textit{Pe}_c$
are no longer constrained to be symmetric in the two monodisperse limits. Plotting
$\phi _m(\alpha )$
at fixed
$\textit{Pe}_c=1$
, figure 3(c), the jamming point coincides with the value at
$\textit{Pe}_{\textit{cg}}=1$
at
$\alpha =0$
(by definition), closer to the frictionless branch, but peels below this curve as
$\alpha$
increases and transitions to the frictional branch in the granular (
$\alpha \to 1$
) limit. Finally, with these tools to compute
$\phi _m(\alpha )$
at fixed
$\textit{Pe}_c$
, we now have a comprehensive framework to understand the differing effects of adding either colloids or grains to an initially monodisperse suspension of the other species.
We show interpolated plots of
$\phi _m(\alpha )$
at fixed
$\textit{Pe}_c=0.01$
,
$1$
and
$100$
in figure 4(a)(i),(b)(i),(c)(i), respectively, capturing cases where
$\phi _m(\alpha )$
either remains along the upper frictionless branch (figure 4
a), remains along the lower frictional branch (figure 4
c) or transitions between the two as
$\alpha$
is varied (figure 4
b). The inaccessible jammed region above these curves is shaded in grey. Solid symbols in the same panels show the volume fraction
$\phi$
when either adding grains to a colloidal suspension (green) or colloids to a granular suspension (purple) starting from differing initial volume fractions (
$\phi =0.55$
,
$0.6$
,
$0.62$
). As in figure 1, this particle addition is done holding both the liquid volume
$V_l$
and the volume of the initial species (either
$V_c$
or
$V_g$
) fixed.
Effect of particle addition on viscosity at three fixed values of
$\textit{Pe}_c$
. Shown in (a)(i), (b)(i) and (c)(i) are theoretical predictions (black lines) together with interpolated plots of
$\phi _m(\alpha )$
at
$\textit{Pe}_c=0.01$
,
$1$
and
$100$
(solid blue lines). Solid green and purple markers in (i) are not jamming points but represent changes in
$\phi$
under particle addition. Shaded regions show parameter values for
$\alpha$
and
$\phi$
for which the system is jammed. For each value of
$\textit{Pe}_c$
, we explore six cases of particle addition, adding either grains (green) or colloids (purple) to suspensions initially at
$\alpha =0$
and
$\alpha =1$
, respectively, initially with
$\phi =0.55$
,
$0.6$
and
$0.62$
(shown by light to dark lines and markers coloured green (grain addition) and purple (colloid addition)). Panel rows (ii) and (iii) show the consequent variations in the viscosity
$\eta$
with volume fraction
$\phi$
for addition of grains (ii) and colloids (iii). The inset of (c)(ii) shows the time-averaged large particle radial distribution function
$g(r^\dagger )$
corresponding to colloid addition at
$\textit{Pe}_c=100$
for an initial
$\phi =0.6$
. The growth of the peak at
$r^\dagger =r(a_i+a_j)\approx 1.1$
is due to colloids disrupting grain contacts as sketched in the inset of (c)(iii).
Adding grains to an initially colloidal (
$\alpha =0)$
suspension moves
$\phi (\alpha )$
closer to the jamming point
$\phi _m(\alpha )$
independent of
$\textit{Pe}_c$
or the initial
$\phi$
, so that the corresponding
$\eta (\phi )$
curves always increase monotonically, figure 4(a)(ii),(b)(ii),(c)(ii), eventually diverging as they cross the
$\phi _m$
boundary. This accounts for the monotonic rise in
$\eta$
shown in figure 1(a). Interestingly, the measured viscosities
$\eta (\phi )$
below jamming are lower for
$\textit{Pe}_c=1$
than at
$\textit{Pe}_c=0.01$
, despite the former being closer to jamming. This reflects the decreasing contribution of the Brownian stress at higher
$\textit{Pe}_c$
, also seen in figure 3(a). As
$\textit{Pe}_c$
increases, the window of accessible volume fractions shrinks as particle friction becomes more important, so that in figure 4(c) the systems starting at
$\alpha =0$
and
$\phi =0.6$
,
$0.62$
are already above their respective jamming points. Thus the scope for manipulating the viscosity of initially colloidal suspensions by adding grains is reduced with increasing
$\textit{Pe}$
as the jamming point shifts to lower volume fractions.
In contrast, the effect of adding smaller colloids to an initially granular (
$\alpha =1)$
suspension varies with
$\textit{Pe}_c$
. The asymmetry in the
$\phi _m(\alpha )$
curves, with maxima around
$\alpha \sim 0.7$
independent of
$\textit{Pe}_c$
, results in
$\phi$
increasing with decreasing
$\alpha$
slower than the growth of
$\phi _m$
. In other words, adding colloids to a granular suspension can move the system out of the inaccessible region, allowing the jammed suspension to flow despite the increase in solids content. This flowable region, in which the viscosity nonetheless remains dominated by contributions from particle contacts, extends until the
$\phi _m(\alpha )$
boundary is again crossed at some point below the maxima in the
$\phi _m(\alpha )$
curve. This suggests in principle that adding colloids while keeping
$0.7\lt \alpha \lt 1$
should always move the suspension further away from jamming. This holds at
$\textit{Pe}_c=1$
and
$\textit{Pe}_c=100$
, where in both cases it is possible to fluidise initially jammed suspensions by adding colloidal particles, figure 4(b)(iii),(c)(iii), going from a divergent to finite viscosity with increasing
$\phi$
. Starting from a granular suspension below jamming (
$\phi = 0.55$
), the addition of colloids reduces the suspension viscosity, with this reduction most pronounced around
$\textit{Pe}_c=1$
where
$\phi _m(\alpha )$
rises most steeply with decreasing
$\alpha$
and then falls more gradually below the maximum.
We rationalise this rheological measurement at the microstructural level by computing the radial distribution function of the grains as colloids are added. From the grain positions, we obtain a list of pairwise distances
$r$
(i.e. the lengths of the centre-to-centre vectors
$\boldsymbol{r}_{i,j}$
between all grain pairs in the system), which we then normalise by the sum of the respective radii
$a_i$
and
$a_j$
of the grains in each pair as
$r^\dagger = r/(a_i + a_j)$
. Shown in the inset of figure 4(c)(ii) are time-averaged
$g(r^\dagger )$
for each of the cases shown in figure 4(c)(i) for colloid addition starting from
$\alpha =0$
and
$\phi =0.6$
. Analogous results are obtained at other
$\textit{Pe}_c$
. As colloids are added, there is a systematic decrease of the peak in
$g(r^\dagger )$
at
$r^\dagger =1$
and growth of a wider peak around
$r^\dagger =1.1$
. This latter value represents the expected distance between two particles of radius 5 or 7, separated by a single small particle of radius 1 or 1.4, indicating that colloids disrupt grain–grain contacts. In the inset of figure 4(c)(iii), we sketch this scenario, in which grain–grain interactions are effectively mediated by single layers of colloidal particles. The ability of the colloids to roll at these contact points potentially provides the degree of freedom necessary to facilitate unjamming.
At
$\textit{Pe}_c = 0.01$
, where Brownian contributions to the viscosity are more significant, the effect of adding colloids is more complicated. We do not observe a viscosity reduction when adding colloids to a granular suspension at a relatively low initial packing density of
$\phi =0.55$
(figure 4
a(iii), see also the dot-dashed orange line in figure 1
d). In this case, given our greater initial distance from jamming, (1.1) likely requires higher-order terms to fully describe
$\eta (\phi )$
, so that we cannot simply map between the viscosity and the distance to jamming. If we instead consider initial packing densities either just above or very close to jamming (
$\phi =0.62$
and
$\phi =0.6$
), we recover the unjamming effect and initial viscosity reduction with the addition of colloids, demonstrating that this framework for understanding the effect of added colloids holds even at low
$\textit{Pe}_c$
provided we are sufficiently close to jamming.
Ultimately our modelling predicts that it is the asymmetry in the jamming plots that allows for the contrasting effect of adding colloids or grains. Adding grains always moves
$\phi$
closer to
$\phi _m$
; adding colloids can move
$\phi$
further from
$\phi _m$
, exaggerated when considering fixed
$\textit{Pe}_c$
so that
$\phi _m$
is both
$\textit{Pe}$
- and
$\alpha$
-dependent. Interestingly, we can capture this phenomenology considering only the frictional and frictionless limits of
$\phi _m$
, largely ignoring the colloidal glass transition at some
$\phi _G\lt \phi _m^{\mu =0}$
even at our lowest
$\textit{Pe}_c=0.01$
. While we would expect this glass transition to play a critical role in the
$\textit{Pe}_c\to 0$
limit, our results indicate that the dominant role of thermal motion in the range of finite shear rates explored here is simply to inhibit frictional particle contacts. The phenomenology reported here is likely sensitive to the polydispersity of each species. For instance, a broad size distribution in the colloidal phase could substantially increase
$\phi _m(\alpha =0)$
which, dependent on details, could enhance or suppress the viscosity reduction associated with colloid addition.
5. Conclusion
We have systematically explored the constitutive curves of suspensions comprising mixed colloids and grains, demonstrating the complex dependence on the composition
$\alpha$
and the shear rate quantified through
$\textit{Pe}_c$
and
$\textit{Pe}_{\textit{cg}}$
. As well as the expected shear thinning and thickening behaviour, the simulation predicts a counterintuitive viscosity drop upon addition of colloids to an initially pure granular suspension, an effect that is most pronounced at shear rates near the colloidal-to-granular crossover. We rationalised all of the observed flow curves by mapping the jamming volume fraction as a function of
$\alpha$
and
$\textit{Pe}$
. By introducing a rescaled Péclet number, we quantified the interplay between shear flow and thermal motion, providing a clear separation of granular, colloidal and crossover regimes for bidisperse Brownian suspensions. Additionally, we identified that
$\phi _m$
increases as thermal energy increases, which aligns with the shift from frictional to frictionless behaviour.
These findings provide important insights for the development and optimisation of materials where viscosity control is critical. The ability to predict and manipulate the viscosity of bidisperse suspensions by adjusting
$\alpha$
,
$\dot \gamma$
and
$k_BT$
has potential applications in industries such as concrete formulation (Roussel et al. Reference Roussel, Lemaître, Flatt and Coussot2010), dip coating (Jeong et al. Reference Jeong, Lee, Thiévenaz, Bazant and Sauret2022) and printing (Thiévenaz et al. Reference Thiévenaz, Rajesh and Sauret2021). Future research could focus on extending this framework to more complex particle shapes and interactions, offering broader applicability to real-world systems.
Funding
C.N. acknowledges support from the Royal Academy of Engineering under the Research Fellowship scheme and from the Leverhulme Trust under Research Project Grant RPG-2022-095. For the purpose of open access, the authors have applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising from this submission.
Declaration of interests
The authors report no conflict of interest.
Data availability statement
Scripts required to reproduce the data reported in this paper are available on request.
Appendix A. Simulation details
A comprehensive description of our model is given in Li et al. (Reference Li, Royer and Ness2024). In short, particle motion is governed by a set of Langevin equations incorporating forces
$\boldsymbol{F}$
originating from direct particle contacts, hydrodynamics and Brownian motion. The translational equation of motion for particle
$i$
with position
$\boldsymbol{x}_i$
and mass
$m_i$
is
\begin{align} m_i \frac {\textrm{d}^2\boldsymbol{x}_i}{{\textrm d}t^2} = \sum _j \boldsymbol{F}_{i,j}^{C} + \boldsymbol{F}_{i}^{H,D} + \sum _j \boldsymbol{F}_{i,j}^{H,L} + \boldsymbol{F}_{i}^{B,D} + \sum _j \boldsymbol{F}_{i,j}^{B,L}\text{.} \end{align}
Forces with subscript
$i$
are single-body whereas those with subscript
$i,j$
act pairwise between particles
$i$
and
$j$
. Superscripts C, H, B, D and L refer to forces arising due to contacts (C), hydrodynamics (H) and Brownian (B) effects, with the latter two having drag (D) and lubrication (L) terms. An equivalent set of torques are computed and used to update the rotational velocities of the particles, but these are omitted here for brevity.
Contact forces are modelled as linear springs, with the normal force between
$i$
and
$j$
proportional to their overlap
$\delta _{i,j}$
:
\begin{align} \boldsymbol{F}_{i,j}^{C} = k_n \delta _{i,j} \boldsymbol{n}_{i,j} - k_t \boldsymbol{\xi }_{i,j}\text{,} \end{align}
with
$k_n$
and
$k_t$
the spring constants and
$\boldsymbol{n}_{i,j}$
the unit vector pointing from
$i$
to
$j$
. The tangential force is linear in
$\boldsymbol{\xi }_{i,j}$
, the accumulated displacement of the particle pair perpendicular to
$\boldsymbol{n}_{i,j}$
since the initiation of the contact. Tangential forces are regulated by a Coulombic friction coefficient set as
$\mu =1$
throughout. The contact stress tensor
${\Sigma }^{C}$
is obtained by summing the outer product
$-\boldsymbol{F}_{i,j}^{C}\otimes \boldsymbol{r}_{i,j}$
(with
$\boldsymbol{r}_{i,j}$
the centre-to-centre vector) over all contacting pairs.
Hydrodynamic interactions comprise Stokes drag and pairwise lubrication. The former is given for particle
$i$
as
\begin{align} \boldsymbol{F}_{i}^{H,D} = 6 \pi \eta a_i \left(\boldsymbol{U}^{\infty }(\boldsymbol{x}_i) - \boldsymbol{U}_{\! i} \right)\!\text{,} \end{align}
and the latter are given for
$i$
and
$j$
as
\begin{align} \begin{split} \boldsymbol{F}_{i,j}^{H,L} = & \big(X_{11}^{A} \mathbb{N}_{i,j} + Y_{11}^{A}\mathbb{ T}_{i,j} \big) \left(\boldsymbol{U}_{\! j} - \boldsymbol{U}_{\! i} \right) \\ &+ Y_{11}^{B}(\boldsymbol{\varOmega }_i \times \boldsymbol{n}_{i,j}) + Y_{21}^{B}(\boldsymbol{\varOmega }_j \times \boldsymbol{n}_{i,j})\text{.} \end{split} \end{align}
Here
$\mathbb{N}$
and
$\mathbb{T}$
are the normal and tangential projection operators,
$\boldsymbol{U}_{\! i}$
is the velocity of particle
$i$
,
$\boldsymbol{U}^{\infty }(\boldsymbol{x}_i)$
is the liquid streaming velocity at the position of particle
$i$
and
$\boldsymbol{\varOmega }_i$
is the rotational velocity of particle
$i$
. The scalar prefactors
$X$
and
$Y$
encode the geometry of the interacting pair and are given by Jeffrey (Reference Jeffrey1992). These are obtained by solving for the inertia-free Stokes flow in the narrow gaps between neighbours, assuming no-slip of the viscous fluid on the particle surfaces. This produces a normal force associated with squeezing motion of the particle pair and tangential forces associated with rotations of the particles. These forces have been written in pairwise, frame-invariant form amenable to particle-based simulation by several authors, most notably Jeffrey & Onishi (Reference Jeffrey and Onishi1984); Jeffrey (Reference Jeffrey1992); Kim & Karrila (Reference Kim and Karrila2013). Though these expressions are limited to modest size ratios, more recent work has provided exact solutions applicable for arbitrary size ratios at arbitrary separations (Goddard, Mills-Williams & Sun Reference Goddard, Mills-Williams and Sun2020). The hydrodynamic stress tensor
${\Sigma }^{H}$
has contributions from drag and lubrication, the latter obtained by summing
$-\boldsymbol{F}_{i,j}^{H}\otimes \boldsymbol{r}_{i,j}$
over all interacting pairs.
Brownian forces are generated to satisfy fluctuation–dissipation theorem, with the single-body and pairwise forces given, respectively, as
\begin{align} \boldsymbol{F}_{i}^{B,D} = \sqrt {\frac {2k_b T}{\Delta t}} \sqrt {6 \pi \eta a_i} \boldsymbol{\psi }_i\text{,} \\[-28pt] \nonumber \end{align}
\begin{align} \boldsymbol{F}_{i,j}^{B,L} = \sqrt {\frac {2k_b T}{\Delta t}} \left (\sqrt {X_{11}^{A}} \mathbb{N}_{i,j} + \sqrt {Y_{11}^{A}} \mathbb{T}_{i,j} \right ) \boldsymbol{\theta }_{i,j}\text{,} \\[3pt] \nonumber \end{align}
where
$k_b T$
is the thermal energy,
$\Delta t$
is the time step and
$\boldsymbol{\psi }_i$
,
$\boldsymbol{\theta }_{i,j}$
are unit vectors of random direction.
The model produces constitutive curves that, for near-monodisperse suspensions, compare well with other numerical methods (Li et al. Reference Li, Royer and Ness2024). The parameters used to generate the data in figure 1 are as follows, with the dimensions of each quantity given in terms of mass (
$M$
), length (
$L$
) and time (
$T$
), and (–) representing dimensionless quantities. Small particles have
$a_S = 1$
$(L)$
,
$k_n = 10\,000$
$(M/T^2)$
,
$k_t = 7\,000$
$(M/T^2)$
and density
$\rho = 0.1$
$(M/L^3)$
. Large particles have
$a_L = 5$
$(L)$
,
$k_n = 50\,000$
$ (M/T^2)$
,
$k_t = 35\,000$
$ (M/T^2)$
and
$\rho = 0.004$
$ (M/L^3)$
. The liquid viscosity is
$\eta = 0.1$
$(M/L T)$
and the time step is
$\Delta t = 0.0001$
$(T)$
. These parameters ensure approximately inertia-free and hard-sphere conditions as described in Li et al. (Reference Li, Royer and Ness2024). The values of
$\dot \gamma$
and
$k_BT$
used to achieve each desired Péclet number are given in table 1. For the range of parameters used, we did not observe evidence of depletion interactions between the grains.
Parameters used to generate the Péclet numbers in figures 1(b) and (c), and their dimensions.
Open access
