Flow and rheology of frictional elongated grains

The rheology of a 3-dimensional granular system consisting of frictional elongated particles was investigated by means of discrete element model (DEM) calculations. A homogenous shear flow of frictional spherocyliders was simulated, and a number of rheological quantities were calculated. In the framework of the $\mu(I)$ rheology, the effective friction was found to be a non-monotonic function of the aspect ratio for interparticle friction coefficient $\mu_p \lesssim 0.4$, while it was an increasing function for larger $\mu_p$. We reveal the microscopic origin of this peculiar non-monotonic behavior. We show the non-trivial dependence of the velocity fluctuations on the dissipation regime, and trace back the behavior of the normal stress differences to particle-level quantities.


Introduction
Dense granular flows, with their rich phenomenology, are of great interest for fundamental questions as well as for their relevance in applied problems. They have been the subject of a large number of studies in the last decades, with experimental, theoretical and numerical approaches, see the general book [1]. Taking advantage of the constant increase of computational power, numerical simulations of granular systems have now reached an impressive degree of realism, allowing them for reliable predictions, even in rather sophisticated configurations [2]. In an effort to go beyond the ideal case of frictionless hard spheres, particles with various shapes [3][4][5][6][7][8][9][10][11][12][13][14][15] and surface or bulk conditions [16][17][18][19][20][21][22][23][24] (e.g., friction, cohesion, stiffness), have been modelled and investigated.
Here, our interest is focused on the rheological behaviour of assemblies of frictional elongated grains close to jamming. The fundamental question is, how large is the resistance (i.e., the effective friction µ) of the material against slow shearing, and how this effective friction changes with grain elongation. In such systems the shear flow induces particle rotation which leads to more intensive collisions between neighbouring particles than for spherical grains. The speed of the shear induced rotation depends on the particle orientation, faster rotation for particles parallel to the shear gradient and slower rotation for particles pointing in the flow direction, which results in orientational order. Both of these phenomenacollisions due to rotation and orientational ordering -affect the flowing and mechanical properties of the system [25][26][27][28]. This problem has previously only been addressed numerically in simplified situations: either with frictionless grains [29] or in a 2D system [30]. Those studies revealed an unexpected, peculiar behaviour: the effective friction was found to be non-monotonic (increasing and decreasing) with increasing aspect ratio α for the 3D frictionless case [29], and this non-monotonic tendency was shown to persist even for frictional grains in a 2D system, although only at small values of the interparticle friction coefficient (up to around µ p = 0.15) [30]. In a real world situation the interparticle friction is significantly larger (µ p is around 0.3). Moreover, a 3D system is substantially different from the 2D case as the rotating particles have extra degrees of freedom to evade and reduce the effect of collisions, which is not possible in 2D. DEM simulations provide both macroscopic and microscopic information about these processes. An effective way is to use a pressure-imposed shear geometry, where the constitutive laws for elongated grains followed [29,30] the general framework of the so-called µ(I) rheology [4,[31][32][33]. In this description the effective friction µ = τ /P as well as the volume fraction φ are functions of the inertial number I =γd/ P/ρ, where τ is the shear stress, P the pressure,γ is the shear rate, d is the grain size and ρ is the grain density. These rheological functions also depend on the aspect ratio α of the grains. In addition to the above mentioned non-monotonic behaviour of the effective friction on α, we found that these flows develop normal stress differences for α > 1.
The value of the interparticle friction has an important role in the rheological behavior of dense granular flows. One example is that numerical simulations with circular or spherical grains have shown [4] that the effective friction in the low shear rate limit µ c is an increasing function of µ p . Another one concerns the exponent of the constitutive law (see Eq. 1 below). In the simulations of Favier et al. [21] it is found that the exponent of the power law term switches from β = 0.5 in a low friction limit (µ p 10 −2 ) to β = 1 in the high friction limit (µ p 10 −1 ). Three regimes have been identified and associated with different dissipation mechanisms [34,35].
In this work we extend the µ(I) rheology to the case of a 3-dimensional system of frictional spherocylinders. We show, that in a realistic 3D frictional system the peculiar non-monotonic behaviour of µ(α) is observed in a much wider range of interparticle friction (up to around µ p = 0.4, thus including common granular materials) than in the previously reported case of a simplified 2D system (up to µ p = 0.15) [30]. We reveal the microscopic origin of this observation and relate the behaviour of the constitutive coefficients to the above mentioned dissipation regimes. The paper is organised as follows: in Section 2 we briefly recall the numerical setup of the simulations, which builds on that of our earlier work [29]. In Section 3 we present and discuss our numerical results, explain a number of observed phenomena, including scaling arguments to understand and interpret the data. We summarize and draw perspectives in Section 4.

Setup
We performed numerical simulations of homogenous shear flow in 3D plane Couette geometry, see Fig. 1. The particles were frictional spherocylinders with length-to-diameter aspect ratio α = /2R. The interparticle force F ij which particle j exerts on particle i consists of a soft core repulsion and a tangential Here we plot only one particle near the center of the box (orange) and particles in its hemispherical neighborhood (blue); their aspect ratio is α = 2. The system is sheared with shear rateγ = dvx/dy, and a constant normal stress σyy = py is maintained by adjusting box side Ly by a feedback loop. component due to friction. The repulsive force points in the direction of the local surface normalĉ ij , its amplitude is proportional to the virtual overlap δ ij between the particles, and it contains a dissipative term proportional to the velocity difference v c,ij at the contact point: The prefactor b was determined by requiring a given coefficient of restitution e for binary collisions (we used e = 0.5 in this paper). The frictional tangential force is based on the Mindlin force law: the force increment between time steps is ∆F fric ij = k∆t c,ij , where ∆t c,ij is the tangential displacement (projected to the plane perpendicular toĉ ij ) during a time step between the touching contact points of the particles. The magnitude of the frictional force is limited by the interparticle friction coefficient: |F fric ij | < µ p |F rep ij | . The length, time and mass units of the simulation were set implicitly by setting the mean particle diameter 2R, density ρ and contact stiffness k (equal for the normal and tangential force) to unity. To prevent crystallization for frictionless particles [36], especially at larger values of α, we used size polydispersity of 10% (standard deviation to mean ratio in a uniform distribution). While crystallization was less critical for frictional particles, we kept the polydispersity fixed for consistency. The rheological measurements were performed under fixed normal stress: the p y := −σ yy component of the stress tensor (where y is the velocity gradient direction) was controlled around a fixed value of 10 −3 by a feedback loop adjusting the L y side of the periodic simulation box.
The rest of the simulation details, including the preparation protocol for the initial conditions, are detailed in Ref. 29.

Rheology
The inertial number dependence of the following rheological parameters are shown on the top row of Fig. 2: effective friction µ = σ xy /p y , normalized first normal stress difference N 1 /p y = (σ xx − σ yy )/p y , and normalized second normal stress difference N 2 /p y = (σ yy − σ zz )/p y . The solid curves are fit (in the range 10 −3 ≤ I ≤ 10 −1 ) to the following empirical form: and similarly for N 1 /p y and N 2 /p y . The values of the exponent β ranged between about 0.4 and slightly more than 1, where we observed the smallest values for the frictionless case, and the largest values for   3.1.1. Effective friction One remarkable finding is that the quasistatic effective friction coefficient µ c is a non-monotonic function of the aspect ratio for negligible to moderate interparticle friction (µ p 0.4), and the aspect ratio where the maximum occurs shifts to larger values for increasing µ p . For large µ p the effective friction is a monotonically increasing function of the aspect ratio, at least in the range 1 ≤ α ≤ 3 we explored. ‡ As discussed above, the nonmonotonicity of µ c (α) has been observed earlier for purely frictionless spherocylinders [29], and for 2D ellipses with low friction coefficient [30]. The data in Fig. 2(d) clearly ‡ We note that µc extrapolated from high inertial number non-homogenous flow down an inclined plane [37] (for µp = 0.5) is remarkably close to our measurements.  Fig. 2, with the blue full circles with dashed lines corresponding to random close packed simulations without particle friction in non-sheared systems [3].
show, that in a 3D system this is observed in an extended friction range, which already includes realistic materials. We now explore the microscopic ingredients leading to this behaviour.
The explanation can be traced back to the shear induced alignment of elongated particles, initially described in Refs. 25 and 38. Let us denote the deviation of the particle axis from the streamlines by θ within the x-y plane, and by ϕ out of this plane. Due to shear the elongated particles develop nematic ordering, where θ , which we call shear alignment angle, is interestingly non-zero [see Fig. 3(a)] §, while ϕ = 0 by symmetry. (With increasing elongation α the distributions of these angles become typically narrower, see Fig. 3(b-c). (Due to the periodicity directional statistics have to be used, and since both θ and ϕ are periodic by π, the relevant quantities are the circular variances Var[2θ] and Var[2ϕ].) For small µ p the circular variances drop sharply with α, resulting in more orientationally ordered configurations. In addition θ decreases as well, which altogether leads to a situation where the particles obstruct each other's motion less, thus despite the more elongated shape the shear resistance µ c decreases. For larger particle friction however, Var[2θ] barely decreases with α, and the drop in Var[2ϕ] is also very small; these packings remain orientationally rather disordered. The disoriented particles with increasing elongation hinder each other's motion more, leading to a monotonically increasing µ c as a function of α. The insets of the bottom row of Fig. 2 present the effective friction and the normal stress differences as a function of the interparticle friction coefficient µ p for α = 2. As it is expected the effective friction of the system is first gradually growing with increasing µ p , but we find an interesting unexpected breakdown § The periodicity of θ by π must be taken into account when calculating its average. This can be done by calculating the nematic order tensor (3/2)ê •ê − 1/2 (whereê is the unit vector in the particle's axis) and considering its largest eigenvalue's eigendirection, which we do to obtain Fig. 3(a), or by averaging on the complex unit disk: (1/2) arg exp(2iθ) . of this growing trend between 0.4 < µ p < 1. We also see, that N 1 , i.e. the normal stress in the gradient direction (with respect to the flow direction) reaches a peak at µ p = 0.1, where µ c has the highest growth rate, and then decreases back to a small value. At the same time N 2 , the normal stress in the neutral direction with respect to the gradient direction also reaches its minimum. This can be explained by considering the strong change in the particle orientational order, and the distribution of the forces on the particles.
The first and second normal stress differences can be rewritten in terms of average particle level quantities, like orientational angles, their variance, and quantities related to the force distribution on the particle surface. In order to obtain an analytical expression to the stress differences, we have to make approximations, most importantly neglect correlations between the in plane angle θ, the out of plane angle φ, and the eigenvalues of the single particle stress tensor. By doing so we get curves that closely resemble the values obtained by the simulation, including their dependence on α and µ p in most cases [ Fig. 2(e,f) and their insets]. In particular, these calculations recover that N 1 = 0 for α = 1 (with the exception of very small µ p ); that N 1 is increasing and N 2 is decreasing function of α, the dependence of N 1 on µ p is non-monotonic (for example for α = 2), and the decreasing trend of N 2 on small to medium values of µ p . These derivations are technical and we have gathered them in the Supplementary Material accompanying this article.

Volume fraction and coordination number
To complete the rheological description, the quasistatic values of the volume fraction φ c and the coordination number z c are plotted as a function of the aspect ratio on Fig. 4. It is interesting to note, that for the packing fraction the random close packed (RCP) values, obtained as simulation of frictionless particles without shear [3], follow closely our µ p = 0 case for α 2, but deviate for larger aspect ratios: our values start to increase while RCP shows decreasing trend for growing α. For moderate interparticle friction the packing fraction of sheared spherocylinders also decreases for large α, but the curves are still shallow. Our explanation is the following. The most significant difference between our measurements and RCP is that in our case the system is sheared, while for RCP it is not. Shear induces orientational order for elongated particles [25], which gets increasingly pronounced for larger aspect ratios, and orientational order increases the packing fraction. Similar effect (sheared spherocylinders and RCP deviates only for α 2) is observed for the coordination number z c as well. A volume controlled simulation has been performed recently by Nath and coworkers [15]; their jamming density agrees with our measurements using stress controll, which shows the robustness of these results.  The flow of granular materials is -like any typical material -dissipative, quantified by the effective friction coefficient µ. For granular materials dissipation has two sources: collisional loss (parameterized by the coefficient of restitution e) and sliding friction (parameterized by the interparticle friction coefficient µ p ). Since the two mechanisms have different nature, it is worth considering which one is dominant as a function of the parameters [34]. When comparing the two mechanisms, we keep the coefficient of restitution at a fixed intermediate value e = 0.5, and vary µ p and I. Figure 5(a) shows the fraction of the collisional loss (calculated on the contacts during force evaluation) to the total loss as a function of I and µ p . The 50% level set divides the parameter space into three regimes: collisional loss dominated for large I and small µ p , sliding friction dominated for intermediate values of µ p , and a third regime, which we call collisional-frictional , where sliding friction is somewhat suppressed as µ p is so large that the contacts are rarely sliding. Since the transitions are smooth, we call these regions as regimes, instead of phases (which would imply sharp transition). These regimes have been identified earlier for spherical particles [34,35]; here we confirm their presence for elongated grains, and investigate how they are affected by changing the particle elongation. Figure 5(b) shows the dependence of the 50% level set on the aspect ratio. The borderline between sliding and collisional-frictional shifts to larger µ p for more elongated particles: the nearly spherical particles stop sliding at around µ p ≈ 2, while for α = 3 this only happens beyond µ p ≈ 4 or 5. The elongated particles While Ref 34 calls this region "rolling", we consider collisional-frictional more appropriate, as the dominant contribution to dissipation is collisional, and not rolling dissipation. are more entangled by their neighbors, their rotational degrees of freedom are suppressed, resulting in a decrease in collisional dissipation, so sliding friction dominated regime extends further. The border between collisional and sliding friction dominated regime also shows slight aspect ratio dependence, but there it is not monotonic on α. Figure 5(c) and (d) show the dependence of the collisional and the sliding friction contribution to the effective friction as a function of the inertial number. Power law scaling can be observed, with the scaling exponents depending on the dissipation regime. This is especially striking for µ p = 0.01, where increasing I switches regimes from sliding to collisional at around I = 10 −2 .

Velocity fluctuations
The velocity fluctuations also depend strongly on the dissipation regime. On Fig. 6(b) the distribution of the random velocity is shown (v * is the excess velocity on top of the linear velocity profile). The velocity distribution displays a heavy tail in both the collisional and the collisionalfrictional regimes. The same phenomenon is displayed by the width (second moment) of the velocity distribution ( Fig. 6(a)). Not only the value, but also the I-dependence of the velocity fluctuations varies by the dissipation regime. Figure 6(c) shows that v * 2 /(γd) 2 ∼ I −1 in the collisional and collisionalfrictional regimes, while practically independent of I (i.e., ∼ I 0 ) in the sliding friction regime. Below we provide scaling arguments based on the underlying microscopic processes, which also explain the behavior of the autocorrelation functions plotted on Fig. 6(d).
The velocity fluctuations can be understood by the following simple microscopic picture. One must distinguish between the regime of low and large values of µ p on the one hand, for which these fluctuations are large, and the regime for intermediate values of µ p on the other hand, for which they are significantly smaller [ Fig. 6(a)]. In the first highly fluctuating case, the grains move in an intermittent way. They experience short phases of typical duration T = d/ p/ρ during which they are suddenly accelerated by the pressure p to a velocity d/T with respect to their neighbours. The average fluctuating kinetic energy per grain then scales as ρd 3 (d/T ) 2 × f , where f is the fraction of time during which this acceleration phase occurs, i.e. f = Tγ = I. This argument [21,31] gives m v * 2 ∼ d 3 pI. Dividing by the average relative velocity, we obtain v * 2 /(γd) 2 ∼ I −1 , as shown in Fig. 6(c).
By contrast in the second case, the particles' motion does not appear intermittent. The grains move continuously (f = 1), at a time scale that follows the overall shear rate: T ∼ 1/γ. The average fluctuating kinetic energy per grain then scales as d 3 pI 2 , leading to v * 2 /(γd) 2 ∼ I 0 . This behaviour is also consistent with the corresponding curves in Fig. 6(c), which are almost flat for intermediate µ p .
This change in the relevant time scale T is supported by the computation of the autocorrelation function, displayed in Fig. 6(d). For µ p = 10 −1 , the curves lie above those for µ p = 10 −3 and µ p = 10 1 , indicating more persistent grain motion. Also, plotted as functions ofγt the curves for intermediate µ p show a collapse when varying I, while the others rather follow the scale t p/ρ/d =γt/I.

Fluctuations of local shear and rotation
The fluctuations display similar trend on the mesoscopic scale as well. To extract local deformation rates, the simple shear can be written as a sum of pure shear and solid body rotation: so˙ loc and ω loc can be obtained as the symmetric and antisymmetric xy component of the local deformation rate tensor. For homogenous simple shear,˙ = ω =γ/2. Figure 7(a) shows the distribution of the local pure shear rate and the local angular velocity of mesoscopic regions. The local strain rate tensor is obtained by linear regression of the relevant components of the matrix, which projects the particle positions onto their velocity space. The particles are sampled from localized regions of linear extent 1/4 of the largest side of the simulation box. As evident from this panel and Fig. 7(b) (the distribution width of the same quantities), both the local pure shear rate and the local angular velocity have moderately narrow distribution around their mean (which is 1/2 for both the normalized˙ loc /γ and ω loc /γ) in the sliding frictional regime, while the distribution is very wide in the other two regimes. This includes non-negligible fraction of mesoscopic regions, which deform and/or rotate with opposite sign compared to the bulk average. Similarly to the velocity fluctuations, due to the intermittency of grain motion at small and large µ p , the width of these distributions around their average values are an order of magnitude larger than that for intermediate interparticle friction. In the latter case, the grain's velocity is then typically affine, following the global shearing dynamics.

Force and dissipation spatial distribution
It is interesting to consider which part of the surface of the particles experience the strongest confining forces, and whether these areas coincide where most of the dissipation takes place. In the left panel's first column of Fig. 8 the normalized forces on an α = 1.3 particle are displayed: the absolute value of the vectorial average of the forces acting on a surface element is normalized by a typical force based on the confining pressure. In the frictionless regime the forces are concentrated on the cylindrical belt, which holds, somewhat less sharply, in the sliding friction regime as well. This concentration of the forces on the cylindrical belt can be understood by looking at the torque. For spherocyliders the torque on particle i can be expressed as: where o i is the orientation vector of the particle, Z (i) is the set of particles that particle i is in contact with, λ ij ∈ [−d (α − 1) , d (α − 1)] is the (signed) distance between the centre of mass and the normal projection of the contact point to the symmetry axis, and ζ = |F fric ij |/µ p |F rep ij | ∈ [0, 1] is the friction mobilization or plasticity index. In the frictionless case (µ p = 0) only one component remains (sincê c ij ×ĉ ij = 0): Based on this expression the torque of a single contact is zero only at the two singular points at the tips of the particle (o i ×ĉ ij = 0) which are unstable, and at the circle along the centre of the cylinder (λ ij = 0) Figure 8.
The spatial distribution of forces and the dissipation on the surface of the particles. On the left panel the particles are only slightly elongated, aspect ratio is α = 1.3; this is the shape where the effective friction µc is maximal for µp = 0.1. The right panel displays more elongated particles with α = 2, corresponding to the shape used on the insets of Fig. 2. On both panels the first columns show the force distribution on the surface of a particle, non-dimensionalized as | F k |/A k σyy, where the numerator is the absolute value of the vectorial average of the forces on the k-th surface element, and A k is the area of the surface element. The second columns show the dissipation distribution non-dimensionalized as N P k A/σxyγV box A k , where P k is the power dissipated on the surface element, N is the number of particles, A is the surface area and V box is the volume of the box. For this figure α = 1.3 and I = 10 −3 . The normalization in the first columns is fixed across the three images (σyy = 10 −3 ), while in the second columns it is different (σxy, like µc, varies by a factor of ≈ 4 across small to large µp).
which is stable. This shows that frictionless particles, regardless their aspect ratio, prefer to align in a way that their contact points (especially those carrying large forces) are at the middle of the cylindrical belt, as simultaneous force and torque balance is easier achieved when many of the torque contributions are small; this is especially apparent when only a few (eg. two) contacts dominate the force balance on a particle. This argument is valid with good approximation for small nonzero µ p , as the Coulomb cone is still restricted close to the normal direction (µ p ζ 1). In the collisional-frictional regime, however, the areas experiencing the largest forces are the edges of the spherical caps.
The second column on the left panel of Fig. 8 shows the distribution of the dissipation, also normalized by the fraction of the total dissipation projected onto a surface element. There is strong correlation between areas of large forces and large dissipation, with two remarks. First, in the collisional regime there is an extra high-dissipation ring on each of the spherical caps, their location is determined by geometrical constraints on the neighboring particles. Second, the noise in the obtained distribution is largest in the collisional-frictional regime. This can be explained by the fact that the contacts are long lasting in this regime, so fewer separate contacts contribute to the average; and unlike in the sliding friction regime, the contacts are more often localized to their original formation point, and thus contribute only to a single surface element.
On the right panel of Fig. 8 we show the same quantities for a more elongated, α = 2 particle shape. We can observe, like for less elongated particles, that the forces are concentrated on the middle of the cylindrical belt, especially for small to intermediate particle friction. For µ p = 10 −3 , the spherical caps receive also significant incoming forces, ¶ while for µ p = 10 −1 they do not. This can be explained by considering that the elongated particles, regardless of µ p , are oriented by shear to an average direction which is close to the streaming (x) direction. We have seen in Section 3.2.2 that in the collisional regime the velocity fluctuations are large; for oriented elongated particles this fluctuation mainly happens in the axial direction of the particles + , causing frequent collisions on the caps. For the sliding friction regime the velocity fluctuations are much smaller, therefore the collisions on the caps are less frequent and weaker, resulting in less incoming force on those surface elements. This in fact explains why N 1 = σ xx − σ yy is small in the collisional regime and large in the sliding friction regime: up to moderate particle friction the Coulomb cone of the forces is narrow, the direction of the contact forces is close to the surface normal, which for particles oriented with their axis close to the x direction means that forces on the spherical caps contribute mostly to σ xx , while those on the top and bottom of the cylindrical belts contribute to σ yy . In the collisional regime the difference is small, while in the sliding friction regime it is large as the caps receive little incoming force. In the high particle friction collisional-frictional regime both the Coulomb cone is very wide and the orientational ordering is weak, thus the forces can point in almost any direction, resulting in a more isotropic force direction distribution, therefore again small N 1 . This completes the explanation of the non-monotonic behavior of N 1 on µ p for elongated particles [ Fig. 2(e) inset].

Summary and perspectives
We performed DEM simulations to investigate the rheology of a realistic 3-dimensional frictional granular material consisting of elongated particles (spherocylinders). Such systems develop orientational ordering when exposed to shear flow. The degree of this ordering depends on the interparticle friction and particle elongation on a nontrivial manner. Namely, the shear induced orientational ordering is in principle increasing with particle elongation, but the characteristics of collisional and frictional interactions between neighbours (which hinder each others rotation) changes with the interparticle friction coefficient. We measured how key rheological quantities, including effective friction and normal stress differences depend on these two key parameters. We found that the aspect ratio dependence of the effective friction is non-monotonic not only for frictionless particles as we saw earlier, but also for frictional particles up to µ p 0.4, -a range already relevant for every day materials. For higher µ p the effective friction is monotonically increasing. We explained the microscopic origins of both the non-monotonic behavior for small and intermediate µ p and the monotonic one for large µ p . These observations are connected to the fact, that for small friction coefficient the increasing particle aspect ratio leads to stronger ordering and smaller average alignment angle -consequently less obstruction between particles -leading to less resistance against shearing. For particles with large surface friction, however, for increasing aspect ratio the stronger entanglement is not counteracted by the ordering -as it is weaker in this case -leading to monotonically increasing shear resistance. We showed that the collisional, sliding frictional, and collisional-frictional dissipation regimes, which have been identified before for spherical particles, are found also for elongated ones, and observed that the boundary between the sliding frictional and the collisional-frictional regimes moves towards higher µ p for increasing aspect ratio, i.e. increasing grain elongation leads to the expansion of the sliding frictional regime to higher values of the interpaticle friction. We explain this by considering the effect of entanglement on motion: for more elongated particles the larger entanglement leads to the suppression of the rotational motion, shifting the balance from collision-dominated to sliding friction dominated dissipation. We observed that the velocity fluctuations behave differently in the dissipation regimes, and explained its microscopic origins based on the different characteristic time scales of the fluctuations. We measured the spatial distribution of the forces on the particle surfaces, and observed that for small µ p the forces are concentrated on the cylindrical belt. We gave the explanation of this phenomenon based on torque balance on the particles and the nature of the contacts. Finally we expressed the first and second normal stress differences in terms of average particle level quantities, which explained some of the properties of the normal stress differences. One particular non-monotonic behavior, i.e., how N 1 depends on the particle friction µ p for elongated particles, can be explained by the interplay between the amount of velocity fluctuations and the orientation of the particles: large velocity fluctuations (collisional regime) or large Coulomb cone with weak orientational order (collisional-frictional regime) increase the isotropy of the force network, resulting in small values of N 1 ; while its value in the intermediate (sliding friction) regime is high.
This work opens towards the rheology of elongated particles with more complicated shapes or fibers with some flexibility, for which entanglement effects are enhanced [39][40][41]. Also, in the context of active matter, it occurs that swimmers or bacteria can present an elongated shape, which matters for their behaviour [42]. Beyond the properties of clustering of self-propelled rods [43][44][45], the extension of the rheology of active dense granular flows [46] for such long particles remains to be studied. Finally, the limitations of the µ(I) rheology have been recently emphasised, especially in the presence of strong gradients with non-local effects coming into play [19,[47][48][49][50][51][52][53][54][55][56][57], or as a source of ill-posedness in time dependent calculations [58][59][60][61][62]. Because elongated particles can develop secondary flows and consequently build gradients over time, these issues become crucial for the description of their flows [63].