Thinning or thickening? Multiple rheological regimes in dense suspensions of soft particles

Understanding the shear rheology of dense colloidal and granular suspensions remains a central challenge at the crossroad between nonequilibrium statistical mechanics and soft matter physics, with a clear technological relevance [1–3]. Simple liquids display simple rheological properties characterized by linear Newtonian behaviour [4]. In a simple shear flow, for instance, the rate of deformation, $\dot{\gamma}$, is proportional to the applied shear stress, σ, such that the viscosity $\eta = \sigma/ \dot{\gamma}$ uniquely characterizes the rheological response.

In dense particle suspensions such as emulsions, colloidal assemblies, or granular materials, the viscosity is usually a nonlinear function of the applied flow rate. To characterise these materials, an entire flow curve $\eta = \eta(\dot{\gamma})$ is thus needed [1–3]. Because the applied deformation now determines the response of the system, understanding nonlinear flow curves obviously requires a more detailed analysis, which must deal with both nonlinear and nonequilibrium effects. When the viscosity varies with the applied shear rate, the system can either flow more easily as $\dot{\gamma}$ increases (shear-thinning), or offer increasing resistance to flow (shear-thickening). We are familiar with both these effects, as most of the complex fluids used for cosmetics or in food products display these amusing nonlinearities, which can be technologically both useful and annoying [3,5].

In practice, most of the experimental flow curves measured even in model suspensions display a complex mixture of both these nonlinear effects [2,3,5]. As two typical examples, we show flow curves measured in a colloidal dispersion of latex particles [5,6] (diameter $a = 250\ \text{nm}$, fig. 1(a)), and in an oil-in-water emulsion [7] (diameter $a=20\ \mu \text{m}$, fig. 1(b)). For a given volume fraction φ, the flow curves may display an initial Newtonian regime at low enough $\dot{\gamma}$ and φ, or a strong shear-thinning regime when φ is larger. This thinning regime is followed, for intermediate φ and larger $\dot{\gamma}$, by a Newtonian plateau regime. At larger $\dot{\gamma}$, shear-thickening sets in, and the magnitude of the viscosity increase clearly depends on the density regime. In some cases, shear-thickening is interrupted and flow curves display a viscosity maximum. Finally shear-thickening is not observed when density is too large, see, for instance, the large density data in fig. 1(b). Books and reviews, of course, offer an even broader range of possible behaviours [1–3,5,8,9], but the data in figs. 1(a), (b) are representative of the typical behaviour of dense suspensions.

Zoom In Zoom Out Reset image size

The primary purpose of this work is to show that a simple model of soft repulsive particles can exhibit a similarly complex rheology, despite the fact that it does not incorporate several of the physical ingredients usually put forward to account for these nonlinearities. We argue that understanding first such a simple model is useful before turning to more complicated physical explanations which can only increase the complexity of the results. In our model of spherical particles, we use a canonical pairwise, purely repulsive particle interaction. Dynamics is controlled by a Langevin dynamics, which we simulate either at finite or zero temperature. The system is homogeneously sheared (no shear-banding) at constant shear rate $\dot{\gamma}$ in a simple shear geometry, without any kind of density or stress gradients. We do not introduce more complex ingredients such as hydrodynamic interactions, attractive forces, lubrication forces, or frictional contacts between the particles. Our model, we believe, represents an ideal starting point to understand the basic physics of suspensions, to which more ingredients could then be added, if needed.

In fig. 1(c), we present a set of flow curves obtained from our simulations, for a specific set of parameters, which we shall discuss shortly. For moderate volume fractions, we obtain a complex succession of Newtonian, shear-thinning, Newtonian, shear-thickening, shear-thinning regimes as the flow rate is varied. The first Newtonian regime disappears as φ is increased, and at even larger φ the intermediate Newtonian and thickening regimes also disappear, such that the flow curve becomes purely shear-thinning at large density. The existence of these multiple rheological regimes and their evolution with volume fraction is in excellent qualitative agreement with the experimental flow curves presented in figs. 1(a), (b). Given the extreme simplicity of our model, we believe that such an agreement is an important achievement of our work.

