## Abstract

The viscosity of suspensions of large (≥10 μm) particles diverges at high solid fractions due to proliferation of frictional particle contacts. Reducing friction, to allow or improve flowability, is usually achieved by tuning the composition, either by changing particle sizes and shapes or by adding lubricating molecules. We present numerical simulations that demonstrate a complementary approach whereby the viscosity divergence is shifted by driven flow tuning, using superimposed shear oscillations in various configurations to facilitate a primary flow. The oscillations drive the suspension toward an out-of-equilibrium, absorbing state phase transition, where frictional particle contacts that dominate the viscosity are reduced in a self-organizing manner. The method can allow otherwise jammed states to flow; even for unjammed states, it can substantially decrease the energy dissipated per unit strain. This creates a practicable route to flow enhancement across a broad range of suspensions where compositional tuning is undesirable or problematic.

## INTRODUCTION

Densely packed suspensions arise widely in industry and manufacturing, where reliable, predictable, and prescribable flow properties are essential (*1*). A major limiting factor in their processability is the very steep increase of viscosity upon increasing the volume fraction of solid material φ toward the jamming transition (*2*). This is particularly evident in the non-Brownian regime (particle size, ≥10 μm), where frictional particle contact interactions reduce the jamming density φ_{m} (*3*–*5*) and increase dissipation in process flows, resulting in high energy costs.

Empirical strategies that reduce the viscosity and/or net dissipation include tuning the physical properties of the particles—for example, their size, shape, and polydispersity—or modifying their interactions through chemical additives known as plasticizers, emulsifiers, or friction modifiers. These lubricate interparticle contacts and reduce the suspension viscosity by raising φ_{m}. Often, though, end-use requirements leave little room for maneuver in the formulation. In calcium phosphate cements for bone injection (*6*, *7*), for example, chemistry and biology both constrain the use of molecular additives. There is, therefore, a practical need for methods of dense suspension flow control that do not require changes to formulation.

Two recent experiments suggest a possible route toward this goal, achieving driven viscosity reduction by superimposing an oscillatory cross shear (OCS) on a primary desired flow. Using OCS, Blanc *et al*. (*8*) demonstrated a twofold increase in the sedimentation velocity of an intruder in a granular suspension of rate-independent rheology, whereas Lin *et al*. (*9*) measured a two-decade viscosity drop in one of shear-thickening rheology. The latter effect was argued to be a consequence of the fragility of shear-induced particle contacts (*9*, *10*), suggesting that good flowability might be achieved only when the OCS is sufficiently fast to keep the microstructure in a load-incompatible state (*10*, *11*). In this limit, the reduction in primary flow viscosity (unless this is infinite) might easily be outweighed by the high energy cost of implementing fast OCS, particularly for rate-independent suspensions (*8*) whose primary viscosity drop is much less than in shear-thickening ones (*9*). More generally, it is not clear how far the benefits of OCS depend on the underlying suspension rheology: The short-range repulsions that prevent frictional particle contacts at low stresses in thickening suspensions (*12*–*14*) may or may not play a major role during OCS-assisted viscosity reduction.

Here, we present numerical simulations showing that the viscosity drop induced by OCS is generic to suspensions with friction-dominated stress. This includes noninertial flows of most dense suspensions of super-micrometer–sized particles (*5*). The transverse flow oscillations directly inhibit particle contacts without requiring short-range repulsions, enhancing lubrication and shifting φ_{m} to higher values. Consequently, the viscosity reduction increases with increasing φ, so that near jamming the saving in primary flow dissipation outweighs the cost of OCS at any primary flow rate, giving a net reduction in the energy expended per unit strain in the primary direction. We then show that the reduction in particle contacts stems from an OCS-induced “random organization” mechanism (*15*–*17*); this leads us to an enhanced version of the flow protocol that can reduce the dissipation further. Guided by these results, we argue that driven viscosity control should extend flowability and reduce the associated energy cost across a broad class of materials, including slurries, muds, cement, and other immersed granular systems.

