Effective friction of granular flows made of non-spherical particles

Understanding the rheology of dense granular matter is a long standing problem and is important both from the fundamental and the applied point of view. As the basic building blocks of granular materials are macroscopic particles, the nature of both the response to deformations and the dissipation is very different from that of molecular materials. In the absence of large gradients, the best approach formulates the constitutive equation as an effective friction: for sheared granular matter the ratio of the off-diagonal and the diagonal elements of the stress tensor depends only on dynamical parameters, in particular the inertial number. In this work we employ numerical simulations to extend this formalism to granular packings made of frictionless elongated particles. We measured how the shape of the particles affects the effective friction, volume fraction and first normal stress difference, and compared it to the spherical particle case. We had to introduce polydispersity in particle size in order to keep the systems of the more elongated particles disordered.


Introduction
After a few decades of research, the physics of granular material has reached today a good level of general understanding [1].In particular, it has been shown that the rheology of dense granular flows takes a simple form when working at controlled pressure P [2].It this case, the constitutive laws have been formalized in terms of the socalled inertial number I = γd/ P/ρ, where γ is the shear rate and d and ρ are the grain size and density [3][4][5][6].This dimensionless quantity compares the macroscopic shearing time scale 1/γ to the microscopic time scale of local plastic events d/ P/ρ.In the limit of rigid grains, the effective friction coefficient μ, which compares the shear stress τ to the pressure P, as well as the volume fraction φ of the system must be functions of I, and one now refers to this framework as the μ(I)-rheology for granular flows.The shape of the function μ(I) has first been determined by means of numerical simulations [4][5][6], also complemented by experimental measurements [7,8].Its behavior has been related to grain contact distribution, whose anisotropy increases with I [9].This rheology is able to recover Bagnold scaling in flows down an incline plane [10].More generally, this approach has been successfully applied in different geometries [2,3], e.g.silo discharge [11], granular chute flows [12][13][14][15] and granular column collapse [16,17], as well as to describe dynamic compressibility effects associated with spontaneous oscillatory motion [18,19].Interestingly, experiments and numerical simulations have recently provided evidence for limie-mail: somfai.ellak@wigner.mta.hutations of the μ(I)-rheology, when some non-local effects come into play (see e.g.[20] and references therein).Many works are now interested in possible extensions of this rheology, when another effect is added, such as cohesion [21], particle softness [22] or activity [23].This framework has also been extended to the case of granular suspensions [24,25].In this case, the dimensionless viscous number J, which compares the shear rate to the viscous frequency P/η f , where η f is the viscosity of the fluid, must be used instead of, or in complement to I. Interestingly, the case of Brownian suspensions can also be described in this framework [26].Finally, these constitutive laws have been incorporated to models for a diphasic description of sediment transport [27][28][29][30].Several recent studies investigated the behevior of elongated grains under shear [31][32][33][34][35][36][37].These studies mostly focused on quasistatic flows [31,32] or on the orientational distribution of particles as a function of grain elongation and shear rate [33][34][35][36].Stresses were analyzed as a function of particle shape for glued-spheres and true cylinders [38].In these systems, the shear induced alignment (average orientation of the grains) is not parallel to the streamlines but encloses a small angle with it.In the present work, we numerically consider the case of elongated grains shaped as spherocylinders, and investigate how the μ(I)-rheology is modified with respect to the aspect ratio Q of the grains.

Numerical system
In the numerical measurements we used the 3-dimensional plane Couette geometry [Fig.1(a)], where the boundary conditions were periodic in the flow (x) and neutral (y) directions, and Lees-Edwards in the velocity gradient (z) direction.The box was filled with elongated particles: we used spherocylinders with aspect ratio Q = /2R = 1.5, 2, 2.5 and 3, where the R is the radius of a particle (both of its spherical cap and of its cylindrical body), and is its total length.For reference we used spheres (Q = 1) as well.The interaction between the particles in this work was frictionless: the repulsive force was proportional to the overlap and contained an additional velocitydependent dissipative term: where F i j is the force between particle i and j, the overlap is δ i j = R i + R j − d i j ; the distance d i j and its unit length direction ĉij between the core line segments of the particles (connecting the centers of the spherical caps) is calculated by an efficient algorithm [39].The velocity difference between particles at contact, v c,i j , is calculated using the velocity difference between the centers of the particles, their rotation angular velocities and the relative position of the contacts.The prefactor b of the dissipative term is set such that the restitution coefficient e takes a specified value, we used e = 0.5.The length unit was the average particle diameter (not length), the mass unit was set by specifying unit density for the particles (they were solid bodies), and fixing the stiffness to k = 1 defined the force scale and consequently the pressure and time scales.
The time evolution was calculated by a velocity-Verlet scheme, with particle rotations represented by quaternions [40].The time step was taken to 1/100 of what would be the duration of a binary collision between two particles.
The initial conditions were prepared by a strongly overdamped dynamics, afterwards the system was sheared by a constant shear rate γ.During shear the L x and L y sides of the box were kept fixed, but L z was adjusted by a feedback loop such that the P z component of the stress tensor, monitored continuously, fluctuated around a predefined value P; we used P = 10 −3 in our units.
To ensure stationary state, before making measurements we sheared the system for shear strain γ = 25 starting from a random initial condition, or for γ = 10 starting from a final state of the same system sheared with a different shear rate; all transients have decayed to below noise level at this deformation.After ensuring no crystallization (see next Section), we checked that in the stationary state the density and particle level quantities like the average alignment angle and nematic order parameter were homogenous in the simulation box.The velocity profile was linear, consistent with the prescribed γ -i.e.no shear banding was observed.
After the stationary state has been reached, we measured the particle volume fraction φ and the stress tensor σ at regular intervals and recorded their time average over strain scale γ = 5.The stress tensor was defined as where v * i is the non-affine velocity of particle i in excess of the average shear velocity at that position; the second, kinetic term is negligible except at high shear rates.The most relevant elements of σ are the shear stress σ zx and the diagonal elements σ xx = −P x , σ yy = −P y and σ zz = −P z (the compressive pressures are positive).