Shear-thinning is the most commonly observed nonlinearity in experimental flow curves [1]. This effect typically results from the competition between the complex structure of the suspension under study that is responsible for the large viscosity of the system at rest, and the flow rate which tends to disrupt this static organization. For dense suspensions, large viscosities emerge when either the glass transition (for thermal systems [10]), or the jamming transition (for athermal ones [11]) is approached. In both cases, shear-thinning occurs when the shear rate competes with the relevant timescale associated to structural relaxation.

By contrast, shear-thickening is typically less common in experiments [8], and its detailed characterization is more difficult because it is often mixed with issues such as flow heterogeneities, shear bands, density gradients, or particle migration [12,13]. However, shear-thickening has been documented for a large number of dense suspensions of various types, for soft and hard particles, for small colloids and large granular particles [3,5,8]. This variety of systems has led to a similar variety of theoretical arguments and models to account for shear-thickening, but a broad consensus has yet to emerge. In some systems, a strong increase of the viscosity can even become a nearly discontinuous jump, which we do not observe in our simple system. Another regime that we do not explore consists in suspensions at moderate densities, where a shear-thickening of modest amplitude likely results from hydrodynamic interactions promoting cluster formation under shear [5,14].

Recent experiments in dense suspensions show that large-amplitude shear-thickening is observed in the vicinity of the jamming transition, suggesting that particle crowding promotes shear-thickening [9,13,15–22]. However, the emergence of a yield stress above the jamming transition seems to suppress shear-thickening [9,15]. These observations suggest that neglecting hydrodynamic interactions to focus instead on the specific competition between steric constraints and the shear flow is an interesting path.

Recent works [13,18,20] explored the idea that a strong continuous increase of the viscosity under shear results from a change of dissipation mechanism from viscous to inertial damping, in the spirit of Bagnold [23]. We extend these ideas to include also the effect of thermal fluctuations and particle softness, which are experimentally relevant. In recent work, we analysed in detail the complex interplay of thermal fluctuations and steric constraints in the overdamped limit [24,25]. Here we add inertia, which introduces another timescale and adds a new level of complexity to flow curves that are already nontrivial.

For completeness, we mention recent studies invoking frictional forces to explain discontinuous shear-thickening [16,17,19,20,22]. Numerical studies seem in good qualitative agreement with experiments, although quantitative understanding and detailed comparison to experimental data are unavailable.

We study a system composed of N particles with equal mass m, interacting through a purely repulsive truncated harmonic potential [26], $V(r) = \frac{\epsilon}{2} (1 - r/a)^2 \Theta (a-r)$, where a is particle diameter, an energy scale, and r the interparticle distance; $\Theta(x)$ is the Heaviside function. In practice, to avoid crystallization which might occur at large density, we use a 50 : 50 binary mixture of spheres with diameter ratio 1.4. The resulting volume fraction is $\varphi = \frac{\pi}{12L^3} N (a^3 + (1.4 a)^3)$, where L is the linear size of the simulation box. We use $N=10^3$ particles in a cubic box in three spatial dimensions, and could not find significant finite-size effects in additional simulations performed with $N=10^4$ particles. We inforce Lees-Edwards periodic boundary conditions [27].

These soft particles evolve with Langevin dynamics,

$$\begin{equation} m{\frac{{\rm d} {\vec{v}_i}}{{\rm d} {t}}} + \xi (\vec{v}_i - \dot{\gamma} y_i\vec{e}_x) + \sum_{j \neq i} {\frac{\partial {V(|\vec{r}_i - \vec{r}_j|)}}{\partial {\vec{r}_i}}} + \vec{f}_i = 0, \end{equation} \tag{ 1 }$$

