Dániel B Nagy

Philippe Claudin

Tamás Börzsönyi

Ellák Somfai
email: somfai.ellak@wigner.hu

Flow and rheology of frictional elongated grains

Keywords: numbers: 45, 70, Mg,

published or not. The documents may come   L'archive ouverte pluridisciplinaire

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 [START_REF] Andreotti | Granular media, between fluid and solid[END_REF]. 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 [START_REF]Discrete numerical modeling of granular materials[END_REF]. In an effort to go beyond the ideal case of frictionless hard spheres, particles with various shapes [START_REF] Donev | [END_REF][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 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:

F rep ij = (-k δ ij + b v c,ij • ĉij ) ĉij .
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.

Results

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:

µ(I) ≈ µ c + µ 1 I β , (1)
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

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 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.

(d) clearly

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 α.

Normal stress differences

The first normal stress difference [Fig. 2(e)] is zero for spherical particles regardless of friction, and increases monotonically with aspect ratio (except for very small µ p ). The second normal stress difference [Fig. 2(f)] is nonzero even for spherical particles.

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 [START_REF] Donev | [END_REF], 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 curve) fall off by γt, so the curves for different I collapse. However, in the collisional and the collisionalfrictional regimes (µp = 10 -3 and 10 1 ) the dependence is on γt/I; the curves would approximately collapse when shifted by the difference in I.

Dissipation regimes

α = 1 α = 1.1 α = 1.5 α = 2 α = 3
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.

Regime diagram 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 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:

γ γ γ =   0 γ 0 0 0 0 0 0 0   =   0 ˙ 0 ˙ 0 0 0 0 0   +   0 ω 0 -ω 0 0 0 0 0   ,
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:

τ i = j∈Z (i) λ ij o i + d(1 - 1 2 δ ij )ĉ ij × |F rep ij | ĉij + µ p ζ tc,ij , (2)
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 (since ĉij × ĉij = 0):

τ i = j∈Z (i) λ ij |F rep ij |o i × ĉij . (3)
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) 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][START_REF] Gravish | Entangled granular media Fluids, Colloids and Soft Materials: An Introduction to Soft Matter Physics[END_REF]. Also, in the context of active + We confirmed that the velocity fluctuations are indeed larger in the x direction compared to the y and z directions.

matter, it occurs that swimmers or bacteria can present an elongated shape, which matters for their behaviour [START_REF] Ilkanaiv | [END_REF]. 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].

Figure 1 .

 1 Figure 1. Setup of the simulation.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.

Figure 2 .

 2 Figure 2. The inertial number dependence of rheological parameters are shown on top row for α = 2: (a) effective friction, (b) normalized first normal stress difference, and (c) normalized second normal stress difference. The colors and symbol shapes indicate different interparticle friction coefficient values, including frictionless (µp = 0) and frictional in the range 10 -3 ≤ µp ≤ 10. The quasistatic limit (I → 0) of the same quantities are plotted in the bottom row for a range of shapes 1 ≤ α ≤ 3; on the insets also the quasistatic limit values are plotted, but for α = 2 and against µp, see text for details.

Figure 3 .

 3 Figure3. The aspect ratio dependence of (a) the average of the shear alignment angle θ, (b) circular variance of 2θ, and (c) circular variance of 2ϕ. The angle θ denotes the deviation angle of the particle orientation from the streaming direction towards the gradient direction, and ϕ is the deviation towards the neutral direction. Circular variance of 1 corresponds to completely uniform angular distribution, while 0 corresponds to zero-width (Dirac-delta) distribution. Symbols and colors represent different µp. The insets, like on Fig.2, show the dependence on µp for aspect ratio value α = 2.

Figure 4 .

 4 Figure 4. The quasistatic value of the (a) volume fraction φc and (b) coordination number zc, plotted against the aspect ratio α for different interparticle friction coefficient values. The symbols are the same as in Fig. 2, with the blue full circles with dashed lines corresponding to random close packed simulations without particle friction in non-sheared systems [3].

Figure 5 .

 5 Figure 5. (a) The fraction of the collisional loss of the total dissipation as a function of the inertial number and interparticle friction coefficient, plotted as level sets. (b) Dependence of the 50% level set on the aspect ratio in the range 1 ≤ α ≤ 3. (c) Collisional and (d) sliding friction contribution to the effective friction coefficient as a function of the inertial number for a range of interparticle friction coefficient values. The power law exponent depends on the dissipation regime. The aspect ratio on panels (a), (c) and (d) is α = 2. The symbols on panels (c) and (d) are the same as on Fig. 2.

