A dilation-driven vortex flow in sheared granular materials explains a rheometric anomaly

Granular flows occur widely in nature and industry, yet a continuum description that captures their important features is yet not at hand. Recent experiments on granular materials sheared in a cylindrical Couette device revealed a puzzling anomaly, wherein all components of the stress rise nearly exponentially with depth. Here we show, using particle dynamics simulations and imaging experiments, that the stress anomaly arises from a remarkable vortex flow. For the entire range of fill heights explored, we observe a single toroidal vortex that spans the entire Couette cell and whose sense is opposite to the uppermost Taylor vortex in a fluid. We show that the vortex is driven by a combination of shear-induced dilation, a phenomenon that has no analogue in fluids, and gravity flow. Dilatancy is an important feature of granular mechanics, but not adequately incorporated in existing models.

F lowing granular materials exhibit superficial similarities with fluids, but the details of their mechanical response are significantly different, as the micromechanics of grain interactions is quite distinct from that of fluid molecules 1,2 . The design of industrial processes and understanding of natural phenomena involving granular flows call for the development of reliable continuum models. Although several advances have been made in this direction [2][3][4][5][6][7][8][9] , most models are confined to narrow-flow regimes, or leave out key features. For critical evaluation and refinement of available models, experimental measurement of the rheology is indispensable, for which the protocols of fluid rheometry may be usefully employed.
The rheological properties of fluids are most conveniently measured using devices that generate a class of simple shear flows 10 . The cylindrical Couette cell is such a device, but when it is used for granular materials, unexpected behaviour emerges. Recent studies that used this device for granular rheometry found a striking anomaly in the stress 11,12 : the vertical shear stress changes sign on shearing, and the magnitudes of all components of the stress increase roughly exponentially with depth. This behaviour is contrary to previous experiments 13,14 and the predictions of plasticity theories 8,13 , which yield a fluid-like stress. After considering several plausible mechanisms, it was speculated 11,12 that an anisotropic microstructure 15,16 is the likely cause of the anomalous stress.
Here we show the cause to be one of the least-expected mechanisms-a single vortex spanning the entire granular column. Vortices arise in fluids too when the rotation rate of the inner cylinder exceeds a critical value, due to the centrifugal Taylor-Couette instability 17 ; this instability has also been observed in a fluidized granular bed sheared at high rates 18 . We show that the vortex in a slowly sheared granular material is fundamentally different in its origin and manifestation.

