Unsteady flow, clusters and bands in a model shear-thickening fluid

TL;DR

This paper investigates the complex unsteady flow behaviors, including clusters and bands, in a 2D model of shear-thickening granular fluids near jamming, revealing large-scale patterns and their dynamics.

Contribution

It uncovers the formation and dynamics of large-scale flow patterns like bands and clusters in shear-thickening granular fluids, highlighting their dependence on strain rate.

Findings

Large-scale patterns depend on strain rate.

Force chains form giant bands at lower stress.

Clusters extend along dilational direction at higher stress.

Abstract

We analyse the flow curves of a two-dimensional assembly of granular particles which are interacting via frictional contact forces. For packing fractions slightly below jamming, the fluid undergoes a large scale instability, implying a range of stress and strainrates where no stationary flow can exist. Whereas small systems were shown previously to exhibit hysteretic jumps between the low and high stress branches, large systems exhibit continuous shear thickening arising from averaging unsteady, spatially heterogeneous flows. The observed large scale patterns as well as their dynamics are found to depend on strainrate: At the lower end of the unstable region, force chains merge to form giant bands that span the system in compressional direction and propagate in dilational direction. At the upper end, we observe large scale clusters which extend along the dilational direction and propagate along the compressional direction. Both patterns, bands and clusters, come in with infinite correlation length similar to the sudden onset of system-spanning plugs in impact experiments.

Equations

$${\mathbf{f}}_{ij}={\mathbf{f}}_{ij}^{(n)}+{\mathbf{f}}_{ij}^{(t)}.$$

$${\mathbf{f}}_{ij}^{(n)}$$

$${\mathbf{f}}_{ij}^{(t)}$$

Full text

I Introduction

The defining feature of non-Newtonian fluids is that the viscosity, i.e. the resistance of the fluid to flow, is not a material constant but depends on the flow itself. A shear thickening fluid, in particular, has the property that the viscosity increases with the speed of the flow. Shear thickening (ST) can be modest, with only a small increase of the viscosity above its zero-strainrate value. It may also be an order of magnitude effect, and even lead to a discontinuous flow arrest upon increasing the driving force beyond a threshold. Such a spectacular phenomenon is quite opposite to “normal” materials that start to flow or break upon increasing the force.

A related phenomenon is the establishment of system-spanning solid structures upon impact of an object on the surface Maharjan et al. (2018); Allen et al. (2018). These structures may even be strong enough to carry persons, at the same time absorbing enough kinetic energy to be used as shock-absorbers Decker et al. (2007).

It has been shown by experiments Brown and Jaeger (2012, ); Clavaud et al. (2017); Comtet et al. (2017); Pan et al. (2015); Hsu et al. (2018) and simulations Heussinger (2013); Seto et al. (2013); Ness and Sun (2016) that solid-solid friction between particles is the relevant force for the ST effect. The onset for ST is governed by an intrinsic force scale, sometimes modeled via a switch, at which frictional forces start to be relevant Seto et al. (2013); Wyart and Cates (2014). The upper limit for ST is equally set by a force-scale, that represents the weakest link in the system, e.g. surface tension between the sample and surrounding air in a rheometer Brown and Jaeger . In simulations this force can be fixed by the stiffness of the particles, which leads to a finite-yield stress even in the arrested state (in contrast to hard-sphere systems) Heussinger (2013). The resulting phase diagram has a re-entrant shape with a fluid-solid transition at low stress and a re-fluidization at higher stress Grob et al. (2014). Several variations of the basic system have been considered, including e.g. inertia Fernandez et al. (2013); Grob et al. (2014) or Brownian forces Mari et al. (2015). Subsequently, it turned out that ST represents an unsteady coexistence of different fluid states Grob et al. (2016); Maiti et al. (2016). Vorticity banding has also been observed Chacko et al. (2018), albeit with dynamic bands that move along the vorticity direction. Experiments also observe unsteady states, e.g. Saint-Michel et   al. Saint-Michel et al. (2018) report propagating bands at the onset of ST, and a proliferation of the dynamics deeper in the ST state. Rathee et al. Rathee et al. (2017) also find localized regions of increased stress that occur intermittently. Interestingly the size of these regions seem to grow with the gap-width of the rheometer, i.e. with the system size.

