Non-spherical granular ﬂows down inclined chutes

. In this work, we numerically examine the steady-state granular ﬂow of 3D non-spherical particles down an inclined plane. We use a hybrid CPU / GPU implementation of the discrete element method of non-spherical elongated particles. Thus, a systematic study of the system response is performed varying the particle aspect ratio and the plane inclination. Similarly to the case of spheres, we observe three well-deﬁned regimes: arresting ﬂows, steady uniform ﬂows and accelerating ﬂows. Both steady and dynamic macroscopic ﬁelds are derived from microscopic data, by time-averaging and spatial smoothing (coarse-graining), including density, velocity, as well as the kinetic and contact stress tensors. The internal morphology of the ﬂow was quanti-ﬁed exploring the solid fraction proﬁles and the particle orientation distribution. Furthermore, the system’s characteristic time and length scales are investigated in detail. Our aim is to achieve a continuum mechanical description of granular ﬂows composed of non-spherical particles based on the micromechanical details. Thus, to evaluate the inﬂuence of particle shape on the constitutive response in granular of those systems. However, to meet the proceeding’s page restrictions here we will only discuss the dependence of some terms of the continuum averaged equations on the coarse-graining scale, speciﬁcally the case of the kinetic part of the stress tensor.


Introduction
Granular flows [1] usually show quite complex behaviors that are also found in diverse systems, such as colloidal suspensions [2], pedestrian dynamics [3] and animal herds [4].Unfortunately, a complete theoretical framework or rheological model explaining the constitutive response of granular systems is not available yet.However, several promising approaches have recently been published [5][6][7], and it is known that several modes of energy dissipation, such as mechanic friction or plastic deformations lead to non-equilibrium steady state situations.
In examining granular flows there are several experimental restrictions, due to the opaque nature of the grains.Thus, full access to the 3D-behavior of the grains is generally not feasible and, hence, there is a real need to perform numerical simulations in this framework.Discrete element modeling (DEM) is widely accepted as an effective method to address engineering problems concerning dense granular media [8].Moreover, in typical applications the formulation of granular macroscopic fields are also necessary.The micro-mechanical details, i.e. velocity and position of individual particles, allow one to find the continuum field profiles using coarse-grain averaging technique [9][10][11][12][13].Furthermore, with this homogenization approach, the static and dynamic parts of the stress tensor e-mail: raulcruz@unav.ese-mail: t.weinhart@utwente.nl are deduced in terms of contact forces and velocity fluctuations, respectively.
Although spherical particles are a special case in nature, the majority of both experimental and numerical studies have considered packings of spherical particles so far.However, the appropriate description of systems composed of non-spherical particles is of significant importance for practical applications such as the handling of rocks, rice, wheat or tablets.Thus, there is currently and increasing interest in the behavior of non-spherical grains both experimentally [14][15][16][17] and numerically [18,19].The goal of this study is to extend existing theories of frictional avalanching flows of spherical particles in steady state situations, to the case of non-spherical grains.The ultimate aim is to achieve a continuum-mechanical description of granular flows composed by non-spherical grains based on micromechanics, This will help us to understand the role that particle shape plays in governing the system's mechanical response.
Our DEM simulations and the coarse-grain averaging technique [9][10][11][12][13] allows us to compute the mean velocity, density and stress fields in detail.However, the paper size limitation precludes us to present a complete description.Here we will only discuss the dependence of some terms of the continuum averaged equations on the size of the averaging domain, specifically the case of the kinetic part of the stress tensor, which is know to be scale-dependent [10][11][12][13] even for spherical particles.