## RESULTS AND DISCUSSION

We study a suspension of nearly hard, athermal spheres subject to short-range hydrodynamic and contact interactions with static friction coefficient μ_{s} as described in Materials and Methods below. This numerical model [and similar ones (*18*, *19*)] is known to yield accurate predictions for the rheology of non-Brownian hard sphere suspensions. The suspension shows rate-independent rheology, well described under steady simple shear by the viscous number formalism [see the study by Boyer *et al*. (*20*) and the Supplementary Materials]. A snapshot of the simulated system is shown in Fig. 1A (i).

### Manipulating suspension viscosity using superimposed oscillations

From the argument that fragility makes contact stresses in suspensions susceptible to driven perturbations (*9*), it follows that the addition of any arbitrary oscillating flow might lead to viscosity reduction. This hypothesis is in line with experimental (*21*, *22*) and theoretical (*23*) works that propose applied and endogenous noise, respectively, as sources of opening and closing granular contacts and consequent unjamming. To test this, we first explore a generalization of OCS comprising primary steady shear with rate and superimposed oscillatory shears in both the primary and cross shear directions, leading to an overall strain in *xy* as and in *zy* as γ^{OCS}(*t*) = γ sin (ω*t*). For simplicity, we keep γ = 1% in each case, which Lin *et al*. (*9*) found to be an optimal amplitude for viscosity reduction. The remaining dimensionless control parameters are then , , and the phase shift δ. This protocol gives strain paths such as those illustrated in Fig. 1A (ii and iii). A characteristic viscosity is computed as averaged over time units, where η is the solvent viscosity and σ_{xy} is the *xy* component of the stress. We find that δ has very little effect on the viscosity (see the Supplementary Materials), and present a contour map of η_{r} in the (, ) plane at δ = 0 and φ = 0.55 in Fig. 1B. At fixed δ, viscosity minima are obtained as and . In this limit, that is, with cross shear oscillations only, we obtain a viscosity drop comparable to that for a thickening suspension (*9*). This suggests that OCS—by keeping frictional particle contacts open—effectively brings the suspension to a low-friction state. We similarly find a maximal rate of viscosity reduction when is close to unity. Contrary to the hypothesis made above, however, our results show that the orientation of the oscillatory flow is crucial: At this strain amplitude, any oscillatory component along the primary flow direction makes no useful contribution to improving flowability. If the shear is constrained to a single direction and the amplitude of the oscillations is very small compared to the primary flow, then the net displacements of the particles over large strains are, for rate-independent flow, the same as for steady shear. This is not the case when the oscillations are applied transverse to the primary flow.

### Simple OCS: Viscosity reduction using transverse oscillations

In what follows, we therefore revert to the purely transverse case with ω^{pri} = 0, leading to , and γ^{OCS}(*t*) = γ sin (ωt) (see Fig. 2A, inset), hereafter called the “simple OCS” (SO) protocol. (This is to distinguish it from an alternative protocol introduced below.) In Fig. 1C, we report the viscosity η_{r} under this protocol at γ = 1% as a function of the reduced frequency , and in Fig. 1D we compare, as a function of volume fraction φ, the steady shear viscosity (obtained when ) to the limiting viscosity under SO (obtained when ). The viscosity drop increases rapidly with φ, reaching a decade at φ = 0.56 and actually diverging between φ = 0. 57 and 0. 58 (Fig. 1D, inset). This reveals that as well as reducing the viscosity, the effect of SO is to slightly delay the jamming transition from φ_{m} ≈ 0.58, for steady shear to φ_{m,SO} ≈ 0.60, at . Although small in absolute terms, shifts of jamming by a couple of percent can have dramatic consequences for formulation and processing (*24*), as discussed further below.