In this study we characterize the unsteady state of a model ST fluid. We build on our previous work published in Refs. Grob et al. (2016, 2014).

II Model

Our starting point are Newton’s equations of motion for a granular mixture of dry frictional particles in two space dimensions. Forces between particles arise at contact and are frictional. The particles are modeled as soft spheres of radius $R_{i}$, interacting with normal and tangential forces:

[TABLE]

Denoting particles’ positions and velocities by $\{{\bf r}_{i}\}$ and $\{{\bf v}_{i}\}$, the visco-elastic normal force can be written as Silbert et al. (2001)

[TABLE]

The unit vector, ${\bf n}_{ij}\equiv{\bf r}_{ij}/r_{ij}$, points from the center of particle $i$ to the center of particle $j$ and the particles only interact, when they overlap, $\delta_{ij}^{(n)}\equiv R_{i}+R_{j}-r_{ij}>0$. The $k^{(n)}$ and $\eta^{(n)}$ are the elastic and damping coefficients along the normal direction and $v_{ij}^{n}\equiv({\bf v}_{i}-{\bf v}_{j})\cdot{\bf n}_{ij}$ denotes the normal component of the relative velocity.

The tangential force between particle $i$ and $j$ is given by Silbert et al. (2001)

[TABLE]

where $k^{(t)}$ and $\eta^{(t)}$ are the elastic and damping coefficients along tangential direction ${\bf t}_{ij}$, defined by ${\bf t}_{ij}\cdot{\bf n}_{ij}=0$. The function $\mathrm{min(m,n)}$ yields the lower value between $m$ and $n$, enforcing the Coulomb criterion $|{\bf f}_{ij}^{(t)}|\leq\mu|{\bf f}_{ij}^{(n)}|$ with the friction coefficient $\mu$. The relative tangential velocity at contact $v_{ij}^{(t)}=({\bf v}_{i}-{\bf v}_{j})\cdot{\bf t}_{ij}+(R_{i}\omega_{i}+R_{j}\omega_{j})$ is the sum of a translational and a rotational contribution. It determines the tangential displacement according to $\delta_{ij}^{(t)}=\int_{t}^{t+\Delta t}dt~{}v_{ij}^{(t)}$, where the integration is over the time ineterval $[t,t+\Delta t]$, when the contact is present and not sliding.

We consider an equi-quaternary mixture with particle sizes: $2R_{A}=0.7$, $2R_{B}=0.8$, $2R_{C}=0.9$ and $2R_{D}=1.0$, all of equal mass $m$. The units of length, time and stress are chosen as $2R_{D}$, $(m/k^{(n)})^{1/2}$ and $k^{(n)}$. The friction coefficient has been set to $\mu=2$ in accordance with previous work Otsuki and Hayakawa (2011); Grob et al. (2016) and the damping constants are set to $\eta^{(n)}=\eta^{(t)}=1/2$.

A constant strainrate $\dot{\gamma}$ along the $x$-axis is imposed with help of Lee-Edwards boundary condition. Molecular dynamics simulations have been performed using the LAMMPS simulation package Plimpton (1995); Silbert et al. (2001) for $N=8000$ up to $N=80000$ particles at various packing fractions $\phi$.

III Results and Discussion

III.1 Global stress-strain relation

We want to understand the heterogeneous, time-dependent shear stress, which develops in response to an applied strainrate in large systems. Before analysing these structures in detail, we briefly recall the global stress-strain relations and discuss the non-trivial dependence on system size. Some representative flow curves are shown in Fig. 1. One observes the well-known Bagnold scaling at small $\dot{\gamma}$: $\sigma=\eta\dot{\gamma}^{2}$ and Herschel-Bulkley (HB) like behaviour $\sigma\propto\gamma^{1/2}$ for large $\dot{\gamma}$. For “small” volume fraction (black curve), these two regimes are connected by a smooth crossover and ST is relatively weak. For increasing volume fraction $\phi$, one observes stronger shear thickening, which is discontinuous for the small system (blue curve) and continuous for the large one (green curve). For even larger $\phi$ a finite yield stress is required for the system to flow (red curve, inset).

