Dynamic vorticity banding in discontinuously shear thickening suspensions
It has recently been argued that steady-state vorticity bands cannot arise in shear thickening suspensions, because the normal stress imbalance across the interface between the bands will set up particle migrations. In this Letter, we develop a simple continuum model that couples shear thickening to particle migration. We show by linear stability analysis that homogeneous flow is unstable towards vorticity banding, as expected, in the regime of negative constitutive slope. In full nonlinear computations, we show however that the resulting vorticity bands are unsteady, with spatiotemporal patterns governed by stress-concentration coupling. We furthermore show that these dynamical bands also arise in direct particle simulations, in good agreement with the continuum model.
Recent years have seen rapid advances in understanding the rheology of dense non-Brownian suspensions, comprising solid particles in a Newtonian fluid at volume fraction close to isotropic jamming. In particular, the phenomenon of shear thickening M&W (); Brown2014 (), in which the viscosity increases with stress , has recently been understood as an evolution from lubricated to frictional particle interactions, as the hydrodynamic forces that push particles together overcome short-ranged repulsive forces keeping them apart Lootens2005 (); Bashkirtseva2009 (); Brown2012 (); Fernandez2013 (); Seto2013 (); Heussinger2013 (); Mari2014 (); Mari2015 (); Lin2015a (); Guy2015 (); Comtet2017 (); Clavaud2017 (); Hsiao2017 (); 2018arXiv180103799H (). When strong, this effect creates, for states of homogeneous shear rate , a constitutive curve that is S-shaped Wyart2014 (); Mari2015a (): a positively sloping hydrodynamic branch of low viscosity at low stresses connects to a positively sloping frictional branch of high viscosity at high stresses via a negatively sloped region at intermediate stresses. A slowly increasing imposed shear rate then provokes a discontinuous jump, between the low and high viscosity branches, in the measured or ‘macroscopic’ flow curve. This is known as discontinuous shear thickening (DST) M&W (); Brown2014 ().
At imposed macroscopic shear stress, when one expects homogeneous steady flow to be unstable, at least for large system sizes Yerushalmi1970 (). (For the system sizes used in particle-based simulations, this expectation is not always met Mari2015a ().) Fluctuations with wave vector along the vorticity axis, , are then predicted to grow with time, leading to the formation of “vorticity bands” of differing local shear stress, coexisting at a common shear rate, with layer normals in Olmsted2008 (). In dense suspensions, however, a steady state coexistence of vorticity bands is argued to be ruled out by the differences in normal stress that generally arise across the interface between bands, leading to a particle migration flux Hermes2016 (). Suggestively, experiments have revealed an unsteady strain rate signal under conditions of a constant imposed macroscopic shear stress in the DST regime, with complicated time dependence Neuville2012 (); Hermes2016 (); Bossis2017 ().
In this Letter, we provide a fourfold advance in understanding dynamic vorticity banding in dense suspensions. First, we propose a scalar continuum constitutive model for the relevant rheology, by combining the Wyart-Cates theory Wyart2014 () (which captures shear thickening but assumes homogeneous flow) with a suspension balance model of particle migration Nott1994 (); Morris1999 (); Yurkovetsky2008 (); Lhuillier2009 (); Nott2011 (). Second, for this model we use linear stability analysis to determine when a homogeneous shear flow is unstable to fluctuations along the vorticity axis, finding instability whenever in the limit of large system size. Third, we elucidate numerically the model’s full nonlinear vorticity-banding dynamics, identifying two distinct spatio-temporal patterns that we shall term “travelling bands” (TB) and “locally oscillating bands” (LOB). The LOB state shows an oscillating bulk shear rate signal, as seen experimentally Hermes2016 (); Bossis2017 (). Finally, we perform particle-based simulations using the so-called Critical Load Model Mari2015a () and show that this also has TB and (at least transiently) LOB states, in close counterpart to the continuum model.
We consider Stokes flow of a dense suspension sheared between hard flat plates at under conditions of constant imposed shear stress. This produces a velocity field with a shear rate that is in general time-dependent. (Here denotes the suspension velocity averaged over the particle and solvent components introduced below.) The dynamical variables that we consider are the component of the particle phase stress tensor, the fraction of frictional contacts (the microstructural order parameter entering the Wyart-Cates theory Wyart2014 ()), and the volume fraction . Note that the component of stress is actually negative in dense suspensions M&W (), but we work throughout with its absolute value and denote this simply by . We assume spatial invariance in the flow direction and flow-gradient direction , considering spatial variations only in the vorticity direction . The shear rate however follows the relative speed of the plates, and so must remain uniform along .
The Wyart-Cates theory Wyart2014 () gives a scalar constitutive model for the steady state homogeneous shear rheology of dense suspensions. While initially presented as a model for the shear stress , this equally describes , because all stress components evolve in a similar way near jamming Boyer2011 () and across DST Mari2014 (). The associated viscosity is taken to diverge as the volume fraction approaches a critical jamming point :
where is, within the range of of interest here, of order the solvent viscosity and effectively constant, so hereafter we set for simplicity; is likewise a constant. Shear thickening is then captured by assuming that at low stresses repulsive forces maintain a lubrication film between particles, with a fraction of frictional contacts , whereas at high stresses frictional contacts dominate the rheology, . The critical volume fraction for jamming correspondingly also varies with stress, interpolating linearly with from (at low stress and ) the critical value for the jamming of frictionless particles to (at high stress and ) that for the jamming of frictional particles :
Here is the particle friction coefficient.
For the dependence of the fraction of frictional contacts on stress , particle simulations suggest a relation
in steady state, where . This depends on the typical repulsive force that must be overcome to create a contact and the typical particle radius ; particle simulations suggest Mari2014 (); Singh2017 (). Departing now from the steady-state assumptions of Wyart2014 (), we assume that does not react infinitely fast to changes in stress Mari2015a (); 2017arXiv171102196H (), following instead
Note that the evolution is on a characteristic strain scale . For the typical volume fractions considered here, particle simulations suggest Mari2015a ().
Next we assume that particle migration between vorticity bands is driven by the difference in normal stress that will in general exist across the interface between them. We model this via a two-fluid approach DoiBook (), known in this context as the suspension balance model Nott1994 (); Morris1999 (); Yurkovetsky2008 (); Lhuillier2009 (); Nott2011 (). The divergence of the particle stress gives a force imbalance on the particle phase, which must be rebalanced by a drag between the particles (p) and fluid (f) due to an interphase relative velocity , using the balance condition (that is, particles moving from regions of large stresses to regions of small stresses). Via mass conservation, , this gives
Particle simulations Hoef2005 () suggest the drag coefficient ranges from for to for , with the solvent viscosity. However, in this work, variations in will turn out to be , so we assume constant.
Eqs. 1–5 define our model. It contains the parameters and , along with the cell length in the vorticity direction , the average volume fraction , and the imposed average particle stress . Also relevant in our linear stability analysis will be the perturbation wave vector . We choose as length unit, as stress unit, and as time unit. Except when explicitly comparing with the simulation data, we also choose to rescale all strains by , so setting . We set rheological parameters compatible with the ones of spherical particles with moderate polydispersity, setting Boyer2011 (); Singh2017 (), (for friction Boyer2011 (); Mari2014 ()) and scott_packing_1960 (); Berryman1983 (). There then remain just three dimensionless parameters: rescaled drag (which for given bulk parameters and , is a measure of the system size ), rescaled stress and volume fraction . We drop tildes and denote these and hereafter. For any states that are imposed to be purely homogeneous, we denote the stress and particle fraction simply and .
The constitutive curves for homogeneous steady shear follow as the stationary solutions of Eqs. 1–4, and coincide directly with those of Wyart2014 (). They are shown as black lines in Fig. 1. At low volume fraction , they are monotonic. For they are S-shaped, with a regime in which , giving discontinuous shear thickening. At even larger , they bend right back to ascend the axis above a -dependent shear jamming stress , with flow only possible for stresses .
For any initial state on such a constitutive curve, the volume fraction and fraction of frictional contacts are defined to be uniform. We now perform a linear analysis to determine whether any such “base state” is stable, by adding to it small amplitude perturbations , . Expanding Eqs. 1–5 to first order in the amplitudes , we find instability (that is, positive real part for the eigenvalue of the resulting linearized dynamics), whenever,
In regimes of high or low , we also find a non-zero imaginary part , giving an oscillatory component to the growing perturbations. For an infinitely large system, that is, , this yields the familiar mechanical instability criterion . When the system size is finite, the unstable region shrinks to an extent prescribed by , with stability boundaries shown as colored lines in Fig. 1.
The mechanism of this instability is as follows. Temporarily ignoring variations in , Eqs. 1–4 effectively reduce to (a) and (b) . Recalling that must remain uniform in while and can vary, we imagine a localised fluctuation in which slightly increases at some : (b) then requires that correspondingly also increases. For large enough , must increase even further to maintain uniform along via (a). This gives positive feedback and instability, irrespective of wave vector . We now relax the assumption of constant , noting that the relaxation depends on , and is always stabilizing: whenever increases locally, there is an outward migration of particles, thereby locally decreasing and hence the viscosity. The competition between these processes restricts the unstable dynamics to large length scales, because the time for particle migration increases with distance.
Having shown a state of initially homogeneous flow to be linearly unstable if Eq. (6) is satisfied, we now numerically integrate the model equations to elucidate the full nonlinear dynamics that prevails at long times. We do so for a representative value of and for two volume fractions: , for which the constitutive curve is S-shaped; and , for which it folds right back to the axis. We take as an initial condition a state of homogeneous shear on the constitutive curve for some imposed , subject to small amplitude perturbations. We then quasi-statically sweep up and (in a separate run) down from this value. By this procedure, we find two different long-time states of dynamical vorticity bands, with associated macroscopic flow curves shown separately by black and red lines in Fig. 2.
The black flow curves pertain to states with a steady state bulk shear rate at any given stress value. Within the sample, these exhibit travelling bands (TB), shown in Fig. 3, top left: localised pulses resembling solitons travel along the vorticity axis at constant speed in one direction. (The direction represents a spontaneously broken symmetry, depending sensitively on initial conditions.) The red flow curves in Fig. 2 pertain to states in which the bulk shear rate shows sustained oscillations in time, at any given imposed stress, reminiscent of experimental observations Hermes2016 (); Bossis2017 (). The average of the oscillation is shown by the solid red line, and the limits by the dotted red lines. These lie within the region where . Within the sample, these states exhibit locally oscillating bands (LOB) comprising two soliton-like excitations that travel in opposite directions and intermittently collide. The time of collision corresponds with the falling phase in the oscillating bulk shear rate signal. The model also predicts TB to move at a speed a little larger than observed in the simulations (by about a factor of two). For the higher volume fraction, homogeneous flow is recovered above some stress ; flow then arrests completely at , where the underlying constitutive curve re-joins the vertical axis. We note that the composite (macroscopic) flow curves are tilted – as opposed to vertical, as would have been expected without concentration coupling Olmsted2008 () – and that negatively sloping macroscopic flow curves were seen experimentally in Pan2015 ().
We now compare these continuum results to direct particle simulations. These use 8000 bidisperse spheres (radii and in equal volume proportions), sheared in a tri-periodic box of size and , using Lees-Edwards boundary conditions Lees1972 (). The spheres interact through both lubrication and contact forces Mari2014 (), with any contact force becoming frictional (with friction coefficient ) once the normal part exceeds a fixed value . This “Critical Load Model” captures both DST and jamming, but previous work Mari2014 () used much smaller values, precluding the vorticity instabilities addressed here. With , the values of , and are consistent with the ones we set for the continuum model, giving a good agreement for the homogeneous steady state flow curves Singh2017 ().
In our large- simulations, we indeed find dynamic vorticity banding, with two distinct dynamical states, in close analogy with those of our continuum model. As seen in the top right of Fig. 3, we recover the TB solutions at lower stresses. We find good qualitative agreement between simulations and our continuum model (top left of Fig. 3) in the overall dynamics of stress and volume fraction fields. To perform the comparison, we used model parameters and , based on their measured value in the simulations ( from Mari2015a () and from fits to the measured relation between the particle phase velocity and stress gradient ). Both continuum model and particle simulations show slight accumulation of particles at the front of the travelling thickened band, albeit with somewhat flatter stress profiles in the simulations than in the model (not shown). However, the model predicts TB to move at a speed roughly twice larger than observed in the simulations. At higher stresses, our simulations also exhibit the LOB states of the continuum model, but we have so far only found these as a transient effect in the particle simulations, with LOB always eventually giving way to TB: see the lower panel of Fig. 3.
In summary, we have proposed a continuum model for the vorticity instabilities of a shear-thickening suspension, held at a constant macroscopic stress in the unstable part of the constitutive curve where . Its predictions compare very well with our particle based simulations, including a regime of oscillating macroscopic shear rates as found experimentally Hermes2016 (); Bossis2017 (). Crucially, the unsteady behavior results from a bulk rheological mechanism, not from coupling with the mechanical response of the rheometer, even if the latter plays a part in some experiments Bossis2017 (). Particle migration is crucial: the banding dynamics relies on small concentration variations that have large rheological effects close to jamming, as reported previously for colloidal glasses in pipe flow Besseling2010a ().
Observing the predicted spatiotemporal bands directly in experiments, by measuring local stress fields and small concentration fluctuations, may prove challenging. The velocity field along the vorticity direction also bears a signature of the bands due to particle migration, which could be more accessible. Very recent experiments do report vorticity bands in cornstarch suspensions under controlled stress, similar to those presented here in shape, size and velocity Saint-Michel2018 (). However, the banding signature involves the flow () velocity component; this may stem from differential wall slip induced by a frictional band moving along the vorticity direction.
While we focused here on a constitutive model involving only the normal stress along the vorticity direction, the physical ingredients in our model likely also admit instabilities along the gradient and/or flow directions. Without volume fraction variations, homogeneous shear flow is stable against gradient perturbations in absence of inertia Mari2015a (), for constitutive curves of the shapes considered here, so that such instability again requires particle migration (and/or inertia). Indeed some experiments in the discontinuous shear thickening regime report a flowing gradient band of lower concentration coexisting with a densely jammed band Fall2015 (). For very large systems, inertia will separately trigger gradient instabilities Bashkirtseva2009 (); Nakanishi2012 (). Our coupled model of shear thickening and particle migration represents a first step towards explaining the full range of unsteady flows close to the jamming transition in dense suspensions.
We thank M. Hermes, I. Peters, B. Saint-Michel, T. Gibaud and S. Manneville for insightful discussions. Work funded in part by SOFI CDT, Durham University, and EPSRC (EP/L015536/1); the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013) and Horizon 2020 Programme / ERC grant agreements 279365 and 740269. MEC is funded by the Royal Society.
- preprint: APS/123-QED
- J. Mewis and N. J. Wagner, Colloidal Suspension Rheology (Cambridge University Press, Cambridge, England, 2011).
- E. Brown and H. M. Jaeger, Rep. Prog. Phys. 77, 046602 (2014).
- D. Lootens, H. Van Damme, Y. Hémar, and P. Hébraud, Phys. Rev. Lett. 95, 268302 (2005).
- I. A. Bashkirtseva, A. Y. Zubarev, L. Y. Iskavova, and L. B. Ryashko, Colloid J. 71, 446 (2009).
- E. Brown and H. M. Jaeger, J. Rheol. 56, 875 (2012).
- N. Fernandez, R. Mani, D. Rinaldi, D. Kadau, M. Mosquet, H. Lombois-Burger, J. Cayer-Barrioz, H. J. Hermann, N. D. Spencer, and L. Isa, Phys. Rev. Let. 111, 108301 (2013).
- R. Seto, R. Mari, J. F. Morris, and M. M. Denn, Phys. Rev. Let. 111, 218301 (2013).
- C. Heussinger, Phys. Rev. E 88, 050201 (2013).
- R. Mari, R. Seto, J. F. Morris, and M. M. Denn, J. Rheol. 58, 1693 (2014).
- R. Mari, R. Seto, J. F. Morris, and M. M. Denn, Proc. Natl. Acad. Sci. U.S.A. 112, 15326 (2015).
- N. Y. C. Lin, B. M. Guy, M. Hermes, C. Ness, J. Sun, W. C. K. Poon, and I. Cohen, Phys. Rev. Lett. 115, 228304 (2015).
- B. M. Guy, M. Hermes, and W. C. K. Poon, Phys. Rev. Lett. 115, 088304 (2015).
- J. Comtet, G. Chatté, A. Niguès, L. Bocquet, A. Siria, and A. Colin, Nat. Commun. 8, 15633 (2017).
- C. Clavaud, A. Bérut, B. Metzger, and Y. Forterre, Proc. Natl. Acad. Sci. U.S.A. 114, 5147 (2017).
- L. C. Hsiao, S. Jamali, E. Glynos, P. F. Green, R. G. Larson, and M. J. Solomon, Phys. Rev. Lett. 119, 158001 (2017).
- C.-P. Hsu, S. N. Ramakrishna, M. Zanini, N. D. Spencer, and L. Isa, arXiv:1801.03799.
- M. Wyart and M. E. Cates, Phys. Rev. Lett. 112, 098302 (2014).
- R. Mari, R. Seto, J. F. Morris, and M. M. Denn, Phys. Rev. E 91, 052302 (2015).
- J. Yerushalmi, S. Katz, and R. Shinnar, Chem. Eng. Sci. 25, 1891 (1970).
- P. D. Olmsted, Rheol. Acta 47, 283 (2008).
- M. Hermes, B. M. Guy, W. C. K. Poon, G. Poy, M. E. Cates, and M. Wyart, J. Rheol. 60, 905 (2016).
- M. Neuville, G. Bossis, J. Persello, O. Volkova, P. Boustingorry, and M. Mosquet, J. Rheol. 56, 435 (2012).
- G. Bossis, P. Boustingorry, Y. Grasselli, A. Meunier, R. Morini, A. Zubarev, and O. Volkova, Rheol. Acta 56, 415 (2017).
- P. R. Nott and J. F. Brady, J. Fluid Mech. 275, 157 (1994).
- J. F. Morris and F. Boulay, J. Rheol. 43, 1213 (1999).
- Y. Yurkovetsky and J. F. Morris, J. Rheol. 52, 141 (2008).
- D. Lhuillier, Phys. Fluids 21, 023302 (2009).
- P. R. Nott, É. Guazzelli, and O. Pouliquen, Phys. Fluids 23, 043304 (2011).
- F. Boyer, É. Guazzelli, and O. Pouliquen, Phys. Rev. Lett. 107, 188301 (2011).
- A. Singh, R. Mari, M. M. Denn, and J. F. Morris, J. Rheol. 62, 457 (2018).
- E. Han, M. Wyart, I. R. Peters, and H. M. Jaeger, arXiv:1711.02196.
- M. Doi, Soft Matter Physics (Oxford University Press, Oxford, England, 2013).
- M. A. van der Hoef, R. Beetstra, and J. A. M. Kuipers, J. Fluid Mech. 528, 233 (2005).
- G. D. Scott, Nature 188, 908 (1960).
- J. G. Berryman, Phys. Rev. A 27, 1053 (1983).
- Z. Pan, H. de Cagny, B. Weber, and D. Bonn, Phys. Rev. E 92, 032202 (2015).
- A. W. Lees and S. F. Edwards, J. Phys. C 5, 1921 (1972).
- R. Besseling, L. Isa, P. Ballesta, G. Petekidis, M. E. Cates, and W. C. K. Poon, Phys. Rev. Lett. 105, 268301 (2010).
- B. Saint-Michel, T. Gibaud, and S. Manneville, Submitted to Phys. Rev. X (2018).
- A. Fall, F. Bertrand, D. Hautemayou, C. Mezière, P. Moucheront, A. Lemaître, and G. Ovarlez, Phys. Rev. Lett. 114, 098301 (2015).
- H. Nakanishi, S.-I. Nagahiro, and N. Mitarai, Phys. Rev. E 85, 011401 (2012).