This shift of jamming under SO naturally raises the question of the sensitivity of the viscosity reduction to particle friction, which may stem from, for example, surface roughness (*25*). It is expected that in the absence of ordering, which we do not observe in our binary system, the random close packing density φ_{RCP} ≈ 0.64 is an upper limit for both φ_{m} and φ_{m,SO}. Furthermore, it is established that the jamming point φ_{m} approaches φ_{RCP} as surface friction μ_{s} is decreased (*26*, *27*). Because φ_{m} < φ_{m,SO}, it follows that φ_{m} < φ_{m,SO} < φ_{RCP} and the window between φ_{m} and φ_{m,SO} consequently vanishes in the limit of low friction (as μ_{s} → 0). The performance of SO, thus, diminishes as friction decreases. We demonstrate this in Fig. 1E in the limits of steady shear and SO with . As a result, suspensions of rough particles, which are typically the most problematic in terms of processing (*18*), are best placed to benefit from driven flow control. In the context of friction-driven shear thickening of colloids, this result thus confirms that SO can be successful in reducing the viscosity of a thickened sample, as demonstrated by Lin *et al*. (*9*), but that it would fail to reduce the viscosity of a nonthickened sample, that is, one at which the applied stress lies below the onset stress (*5*, *12*, *28*, *29*).

### SO-enabled reductions in energy dissipation

The ability to control suspension viscosity during flow is itself desirable for mitigating instabilities (*30*) and, for example, when pumps are desired to operate within narrow bounds. Often, though, rheological tuning has a somewhat different objective: to minimize the energy cost of processing. For 0.58 < φ < 0.60, OCS triumphs: It permits flow at finite dissipation rates not otherwise possible. Below jamming (φ < φ_{m}) however, its benefits are less obvious. The energy dissipation is given per unit volume and per unit primary strain as . Figure 1F (i to iii) shows this quantity (rescaled by the φ-dependent steady shear dissipation *W*_{SS}) as a function of oscillation rate for our SO protocol, separating out the primary () and cross flow () contributions. The primary dissipation decreases in line with the viscosity, but the direct cost of the cross shear increases as . Summing these, we identify oscillation rates for which *W* is usefully decreased, highlighted green in Fig. 1F (i and ii) and outlined in white in Fig. 1F (iii). This operating window, although it grows as the density approaches φ_{m}, remains narrow at lower densities. We show below that it can be extended significantly by a simple modification to the oscillatory protocol.

### Random organization drives the viscosity reduction

The proposed modification exploits mechanistic insights, gleaned from our simulations, into how OCS promotes flowability. To gain these insights, we start by decomposing the viscosity into its hydrodynamic and frictional particle contact contributions, revealing that at φ = 0.54, the stress is dominated by friction for any oscillatory frequency (Fig. 2A). Significantly, the effect of the cross shear oscillations is to decrease this frictional part while leaving the hydrodynamic part unchanged. The loss of friction parallels the shift of jamming to higher φ (Fig. 1D), indicative of a shift from rolling to sliding contacts as is increased (*5*, *29*). Defining interparticle gaps *h*_{ij} = 2(*r*_{ij} − *a*_{i} − *a*_{j})/(*a*_{i} + *a*_{j}) with *r*_{ij} as the center-to-center distance between particles *i* and *j* with radii *a*_{i} and *a*_{j}, respectively, we compute *G*(*h*), the average number of neighbors around a particle separated at most by *h* (Fig. 2B). The loss of frictional particle contacts occurs by a “room-making” process, whereby the mean distance between nearest neighbors increases. Consequently, starting with a presheared sample, there is a gradual decrease of particle contacts over cycles after SO startup (Fig. 2B, inset). Strikingly, room-making does just enough to hinder the stress-generating contacts. This is strongly reminiscent of “random organization,” whereby application of oscillatory shear to a suspension of hard particles drives collective self-organization, leading to configurations that minimize the number of particle contacts generated per cycle (*15*–*17*). Below a critical volume fraction, the system evolves to an “absorbing state” for which configuration invariance is ensured under further oscillations. Although first elaborated for dilute suspensions, a similar scenario applies at higher density where the absorbing state is defined not by absence of collisions but by absence of plastic rearrangements (*31*, *32*).