where $\vec{r}_i$ and $\vec{v}_i$ respectively represent the position and velocity of particle i. The first term is the particle acceleration, and the second represents the viscous damping, controlled by the damping constant ξ. The system is sheared in the xy plane, and advection occurs along the x-axis. The shear rate is $\dot{\gamma}$, y_i is the y-coordinate of particle i, and $\vec{e}_x$ is the unit vector along the x-axis. The third term in eq. (1) incorporates pairwise harmonic repulsion between the particles, and the final term is the Brownian random force acting on particle i, which we draw from a Gaussian distribution with zero mean and variance obeying the fluctuation-dissipation relation, $\langle \vec{f}_i(t)\vec{f}_j(t^{\prime})^{T}\rangle=2k_BT\xi\delta_{ij}\textbf{1} \delta(t-t^{\prime})$, where $k_{\text{B}}$ is the Boltzmann constant and T the temperature.

We impose the shear rate $\dot{\gamma}$ and measure the shear stress σ as the xy component of stress tensor using the Irving-Kirkwood formula [27]. We deduce the shear viscosity, $\eta = \sigma/\dot{\gamma}$.

While very simple, eq. (1) contains a number of distinct ingredients which typically control the physics of dense suspensions: inertial forces, viscous damping, particle interactions, and thermal noise. Their competition is more easily understood by introducing characteristic timescales to each term. We can define three independent timescales. The damping time of particle velocity is given by $\tau_v = m/\xi$. Energy is dissipated through viscous damping over a timescale ${\tau_0 = \xi a^2 /\epsilon}$, while thermal fluctuations occur over a Brownian timescale $\tau_T = \xi a^2 /(k_BT)$.

We can construct two dimensionless parameters out of these three timescales, which we will use to quantify the relative effects of inertia, thermal fluctuations and viscous damping. Comparing $\tau_v$ to $\tau_0$, we can create a dimensionless particle mass, $m^{\star} = \tau_v/\tau_0 = m \epsilon/(\xi a)^2$. In the same spirit, comparing $\tau_T$ to $\tau_0$ allows the definition of a dimensionless temperature, $T^{\star} = \tau_0/\tau_T = k_BT/\epsilon$. We use both $T^{\star}$ and $m^{\star}$ to specify the values of the control parameters employed in a given study, which amounts to choosing $\tau_0$ as our microscopic time unit. Therefore, a set of simulations is fully specified by the values of $(m^{\star}, T^{\star})$, for which we can then vary both the packing fraction φ and the flow rate $\dot{\gamma}$.

We use the "solvent" viscosity as a unit of viscosity, which is defined through the Stokes law, $\eta_s = \xi/ (3 \pi a)$. Accordingly, the natural stress scale is $\sigma_s = \epsilon / a^3$.

Having defined the relevant timescales, we are now in a position to properly identify and interpret the multiple rheological regimes observed in the simulation results presented in fig. 1(c), which were obtained by fixing $m^\star = 10^3$ and $T^\star = 10^{-7}$. As should now be obvious, these values imply a clear separation of the relevant timescales, namely $\tau_T \gg \tau_v \gg \tau_0$. As a result, the imposed shear rate successively competes with all three timescales, which directly impacts the measured flow curves at various densities shown in fig. 1(c), as we now explain.

At $\varphi \leq 0.56$ a Newtonian behaviour is observed in the limit $\dot{\gamma} \to 0$, which gives way to shear-thinning behaviour as $\dot{\gamma}$ is increased. Both these behaviours occur in the regime $\dot{\gamma} \tau_T \lesssim 1$, which implies that thermal fluctuations control this regime. ($\dot{\gamma} \tau_T$ is the "bare" Péclet number.) The corresponding Newtonian viscosity, $\eta_T$, is therefore directly related to the equilibrium structural relaxation timescale, $\tau_\alpha(T,\varphi)$, of the thermalized suspensions of harmonic spheres [28], which can be calculated through equilibrium linear response theory. Instead, shear-thinning is characteristic of the nonlinear rheology of viscous glassy fluids [10]. It occurs when structural relaxation is provoked by the imposed shear flow (i.e. when the "dressed" Péclet number $\dot{\gamma} \tau_{\alpha}$ is not small).