Discontinuous ST is visible as a sudden jump of the stress from the fluid to the HB branch, for $\phi=0.7975$ at a strainrate $\dot{\gamma}_{l}\approx 10^{-4}$. There is clearly a range of forbidden values of $\sigma$ and when observed with temporal resolution, one finds that the small system jumps frequently between the two branches. In Fig. 2 we show such a trajectory together with the distribution of $\sigma$-values for this run. Clearly the distribution is bimodal with the two peaks corresponding to the low $\sigma$ Bagnold regime and the high $\sigma$ HB regime. The jump is associated with hysteresis within a finite range of strainrates, $\dot{\gamma}_{l}\leq\dot{\gamma}\leq\dot{\gamma}_{u}$, when comparing a ramping simulation with slowly increasing strainrate (jump at $\dot{\gamma}_{u}$) with simulations with decreasing strainrate (jump at $\dot{\gamma}_{l}$).

These features bear some resemblance with the phenomenon of phase coexistence, where the jump represents the forbidden (unstable) region in the $\sigma-\dot{\gamma}$ phase diagram, and the “critical point” lies in the flowcurve with diverging slope $d\sigma/d\dot{\gamma}\to\infty$. However, as it turns out there is a non-trivial dependence on system-size which is not expected for that type of phenomena. As can be seen in Fig. 1, the discontinuity in the flow-curve is a finite-size effect and vanishes if the system is large enough. In thermodynamic systems with phase coexistence on the other hand, larger systems imply longer time scales for the nucleation of domain walls, and therefore a more pronounced discontinuity or hysteresis in finite-time simulations.

What happens in large systems? There is still a range of forbidden $\sigma$-values in the sense that in this range no stationary homogenous flow exists. Instead we observe heterogeneous time-dependent flow, whose properties will be analysed in the following paragraphs. The smooth curve in Fig. 1 does not represent stationary flow, but is obtained by averaging the flow over space and time.

The range of forbidden $\sigma$-values, depends on volume fraction and so does the the critical system-size for the vanishing of the discontinuity. While for $\phi=0.78$ the flowcurve is continuous already for $N=3600$, a much larger system ($N=32000$) is needed for the volume fraction $\phi=0.7975$. At even higher $\phi=0.801$ (inset) a discontinuity persists up to $N=80000$.

Interestingly, at this volume fraction one observes a discontinuity that does not lead upwards in stress when increasing the strainrate, but downwards (red curve, inset). Here, a metastable yield-stress branch becomes unstable at $\dot{\gamma}\approx 5\cdot 10^{-5}$. After the instability, at higher strainrates, the system follows a flowcurve that resembles the heterogeneous time-dependent flows that are displayed in the main panel. This downward jump is remarkable for two reasons. First, its presence depends on initial conditions. Choosing different starting configurations, one also observes a continuous branch (black curve, inset) which, instead of a yield-stress, shows Bagnold behavior at small strainrates. We expect the time-scale that governs this memory to diverge with increasing system size and also when crossing the jamming transition. Secondly, a downward jump opens room for additional instabilities in the form of shear-bands. Such a phenomenology of “discontinuous shear-thinning” has recently been discussed by one of us in Ref. Irani et al. (2018).

As a result of this section we conclude that ST in infinite systems always seems to be continuous, but flow is spatially heterogeneous and time-dependent. The actual system size that is necessary to reach this limit depends on volume fraction. For smaller volume fractions, like for the black curve in Fig. 1, the relatively small system ($N=8000$) is already large enough. Thus, the unsteadyness might already be present at the lowest volume-fractions that show only mild ST and could therefore be an inherent feature of any ST flow.

III.2 Heterogenous stress states

To better understand the heterogenous stress states in large systems, we analyse representative snapshots, such as those displayed in Fig. 3. The local shear stress of each particle is calculated as $\sigma_{i}=\sum_{j}F_{ij}^{x}y_{ij}+v_{i}^{y}(v_{i}^{x}-y_{i}\dot{\gamma})$, where $F_{ij}^{x}$ is the x-component of the interaction force between particles i and j, and $y_{ij}=y_{i}-y_{j}$ the distance in y-direction. Except for the largest strainrates, the kinetic contribution to the stress is negligible. Particles are colored black if their local stress exceeds the average value, and gray otherwise.