DEM and Coarse-Graining Formulation
We have developed a hybrid GPU-CPU discrete element algorithm to investigate granular chute flows of rod-shaped particles.In the model, the rods are considered spherocylinders characterized by their length l and sphero-radius r.Thus, their aspect ratio is defined by ζ = (l + 2r)/2r (see figure 1).For calculating the interaction force between two particles i and j, F i j , we use an algorithm based on the concept of spherocylinder [20], which is defined by two vertices and the sphero-radius r; the surface of a spherocylinder is delineated by all points at distance r from the edge between the two vertices.Thus, the contact detection between to spherocylinders is reduced to find the closest point between two edges, resulting an overlap distance δ, defined by the overlap of two spheres of radius r.Fig. 1 illustrates the sketch of the interaction between two spherocylinders.Thus, the force F i j exerted on particle i by the particle j is defined by: F i j = − F ji .The force F i j can be decomposed as where F n is the component normal to the contact plane n.Additionally, F t acts on the tangential direction t.To define the normal interaction F n , we use a linear elastic force, which is governed by the overlap distance δ between two spherocylinders.To account for dissipation, a velocity dependent viscous damping is assumed.Hence, the total normal force reads as where k n is the spring constant in the normal direction, m r = m p 2 stands for the pair's reduced mass, γ n is the damping co- rel is the normal relative velocity between i and j.The tangential force F t also contains an elastic term and a tangential frictional term accounting for static friction between the grains.Taking into account Coulomb's friction constrain, which reads as, , where γ t is the damping coefficient in tangential direction, v T rel is the tangential component of the relative contact velocity of the overlapping pair.ξ represents the elastic elongation of an imaginary spring with spring constant k t at the contact, which increases as dξ(t) dt = v t rel as long as there is an overlap between the interacting particles.The elastic tangential elongation ξ is kept orthogonal to the normal vector (truncated if necessary) [11].μ is the friction coefficient of the particles.
We use a coordinate system where x denotes the flow direction, z the in-plane vorticity direction, and y the depth direction normal to the base.The chute is inclined at an angle θ such that gravity acts in the direction (sin(θ), −cos(θ), 0).The size of the simulation cell was L x = L z = 16 l and L y = 32 l with periodic boundaries conditions on the x− and z− directions.The base of the system is a rough surface consisting of N b = 385 fixed particles spherical particles with R = l/2.Note those conditions resemble one particular surface roughness of Ref. [12].Moreover, similar to Ref. [12] we have use k t = (2/7)k n , γ t = γ n with e n = 0.88 and ρ p = 2500.The microscopic friction coefficient is set to μ = 0.5, the gravity to g = 1 and k n = 2 × 10 5 m p g/l.
In the simulations presented here, all particles are mono-dispersed with the same length l = 1/8.Simulations are computed using rods of different aspect ratios, from ξ = 1.01 to ξ = 3.0, but always keeping the total mass of particles of M T ≈ 1.4 × 10 4 .In each case, the values of r and N p are adjusted to the choice of ζ (see Table 1).The N p flowing particles are introduced to the system at random non-overlapping positions well above the base.The gravity field induces the particle motion and they fall and accelerate down the plane until they reach a steady state, which is then analyzed.
In order to explore the dynamical and mechanical properties of the particle flow, a coarse graining methodology is used to analyze the results [9][10][11][12][13] The numerical data provide the position and velocities of every particle as well as the particle interaction forces.According to [10][11][12], the macroscopic mass density of a granular flow at time t is defined by ρ r, t = N p i=1 m i φ r − r i (t) where the sum runs over all the particles within the system and the coarse-grained (CG) function, φ( R).In our case, we use a truncated Gaussian coarse-graining function φ( R) = A w e −(|R|/2w) 2 with cutoff r c = 6w where the value of w defines the coarse-grained scale.A w is calculated in or- der to guarantee the normalization condition.Thus, the flow solid fraction can be found ϕ r, t = ρ r, t /ρ p , where ρ p is the material density.In the same way, the coarse grained momentum density function P( r, t) is defined by P( r, t) = N p i=1 m i v i φ r − r i (t) , where the v i represent the velocity of particle i.The macroscopic velocity field V( r, t) is then defined as the ratio of momentum and density fields, V( r, t) = P( r, t)/ρ( r, t).
On the other hand, the mean stress field σαβ is composed by the mean contact stress field σ c αβ ( r) and the mean kinetic stress field σ k αβ ( r).The spatial behavior of σ c αβ ( r) can be deduced in terms of the microscopic entities that characterize each contact, i.e.F i j and the position of the contacting particles r i and r j [10][11][12].In general the outcomes obtained for σ c αβ ( r) are independent of the coarse-graining scale.
Following Refs.[9][10][11][12], the kinetic stress reads as, where the sum runs over all the particles, v i is the velocity fluctuation of particle i, respect to the mean field.v i ( r, t) = v i (t) − V( r, t).Note that the velocity fluctuation are defined with respect to the center of the averaging volume at r. Very recently, Artoni and Richard [13] has proposed to define the velocity fluctuation of the particles respect to the particle location r i .As a result they found a decomposition of the kinetic stress tensor in two terms.The first, T k αβ ( r), that is independent of the averaging domain size and a new term Remarkably, the correction T γ αβ ( r) accounts for the dependence of averaging domain size, the used coarse-graining function, particle shape and polidespersity.Following their scheme the vector D read as, Hence, using Eqs( 1) and ( 2) one can find