Crystallization
With our numerical simulation we aimed to model the behavior of real, disordered granular systems.However we observed that our numerical model produces configurations ordered spatially in the y and z directions, when the shear rate was not too large for sufficiently elongated particles.
According to common wisdom 2D packings of monodisperse disks crystallize easily, while in 3D this phenomenon is less pronounced.It has been observed, however, that in 3D numerical systems of identical spheres crystallization may happen after sufficiently long shear [41] or shaking [42].Experiments of sheared identical spheres also displayed ordering, where chains of particles formed along the streamlines, and neighboring chains were arranged on a triangular lattice in the perpendicular plane [36].
We also observed in our simulations that small periodic packings (around 500 particles) crystallized easily unders shear, though that order was quickly destroyed by a few percent of polydispersity, slight deviation from spherical shape, and even by increasing the system size.
For the typical system size used in this work (1000 particles) no crystallization was observed for spheres and Q = 2 spherocylinders, but those with Q = 3 displayed strong spatial order.To quantify this we calculated the pair correlation function in the neutral y direction, the results are shown in Fig. 2. In order to destroy this order we introduced polydispersity: the particles were of identical aspect ratio, but their radii were distributed according to a uniform distribution, characterized by the variation coefficient δ R (the ratio of standard deviation of radii to their average).At polydispersity δ R = 0.1 no spatial order has been observed for the Q = 3 spherocylinders; for consistency we used this value of polydispersity for all particle shapes in the following rheological measurements.

Rheological laws
Our main goal in this work was the measurement of rheological quantities as a function of the inertial number I for different particle shapes; the results are shown on Fig. 3.
The effective friction, defined as the ratio of the shear stress and the control pressure, μ = σ zx /P z , is shown on Fig. 3(a).The general trend of the elongated particles is similar to that of spheres: μ tends to a nonzero value in the quasistatic (I → 0) limit even for frictionless particles, and gradually increases with the inertial number.It is interesting to note that the dependence appears to be non-monotonic with Q, as in the quasistatic limit the effective friction for Q = 1.5 appears to be larger than for either spheres or more elongated particles; though this value needs to be extracted systematically for a sharp statement.
The particle volume fraction φ is plotted in Fig. 3(b).For spheres in the quasistatic limit it tends to a value close to that quoted for random close packing, though for comparing actual values one needs to take into account the polydispersity we used.For elongated particles the quasistatic volume fraction is higher, similarly to what has been observed for ellipsoids [43].For increasing inertial numbers the constant pressure volume starts to increase as the collisions become increasingly violent, consequently the volume fraction starts to drop.
Finally the first normal stress difference was measured as well, which we normalize by the control pressure as N 1 /P z = (P x − P z )/P z ; see Fig. 3(c).As expected, for spheres this quantity is zero up to moderate inertial numbers, and starts to deviate from zero only in the collisional regime.For elongated particles the normal stress difference is nonzero even for small inertial numbers.The quasistatic limit needs to be determined carefully before a conclusive statement can be made about its value for spherocylinders.The second normal stress difference (not shown) is smaller in magnitude, its sign is positive, and the values for Q > 1.5 seem to collapse on each other.
In summary, we performed numerical simulations to determine the rheological properties of homogenously sheared granular packings of frictionless spherocylinders.We found that the qualitative behavior of the effective friction is similar to that of spheres: takes a nonzero value in the quasistatic limit, and increases with the inertial number; furthermore its aspect ratio dependence is nonmonotonic.The volume fraction of sheared spherocylinders is higher than spheres for all aspect ratios we considered.The first normal stress difference for elongated particles deviates from zero even for small inertial numbers, unlike for spheres.
Our short term plans include extracting systematically the quasistatic limit of μ, φ and the normal stress differences, determining their aspect ratio dependence, and trying to establish a link between these macroscopic rheological properties and the particle level behavior like shear induced alignment and coordination of elongated particles.

Figure 1 .
Figure 1.Numerical setup.(a) The simulation box.The velocity difference between the top and bottom sides is v = γL z , where L z is adjusted continuously as if a contstant pressure was applied to the top and bottom sides.(b) Particles used: sphere (Q = 1) and spherocylinders with aspect ratio Q = 1.5, 2, 2.5 and 3.

Figure 2 .
Figure 2. Crystallization.The pair correlation function is plotted in the neutral y direction for (a) spheres, (b) Q = 2 and (c) Q = 3 spherocylinders.Monodisperse spheres and Q = 2 spherocylinders show no crystallization: the pair correlation function rapidly converges to its asymptotic value of 1 after a few oscillations.However, for Q = 3 monodisperse spherocylinders the peaks are very sharp, and the oscillation is still very strong even at the 7th neighbor.Introducing polydispersity at δ R = 0.1 (red curve) restores disorder.The insets are side views of the simulation box taken from the x direction.The ordering for Q = 3, δ R = 0 is clearly visible.

Figure 3 .
Figure 3. Rheological laws: (a) effective friction coefficient, (b) volume fraction and (c) normalized first normal stress difference as a function of I for five different particle shapes.