In the flowing state (panel a) linear structures that bear high stresses can be seen. These forcechains are oriented along the compressive direction of the flow (flow in positive x-axis) and seem to have a typical size of 5-10 particles. This length presumably depends on volume fraction and diverges at jamming Lois et al. (2007). The forcechains are distributed homogeneously throughout the system and seem to exist independently from each other. Similarly, for the highest strainrates (panel d), the flow is approximately homogeneous and time independent

In the continuous shear thickening (CST) region, on the other hand, the distribution of forcechains is rather inhomogeneous (b and c) and we observe large patches of high-stresses coexisting with regions of small stresses. Thus, the emerging CST in large systems is reminiscent of spatial coexistence of an inertial flow state (small-stress) and a plastic flow state (large-stress). However, coexistence cannot represent stationary two-dimensional flow. Hence the observed patterns are inherently time-dependent, as will be discussed in the next section. In small systems, spatial coexistence is not possible. Instead the whole system corresponds to a patch of either inertial or plastic flow and switches as a function of time between these two homogeneous states.

If we look in more detail then we can distinguish different features. Panel b is taken at a strainrate $\dot{\gamma}_{l}\approx 10^{-4}$ at the lower end of the CST region. Here, the force-chains seem to merge together to form giant bands that span the system in compressional direction, i.e. parallel to the chains themselves. Panel c on the other hand corresponds to a strainrate $\dot{\gamma}_{u}\approx 10^{-3}$ which is one order of magnitude larger than $\dot{\gamma}_{l}$ and marks the upper end of CST. Here, the structures do not seem to percolate in the compression direction but rather form clusters that extend along the dilational direction. Thus, the type of structures that form (clusters, bands) depends on and changes with the strainrate Cates et al. (1999).

To put these observations on a quantitative basis, we calculate the stress-stress spatial correlation function $C_{\sigma}({\mathbf{r}})=\frac{1}{N}\sum_{i\neq j}\langle\sigma_{i}\sigma_{j}\delta({\mathbf{r}}-{\mathbf{r}}_{ij})\rangle$. It is obvious from Fig. 3 that the correlations in the CST regime are anistropic as has been observed previously in experiment Majmudar and Behringer (2005). To capture these anisotropies we consider correlations in the compressive and dilational directions separately by restricting the separation vector ${\mathbf{r}}$ to the compressive and dilational direction, respectively. A length-scale $\xi$ can be extracted by monitoring the distance, at which the correlation function drops to a certain fraction ($C(\xi)=0.07\cdot C(0)$)). This value needs to be chosen small enough to not interfere with short distance effects, and large enough to avoid noise and finite-size effects due to the large distances involved. Within these bounds the results presented are independent of the value chosen.

The resulting correlation lengths in these directions are plotted in Fig. 4. Different system sizes are included. For the smallest system $N=8000$ no length scale is observed in the range of intermediate strainrates, where the CST regime resides. With an eye on Fig. 1 we confirm that this system does not show CST, but rather the discontinuous jump from Bagnold to HB. Once the system is large enough to show CST a large length-scale builds up that ever increases with system size. In fact, we find $\xi\propto L$, with $L$ the linear dimension of the simulation box.

The peak value of the correlation length is different in the two directions. In the compressive direction the maximum length is obtained at the lower end of CST, at $\dot{\gamma}_{l}$. In dilational direction, maximal correlation is achieved deep within CST and towards $\dot{\gamma}_{u}$. This proofs the visual impression from Fig. 3 that structures change their orientation.

To understand these local correlations in more detail, we also looked at other fields, such as the local pressure, the local density and the local connectivity. In Fig. 4 we show a profile of the local shear stress and the local connectivity along the compressive direction, averaged over the dilational direction. Stress $\sigma$ and particle connectivity $z$ are seen to be strongly correlated, as one might expect, because stresses are carried by contacts between the particles. The figure suggests a relation $\sigma=\sigma_{0}\exp(z/z_{0})$ with constants $\sigma_{0}\approx 5\cdot 10^{-7}$ and $z_{0}\approx 0.4$, thus stress is exponentially amplified if local connectivity is large enough.