At higher shear rate, $\dot{\gamma} \tau_T \gtrsim 1$, thermal fluctuations cannot affect the physics which thus becomes equivalent to zero-temperature (or "athermal") rheology [24]. In this regime, we first observe a Newtonian behaviour, with an apparent viscosity $\eta_0$ which has a different value than in the first Newtonian regime, $\eta_0 \neq \eta_T$. The difference between these thermal and athermal Newtonian regimes for the overdamped case with $m^{\star} =0$ is discussed in ref. [24]. Comparison with these earlier results reveals perfect agreement, which is expected as long as inertial effects do not affect the physics, i.e. when $\dot{\gamma} \tau_v$ is small.

Increasing further the shear rate, deviations from the overdamped limit start to appear when $\dot{\gamma} \tau_v$ becomes large [13,18]. In this regime, we observe a succession of continuous shear-thickening followed by a shear-thinning regime. As a result, the viscosity exhibits a maximum at a well-defined $\dot{\gamma}$ value. The shear-thickening and the viscosity maximum are not present in the overdamped simulations [24], and are analysed in more detail below.

Having properly identified the various regimes, we turn to the evolution of the flow curves with volume fraction. The value of the viscosity $\eta_T$ in the thermal Newtonian regime increases rapidly with φ as the glass transition density is approached, with $\varphi_g \approx 0.59$ [28], and the thermal regime is characterized by a strong shear-thinning behaviour, $\eta \propto \dot{\gamma}^{-1}$, as is typical for a glassy material with a finite yield stress.

The value of the second Newtonian viscosity $\eta_0$ in the athermal regime at larger $\dot{\gamma}$ also displays a strong increase with volume fraction, but with a density dependence distinct from the one of $\eta_T$, diverging near the jamming transition at $\varphi_J \approx 0.64$ [29].

An interesting behaviour is observed when the jamming transition is approached because $\eta_0$ increases more rapidly than the value of the viscosity maximum observed at the end of the shear-thickening regime. At a result, close enough to the jamming density these two values become equal, which implies that the shear-thickening regime eventually disappears as the volume fraction is larger than the jamming density, see fig. 1(c). Indeed, flow curves simplify above jamming and we simply observe shear-thinning behaviour for $\varphi \gtrsim 0.64$. In this density regime, the system behaves as a jammed athermal assembly of soft particles, and is therefore characterized by a finite yield stress [24,29] which entirely dominates the flow curve, such that again $\eta \propto \dot{\gamma}^{-1}$.

A comparison to experimental flow curves in figs. 1(a), (b) shows that all features described in the simulations are also present in experiments. The thermal Newtonian and shear-thinning regimes, and the glass transition physics are observed in the colloidal latex dispersion at small shear rates (small Péclet number). The athermal Newtonian viscosity is observed in both the colloidal dispersion at large enough $\dot{\gamma}$ (large Péclet number), and in the emulsion where the large droplet size ensures that thermal fluctuations are irrelevant. The strong continuous shear-thickening is also observed in both systems at larger shear rates, with a viscosity maximum also observed in some cases. Finally, when the soft emulsion is compressed above jamming, shear-thickening is not observed anymore. All these features are in excellent agreement with our numerical observations, as claimed in the introduction.

We now focus on the specific deviations brought about by the inclusion of inertia into the equation of motion, eq. (1). To simplify the discussion, it is useful to analyse the case where thermal fluctuations are completely absent, $T^{\star} = 0$, such that the competition is between viscous dissipation, inertial effects, and particle softness in the vicinity of the jamming transition. In other words, we consider the physics at large Péclet number.

We first concentrate on a fixed volume fraction below the jamming transition, $\varphi = 0.6 < \varphi_J$, and study the evolution of the flow curves as $m^{\star}$ is varied over a broad range, see fig. 2(a). For the smallest value, $m^{\star} = 10$, inertia has a negligible effect on the flow curves, which thus resemble the result obtained in the overdamped limit $m^{\star} = 0$, namely a Newtonian plateau with viscosity $\eta_0$, followed by shear-thinning at large $\dot{\gamma}$. When $m^{\star}$ is increased, inertial effects set in and the shear-thickening behaviour and viscosity maximum become apparent. Clearly, the onset of shear-thickening, the magnitude of the viscosity increase, and the viscosity maximum are strongly dependent on $m^{\star}$, showing that they directly result from inertial effects.

