Flow and rheology of frictional elongated grains
Dániel B. Nagy, Philippe Claudin, Tamás Börzsönyi, Ellák Somfai
FFlow and rheology of frictional elongated grains
D´aniel B. Nagy, , Philippe Claudin, Tam´as B¨orzs¨onyi, and Ell´ak Somfai ∗ Institute for Solid State Physics and Optics, Wigner Research Centre for Physics, P.O. Box 49,H-1525 Budapest, Physique et M´ecanique des Milieux H´et´erog`enes, PMMH UMR 7636 CNRS, ESPCI PSL ResearchUniversity, Sorbonne Universit´e, Universit´e de Paris, 10 rue Vauquelin, 75005 Paris, FranceE-mail: [email protected]
Abstract.
The rheology of a 3-dimensional granular system consisting of frictional elongated particleswas investigated by means of discrete element model (DEM) calculations. A homogenous shear flow offrictional spherocyliders was simulated, and a number of rheological quantities were calculated. In theframework of the µ ( I ) rheology, the effective friction was found to be a non-monotonic function of theaspect ratio for interparticle friction coefficient µ p (cid:46) .
4, while it was an increasing function for larger µ p . We reveal the microscopic origin of this peculiar non-monotonic behavior. We show the non-trivialdependence of the velocity fluctuations on the dissipation regime, and trace back the behavior of thenormal stress differences to particle-level quantities.PACS numbers: 45.70.Mg, Published in:New J. Phys. https://doi.org/10.1088/1367-2630/ab91fe (open access) a r X i v : . [ c ond - m a t . s o f t ] N ov low and rheology of frictional elongated grains
1. Introduction
Dense granular flows, with their rich phenomenology, are of great interest for fundamental questions aswell as for their relevance in applied problems. They have been the subject of a large number of studiesin 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 granularsystems have now reached an impressive degree of realism, allowing them for reliable predictions, evenin rather sophisticated configurations [2]. In an effort to go beyond the ideal case of frictionless hardspheres, particles with various shapes [3–15] and surface or bulk conditions [16–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 grainsclose 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. Insuch systems the shear flow induces particle rotation which leads to more intensive collisions betweenneighbouring particles than for spherical grains. The speed of the shear induced rotation depends onthe particle orientation, faster rotation for particles parallel to the shear gradient and slower rotation forparticles pointing in the flow direction, which results in orientational order. Both of these phenomena –collisions due to rotation and orientational ordering – affect the flowing and mechanical properties of thesystem [25–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) withincreasing aspect ratio α for the 3D frictionless case [29], and this non-monotonic tendency was shownto persist even for frictional grains in a 2D system, although only at small values of the interparticlefriction coefficient (up to around µ p = 0 .
15) [30]. In a real world situation the interparticle friction issignificantly larger ( µ p is around 0.3). Moreover, a 3D system is substantially different from the 2D caseas the rotating particles have extra degrees of freedom to evade and reduce the effect of collisions, whichis not possible in 2D. DEM simulations provide both macroscopic and microscopic information aboutthese processes. An effective way is to use a pressure-imposed shear geometry, where the constitutivelaws for elongated grains followed [29, 30] the general framework of the so-called µ ( I ) rheology [4, 31–33].In this description the effective friction µ = τ /P as well as the volume fraction φ are functions of theinertial number I = ˙ γd/ (cid:112) P/ρ , where τ is the shear stress, P the pressure, ˙ γ is the shear rate, d is thegrain size and ρ is the grain density. These rheological functions also depend on the aspect ratio α of thegrains. In addition to the above mentioned non-monotonic behaviour of the effective friction on α , wefound that these flows develop normal stress differences for α > µ c is an increasing function of µ p . Another oneconcerns the exponent of the constitutive law (see Eq. 1 below). In the simulations of Favier et al. [21] itis found that the exponent of the power law term switches from β = 0 . µ p (cid:46) − )to β = 1 in the high friction limit ( µ p (cid:38) − ). Three regimes have been identified and associated withdifferent dissipation mechanisms [34, 35].In this work we extend the µ ( I ) rheology to the case of a 3-dimensional system of frictionalspherocylinders. We show, that in a realistic 3D frictional system the peculiar non-monotonic behaviourof µ ( α ) is observed in a much wider range of interparticle friction (up to around µ p = 0 .
4, thus includingcommon 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 theconstitutive 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 earlierwork [29]. In Section 3 we present and discuss our numerical results, explain a number of observedphenomena, including scaling arguments to understand and interpret the data. We summarize and drawperspectives in Section 4.
2. Setup
We performed numerical simulations of homogenous shear flow in 3D plane Couette geometry, seeFig. 1. The particles were frictional spherocylinders with length-to-diameter aspect ratio α = (cid:96)/ R . Theinterparticle force F ij which particle j exerts on particle i consists of a soft core repulsion and a tangential low and rheology of frictional elongated grains 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 shearedwith shear rate ˙ γ = dv x /dy , and a constant normal stress σ yy = p y is maintained by adjusting box side L y by a feedback loop. component due to friction. The repulsive force points in the direction of the local surface normal ˆ c ij , itsamplitude is proportional to the virtual overlap δ ij between the particles, and it contains a dissipativeterm proportional to the velocity difference v c,ij at the contact point: F rep ij = ( − k δ ij + b v c,ij · ˆ c ij ) ˆ c ij . The prefactor b was determined by requiring a given coefficient of restitution e for binary collisions (weused e = 0 . F fric ij = k ∆ t c,ij , where ∆ t c,ij is the tangential displacement (projectedto the plane perpendicular to ˆ c 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 particlediameter 2 R , 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 sizepolydispersity of 10% (standard deviation to mean ratio in a uniform distribution). While crystallizationwas less critical for frictional particles, we kept the polydispersity fixed for consistency. The rheologicalmeasurements 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 − by a feedbackloop adjusting the L y side of the periodic simulation box.The rest of the simulation details, including the preparation protocol for the initial conditions, aredetailed in Ref. 29.
3. Results
The inertial number dependence of the following rheological parameters are shown on the top row ofFig. 2: effective friction µ = σ xy /p y , normalized first normal stress difference N /p y = ( σ xx − σ yy ) /p y ,and normalized second normal stress difference N /p y = ( σ yy − σ zz ) /p y . The solid curves are fit (in therange 10 − ≤ I ≤ − ) to the following empirical form: µ ( I ) ≈ µ c + µ I β , (1)and similarly for N /p y and N /p y . The values of the exponent β ranged between about 0.4 and slightlymore than 1, where we observed the smallest values for the frictionless case, and the largest values for low and rheology of frictional elongated grains − − − I . . . . . . µ µ p = 0 µ p = 0 . µ p = 0 . µ p = 0 . µ p = 0 . µ p = 0 . µ p = 1 µ p = 101 . . . . . α . . . . µ c − − µ p . . − − − I . . . . . . . N / p y . . . . . α . . . . . . N c / p y − − µ p . . − − − I − . − . − . − . − . N / p y . . . . . α − . − . . N c / p y − − µ p − . . ( a )( d ) ( b )( e ) ( c )( f ) 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 normalstress difference. The colors and symbol shapes indicate different interparticle friction coefficient values,including frictionless ( µ p = 0) and frictional in the range 10 − ≤ µ p ≤
10. The quasistatic limit ( I → ≤ α ≤
3; on the insets alsothe quasistatic limit values are plotted, but for α = 2 and against µ p , see text for details. . . . . . α h θ i [ ◦ ] ( a ) − − µ p . . . . . α . . . . . . V a r [ θ ] ( b ) − − µ p . . . . . . . . α . . . . . . V a r [ ϕ ] ( c ) . . . . . µ p = − − µ p . . Figure 3.
The aspect ratio dependence of (a) the average of the shear alignment angle θ , (b) circularvariance of 2 θ , and (c) circular variance of 2 ϕ . The angle θ denotes the deviation angle of the particleorientation from the streaming direction towards the gradient direction, and ϕ is the deviation towardsthe 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. moderate friction 0 . (cid:46) µ p (cid:46) .
4. Equation (1) enables the reliable extraction of the quasistatic ( I → One remarkable finding is that the quasistatic effective friction coefficient µ c isa non-monotonic function of the aspect ratio for negligible to moderate interparticle friction ( µ p (cid:46) . µ p . For large µ p theeffective friction is a monotonically increasing function of the aspect ratio, at least in the range 1 ≤ α ≤ ‡ As discussed above, the nonmonotonicity of µ c ( α ) has been observed earlier for purely frictionlessspherocylinders [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 . low and rheology of frictional elongated grains Figure 4.
The quasistatic value of the (a) volume fraction φ c and (b) coordination number z c , plottedagainst the aspect ratio α for different interparticle friction coefficient values. The symbols are the sameas in Fig. 2, with the blue full circles with dashed lines corresponding to random close packed simulationswithout particle friction in non-sheared systems [3]. show, that in a 3D system this is observed in an extended friction range, which already includes realisticmaterials. We now explore the microscopic ingredients leading to this behaviour.The explanation can be traced back to the shear induced alignment of elongated particles, initiallydescribed 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 nematicordering, where (cid:104) θ (cid:105) , which we call shear alignment angle, is interestingly non-zero [see Fig. 3(a)] § , while (cid:104) ϕ (cid:105) = 0 by symmetry. (With increasing elongation α the distributions of these angles become typicallynarrower, 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 (cid:104) θ (cid:105) decreases as well, which altogether leads to a situation where the particles obstruct eachother’s motion less, thus despite the more elongated shape the shear resistance µ c decreases. For largerparticle friction however, Var[2 θ ] barely decreases with α , and the drop in Var[2 ϕ ] is also very small; thesepackings remain orientationally rather disordered. The disoriented particles with increasing elongationhinder each other’s motion more, leading to a monotonically increasing µ c as a function of α . The first normal stress difference [Fig. 2(e)] is zero for spherical particlesregardless of friction, and increases monotonically with aspect ratio (except for very small µ p ). The secondnormal 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 differencesas a function of the interparticle friction coefficient µ p for α = 2. As it is expected the effective friction ofthe 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 calculatingthe nematic order tensor (cid:104) (3 / e ◦ ˆ e − / (cid:105) (where ˆ e is the unit vector in the particle’s axis) and considering its largesteigenvalue’s eigendirection, which we do to obtain Fig. 3(a), or by averaging on the complex unit disk: (1 /
2) arg (cid:104) exp(2 iθ ) (cid:105) . low and rheology of frictional elongated grains . < µ p <
1. We also see, that N , i.e. the normal stress in the gradientdirection (with respect to the flow direction) reaches a peak at µ p = 0 .
1, where µ c has the highest growthrate, and then decreases back to a small value. At the same time N , the normal stress in the neutraldirection with respect to the gradient direction also reaches its minimum. This can be explained byconsidering the strong change in the particle orientational order, and the distribution of the forces on theparticles.The first and second normal stress differences can be rewritten in terms of average particle levelquantities, like orientational angles, their variance, and quantities related to the force distribution on theparticle surface. In order to obtain an analytical expression to the stress differences, we have to makeapproximations, most importantly neglect correlations between the in plane angle θ , the out of planeangle φ , and the eigenvalues of the single particle stress tensor. By doing so we get curves that closelyresemble 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 = 0 for α = 1 (with theexception of very small µ p ); that N is increasing and N is decreasing function of α , the dependence of N on µ p is non-monotonic (for example for α = 2), and the decreasing trend of N on small to mediumvalues of µ p . These derivations are technical and we have gathered them in the Supplementary Materialaccompanying this article. To complete the rheological description, the quasistaticvalues of the volume fraction φ c and the coordination number z c are plotted as a function of the aspectratio 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 casefor α (cid:46)
2, but deviate for larger aspect ratios: our values start to increase while RCP shows decreasingtrend for growing α . For moderate interparticle friction the packing fraction of sheared spherocylindersalso decreases for large α , but the curves are still shallow. Our explanation is the following. The mostsignificant difference between our measurements and RCP is that in our case the system is sheared, whilefor RCP it is not. Shear induces orientational order for elongated particles [25], which gets increasinglypronounced for larger aspect ratios, and orientational order increases the packing fraction. Similar effect(sheared spherocylinders and RCP deviates only for α (cid:38)
2) is observed for the coordination number z c as well. A volume controlled simulation has been performed recently by Nath and coworkers [15]; theirjamming density agrees with our measurements using stress controll, which shows the robustness of theseresults. − − µ c o ll ( c ) − − − I − − µ s li d ( d ) − . − . − . − . − . − . . . . ( µ p ) − . − . − . − . − . − . l og ( I ) ( a ) Collisional Sliding frictional Collisional-frictional . . . . . . . ( b ) α = 1 α = 1 . α = 1 . α = 2 α = 3 Figure 5. (a) The fraction of the collisional loss of the total dissipation as a function of the inertialnumber and interparticle friction coefficient, plotted as level sets. (b) Dependence of the 50% level seton the aspect ratio in the range 1 ≤ α ≤
3. (c) Collisional and (d) sliding friction contribution tothe effective friction coefficient as a function of the inertial number for a range of interparticle frictioncoefficient values. The power law exponent depends on the dissipation regime. The aspect ratio onpanels (a), (c) and (d) is α = 2. The symbols on panels (c) and (d) are the same as on Fig. 2. The flow of granular materials is – like any typical material – dissipative, quantified by the effectivefriction coefficient µ . For granular materials dissipation has two sources: collisional loss (parameterized by low and rheology of frictional elongated grains − − − µ p h v ∗ i ( ˙ γ d ) ( a ) Collisional Sliding frictional Collisional-frictional v ∗ / ( ˙ γd )10 − − p ( v ∗ / ( ˙ γ d )) ( b ) − − − I h v ∗ i ( ˙ γ d ) ( c ) . . . . . µ p = 1 110 − − − − ˙ γt − − C ( ˙ γ t ) / C ( ) ( d ) Collisional Sliding frictionalCollisional-frictional Solid: I = 10 − Dashed: I = 10 − . µ p = 10 − µ p = 10 − µ p = 10 Figure 6. (a) Variance of the random velocity v ∗ , non-dimensionalized by the shear rate, as a function of µ p , for I = 10 − . Line colors represent different α values, as in the legend of Fig. 5(b). (b) Distributionof v ∗ for α = 2. (c) Variance of the random velocity v ∗ , as a function of the inertial number. The symbolson panels (b) and (c) represent the same values of µ p as in Figs. 2-4. (d) Normalized autocorrelation ofthe 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 − curve) fall off by ˙ γt , so the curves for different I collapse. However, in the collisional and the collisional-frictional regimes ( µ p = 10 − and 10 ) the dependence is on ˙ γt/I ; the curves would approximatelycollapse when shifted by the difference in I . 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 dominantas a function of the parameters [34]. When comparing the two mechanisms, we keep the coefficient ofrestitution 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 contactsduring force evaluation) to the total loss as a function of I and µ p . The 50% level set divides the parameterspace into three regimes: collisional loss dominated for large I and small µ p , sliding friction dominated forintermediate values of µ p , and a third regime, which we call collisional-frictional (cid:107) , where sliding frictionis somewhat suppressed as µ p is so large that the contacts are rarely sliding. Since the transitions aresmooth, we call these regions as regimes , instead of phases (which would imply sharp transition). Theseregimes have been identified earlier for spherical particles [34, 35]; here we confirm their presence forelongated 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 andcollisional-frictional shifts to larger µ p for more elongated particles: the nearly spherical particles stopsliding at around µ p ≈
2, while for α = 3 this only happens beyond µ p ≈ (cid:107) While Ref 34 calls this region “rolling”, we consider collisional-frictional more appropriate, as the dominant contributionto dissipation is collisional, and not rolling dissipation. low and rheology of frictional elongated grains α .Figure 5(c) and (d) show the dependence of the collisional and the sliding friction contribution tothe effective friction as a function of the inertial number. Power law scaling can be observed, with thescaling exponents depending on the dissipation regime. This is especially striking for µ p = 0 .
01, whereincreasing I switches regimes from sliding to collisional at around I = 10 − . 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 linearvelocity profile). The velocity distribution displays a heavy tail in both the collisional and the collisional-frictional regimes. The same phenomenon is displayed by the width (second moment) of the velocitydistribution (Fig. 6(a)). Not only the value, but also the I -dependence of the velocity fluctuations variesby the dissipation regime. Figure 6(c) shows that (cid:10) v ∗ (cid:11) / ( ˙ γd ) ∼ I − in the collisional and collisional-frictional regimes, while practically independent of I (i.e., ∼ I ) in the sliding friction regime. Below weprovide scaling arguments based on the underlying microscopic processes, which also explain the behaviorof the autocorrelation functions plotted on Fig. 6(d).The velocity fluctuations can be understood by the following simple microscopic picture. One mustdistinguish between the regime of low and large values of µ p on the one hand, for which these fluctuationsare large, and the regime for intermediate values of µ p on the other hand, for which they are significantlysmaller [Fig. 6(a)]. In the first highly fluctuating case, the grains move in an intermittent way. Theyexperience short phases of typical duration T = d/ (cid:112) p/ρ during which they are suddenly accelerated bythe pressure p to a velocity d/T with respect to their neighbours. The average fluctuating kinetic energyper grain then scales as ρd ( d/T ) × f , where f is the fraction of time during which this accelerationphase occurs, i.e. f = T ˙ γ = I . This argument [21, 31] gives m (cid:10) v ∗ (cid:11) ∼ d pI . Dividing by the averagerelative velocity, we obtain (cid:10) v ∗ (cid:11) / ( ˙ γd ) ∼ I − , as shown in Fig. 6(c).By contrast in the second case, the particles’ motion does not appear intermittent. The grainsmove continuously ( f = 1), at a time scale that follows the overall shear rate: T ∼ / ˙ γ . The averagefluctuating kinetic energy per grain then scales as d pI , leading to (cid:10) v ∗ (cid:11) / ( ˙ γd ) ∼ I . This behaviour isalso 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 autocorrelationfunction, displayed in Fig. 6(d). For µ p = 10 − , the curves lie above those for µ p = 10 − and µ p = 10 ,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 (cid:112) p/ρ/d = ˙ γt/I . The fluctuations display similar trend on the mesoscopicscale as well. To extract local deformation rates, the simple shear can be written as a sum of pure shearand solid body rotation:˙ γγγ = γ
00 0 00 0 0 = (cid:15) (cid:15) + ω − ω , so ˙ (cid:15) loc and ω loc can be obtained as the symmetric and antisymmetric xy component of the localdeformation rate tensor. For homogenous simple shear, ˙ (cid:15) = ω = ˙ γ/
2. Figure 7(a) shows the distributionof the local pure shear rate and the local angular velocity of mesoscopic regions. The local strain ratetensor is obtained by linear regression of the relevant components of the matrix, which projects theparticle positions onto their velocity space. The particles are sampled from localized regions of linearextent 1 / / (cid:15) loc / ˙ γ and ω loc / ˙ γ ) in the sliding frictional regime, while the distribution is very wide in the other two regimes. Thisincludes non-negligible fraction of mesoscopic regions, which deform and/or rotate with opposite signcompared 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 low and rheology of frictional elongated grains − . − . − . . . . . . ˙ (cid:15) loc ˙ γ , ω loc ˙ γ . . . . . . p ( ˙ (cid:15) l o c ) , p ( ω l o c ) ( a ) µ p = 0 .
001 : ˙ (cid:15) loc µ p = 0 .
001 : ω loc µ p = 0 . (cid:15) loc µ p = 0 . ω loc µ p = 1 : ˙ (cid:15) loc µ p = 1 : ω loc − − − µ p σ ˙ (cid:15) l o c ˙ γ , σ ω l o c ˙ γ ( b ) Collisional Sliding frictional Collisional-frictional α = 1 : ˙ (cid:15) loc α = 1 : ω loc α = 1 . (cid:15) loc α = 1 . ω loc α = 2 : ˙ (cid:15) loc α = 2 : ω loc α = 3 : ˙ (cid:15) loc α = 3 : ω loc Figure 7. (a) Distribution of the local pure shear rate ˙ (cid:15) loc (solid symbols) and local angular velocity ω loc (open symbols) on mesoscopic scales. (b) The distribution width (standard deviation) of themesoscopic pure shear rate and angular velocity as a function of µ p . Both quantities have relativelynarrow distribution in the sliding frictional regime, and have very wide distribution in the other tworegimes. that for intermediate interparticle friction. In the latter case, the grain’s velocity is then typically affine,following the global shearing dynamics. It is interesting to consider which part of the surface ofthe particles experience the strongest confining forces, and whether these areas coincide where most ofthe dissipation takes place. In the left panel’s first column of Fig. 8 the normalized forces on an α = 1 . i can be expressed as: τ i = (cid:88) j ∈ Z ( i ) (cid:18) λ ij o i + d (1 − δ ij )ˆ c ij (cid:19) × | F rep ij | (cid:0) ˆ c ij + µ p ζ ˆ t c ,ij (cid:1) , (2)where o i is the orientation vector of the particle, Z ( i ) is the set of particles that particle i is in contactwith, λ ij ∈ [ − d ( α − , d ( α − ζ = | F fric ij | /µ p | F rep ij | ∈ [0 ,
1] is the frictionmobilization or plasticity index. In the frictionless case ( µ p = 0) only one component remains (sinceˆ c ij × ˆ c ij = 0): τ i = (cid:88) j ∈ Z ( i ) λ ij | F rep ij | o i × ˆ c ij . (3)Based on this expression the torque of a single contact is zero only at the two singular points at the tips ofthe particle ( o i × ˆ c ij = 0) which are unstable, and at the circle along the centre of the cylinder ( λ ij = 0) low and rheology of frictional elongated grains Figure 8.
The spatial distribution of forces and the dissipation on the surface of the particles. On the left panelthe particles are only slightly elongated, aspect ratio is α = 1 .
3; this is the shape where the effectivefriction µ 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 forcedistribution on the surface of a particle, non-dimensionalized as | (cid:104) F k (cid:105) | /A k σ yy , where the numeratoris the absolute value of the vectorial average of the forces on the k -th surface element, and A k is thearea of the surface element. The second columns show the dissipation distribution non-dimensionalizedas N (cid:104) P k (cid:105) A/σ xy ˙ γV box A k , where P k is the power dissipated on the surface element, N is the number ofparticles, A is the surface area and V box is the volume of the box. For this figure α = 1 . I = 10 − .The normalization in the first columns is fixed across the three images ( σ yy = 10 − ), while in the secondcolumns it is different ( σ xy , like µ c , varies by a factor of ≈ µ p ). which is stable. This shows that frictionless particles, regardless their aspect ratio, prefer to align in away that their contact points (especially those carrying large forces) are at the middle of the cylindricalbelt, as simultaneous force and torque balance is easier achieved when many of the torque contributionsare small; this is especially apparent when only a few (eg. two) contacts dominate the force balance on aparticle. This argument is valid with good approximation for small nonzero µ p , as the Coulomb cone isstill restricted close to the normal direction ( µ p ζ (cid:28) α = 2 particle shape.We can observe, like for less elongated particles, that the forces are concentrated on the middle of thecylindrical belt, especially for small to intermediate particle friction. For µ p = 10 − , the spherical capsreceive also significant incoming forces, ¶ while for µ p = 10 − they do not. This can be explained by ¶ Similar increased concentration of contacts on the central region at small aspect ratios, and on the tips at large aspectratios are also observed on frictionless 2D elongated particles [14]. low and rheology of frictional elongated grains µ p , are oriented by shear to an average directionwhich is close to the streaming ( x ) direction. We have seen in Section 3.2.2 that in the collisional regimethe velocity fluctuations are large; for oriented elongated particles this fluctuation mainly happens in theaxial direction of the particles + , causing frequent collisions on the caps. For the sliding friction regime thevelocity 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 = σ xx − σ yy issmall in the collisional regime and large in the sliding friction regime: up to moderate particle friction theCoulomb 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 sphericalcaps 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 asthe caps receive little incoming force. In the high particle friction collisional-frictional regime both theCoulomb cone is very wide and the orientational ordering is weak, thus the forces can point in almostany direction, resulting in a more isotropic force direction distribution, therefore again small N . Thiscompletes the explanation of the non-monotonic behavior of N on µ p for elongated particles [Fig. 2(e)inset].
4. Summary and perspectives
We performed DEM simulations to investigate the rheology of a realistic 3-dimensional frictional granularmaterial consisting of elongated particles (spherocylinders). Such systems develop orientational orderingwhen exposed to shear flow. The degree of this ordering depends on the interparticle friction and particleelongation on a nontrivial manner. Namely, the shear induced orientational ordering is in principleincreasing with particle elongation, but the characteristics of collisional and frictional interactions betweenneighbours (which hinder each others rotation) changes with the interparticle friction coefficient. Wemeasured how key rheological quantities, including effective friction and normal stress differences dependon these two key parameters. We found that the aspect ratio dependence of the effective friction isnon-monotonic not only for frictionless particles as we saw earlier, but also for frictional particles upto µ p (cid:46) .
4, – a range already relevant for every day materials. For higher µ p the effective friction ismonotonically increasing. We explained the microscopic origins of both the non-monotonic behavior forsmall and intermediate µ p and the monotonic one for large µ p . These observations are connected tothe fact, that for small friction coefficient the increasing particle aspect ratio leads to stronger orderingand smaller average alignment angle – consequently less obstruction between particles – leading to lessresistance against shearing. For particles with large surface friction, however, for increasing aspect ratiothe stronger entanglement is not counteracted by the ordering – as it is weaker in this case – leadingto monotonically increasing shear resistance. We showed that the collisional, sliding frictional, andcollisional-frictional dissipation regimes, which have been identified before for spherical particles, arefound also for elongated ones, and observed that the boundary between the sliding frictional and thecollisional-frictional regimes moves towards higher µ p for increasing aspect ratio, i.e. increasing grainelongation leads to the expansion of the sliding frictional regime to higher values of the interpaticlefriction. We explain this by considering the effect of entanglement on motion: for more elongated particlesthe larger entanglement leads to the suppression of the rotational motion, shifting the balance fromcollision-dominated to sliding friction dominated dissipation. We observed that the velocity fluctuationsbehave differently in the dissipation regimes, and explained its microscopic origins based on the differentcharacteristic time scales of the fluctuations. We measured the spatial distribution of the forces on theparticle surfaces, and observed that for small µ p the forces are concentrated on the cylindrical belt. Wegave the explanation of this phenomenon based on torque balance on the particles and the nature of thecontacts. Finally we expressed the first and second normal stress differences in terms of average particlelevel quantities, which explained some of the properties of the normal stress differences. One particularnon-monotonic behavior, i.e., how N depends on the particle friction µ p for elongated particles, canbe explained by the interplay between the amount of velocity fluctuations and the orientation of theparticles: large velocity fluctuations (collisional regime) or large Coulomb cone with weak orientationalorder (collisional-frictional regime) increase the isotropy of the force network, resulting in small values of N ; 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 fiberswith some flexibility, for which entanglement effects are enhanced [39–41]. 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. low and rheology of frictional elongated grains µ ( I ) rheology have been recently emphasised, especially in the presence of stronggradients with non-local effects coming into play [19, 47–57], or as a source of ill-posedness in timedependent calculations [58–62]. Because elongated particles can develop secondary flows and consequentlybuild gradients over time, these issues become crucial for the description of their flows [63]. Acknowledgements
This work was supported in part by the Hungarian National Research, Development and InnovationOffice NKFIH under grant OTKA K 116036. We acknowledge funding from the CNRS with the PICSGrant No. 08187 ‘Flow properties of granular materials consisting of elongated grains’ (2019-2021), andthe Hungarian Academy of Sciences (Grant No. NKM-102/2019). DBN was supported by the studentscholarship of the French Embassy in Budapest and Campus France. We acknowledge KIF ¨U for awardingus access to computational resources based in Hungary at Debrecen.
References [1] Andreotti B, Forterre Y and Pouliquen O 2013
Granular media, between fluid and solid (Cambridge University Press)ISBN 978-1-107-03479-2[2] Radjai F and Dubois F (eds) 2011
Discrete numerical modeling of granular materials (Wiley-ISTE) ISBN 978-1-84821-260-2[3] Donev A, Cisse I, Sachs D, Variano E, Stillinger F H, Connelly R, Torquato S and Chaikin P M 2004
Science
Phys. Rev. E J. Fluid Mech.
Soft Matter Phys. Rev. E Soft Matter EPL Soft Matter Soft Matter Phys. Rev. E (1) 010202[13] Mandal S and Khakhar D V 2016 Phys. Fluids Phys. Rev. E (1) 012905[15] Nath T and Heussinger C 2019 Eur. Phys. J. E Europhys. Lett. J. Fluid Mech.
Phys. Rev. E (2) 021305[19] Kamrin K and Koval G 2014 Comput. Part. Mech. New J. Phys. Phys. Rev. Fluids (10) 102301[22] Koivisto J, Korhonen M, Alava M, Ortiz C P, Durian D J and Puisto A 2017 Soft Matter New J. Phys. Soft Matter Phys. Rev. Lett.
Soft Matter Soft Matter J. Fluid Mech.
R5[29] Nagy D B, Claudin P, B¨orzs¨onyi T and Somfai E 2017
Phys. Rev. E (6) 062903[30] Trulsson M 2018 J. Fluid Mech.
Eur. Phys. J. E Nature
Phys. Rev. E Phys. Rev. E (1) 012904[35] Trulsson M, DeGiuli E and Wyart M 2017 Phys. Rev. E (1) 012605[36] Somfai E, Nagy D B, Claudin P, Favier A, K´alm´an D and B¨orzs¨onyi T 2017 EPJ Web Conf.
Phys. Rev. Fluids (7) 074301[38] B¨orzs¨onyi T, Szab´o B, Wegner S, Harth K, T¨or¨ok J, Somfai E, Bien T and Stannarius R 2012 Phys. Rev. E Phys. Rev. Lett. (10) 108004[40] Bertails-Descoubes F, Cadoux F, Daviet G and Acary V 2011 ACM Trans. Graph. Fluids, Colloids and Soft Materials: An Introduction toSoft Matter Physics ed Nieves A F and Puertas A M (John Wiley & Sons, Inc.) chap 17, pp 341–354[42] Ilkanaiv B, Kearns D B, Ariel G and Be’er A 2017
Phys. Rev. Lett.
Phys. Rev. E Phys. Rev. Lett. low and rheology of frictional elongated grains [45] Yang Y, Marceau V and Gompper G 2010 Phys. Rev. E EPL
Philos. Trans. R. Soc. A
Phys. Rev. Lett. (17) 178301[49] Henann D L and Kamrin K 2013
Proc. Natl. Acad. Sci. U. S. A.
Phys. Rev. Lett. (23) 238301[51] Henann D L and Kamrin K 2014
Phys. Rev. Lett. (17) 178001[52] Kamrin K and Henann D L 2015
Soft Matter EPL
Eur. Phys. J. E J. Fluid Mech.
EPL
Soft Matter J. Fluid Mech.
J. Fluid Mech.
Proc. R. Soc. A
J. Fluid Mech.
Phys. Fluids Soft Matter upplementary Material for Flow and rheology of frictional elongated grains
Dániel B. Nagy, a , b Philippe Claudin, b Tamás Börzsönyi, a and Ellák Somfai ∗ aa Institute for Solid State Physics and Optics, Wigner Research Center for Physics, P.O. Box 49, H-1525 Budapest, Hungary.E-mail: [email protected] b Physique et Mécanique des Milieux Hétérogènes, PMMH UMR 7636 CNRS, ESPCI Paris, PSL University, Sorbonne Université,Université de Paris, 10 rue Vauquelin, 75005 Paris, France
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 dif-ferences
We calculate the total stress on the system using the ex-pression: σσσ = − V box ∑ i ∑ j ∈ Z ( i ) r c , i j ◦ F i j + m v ∗ i ◦ v ∗ i = φ h σσσ i i , (1)where r c , i j is the vector between the center of mass of par-ticle i and the contact point, v ∗ is the random velocity (ex-cess velocity on top of the linear velocity profile), φ is thevolume fraction of the packing, and σσσ i = − V part ∑ j ∈ Z ( i ) r c , i j ◦ F i j + m v ∗ i ◦ v ∗ i , (2)is the single particle stress tensor where V part is the meanvolume of a particle. For the normal stress differencesonly the symmetric part matters, so we can define σσσ ∗ i = (cid:0) σσσ i + σσσ Ti (cid:1) , that has orthonormal eigenvectors. From thedistribution of forces seen on Fig. 7 of our paper, on av-erage the eigenvectors for the smallest eigenvalue w i andfor the largest eigenvalue q i are in the plane defined by o i (particle axis unit vector) and ˆ y (unit vector in the y i.e.,gradient direction). And the eigenvector u i correspondingto the intermediate eigenvalue will be close to parallel to ˆ y × o i . Introducing the average angle η between the parti-cle axis o and the eigenvector w : w i = cos η o i − sin η ˆ y ∗ i q i = sin η o i + cos η ˆ y ∗ i , (3) with ˆ y ∗ i k o i × ( ˆ y × o i ) . Using the mathematical properties oforthonormal eigenvectors: σσσ ∗ i = σ w , i + ( σ q , i − σ w , i ) q i ◦ q i + ( σ u , i − σ w , i ) u i ◦ u i . (4)We introduce spherical coordinates for the particle axisvector: 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 par-ticle orientation from the streaming direction towards thegradient direction, and ϕ (out of plane angle) is the devi-ation towards the neutral direction. After trigonometricalcalculations and neglecting the correlation between the av-eraged quantities: N ≈ φ h ( + R ϕ ) R θ h ∆ σ qw i sin 2 ( π − h θ i + h η i ) − ( − R ϕ ) ( h ∆ σ qw − ∆ σ qu i ) i , (5)if h ϕ i ≈ , and R ψ = |h e i ψ i| ∼ − Var [ ψ ] defines the circularvariance of any angle ψ . Very similarly: N ≈ φ h − ( − R ϕ ) R θ h ∆ σ qw i sin 2 ( π − h θ i + h η i )+ ( + R ϕ ) ( h ∆ σ qw − ∆ σ qu i ) i . (6)The quantities in these expressions were measured numer-ically, and plotted together with the resulting normal stressdifferences on Figs. S1 and S2.In the calculation above we made several approxima-tions, most notably neglected the correlations. It is re-markable, that the resulting normal stress differences arein most cases closely resemble the true (numerically mea-sured) values. Thus for the first normal stress difference S1 ig. S1(a) is very similar to Fig. 2(e) for the α -dependence,with the exception of very small µ p . Also the blue ( α = )curve of Fig. S2(a) is similar to Fig. 2(e) inset for the µ p de-pendence: high in the sliding friction regime and low in theother two regimes. Similarly for the second normal stressdifference Fig. S1(b) is closely resembles Fig. 2(f), with aslight shift in values except again for very small µ p . Andagain the α = curve (blue) on Fig. S2(b) resembles theinset of Fig. 2(f): close to zero in the collisional regime,decreases in the sliding friction regime and stays roughlyconstant in the collisional-frictional regime. B – Supplied videos The three supplied videos show a selected particle (or-ange) and its neighborhood (translucent blue particles)during shear. In all three cases the total number of parti-cles were (most of them not shown), with aspect ratio α = , sheared with a moderate shear rate (inertial number I = − . The simulation boxes are oriented such that wesee the x - y (streamlines and gradient direction) plane. Thethree videos differ in the value of the particle friction coeffi-cient: µ p = . , . , and , which correspond to the col-lisional, sliding frictional and collisional-frictional regimesrespecively.As expected [eg., see Fig. 6(a-b)], the motion of the par-ticles were more noisy in the collisional and collisional-frictional regimes, and less noisy in the sliding frictionalregime. The supplied videos can be downloaded from the journal website: https://doi.org/10.1088/1367-2630/ab91fehttps://doi.org/10.1088/1367-2630/ab91fe