Flow-induced random organization offers a natural explanation for the viscosity decrease upon increasing oscillation frequency. At high frequency, the “primary” and “secondary” labels respectively assigned to steady shear and OCS are misnomers. In fact, we have a steady transverse flow that weakly perturbs an oscillatory flow, for which the random organization effect is well established. This makes room around particles, thus decreasing the contact stress, whereas the steady shear slowly consumes this room and simultaneously initiates new particle contacts. At finite primary flow, the absorbing state can never be reached, but its proximity allows particles to avoid frictional particle contacts at densities where these would otherwise cause large viscosities or jamming.

The largest φ at which absorbing states are obtained locates a nonequilibrium phase transition, where the self-organization process is maximized (*15*, *16*, *33*). To confirm the role of random organization, we determine the location of this transition in our system. Starting from a presheared configuration, we apply an oscillatory shear at γ = 1% with no primary flow and measure the fraction of particles following irreversible, “active” trajectories, which, for these purposes, we define as those whose net displacement after a cycle stays below a threshold of 10^{−5}*a*_{i}. For φ ≤ 0.58, this quantity approaches zero, indicating that the system evolves toward absorbing states (Fig. 2C). For an amplitude γ = 1%, the absorbing state transition density is thus estimated as 0.58 < φ < 0.59. Because in practice the transition is cut off by the primary flow, our precise definition of activity is not crucial here, although a more inclusive one (for example, counting all particles that make frictional contact at any point during the cycle) would give a lower estimate for the transition. Nonetheless, our results indicate that the random organization effect is a strong one throughout the density range where our SO protocol is effective.

### Alternating OCS: Separated flow phases reduce dissipation further

In the SO protocol above, particle contacts are eliminated by applying OCS concurrently with the desired primary shear. A relatively high energy cost arises from the need to have sufficiently fast oscillations to ensure that random organization can compete with the restoration of frictional particle contacts caused by the primary shear. If this is indeed the mechanism, though, there is no strict requirement that we perform these flows concurrently. Instead, we can use alternating intervals of OCS without primary shear and of primary shear without OCS. The former eliminates particle contacts; the latter restores them, but not before a finite strain has been achieved. The cross shear dissipation can, in principle, be reduced to zero by having long intervals of very slow oscillations, creating an optimization scenario different to that of SO.

We therefore now test a new flow protocol [“alternating OCS” (AO)] that alternates an interval of *n* periods of oscillation during a time α*T* with and γ^{OCS}(*t*) = γsin (ω*t*) for ω = 2π*n*/(α*T*), with an interval of primary shear during a time (1 − α)*T* with and . Guided by our result in Fig. 1B, the oscillations are applied transverse to the primary flow, where we anticipate that their efficacy will be maximized. Averaged over one cycle *T*, the shear rate in the primary direction is . The primary shear strain during each cycle is Γ, that is, .

In the case of rate-independent dynamics as simulated above, the microstructure depends on the strain path only (sketched in Fig. 3A, inset), not the rate at which it is followed. As a consequence, the viscosity depends on *n*, γ, and Γ, but not on α (as long as 0 < α < 1). We define the relative viscosity as . This viscosity is averaged over the intervals of pure primary shear during the AO protocol, measured over a time period covering 30 strain units in the primary direction. It is reported as a function of *n* in Fig. 3A for φ = 0.56 and γ = Γ = 1%. The viscosity drops rapidly with *n*, and, remarkably, even *n* = 1 is sufficient to achieve a viscosity reduction of ≈ 96% in this unjammed system. This viscosity reduction is already larger than that achieved with SO, as seen by comparing the relative reductions in Fig. 3B (AO) and the inset of Fig. 1D (SO). Finally, we find that the viscosity drop is maximized with AO [as with SO (*9*)] for γ, Γ ≈ 1 to 5%, whereas for larger γ, interparticle gaps close during cycles, allowing frictional particle contacts and a rapid viscosity increase. A viscosity transient for the AO protocol is given in the Supplementary Materials.