Zoom In Zoom Out Reset image size

In fig. 2(b), we show the same data using a rescaled shear rate $\dot{\gamma} \tau_v$, using the inertial timescale $\tau_v = m/\xi$ defined above. Because $\tau_v$ is the typical relaxation time of particle velocities due to the viscous damping, the condition $\dot{\gamma} \tau_v \gg 1$ means that shear deformation occurs faster than velocity relaxation. In this rescaled plot, the onset of shear-thickening in the flow curves at various $m^{\star}$ collapses very nicely, such that $\eta \approx \eta_0$ when $\dot{\gamma} \tau_v \lesssim 10^{-2}$, while the data suggest $\eta \propto \dot{\gamma}^2$ when $\dot{\gamma} \tau_v \gtrsim 10^{-2}$. While the scaling variable $\dot{\gamma} \tau_v$ is a straightforward choice (justified below), the very small value it takes at the crossover, $\dot{\gamma} \tau_v \approx 10^{-2}$, is less intuitive.

We now rationalize the observed behaviours by developing theoretical arguments in the spirit of kinetic theory. The sole source of dissipation in the Langevin dynamics, eq. (1), is the viscous damping. Energy balance between dissipation and energy injection by the flow per unit volume and unit time yields

$$\begin{equation} \sigma \dot{\gamma} = \rho L^3 \times \xi \bar{v}^2 / L^3, \end{equation} \tag{ 2 }$$

where $\bar{v}$ is the typical amplitude of the particle velocity, and $\rho=N/L^{3}$ is the number density. Using a kinetic theory argument, the shear stress can also be expressed as the product of a number of collisions with a typical momentum transfer at collisions. Neglecting the density dependence of the mean free path, this gives us a second (approximate) relation:

$$\begin{equation} \sigma \approx \rho L^2 \bar{v} \times m \dot{\gamma} a / L^2. \end{equation} \tag{ 3 }$$

By combining eqs. (2), (3) we obtain

$$\begin{equation} \sigma \approx (\rho m^2 a^2/\xi) \dot{\gamma}^3, \qquad \bar{v} \approx (m a / \xi ) \dot{\gamma}^2. \end{equation} \tag{ 4 }$$

Remarkably, the obtained constitutive equation between stress and shear rate fully agrees with the observation $\eta/\eta_s \propto (\tau_v \dot{\gamma})^2$ in fig. 2(b). Note that this result is different from Bagnold scaling $\eta \propto \dot{\gamma}$ [23], because the shear rheology in this regime strongly depends on the details of the energy dissipation [30]. Our results suggest therefore that qualitatively distinct shear-thickening behaviours might arise in different systems, depending on the dissipation mechanism, an issue which is of experimental relevance. Bagnold scaling is obtained when energy dissipation in eq. (1) is introduced via collisions, such that the energy balance equation in eq. (2) is modified into $\sigma \dot{\gamma} \propto \bar{v}^3$. Together with eq. (3), Bagnold scaling $\eta \propto \bar{v}^2$ would be recovered [31].

We can extend our kinetic argument to account for the effect of particle softness and explain the origin of the viscosity maximum [32]. The above analysis suggests that particles move faster as $\dot{\gamma}$ increases, see eq. (4). This implies that "collisions" become ill-defined when particles have enough kinetic energy to overcome their repulsive interactions, which happens when $\bar{v} > \bar{v}_c$ with the crossover velocity $\bar{v}_c$ given by $m \bar{v}_c^2/2 \approx \epsilon$. Combined with eq. (4) this argument provides a rough estimate of the height and location of the viscosity maximum, namely

$$\begin{equation} \eta_c = \rho a (m \epsilon)^{1/2},\qquad \dot{\gamma}_c = [ \epsilon\xi^2 / (m^3 a^2) ]^{1/4}. \end{equation} \tag{ 5 }$$