5 µ p = 10 - 3 µ p = 10 - 1 µ p = 10 1 Figure 6 .

 510310116 Figure 6. (a) Variance of the random velocity v * , non-dimensionalized by the shear rate, as a function of µp, for I = 10 -3 . Line colors represent different α values, as in the legend of Fig. (b) Distribution of v * for α = 2. (c) Variance of the random velocity v * , as a function of the inertial number. The symbols on panels (b) and (c) represent the same values of µp as in Figs. 2-4. (d) Normalized autocorrelation of the random velocity for µp values in the three regimes (colours) and two different inertial number values (continuous and dashed lines). The autocorrelation function in the sliding friction regime (µp = 10 -1curve) fall off by γt, so the curves for different I collapse. However, in the collisional and the collisionalfrictional regimes (µp = 10 -3 and 10 1 ) the dependence is on γt/I; the curves would approximately collapse when shifted by the difference in I.

  ), which are almost flat for intermediate µ p .

Figure 7 .

 7 Figure 7. (a) Distribution of the local pure shear rate ˙ loc (solid symbols) and local angular velocity ω loc (open symbols) on mesoscopic scales. (b) The distribution width (standard deviation) of the mesoscopic pure shear rate and angular velocity as a function of µp. Both quantities have relatively narrow distribution in the sliding frictional regime, and have very wide distribution in the other two regimes.

Figure 8 .

 8 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).

‡ 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.

¶ Similar increased concentration of contacts on the central region at small aspect ratios, and on the tips at large aspect ratios are also observed on frictionless 2D elongated particles[14].

The supplied videos can be downloaded from the journal website: https://doi.org/10.1088/1367-2630/ab91fe S3

Acknowledgements

This work was supported in part by the Hungarian National Research, Development and Innovation Office NKFIH under grant OTKA K 116036. We acknowledge funding from the CNRS with the PICS Grant No. 08187 'Flow properties of granular materials consisting of elongated grains' (2019-2021), and the Hungarian Academy of Sciences (Grant No. NKM-102/2019). DBN was supported by the student scholarship of the French Embassy in Budapest and Campus France. We acknowledge KIF Ü for awarding us access to computational resources based in Hungary at Debrecen.

Supplementary Material for

Flow and rheology of frictional elongated grains

Below we (A) show how the first and second normal stress differences can be expressed in terms of simpler quantities, and (B) describe the three supplied videos.

A -First and second normal stress differences

We calculate the total stress on the system using the expression:

where r c,i j is the vector between the center of mass of particle i and the contact point, v * is the random velocity (excess velocity on top of the linear velocity profile), φ is the volume fraction of the packing, and

is the single particle stress tensor where V part is the mean volume of a particle. For the normal stress differences only the symmetric part matters, so we can define σ σ σ * i = 1 2 σ σ σ i + σ σ σ T i , that has orthonormal eigenvectors. From the distribution of forces seen on Fig. 7 of our paper, on average the eigenvectors for the smallest eigenvalue w i and for the largest eigenvalue q i are in the plane defined by o i (particle axis unit vector) and ŷ (unit vector in the y i.e., gradient direction). And the eigenvector u i corresponding to the intermediate eigenvalue will be close to parallel to ŷ × o i . Introducing the average angle η between the particle axis o and the eigenvector w:

Using the mathematical properties of orthonormal eigenvectors:

We introduce spherical coordinates for the particle axis vector: o x,i = cos θ i cos ϕ i , o y,i = sin θ i and o z,i = cos θ i sin ϕ i , where θ (in plane angle) is the deviation angle of the particle orientation from the streaming direction towards the gradient direction, and ϕ (out of plane angle) is the deviation towards the neutral direction. After trigonometrical calculations and neglecting the correlation between the averaged quantities:

if ϕ ≈ 0, and R ψ = | e iψ | ∼ 1 -Var[ψ] defines the circular variance of any angle ψ. Very similarly:

The quantities in these expressions were measured numerically, and plotted together with the resulting normal stress differences on Figs. S1 and S2.

In the calculation above we made several approximations, most notably neglected the correlations. It is remarkable, that the resulting normal stress differences are in most cases closely resemble the true (numerically measured) values. Thus for the first normal stress difference

B -Supplied videos 1

The three supplied videos show a selected particle (orange) and its neighborhood (translucent blue particles) during shear. In all three cases the total number of particles were 1000 (most of them not shown), with aspect ratio α = 2, sheared with a moderate shear rate (inertial number I = 10 -2 . The simulation boxes are oriented such that we see the x-y (streamlines and gradient direction) plane. The three videos differ in the value of the particle friction coefficient: µ p = 0.001, 0.1, and 10, which correspond to the collisional, sliding frictional and collisional-frictional regimes respecively.

As expected [eg., see Fig. 6(a-b)], the motion of the particles were more noisy in the collisional and collisionalfrictional regimes, and less noisy in the sliding frictional regime.