Correlations of stress and density are also present, but less pronounced. A quantitative calculation gives a positive correlation coefficient of 0.07 for stress and density, when coarse-graining the fields on a grid with size 2. However, a small correlation might already induce strong effects. At volume-fractions close to the singular jamming point only small variations in density are necessary to change the character of the flow substantially. This is sometimes used as argument to favor pressure-controlled over volume-controlled flows, as such a singularity is not present when pressure is controlled Guazzelli and Pouliquen (2018).

We have also investigated different aspect ratios $L_{x}/L_{y}$ of the simulation box. Remarkably CST can also be observed for systems with only few particles, provided the system is large in the gradient direction ($L_{y}\gg L_{x}$). In the right part of Fig. 3e we show the local shear stress in a system of only 3000 particles with $L_{x}=10$ and $L_{y}=216$. The temporary shear banding is clearly seen and the global stress follows the same continuous flow curve (CST), as displayed in Fig.1 for the large system (green curve). At the same, a system with an equivalent number of particles but in a square simulation box ($N_{x}=N_{y}$) would feature a discontinuous flow-curve (as with the blue curve).

To conclude this section, CST is due to a large-scale coexistence of high- and low-stress regions. Most remarkable is the seemingly discontinuous onset of the system-spanning bands in compressive direction at the lower end of CST. This discontinuity exists whenever systems are large enough in gradient direction. Still the flowcurves are continuous. It is tempting to relate this to the system-spanning plugs that occur in impact experiments Allen et al. (2018). These plugs only occur above an onset stress and represent a solidified state of the shear-thickening suspension.

III.3 Dynamics of shear bands

The structures visible in Fig. 3 are not static but dynamically evolving. Interestingly, we can observe propagating modes with structures moving with a certain velocity. The clusters in panel c move in the direction parallel to the forcechains they are made of, i.e. in compressive direction. They do not stay coherent over a long time, rather they constantly split and join together over the course of time.

The system-spanning bands visible in panel b move perpendicular to their main axis, i.e. in dilational direction. Bands sometimes seem to evolve out of clusters that join together, percolating the entire system in compressive direction. Motion in this direction then stops and motion in dilational direction sets in. During propagation, the bands do not rotate but keep their orientation. This implies that the Lees-Edwards periodic boundary conditions in y-direction impose some offset that tends to destabilize the bands. Indeed, some time after its establishment we can observe deformations of the band until it splits in two or more parts. These parts then behave similar to the clusters described above.

The creation/destruction of the bands is also visible in the time-dependent global stress and pressure signals, as shown in Fig. 6. Both stress and pressure are highly correlated, in fact $\sigma=\mu_{\rm eff}p$ with an effective friction coefficient $\mu_{\rm eff}\approx 0.3$, quite similar to other studies Peyneau and Roux (2008). Whenever a band is present the global stress is high. Apparently, the signal displays some regularity, periodic oscillations with period of $\Gamma\approx 0.5\gamma$, i.e. a strain of $50\%$. We have checked that also the auto-correlation of the signal displays regular oscillations. Simulations performed without periodic boundary conditions in y-direction (with moving walls) also yield both propagating clusters and bands. The regularity of the oscillations, however, is strongly reduced.

We can characterize the time evolution of the bands by plotting the stress along the propagation (dilational) direction as a function of time. Such a kymograph is shown in Fig. 7 in the first panel. Each horizontal layer gives the instantaneous stress along the dilational direction and averaged over the compressional direction. Moving upwards layer by layer increases time. Motion of a band is therefore observed as a shift to the right or to the left. Four bands can be seen, two of which move to the left, two move to the right. In the X-shaped regions two bands collide and subsequently separate again without delay. From such kymographs band velocities can be estimated to be $|v^{band}|\approx 0.265$, with bands moving in both directions 111Band velocities of four bands are found to be   0.239, 0.269, -0.239 and -0.299. Negative sign indicates that band   moves away from the top layer..