Results
DEM simulations and experiments. To understand the cause of the anomalous stress, we conducted simulations using the discrete element method (DEM) 19,20 (see Methods). The simulations were complemented by video imaging the free surface of an experimental Couette device. In the simulations, the annular gap of a cylindrical Couette cell ( Fig. 1) with smooth or rough walls is filled with spherical grains of density r p and mean diameter d p to a fill height H by two methods, raining and dense packing. The inner cylinder is then rotated at a constant angular speed O until a state of steady shear is reached. We varied O over a range that spans the quasistatic and lower end of the intermediate regimes, which are delineated by the Savage number Sa, defined as the ratio of the stress due to impulsive grain collisions to the total stress (Supplementary Note 1). Unless stated otherwise, the DEM results shown in the paper are for smooth walls, filling by raining and Sa ¼ 2 Â 10 À 6 ; however, the qualitative features of the kinematics and stress are independent of the wall roughness, the Savage number, the parameters in the grain contact model and the method of filling. In the experiments, the steady-state velocity profile at the free surface is determined by video imaging and particle image velocimetry (see Methods).
Validation of the DEM simulations. The results of our DEM simulations are first validated by comparison with experimental data. The kinematics of the azimuthal flow is compared with the data obtained from video imaging (Fig. 2a). The simulations reproduce the exponential decay of v y with distance from the inner cylinder, a feature that is characteristic of slow granular flows [21][22][23] . The mechanics is validated by comparing the stress at the outer cylinder with the data of Gutam et al. 12 (Fig. 2b,c) anomalous stress observed in the experiments: in a static column, the shear stress s rz is positive (that is, the vertical component on the traction on the wall is downwards) and all components of the stress saturate asymptotically with the depth z, in accord with the classical Janssen solution 24,25 (Supplementary Note 2); on shearing, s rz changes sign and the normal stress s rr acquires a positive curvature. The quantitative differences between the simulation and experimental data may be partly due to the Couette cell in the simulations being much smaller, for computational tractability; moreover, we do not attempt to tune the parameters of the contact model to achieve agreement, as our primary aim is to gain a qualitative understanding of the cause of the anomalous stress. Nevertheless, our qualitative validation is useful because the features we verify are significant: the exponential decay of v y (r), a characteristic feature of slow granular flows, and the reversal in the sign of s rz on shearing, a key finding of our previous experimental studies.
Starting from an initial state generated by raining, shearing results in overall compaction, as seen from the fall of the free surface ( Fig. 3a,b) and the profiles of the solids fraction (Fig. 3d), because raining yields a loosely packed initial state. The low-density layers adjacent to the two cylinders are due to the wall-imposed constraints on packing. Nevertheless, significant dilation in the shear layer is evident from the broadening of the low-density layer near the inner cylinder, in accord with numerous experimental observations 2,26 . The volume-fraction distributions at steady state for the two filling methods are almost identical ( Fig. 3b-d), leading us to the conclusion that the initial preparation of the granular bed has little bearing on the properties at steady state.
The secondary vortex. Further examination of the kinematics reveals a secondary flow in the r-z plane, superimposed over the azimuthal flow. Although its velocity scale is small compared to that of the azimuthal flow, the radial and vertical velocities v r and v z are easily measurable, and show a systematic trend (Fig. 4a,b). The shapes of the v r profiles at the free surface determined from experiment and DEM simulations are in good agreement. When the velocities at all the locations are combined to construct the streamlines, a single vortex that extends over the entire width and height of the Couette gap emerges ( We note that the free surface slopes downward from the outer to inner cylinders (Fig. 3b,c), a feature that is also observed in the experiments; while the time taken for the free surface slope and the secondary flow to reach steady state roughly coincide, our evidence is insufficient to infer whether the slope is caused by the secondary flow, or another aspect of the mechanics (such as normal stress differences 10 ).
The vortex explains the anomalous stress. To show that the vortex flow explains the anomalous stress, we consider the simplest plasticity model for the stress tensor r 28,29 , where v is the velocity vector, D the deviatoric part of the deformation rate tensor (with scalar norm _ g 2D : D ½ 1=2 ), d the identity tensor, and p c (f) is the pressure at the critical state 2,5 . The parameters in the model are the bulk and shear plastic moduli, m b and m s , and the exponent m. The above is a model for rate-independent plasticity if m b and m s are  independent of _ g, and rate-dependent plasticity otherwise 9,30 . Applying the model to the problem at hand, at the two cylinders impenetrability requires v r ¼ 0, whence D rz ¼ qv z /qr. From the radial variation of v z , we see that qv z /qr40 at the outer cylinder (inset of Fig. 4b). Equation (1) then implies that s rz o0, which is in accord with experimental observations 11,12 (Fig. 2c). Thus, the vortex explains the reversal in the sign of s rz when the granular column is sheared. This result, in conjunction with Coulomb friction at the boundaries, explains the sharp rise of all components of the stress with z 11,12 (Supplementary Note 2; Supplementary Fig. 2), and thereby all aspects of the anomalous stress. Although we have used a simple plasticity model to illustrate the link between the vortex and the anomalous stress, more elaborate models 2,29 lead to the same conclusion.
Mechanism for the secondary vortex. We now address the cause of the secondary flow. The sense of the vortex indicates that it is not centrifugal in origin, but for a definitive confirmation, we conduct a simulation of shear between two plane vertical walls (Fig. 5a). Figure 5b shows the secondary flow that results-two counter-rotating vortices placed symmetrically about the mid-plane of the Couette gap (see also Supplementary Movie 3). Symmetry arises here because the two walls are indistinguishable, unlike in cylindrical Couette flow. Thus, the vortices closely resembling those in Fig. 4c-e, arise in the complete absence of the centrifugal force.
What then drives the vortex? This question is answered by considering the transient evolution of the secondary flow soon after initiation of shear. The streamlines in Fig. 5c-e show that the flow is initially radially outwards and upwards, and concurrently, there is downward flow close to the inner cylinder.
These two flows combine to curl the streamlines towards the inner cylinder at later time, eventually establishing a steady vortex (Fig. 4c) within a 90°rotation of the inner cylinder. When gravity is turned off and a frictionless wall used to confine the material from the top, we find no secondary vortex. It is now clear that the vortex is driven by the combination of shear-induced dilation and gravity flow-dilation causes the material to flow away from the shear layer and ultimately to the free surface, and the downward flow of grains near the inner cylinder sustains the vortex. The sense of the vortex is determined by the asymmetry in dilation between the inner and outer cylinders (Fig. 3), as a result of the shear rate of the primary flow decaying exponentially with radial distance from the inner cylinder (Fig. 2a).
We have used the transient evolution of the vortex to elucidate the mechanism, but it applies equally well at steady state. In the absence of gravity flow, dilation would cease after the initial transient, and a steady state would be reached wherein the flow is purely azimuthal and there is a radial density gradient. While gravity flow is essential to sustain the vortex, we surmise that it is initiated by dilation, as we cannot think of another physical mechanism that would cause a radial flow.
Importantly, we find the secondary vortex to be present even when the upper surface is not traction-free. Supplementary Figure 3 shows the streamlines at steady state for plane Couette flow when the material is confined at the top by a horizontal rigid plate of fixed weight, which is allowed to move vertically. The presence of two counter-rotating vortices is evident, as in the case of a traction-free surface (Fig. 5b), though their strength and symmetry decrease as the weight of the confining plate increases. We find the vortex to persist in cylindrical Couette flow too when confined at the top, and its basic structure is largely independent of the confinement condition and plate roughness (K.P.K., P. V. Dsouza, T. Murthy and P.R.N., manuscript in preparation). Thus, a traction-free surface is not essential for the formation of the vortex-the combined effects of dilation and gravity break the up-down symmetry, and determine the sense of the vortex.

Discussion
The dilation-driven secondary flow that we report is a novel phenomenon that has no analogue in fluids. Although the scale of the secondary flow is smaller than the primary azimuthal flow, its rheological signature is large, owing to the dependence of the shear stresses on the pressure (equation (1)). Evidence of a secondary flow was provided by two previous studies 31, 32 , the former hypothesizing dilatancy as a possible cause and the latter finding convection to vanish under microgravity conditions; however, they could not discern the detailed form of the secondary flow and its origin. Our results establish that there is a single vortex (for the range of fill heights we have explored), resulting from the combined effects of dilatancy and gravity flow.
It has been known for long that dilation accompanies deformation in dense granular materials, and, conversely, compaction when the material is loose 2,33 . This key feature distinguishes plastic deformation of granular materials from that of metals, and of course fluids. Existing rheological models for granular flows are inadequate for capturing compressibility; while density change along a streamline is predicted (such as in flow through a hopper), to our knowledge no model captures the shear-induced density gradient across streamlines 8 . It is usually assumed that the effect of dilatancy is confined to the (typically) narrow shear layer, but our study establishes that it can act as a driving force for a large-scale flow. Thus, despite some advances, much remains to be done towards formulating a robust and effective rheological model for granular materials-we expect that our study will provide a useful impetus in this direction.
The Taylor vortex in fluids arises from an instability of the base state of azimuthal flow, due to an imbalance between the centrifugal force and the radial pressure gradient 17 . It is quite possible that the granular vortex arises from a similar hydrodynamic instability, with dilation as the driving force. However, a stability analysis must await the formulation of a rheological model that captures cross-streamline dilation. We finally note that the anomalous stress reported by our group 11,12 was then ascribed to microstructural anisotropy; while we now offer a more compelling explanation, there is sufficient evidence in the literature to motivate the incorporation of an anisotropic fabric in a comprehensive rheological model.

Methods
DEM simulations. DEM is a widely used computational tool for granular mechanics, where the positions and interactions of all the particles are tracked in time using simple models for grain interactions. We used the soft-particle contact model 19,20 , wherein the grains are treated as deformable spheres, but rather than calculate their deformation in detail, they are allowed to overlap. The interaction forces are written in terms of the overlap and its time rate of change. The normal and tangential forces have elastic and viscous components, and the tangential force incorporates an additional Coulomb slider to capture rate-independent friction ( Supplementary Fig. 4a). Considering the contact of grains i and j centred at position vectors r i and r j , the normal and tangential forces on particle i are 34 respectively. Here dR i þ R j À |r i À r j | is the overlap, n is the unit normal at the point of contact pointing towards the centre of i, v n and v t are the relative velocities at the point of contact in the normal and tangential directions, respectively, Ds is the tangential displacement and m eff (1/m i þ 1/m j ) À 1 is the effective mass that determines the damping force. The spring stiffness constants k n , k t , the damping constants g n , g t and the Coulomb friction coefficient m are the parameters that define the interaction. The motion of each particle is determined by integrating Newton's second law, assuming pairwise additivity of the interaction forces. The computationally intensive simulations, involving the tracking of B5 Â 10 5 particles for durations corresponding to several rotations of the inner cylinder, were conducted using the LAMMPS package 35 .
To avoid crystalline order, a mixture of grains of diameter 0.9 d p , d p and 1.1 d p (of mass fractions 0.3, 0.4 and 0.3, respectively) was used. Simulations were conducted for two types of boundaries (cylinders and bottom wall): the surfaces were either smooth, or coated with a rigid, close-packed triangular lattice of grains of diameter 0.9 d p ; the two cases are referred to as smooth and rough walls, respectively. For both types of walls, the grain-wall interactions were computed by treating the walls as bodies of infinite mass, but with the same stiffness, damping and friction constants as for grain-grain interactions.
For values of the spring constants k n and k t that are appropriate for hard particles like sand and glass beads, accurately resolving each contact requires so small a time step of integration that the computation time becomes prohibitively high. It is therefore standard practice to optimize their values such that they are low enough for the computations to be tractable, but high enough that the macroscopic behaviour mimics that of hard particles. The values chosen for the simulations are k n ¼ 10 6 m p g/d p , k t ¼ 2/7k n , as per previous studies that have modelled hard grains 16,34 . Here m p is the mass of the particle of diameter d p and g is the gravitational acceleration. The damping coefficients were set to g n ¼ 317 m p ffiffiffiffiffiffiffiffiffi ffi g=d p p , g t ¼ 1/2g n and the friction coefficient to m ¼ 0.5. The chosen value of g n yields a normal coefficient of restitution for a collision of 0.7 (ref. 16). To verify that the results do not depend qualitatively on the parameter values 36 , simulations were conducted for three values of k n varying over two decades. The profiles of the normal and vertical shear stress at the outer cylinder are shown in Supplementary Fig. 4b,c. Increasing k n results in larger magnitudes of the stress components, but the qualitative features that characterize the stress anomaly remain independent of k n .
The Couette cell was filled by two methods, raining and dense packing. In the former, the mixture of grains (of sizes mentioned above) was poured uniformly over the annular gap from a reservoir at the top of the Couette cell, until the fill height H was reached. In the latter, grains of uniform size 1.1 d p were placed in a body-centred cubic lattice within the Couette cell up to a fill height H, in the absence of gravity; the grains were then randomly shrunk to achieve the required size distribution, after which gravity was turned on. The two methods gave an initial average volume fraction of hfi ¼ 0.576 and 0.602, respectively.
The continuum variables were obtained by averaging over space and time. The velocity and volume fraction at coordinates (r, z) were obtained by averaging over all the particles in an annular cylinder of thickness Dr ¼ d p (centred at r) and depth Dz ¼ 3 d p (centred at z); the wall stresses were obtained by averaging over all particle-wall contacts in such a cylinder adjacent to the outer cylinder. The steady-state variables (Figs 2-4 and 5b) were time-averaged for half a period of rotation of the inner cylinder, but the transient velocity fields were averaged over much shorter durations (see caption of Fig. 5).
Experiments. The experiments used spherical beads of soda-lime glass of density r p ¼ 2,500 kg m À 3 and a narrow size distribution with mean diameter d p ¼ 0.83 mm. The Couette apparatus (dimensions in caption of Fig. 1) was filled with the glass beads, and the inner cylinder rotated at 2 revolutions per minute (corresponding to a Savage number Sa ¼ 3.3 Â 10 À 6 ) for a period of about 30 min to reach steady state. Video images of the free surface were then acquired using a digital video camera (Prosilica GE680) at a rate of 200 frames per second and resolution 640 Â 480 pixels (scale factor 0.112 mm per pixel). The instantaneous velocity field was determined by cross-correlation of the intensity maps (averaged over bins of 12 Â 12 pixels) in successive frames using the PIVlab 37 toolbox in MATLAB. The steady-state velocity profile was determined by time-averaging the instantaneous velocity field for 30 s.