In fig. 2(c), the renormalized flow curves $\eta/\eta_c$ vs. $\dot{\gamma}/\dot{\gamma}_c$ are shown, where the height and location of the viscosity maxima are nearly collapsed, suggesting that our crude argument captures its origin, namely, the competition between inertial effects and particle softness. For collisional dissipation and soft particles, we predict again a viscosity maximum [33] but with different scaling properties controlled by $\eta_c = \rho a (m \epsilon)^{1/2}$ and $\dot{\gamma}_c = [\epsilon /(m a^2)]^{1/2}$.

So far, we analysed the shear-thickening at constant density. To obtain a dynamic state diagram of the multiple rheological regimes as determined experimentally [15], we turn to the influence of the volume fraction. We concentrate again on the zero-temperature case, for convenience. In fig. 3(a), we show flow curves with $T^{\star} = 0$, $m^{\star} = 10^5$, and various volume fractions from $\varphi=0.58$ below jamming up to $\varphi=0.66$ above, while fig. 3(b) summarizes these results in a stress/volume fraction dynamic state diagram [15]. Since we include densities above the jamming transition, this diagram is mostly useful for soft particle systems.

Zoom In Zoom Out Reset image size

The data in fig. 3(a) suggest that the onset of shear-thickening is only very weakly dependent on the volume fraction, and occurs when $\dot{\gamma} \tau_v \approx 10^{-2}$, almost independently of φ (the remaining φ-dependence in fig. 1(c) goes beyond the present scaling analysis, and might be useful to explain experimental flow curves). This indicates that the volume fraction dependences of the Newtonian viscosity $\eta_0$ and of the inertial regime $\eta \propto \dot{\gamma}^2$ are actually identical, suggesting that both regimes actually reflect the fundamental flow properties of an athermal assembly of hard spherical particles below jamming. Accordingly, the stress scale controlling the Newtonian-to-thickening transition in fig. 3(b) is given by $\sigma \approx \eta_0(\varphi) \times 10^{-2}/\tau_v$, and its evolution with φ is thus essentially controlled by the one of the viscosity $\eta_0(\varphi)$, which diverges at $\varphi_J$. Our simulations also indicate that the density dependence of the stress maximum is nonsingular, see fig. 3(a). The stress scale delimiting the upper boundary of the shear-thickening region in fig. 3(b) is thus nearly constant.

When the jamming transition is approached, both Newtonian and shear-thickening regimes disappear. The former because it becomes easier to enter the shear-thinning regime as the viscosity $\eta_0$ increases [29]. The latter regime disappears when the viscosity in the Newtonian plateau becomes larger than the viscosity maximum at the end of the thickening regime, which happens at a density slightly below jamming. These two stress boundaries are reported in the diagram of fig. 3(b), which shows that the jammed system at $\varphi > \varphi_J$ either does not flow when σ is below the yield stress (which increases continuously with density above $\varphi_J$ [24]), or displays shear-thinning. The overall structure of the state diagram is similar to the experimental results [15].

We have shown that a simple Langevin model of soft repulsive particles displays a rheology in good agreement with the complex rheology observed in dense suspensions, due to the timescale competition among $\tau_0$, $\tau_T$ and $\tau_v$. While we described the rheology of soft particles, we note that both the thermal regime at low shear rate and the shear-thickening regime are not sensitive to particle softness, while the viscosity maximum and the transition to thinning rheology above the jamming density are specific to soft systems. Our simulations indicate in particular that the onset of shear-thickening occurs when $\dot{\gamma} \tau_v \approx 10^{-2}$, which agrees excellently with the value $\dot{\gamma} \tau_v \sim 0.7 \cdot 10^{-2}$ obtained in fig. 1(b) for the emulsion, while a somewhat smaller value $\dot{\gamma} \tau_v \sim 0.4 \cdot 10^{-3}$ is found for the latex dispersion shown in fig. 1(a). We note that the addition of frictional forces might affect the nature of the shear-thickening onset [19], but not necessarily the shear rate where it occurs [20].