We also display particle accelerations of the same sequence (Fig. 7 middle panel). It is clear that accelerations are large at the edge of the bands, where the stress gradient is high. This is a consequence of Newton’s equation of motion, or in continuum form, Navier-Stokes equation, $\partial\vec{v}/\partial t\propto\nabla\sigma$.

In the right panel we display the propagation of the clusters visible at larger strainrates $\dot{\gamma}=\dot{\gamma}_{u}$ (panel (c) of Fig. 3). Now the propagation is along the compressive direction, which is therefore taken as the spatial axis in the plot. Velocities of these clusters are a factor of two to three higher than for bands, $|v^{cl}|\approx 0.70$  222For three   bands the velocities come out to be -0.67, -0.62 and -0.78.. This velocity is on the order of the sound velocity $c_{s}=\sqrt{kd^{2}/m}=0.85$, with the relevant modulus taken as the spring constant $k^{(n)}$.

Thus, we propose to categorize the dynamics as follows. In the Bagnold regime at small strainrates, individual force-chains of a typical size form. They are oriented along the compressive direction, and move along this axis with the speed of sound. In the ST regime these force-chains join together to form larger clusters that still move in the compressive direction. The clusters themselves can join together, percolating along the compressive direction. Further motion in this direction is now blocked and the system could jam. However, associated stress/pressure gradients at the edges of the bands induce a new instability that leads to motion perpendicular to the band axis, in direction of the gradient. On a local scale this motion may be initiated via buckling-like events. Once started this motion is sustained, because the stress gradient drives particle accelerations. The resulting particle motion slightly densifies surrounding region, and allows more contacts to form. These new contacts are the foundation on which to build the large stress values, which are encountered in the center of the band.

IV Conclusion

A dense two-dimensional fluid of frictional grains undergoes a large scale instability for strainrates which are intermediate between the Bagnold and the HB regime. The unstable region is seen in $s-$shaped flow curves in experiment Pan et al. (2015) and simulations and has been predicted within a simple hydrodynamic analysis Grob et al. (2014). In the unstable region of strain and stress no stationary, homogenous flow is possible. What is observed instead, depends crucially on system size. In the small system, we observe hysteretic jumps between the two branches, the lower stress value corresponding approximately to Bagnold flow and the larger one to HB flow. Shear thickening is thus discontinuous. In a large system both stress values are present simultaneously in large bands or clusters. Such a state cannot be stationary, but is necessarily time-dependent. At first sight, one might expect diffusive behaviour; however, diffusion tends to equalize the stress profiles, implying intermediate values of stress which are not allowed. Hence the bands propagate in a wavelike fashion with a speed comparable to the speed of sound. A spatial average over such heterogeneous stress states gives rise to continuous shear thickening, so that in sufficiently large and presumably infinite systems ST is always continuous. What system size is necessary, depends on the volume-fraction and therefore on the distance to the jamming transition. Away from jamming, where only mild ST occurs, relatively small systems suffice, while closer to jamming, larger and larger systems are necessary to obtain a continuous flowcurve with a continuous shear thickening regime.

By looking closer at the unstable regime, we find that the flow is is qualitatively different at the upper and lower end of the unstable region. At the upper end of the ST regime, at large strainrates, we find large clusters, in which stress (and pressure) is by orders of magnitude larger than outside of the clusters. On the local scale these clusters consist of force-chains oriented along the compressive direction of the flow. Clusters are not stationary objects but move along the force-chain direction. They form and dissolve frequently. At smaller strainrates, towards the lower end of the ST regime, we find another large-scale structure: system-spanning stress bands. These bands span the system along the compressive direction. They seem to form out of percolating clusters, but their motion is perpendicular to the force-chains, along the dilational direction.

We conclude that large-scale instabilities give rise to continuous flowcurves, if spatially averaged. The heterogeneous time-dependent patterns come in with infinite correlation length, in the compressive direction at the lower end and in the tensile direction at the upper end. This sudden onset is reminiscent of the appearance of system-spanning plugs in impact experiments.

Acknowledgements.

We acknowledge financial support by the German Science Foundation (DFG) via the Heisenberg program (CH: HE-6322/2) and the SFB 937 (SS); AZ acknowledges support by the Heraeus-Stiftung.