Although α has no role in setting the viscosity, it is a crucial parameter when it comes to the dissipation, which depends on the deformation rate. In Fig. 3C, we show the work per unit primary strain *W* (rescaled by the φ-dependent steady shear dissipation *W*_{SS}) as a function of α for several values of *n*, for φ = 0.56. For *n* = 1 and for α values between 0.58 and 0.9, the overall dissipation is reduced compared to steady shear, reaching a reduction of around 45% at α ≈ 0.8. Dissipation is minimized when *n* = 1 for all φ, as the shear rate shoots up quickly with increasing *n*, swamping any further viscosity reduction achieved for *n* > 1. To ease the comparison with the SO protocol, we define an oscillatory frequency for AO as ω = 2π*n*γ/Γ (that is, a given frequency corresponds to the same number of cross shear oscillations per unit strain in the primary direction for SO and AO). In Fig. 3D, we show a map of the relative dissipation *W*/*W*_{SS} in the (φ, ω) plane, showing a wider area of reduced dissipation with AO compared to SO [Fig. 1F (iii)]. The AO protocol therefore has clear advantages particularly in avoiding the need to precisely tune the frequency of the driving oscillations. We finally present in Fig. 3E a comparison of the reduction in energy dissipation achieved by the AO and SO protocols as a function of the volume fraction φ. For φ ≲ 0.54, there is no further gain with AO compared to SO, but for larger volume fractions, where the viscosity reduction performance of AO is markedly superior (Fig. 3B), there is an improvement in AO over SO. We quantify this improvement in Fig. 3E (inset), giving the ratio of the minimal dissipation for AO and SO (). This shows that close to the steady shear jamming volume fraction, the improvement of AO over SO can reach almost 40%, suggesting that AO is likely the protocol of choice for dissipation reduction in very concentrated frictional suspensions.

### Concluding remarks

Our results show that nonsteady deformation protocols can lead to substantial viscosity and energy dissipation reductions in any friction-dominated suspension flow. The strategy is applicable for most flows involving granular suspensions and related systems in which frictional particle contacts bear most of the stress in steady shear, including Brownian suspensions under very large stresses (*34*). Because of their simplicity, our protocols, or ones like them, might be readily implemented as precision unblockers and flow controllers in industrial devices, such as extruders or mixers, or as dissipation regulators for active granular damping (*35*). In particular, an extrusion nozzle might be fitted with an internal coaxial cylindrical actuator that oscillates about its axis with a protocol specified to maximize flowability according to our present results. Moreover, these implementations might be applied not only to minimize viscosities but also to regulate them against a desired set point. We tested such a protocol numerically with good success (see the Supplementary Materials). From a fundamental point of view, the relation to random organization opens new research directions. For example, it suggests that protocols other than oscillatory flow that lead to a similar absorbing phase transition (*36*) might also be good candidates for driven flow enhancement in complex fluids. It also suggests an unexpected link between rheological properties and hyperuniformity (*37*, *38*).

## MATERIALS AND METHODS

We simulated the trajectories of athermal, noninertial particles using a minimal model that comprises short-ranged hydrodynamic lubrication and frictional surface contacts. Our simulations comprise particles with size ratio 1:1.4 in a periodic box. For a particle pair with positions **x**_{1} and **x**_{2} and translational and rotational velocities **U**_{1}, **U**_{2} and **Ω**_{1}, **Ω**_{2}, respectively, in a background flow described at **x**_{1} by **U**^{∞}(**x**_{1}) = **E**^{∞}**x**_{1} + **Ω**^{∞} × **x**_{1}, the hydrodynamic forces , and torques , are given by (*39*–*41*)(1)The matrices **R**_{Lub} and **R**_{Stokes} follow our earlier description (*28*), with the scalar resistances therein comprising only the leading short-ranged diverging contributions, following the study by Ball and Melrose (*42*). The hydrodynamic stress contribution for particle 1 resulting from its pairwise interaction with particle 2, with force and particle-particle vector **r** is given by .

The leading terms of **R**_{Lub} diverge according to 1/*h* as particles 1 and 2 approach, with *h* as the surface-surface distance. Following experimental evidence that lubrication layers break down in suspensions under large stress (*43*), and, equivalently, for large particles (*5*), we used a minimum *h*_{min} = 0.001*a* (with *a* the smaller particle radius), below which hydrodynamic forces are regularized and particles may come into contact. For a particle pair with contact overlap δ and center-center unit vector **n**, we computed the contact force and torque according to (*44*)(2a)(2b)
where **u** represents the incremental tangential displacement, reset at the initiation of each contact. *k*_{n} and *k*_{t} are stiffnesses and *a*_{1} is the radius of particle 1. The tangential force component is restricted by a Coulomb friction coefficient μ_{s}, such that |*k*_{t}**u**| ≤ μ_{s}*k*_{n}δ. For larger values of |*k*_{t}**u**|, contacts enter a sliding regime. The contact stress contribution is given by for particle-particle vector **r** and pairwise force . The stress tensor is , where η is the suspending fluid viscosity and the sums are over all relevant pairwise interactions. Throughout the main text, we focused on the shear component in the primary flow direction, σ_{xy}.

Trajectories were computed from the above forces using two equivalent schemes. In the first, contact and hydrodynamic forces and torques were summed on each particle [according to the study by Radhakrishnan (*45*)] and the trajectory was updated according to Newtonian dynamics [using LAMMPS (*46*)], ensuring the Stokes number ( for particle density ρ, suspending fluid viscosity η, and shear rate ) remained ≪ 1 to approximate overdamped conditions. We also set to approximate hard spheres. In the second, per-particle forces were explicitly set to zero and the velocities were computed to balance contact and hydrodynamic forces and torques, ensuring strictly inertia-free flow (*12*, *28*). The numerical model generates results consistent with μ(*J*)-rheology as predicted by the experimental work of Boyer *et al*. (*20*) (see the Supplementary Materials).

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/4/3/eaar3296/DC1

section S1. Demonstration that our numerical model predicts μ(*J*) rheology

section S2. Influence of primary-cross shear phase difference

section S3. Time evolution of the relative viscosity for the SO and AO protocols

section S4. Setting the viscosity using a proportional flow controller

fig. S1. Comparison of the results from the present model to viscous number rheology.

fig. S2. Viscosity maps analogous to Fig. 1B but with varying phase shift δ.

fig. S3. Time evolution of the viscosity for the SO and AO protocols.

fig. S4. Time evolution of the relative viscosity under proportional control in the SO protocol.

This is an open-access article distributed under the terms of the Creative Commons Attribution license, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**C.N. thanks R. Radhakrishnan for sharing his derivation of the hydrodynamic forces.

**Funding:**C.N. acknowledges an Engineering and Physical Sciences Research Council Impact Acceleration Account and grant EP/N025318/1 and, subsequently, the Maudslay-Butler Research Fellowship at Pembroke College, Cambridge, for financial support. M.E.C. is funded by the Royal Society.

**Author contributions:**All authors planned the work, interpreted the results, and wrote the manuscript. C.N. and R.M. conducted the simulations.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**Codes required to generate the results in this article are given at https://doi.org/10.17863/CAM.13416 and https://doi.org/10.5281/zenodo.999197). All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors

- Copyright © 2018 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution License 4.0 (CC BY).