Results and Discussion
The particles are introduced in the system at random positions.Thus gravity induces their motion and they fall and accelerate down the plane.Here, we will only discuss results corresponding to certain domain of angles in which the system reaches a steady state of motion on the x directions.To obtain detailed information about steady flows, we used the expressions defined above and [10][11][12].Since the flows are uniform in x and z, we further averaged over those directions.Thus, we obtained the time-, width-and length-averaged density ρ(y), momentum density P(y) and velocity V(y) = P(y)/ρ(y) fields (data not shown).The spatial profiles of the coarse contact stress, σ c αβ ( r), was   also computed in all cases (data not shown).Our outcomes are totally consistent with earlier findings [11,12].
It is important to remark that in dense granular flows, the values of σ k αβ ( r) generally result several orders of magnitude smaller than σ c αβ ( r).However, the spatial pattern of σ k αβ ( r) has been successfully used to distinguish between the role of force chains and velocity fluctuations in applications [21].That's why a precise calculation of σ k αβ ( r) is necessary.Next we discuss the dependence of some components of the tensor σ k αβ ( r) on the coarsegraining scale, as well as the role that the velocity gradients play.
In figure 2, the profiles of the xx component of the kinetic stress σ k xx ( r) are illustrated.We present systematic study, with results obtained using particles with four different elongations are shown.The fields were deduced using using a Gaussian coarse-graining function φ( R) with different values of w.As it noticeable, the values of σ k xx in general depends on w.Those results correlate with the fact that ∂V x ∂y 0, which causes to the dependency of the results on the averaging domain size.However, as it was obtained earlier [10][11][12], there is a length scale much less than a particle diameter d, where the kinetic stress can be accurately computed using Eq.( 1).The consistency of our outcomes is proven by the values of In figure 3, the profiles of the zz component of the kinetic stress σ k zz (y) are illustrated.As it is noticeable, the component σ k zz (y) in general do not depend on w and all the curves collapse.The fact that the only non zero component of the velocity gradient is ∂V x ∂y , leads to T γ zz (y) ≈ 0 and, accommodatingly, the outcomes of σ k zz (y) are independent of the averaging domain size.Similar outcomes are obtained for the profiles of the yy component of the kinetic stress σ k yy (y), which is consistent for symmetry reasons.Summarizing, we used a hybrid CPU/GPU implementation of DEM, examining the steady-state flow of 3D nonspherical particles down an inclined plane.A systematic study of the system response was performed varying the particle aspect ratio and the inclination angle.Similarly to the case of spheres, we observed three welldefined regimes: arresting flows, steady uniform flows and accelerating flows.Both steady and dynamic macroscopic fields were derived from microscopic data, by timeaveraging and spatial smoothing (coarse-graining).Here, we have discussed that the values of σ k xx (y), in general depends on w.This results correlates with the fact that ∂V x ∂y 0, with the definition of the velocity fluctuation with respect to the center of the averaging volume.However, we also found there is a length scale much less than a particle diameter d, where the dependence on the coarse graining scale is diminished.The consistency of our outcomes was tested comparing with the scale free kinetic stress T k xx (y) (see Eq.( 2) and Ref. [13]).In our case, we found that the relevant component of D k (ζ) depends on the coarse graining scale w, and the particle elongation ζ.Our ultimate aim is to achieve a continuum mechanical description of granular flows composed by non-sherical particles based on the micromechanical details and evaluate the influence of particle shape on the constitive response of those systems.

Figure 1 .
Figure 1.Snapshot of a system of 32768 spherocylinders with elongation of ξ = 3. Sketch of the interaction among two particles is also shown.
normal direction and v n

Figure 2 .
Figure 2. Vertical profiles of the xx component of the kinetic stress σ k xx (y) calculated using Eq.(1) as a function of y.The fields are deduced using using a Gaussian coarse-graining function φ( R) with different values of w. Results for several particle elongation are shown: a) ζ = 1.01, b) ζ = 1.5, c) ζ = 2.0, d) ζ = 3.0.

Figure 3 .
Figure 3. Vertical profiles of the zz component of the kinetic stress σ k zz (y) calculated using Eq.(1) as a function of y.The fields are deduced using using a Gaussian coarse-graining function φ( R) with different values of w. Results for several particle elongation are shown: a) ζ = 1.01, b) ζ = 1.5, c) ζ = 2.0, d) ζ = 3.0.
y), which are also shown for comparison in figure 3. Note that computing T γ xx (y) involves the calculation of the components of the vector D k (ζ), which depends on the coarse graining function φ( R), the coarse graining scale w 2 , and the particle elongation ζ.Once the system reaches a steady state of motion on the x directions we can find T γ xx (y) = D k (ζ) ∂V x ∂y 2 and the other components of T γ (y) diminishes.For the case of spherical particles and a Gaussian coarse-graining function φ( R) with w, result D k (ζ = 1) = w 2 .The values D k (ζ) in terms of D k (1) are shown for clarifying purposes.

Table 1 .
Particle's elongation ζ and number of particles used in each case.In all cases, we kept total mass M T ≈ 1.4 × 10 4 .