Emilien Azéma
email: emilien.azema@univ-montp2.fr\\franck.radjai@univ-montp2.fr

Farhang Radjaï

Force chains and contact network topology in packings of elongated particles

Keywords: numbers: 45.70.-n, 83.80.Fg, 61.43.-j

By means of contact dynamic simulations, we investigate the contact network topology and force chains in two-dimensional packings of elongated particles modeled by rounded-cap rectangles. The morphology of large packings of elongated particles in quasistatic equilibrium is complex due to the combined effects of local nematic ordering of the particles and orientations of contacts between particles. We show that particle elongation affects force distributions and force/fabric anisotropy via various local structures allowed by steric exclusions and the requirement of force balance. As a result, the force distributions become increasingly broader as particles become more elongated. Interestingly, the weak force network transforms from a passive stabilizing agent with respect to strong force chains to an active force-transmitting network for the whole system. The strongest force chains are carried by side/side contacts oriented along the principal stress direction.

I. INTRODUCTION

Most remarkable properties of granular materials are closely related to their specific disorder induced essentially by steric exclusions and the force balance condition for each particle. The broad and strongly inhomogeneous distribution of contact forces, as a hallmark of granular disorder, has been a subject of extensive investigation [1][2][3][4][5][6][7][8][9][10]. In spite of particle mobility and disorder, granular materials exhibit a finite shear strength due to a genuine anisotropic two-phase organization of the contact network involving strong force chains propped by weak forces [11].

The robustness of these micro-structural features with respect to particle geometry and interactions has been addressed only recently by discrete-element numerical simulations. For example, it is found that in highly polydisperse systems the force chains are mainly captured by large particles so that the shear strength of a noncohesive granular material is practically independent of particle size distribution [12]. As another important example, a parametric study shows that when the particles interact by both sliding friction and high rolling resistance at their contacts, the nature of the weak network is affected by the formation of columnar structures which do not need to be propped by a particular class of weak contacts [13].

The particle shape is another major characteristic of granular material. Most applications of real granular materials involve some degree of deviation with respect to simple circular or spherical shapes often used in simulations by the discrete-element method. While the numerical treatment of large packings of complex particle shapes was until very recently out of reach due to demanding computational resources, there is presently considerable scope for the numerical investigation of complex granular packings. This is not only due to the enhanced computer power and memory but also because during more than two decades of research in this field, many properties of granular media have been investigated for model packings composed of circular and spherical particle shapes. Such properties provide thus a rich guideline for the analysis of specific behaviors arising from particle geometry.

Systematic studies of particle shape dependence in granular materials have been recently reported for polygonal/polyhedral [14][15][16][17][18], elliptical/ellipsoidal [19][20][21][22] and non-convex shapes [23]. The force chains are found to be reinforced in packings of polygonal and polyhedral particles leading to enhanced shear strength [15,22,24]. The effect of shape elongation was investigated for packings of rectangular-shaped particles deposited under gravity [25]. The preparation under gravity has strong influence on the particle orientations and thus on the force distributions. On the other hand, a systematic study of the shear behavior of 2D packings of rounded-cap rectangles (RCR) under homogeneous boundary conditions indicates that the shear strength increases with elongation whereas the packing fraction varies unmonotonically [22], as also found for packings of ellipsoidal shapes [20,21]. In all reported cases, the networks resulting from various shapes appear to be highly complex and hardly amenable to simple statistical modeling.

In this paper, we use contact dynamics simulations to investigate the contact and force networks in sheared granular packings of RCR particles with increasing aspect ratio in 2D. We focus more specifically on the organization of the contact force network in correlation with the fabric anisotropy described in terms of branch vectors joining particle centers. Our data reveal a bimodal force network as in disk packings but with qualitatively different roles of fabric and force anisotropies. This behavior involves a short-range nematic ordering of the particles with side/side contacts that capture stronger force chains. On the other hand, the friction mobilization is shown to be anisotropic and it plays a major role for the stability of elongated particles.

In the following, we first briefly describe the numerical procedures, which are essentially the same as those reported in [22]. Then, we analyze the branch vectors and their correlations with the contact forces. Finally, we present a detailed analysis of the partial stresses and fabric anisotropies sustained by force sub-networks. We conclude with salient results of this work its possible prospectives.

II. MODEL DESCRIPTION AND NUMERICAL SIMULATIONS

The simulations were carried out by means of the contact dynamics (CD) method with irregular polyhedral particles. The CD method is a discrete element approach for the simulation of nonsmooth granular dynamics with contact laws expressing mutual exclusion and dry friction between particles without elastic or viscous regularization [10,[26][START_REF] Radjai | Physics of dry granular media[END_REF][START_REF] Jean | [END_REF][START_REF] Moreau | Novel approaches in civil engineering[END_REF][START_REF] Dubois | Analysis and Simulation of Contact Problems[END_REF][START_REF] Radjaï | [END_REF][START_REF]Discrete Numerical Modeling of Granular Materials[END_REF]. Hence, this method is particularly adapted for the simulation of perfectly rigid particles. Nonsmoothness refers to various degrees of discontinuity in velocities arising in a system of rigid particles. In this method, the equations of motion for each particle are formulated as differential inclusions in which velocity jumps replace accelerations [26]. The unilateral contact interactions and Coulomb friction law are treated as complementarity relations or set-valued contact laws. The time-stepping scheme is implicit but requires explicit determination of the contact network. Due to implicit time integration, inherent in the CD method, this scheme is unconditionally stable.

At a given step of evolution, all kinematic constraints implied by lasting contacts and the possible rolling of some particles over others are simultaneously taken into account, together with the equations of dynamics, in order to determine all velocities and contact forces in the system. This problem is solved by an iterative process pertaining to the non-linear Gauss-Seidel method which consists of solving a single contact problem, with other contact forces being treated as known, and iteratively updating the forces and velocities until a convergence criterion is fulfilled. The iterations in a time step are stopped when the calculated contact forces are stable with respect to the update procedure. To check convergence we thus use the relative variation of the mean contact force between two successive iterations. We require this relative variation to be below a given value which sets the precision of the calculation. In this process, no distinction is made between smooth evolution of a system of rigid particles during one time step and nonsmooth evolutions in time due to collisions or dry friction effects. The uniqueness of the solution at each time step is not guaranteed by CD method for perfectly rigid particles. However, by initializing each step of calculation with the forces calculated in the preceding step, the set of accessible solutions shrinks to fluctuations which are basically below the nu- merical resolution. In this way, the solution remains close to the present state of forces.

For our simulations, we used the LMGC90 which is a multipurpose software developed in Montpellier, capable of modeling a collection of deformable or undeformable particles of various shapes (spherical, polyhedral, or polygonal) by different algorithms [START_REF] Dubois | Analysis and Simulation of Contact Problems[END_REF][START_REF]Discrete Numerical Modeling of Granular Materials[END_REF].

A. Simulation of RCR particles

We model the RCR particle as a juxtaposition of two half-disks of radius R ′ with one rectangle of length L and width 2R ′ ; see Fig. 1. The shape of a RCR particle is a circle of radius R ′ for L = 0. The aspect ratio α = (L + 2R ′ )/(2R ′ ) is 1 in this limit and increases with L for a fixed value of R ′ . In this paper, we use an alternative parameter describing the deviation of the particle shape from a circle. Let R be the radius of the circle circumscribing the particle. We have R = L/2 + R ′ . The radius R ′ is also that of the inscribed circle. Hence, the deviation from a circular shape can be characterized by ∆R = R -R ′ = L/2. We use the dimensionless parameter η defined by

η = ∆R R = α -1 α . (1)
It varies from η = 0, for a circle, to 1 corresponding to a line. We will refer to η as the elongation parameter as in rock mechanics [START_REF] Folk | Petrology of Sedimentary Rocks[END_REF].

The contacts between RCR particles belong to different categories, namely cap-to-cap (cc), cap-to-side (cs) and side-to-side (ss); see Fig. 2. Side-to-side contacts results from contacts between two rectangles as well as two contacts resulting from cap-to-side. In the CD method the case of side-to-side contacts for rectangular particle is represented by two points. Hence, for RCR particles, ss contact is composed of four point contacts : two points due to the rectangle-rectange interface and two points due to the cs contacts. In the iterative procedure of determination of the contact forces and velocities, the points representing the contact between two particles are treated as independent points but the resultant of the calculated forces are attributed to the contact with its application point located on the contact plane.

The detection of line contacts between rectangles was implemented through the so-called shadow overlap FIG. 2: Representation of cap-to-cap, cap-to-side and sideto-side contact and they will be referred as cc contacts, cs contacts and ss contact, respectiveley. method devised initially by Moreau [START_REF]Discrete Numerical Modeling of Granular Materials[END_REF][START_REF] Saussine | [END_REF] for polygons. The reliability and robustness of this method have been tested in several years of previous applications to granular materials [15,16,[START_REF]Discrete Numerical Modeling of Granular Materials[END_REF][START_REF] Saussine | [END_REF][35][36][37]. This detection procedure is fairly rapid and allows us to simulate large samples composed of RCR particles.

B. Packing preparation and bi-axial test

We prepared 8 different packings of 13000 RCR particles with η varying from 0 to 0.7 by steps of 0.1. The radius R of the circumscribing circle defines the size of a RCR particle. In order to avoid long-range ordering in the limit of small values of η, we introduce a size polydispersity by taking R in the range [R min , R max ] with R max = 2R min with a uniform distribution in particle volume fractions.

All samples are prepared according to the same protocol. A dense packing composed of disks (η = 0) is first constructed by means of a layer-by-layer deposition model based on simple geometrical rules [38][39][40]. The particles are deposited sequentially on a substrate. Each new particle is placed at the lowest possible position at the free surface as a function of its diameter. This procedure leads to a random close packing in which each particle is supported by two underlying particles and supports one or two other particles. For η > 0, the same packing is used with each disk serving as the circumscribing circle of a RCR particle. The RCR particle is inscribed with the given value of η and random orientation in the disk.

Following this geometrical process, the packing is compacted by isotropic compression inside a rectangular frame of dimensions l 0 × h 0 in which the left and bottom walls are fixed, and the right and top walls are subjected to a compressive stress σ 0 . The gravity g and friction coefficients µ between particles and with the walls are set to zero during the compression in order to avoid force gradients and obtain isotropic dense packings. Fig. 3 displays snapshots of the packings for several values of η at the end of isotropic compaction.

The isotropic samples were sheared by applying a downward displacement on the top wall at constant velocity for a constant confining stress acting on the lateral walls; see http://cgp-gateway.org/ref010 for video samples. During shear, the friction coefficient µ between particles was set to 0.5 and to zero with the walls. The strain rate was low so that the shearing is basically of quasi-static nature. The internal angle of friction ϕ at every stage of shearing is given by

sin ϕ = q p = σ 1 -σ 2 σ 1 + σ 2 , (2)
with σ 1 > σ 2 are the principals values of the stress tensor σ. sin ϕ increases with shear strain and saturates to a constant value corresponding to the critical state which

f n t ℓ FIG. 4: Local contact frame (n, t)
is a state independent of the initial configuration of the packing. The critical-state value of sin ϕ represents the shear strength of the packing and it increases linearly with η from ∼ 0.3 (for η = 0.0) to ∼ 0.51 (for η = 0.7) [22].

In granular media, the expression of stress tensor σ in the volume V is an arithmetic mean involving the branch vectors ℓ c (joining the centers of the two touching particles) and contact force vectors f c at contact c, and it is given by [START_REF] Moreau | Friction, Arching, Contact Dynamics[END_REF][START_REF] Staron | [END_REF]:

σ = 1 V c∈V f c α ℓ c β , (3)
For the analysis of stress transmission from a particlescale viewpoint we need a statistical description of these quantities.

A common approach used by various authors is to express branch vectors and contact force orientations in terms of the contact direction, i.e. in the local contact frame (n, t), where n is the unit vector perpendicular to the contact plane, and t is an orthonormal unit vector oriented along the tangential force; see figure 4. The components of the branch vector and contact force are expressed in the following frame:

ℓ = ℓ n n + ℓ t t, f = f n n + f t t, (4)
where ℓ n and ℓ t are the normal and tangential components of the branch vectors, and f n and f t the normal and tangential components of the contact force. Remark that only for disks or spherical particles we have ℓ = ℓn where ℓ is the length of the branch vector.

In the following we study the shapes of the normal and tangential force and branch distributions in the residual state.

III. DISTRIBUTIONS OF CONTACT FORCES AND BRANCH VECTORS

A specific feature of the contact network of a packing of elongated particles is that the length ℓ of branch vectors strongly varies throughout the network depending on the relative particle orientations. From the definition of η (Eq. ( 1)) and for given values of R min and R max , it is easy to see that

ℓ R max ∈ 2 R min R max (1 -η), 2 (5)
In our simulations, since R min /R max = 0.5, we have

(1 -η)R max ≤ ℓ ≤ 2R max .
With increasing elongation η, the range of ℓ becomes significant and its statistics can be used as a meaningful characterization of the texture as a function of η. On the other hand, the correlation of ℓ with the total reaction force f between neighboring particles seems to be a good descriptor of the organization of forces for particles of non circular shape. The branch vectors are also important as they enter the expression of the stress tensor given by Eq. ( 3). In Ref. [22], a different point of view was adopted: the contact forces were projected along and perpendicular to the branch vectors and their statistics were investigated. The same framework was used for the decomposition of the total stress tensor. Here, we focus on the distribution of contact forces and their correlation with the branch vector as η is increased.

A. Contact forces and friction mobilization

The probability density function (pdf) of normal forces normalized by the mean normal force f n is shown in Fig. 5 in log-linear and log-log scales at large strains (the data are cumulated from several snapshots in the critical state) for all simulated values of η. As usually observed, in all packings the number of forces above the mean f n falls off exponentially whereas the number of forces below the mean varies as a power-law:

P (f n ) ∝ e -αn(η)(fn/ fn ) , f n > f n , fn fn βn(η) , f n < f n , (6)
where α n (η) and β n (η) whose variations are shown in the insets as a function of η. We see that α n decreases with increasing η, implying that the inhomogeneity of normal forces becomes higher as the particles become more elongated. On the other hand, β n declines from 0.1 to -0.4 with η which means that the proportion of weak contacts (carrying a normal force below the mean) increases with elongation. The proportion of weak forces grows from 60% for η = 0 to 70% for η = 0.7. In other words, while the proportion of strong contacts declines with increasing η, stronger force chains occur at the same time.

Figure 6 shows the pdf P (f t ) of tangential forces normalized by the mean tangential force |f t | in each packing. These distributions show also an exponential falloff for the forces above the average force |f t | and a power law for the forces below |f t | :

P (f t ) ∝ e -αt(η)(|ft|/ |ft| ) , |f t | > |f t | , |ft| |ft| βt(η) , |f t | < |f t | , (7)
the corresponding exponents α t (η) and β t (η) decreasing with η. We observe that, in contrast to α n and β n , the exponents α t and β t saturate beyond η = 0.4. This means that the friction forces do not follow the normal forces as η increases. In other words, the most mobilized (largest) friction forces do not occur necessarily at the contacts where the normal forces are higher.

In order to investigate the properties of friction mobilization, we consider the friction mobilization index

I m = |f t |/µf n . Its average I M = |ft|
µfn increases from 0.4 for η = 0 to 0.6 for η = 0.7 as we see in Fig. 7. This increase underlies to a large extent the increase of the shear strength with η, as we shall see below in Sec. IV. However, the friction force is not uniformly mobilized at all contacts. Fig. 8 shows a map of weak (f n < f n ) and strong (f n > f n ) normal forces, represented by the thickness of vectors joining the particle centers to the contact points, and the corresponding values of I m , represented by circles of diameter proportional to I m for η = 0.1 andη = 0.7. Visual inspection reveals that most mobilized contacts belong to the weak force network. In fact, the average friction mobilization I mf defined as the average by force class, plotted as a function of f n in Fig. -0,1 0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0, 9 for all values of η, declines as f n increases. We also see that the friction mobilization increases with η at all force levels.

Figure 10 displays the pdf of I m for different values of η in the critical state. For the disks, the pdf is a nearly decreasing linear function of I m , which means that the proportion of weakly mobilized contacts is larger than that of strongly mobilized contacts. As η is increased, the distribution becomes more uniform, and at even larger η a class of highly mobilized contacts (with I m close to 1) appears whereas the distribution is nearly uniform for all other contacts. This class belongs to weak force network as was shown previously, so that not only the friction mobilization I m but also the number of highly mobilized contacts are larger in the weak force network. A class of very weak forces was also evidenced in [43] in a packing of disks deposited under gravity and tilted towards its angle of stability. This subclass of the weak network can be defined as the class of contacts where the normal force is below the mean but the friction is highly or fully mobilized.

This enhanced friction mobilization implies that the equilibration of the particles is more complex than in disk packings. In particular, the nematic ordering due to the "geometrical" chains of side/side contacts between particles means that the statistics of forces and the mobilization of friction are closely related to the equilibrium of such chains rather than that of individual particles. These chains are evidenced in Fig. 11 for η = 0.7 where the force bearing particles belonging to the chains are represented by a color level proportionally to their orientations. The friction needs to be highly mobilized inside the chains in order to ensure their stability.

B. Branch vectors

The branch vectors in a packing of elongated particles reflect both the relative orientations of the particles in contact and their size distribution. The latter may be in- Figure 13 displays the pdf of reduced branch-vector lengths for all values of η in the critical state. These pdf's are nothing but normalized radial functions with ℓ r varying in a limited range as only the touching particles are considered. They are nearly similar for all values of η. The first mode, centered on ℓ r = 1 -η reveals the presence of a broad population of side/side contacts with a peak increasing in amplitude with η as displayed in Fig. 14 (b). We also observe a less pronounced mode, centered on ℓ r ≃ 1, corresponding to a distinct population of aligned cap/cap contacts, also marked in Fig. 14 (a). The intermediate mode occurs approximately at ℓ r ≃ 1-η/2 which is the midpoint of the interval [1-η, 1]. This length corresponds to an orthogonal side/cap contact as shown in Fig. 12(c). The presence of such a distinct mode, through decreasing in amplitude as η increases, is a clear proof of the occurrence of orthogonal layers some of which are observed in Fig. 14. This mode is also characterized by a broad extension reflecting the intermediate angles between the orientations of touching FIG. 14: A snapshot of cap-cap modes contact (blue), side-side modes contacts (red) and cap-side modes contacts (green).

particles.

We expect the branch vectors lengths to be correlated with contact forces because of either the contact configurations they represent or simply the fact that force chains tend to be captured by larger particles (hence, longer branch vectors) [12]. This correlation can be estimated with the Pearson coefficient, which for two random variables x and y is defined by the scalar

C xy = (x -x )(y -y ) (x -x ) 2 (y -y ) 2 , (8)
Note that C = 1 corresponds to a full inter-dependence whereas C = 0 means full statistical independence of the two variables. Fig. 15 shows the Pearson coefficients C f ℓr , between the force amplitude f and ℓ r , as well as C f ℓ , between f and ℓ, as a function of η. Both coefficients decrease with η from positive values for η ≤ 0.3 to negative values down to -0.22. The positive correlation (larger forces at longer branch vectors) is a consequence of the fact that the distribution of branch lengths at low values of η is governed by particles sizes. On the other hand, the negative correlation (larger forces at shorter branch vectors) reflects the effect of the increasing number of side/side contacts as the particles become more elongated. Further insight into this force/branch-length correlation can be obtained from the average force amplitude f ℓr , calculated by taking the average force in a class of contacts in the interval [ℓ r -∆ℓ r /2, ℓ r + ∆ℓ r /2], as a function of ℓ r , shown in Fig. 16(a) for all values of η. This plot shows that for all contact classes the associated mean force f ℓr is nearly equal to the global mean force f except to the class of the shortest branch vectors (side/side mode), which concentrates a mean force above f , and the class of the longest branch vectors (cap/cap mode), which seem to carry a considerably lower force on the average. In his way, the rather weak correlation between the reduced branch length and force appears here to be governed by the two afore-mentioned modes. In order to evidence the effect of particle size distribution, let us consider the average force amplitude f ℓ as a function of ℓ as shown in Fig. 16(b) for all values of η. For η ≤ 0.3, the contact force is on the average an increasing function of ℓ. For disks (η = 0), the variation of ℓ is a consequence only of the particle size distribution and, therefore, the increase of the mean force with ℓ means that the larger particles, involved in the longer branch vectors, capture higher forces. The same effect seems thus to underly also the increasing mean force with η for elongated particles with η ≤ 0.3. But at larger elongations, the trend is reversed and we see that the mean force declines as ℓ increases, reflecting thus the effect of the side/side contact mode as discussed previously.

IV. WEAK AND STRONG FORCE NETWORKS

The complex network of contact forces in a packing of elongated particles can also be analyzed by considering the contribution of various classes of forces and/or branch vectors to stress transmission. Indeed, according to equation (3), the stress tensor is expressed as an average involving branch vectors and contact forces, so that partial summations allow one to define partial stress tensors that have been applied in the past to investigate the scale-up of local quanties [11]. For example, the subset of contacts carrying a force below a threshold, reveals the respective roles of weak and strong force chains with respect to the overall shear strength of granular materials [11]. In this section, we apply this methodology to analyze the stress and other texture-dependent quantities in view of elucidating the effect of particle elongation. In what follows, we consider various fabric and force parameters for the "ξ-networks" defined as the subsets S(ξ) of contacts which carry a force below a cutoff force ξ normalized by the mean force (ie f / f n ∈ [0, ξ]), where ξ is varied from 0 to the maximal force in the system. The weak network corresponds to S(1) whereas the strong network is its complement. In section III, we focused on scalar descriptors of granular texture such as the distributions and correlations of force magnitudes and branch lengths. Beyond these low-order quantities, the granular texture is characterized by a disordered but anisotropic structure of both the contact and force networks, which require higher-order description in terms of various fabric and force tensors. We analyze below different parameters pertaining to this tensorial organization of our packings as a function of ξ and for increasing elongation η.

A. Granular texture

A relevant description of granular texture is given by the probability distribution P (n) of the contact normals n ; see Fig. 4. In two dimensions, the unit vector n is described by a single angle θ ∈ [0, π]. The distribution P θ (θ) of contact orientations can be evaluated from the numerical data at different stages of its evolution. In our simulations, all numerical samples are prepared in an isotropic state so that P θ = 1/π in the initial state. This distribution evolves with shear strain and becomes increasingly more anisotropic as the critical state is approached. By restricting the data to those belonging to the ξ-networks, we obtain a continuous family of distributions P θ (θ, ξ) that describe the geometrical state of the system. In practice, however, such functions can be estimated with meaningful statistics only in the critical state where the data can be cumulated from independent configurations representing all the same state.

Figure 17 shows the distributions P θ (θ, ξ) in polar coordinates for η = 0.5 and for several values of ξ. The distributions are similar with nearly the same privileged direction aligned with the principal stress direction θ σ but with increasing anisotropy as a function of ξ. They all can be approximated by their truncated Fourier expansion [11,13,15]:

P θ (θ, ξ) = 1 2π {1 + a c (ξ) cos 2(θ -θ σ )}, ( 9
)
where a c (ξ) is the amplitude of contact anisotropy in the ξ-network. In practice, it is more convenient to estimate a c (ξ) through the partial fabric tensors F (ξ) defined by [START_REF] Satake | Proceedings of the IUTAM symposium on deformation and failure of granular materials[END_REF]:

F αβ (ξ) = 1 π π 0 n α (θ, ξ)n β (θ, ξ)P θ (θ, ξ)dθ, (10)
where α and β design the cartesian components. By definition, we have tr(F (ξ)) = 1. Introducing the harmonic expression ( 9) in ( 10), we get

a c (ξ) = 2(F 1 (ξ) -F 2 (ξ)) cos 2[θ c (ξ) -θ σ ], (11)
where the subscripts 1 and 2 refer to the principal values of F (ξ) and θ c (ξ) represents the privileged direction of the partial fabric tensors F (ξ). Note that, up to statistical fluctuations, the principal directions of the fabric and stress tensors coincide in the critical state for each ξ-network, so that the phase factor cos 2[θ c (ξ) -θ σ ] is either equal to 1 when θ c (ξ) = θ σ or equal to -1 when θ(ξ) = θ σ + π/2. Fig. 18 displays a c as a function of ξ for all values of η. For the disk packings (η = 0), the anisotropy of weak contacts is negative but increases in absolute value and reaches its peak value at ξ ∼ 1. This negative value indicates that in disk packings the weak contacts are orientated preferentially perpendicular to the major principal stress direction [11]. As more contacts come into play with increasing ξ, the partial anisotropy a c (ξ) becomes less negative and finally changes sign, showing that the strong contacts are mainly along the major principal stress direction. This bimodal behavior of stress transmission is a nontrivial organization of the force network and holds also in 3D in the case of sphere packings [16]. However, it is remarkable that for elongated particles (η > 0), the partial anisotropies of both weak and strong networks are positive, as observed in Fig. 18. This means that, in contrast to the disk packings, the weak and strong contacts in packings of elongated particles can not be differentiated on the basis of their roles in the ξ-networks. Physically, this behavior may be interpreted by stating that the static equilibrium of the chains of elongated particles does not require the stabilizing effect of the weak contacts. A similar result was observed by Estrada et al. for disk packings at large values of rolling resistance, which allows for the equilibrium of long chains of particles inter-connected by only two contacts [13]. But, as we shall see below, for our elongated particles the differentiation between the two networks operates via the forces carried by the ξ-networks.

The information involved in the angular distribution P θ may be enriched by accounting for the branch vectors ℓ which, as seen in section III, reflects both the particle size distribution and local contact modes. We thus consider here the average normal and tangential branch vector components ℓ n (θ, ξ) and ℓ t (θ, ξ) defined in (4), obtained by averaging ℓ n and ℓ t over the contacts oriented along θ within a centered angular interval ∆θ. As for P θ , we evaluate these functions in the critical state, for different values of η and as ξ. Fig. 19 shows the functions ℓ n (θ, ξ) and ℓ t (θ, ξ) in polar coordinates for η = 0.5 and for several values of ξ. These functions are anisotropic with an anisotropy which depends on ξ. We introduce here their truncated expansion on an orthonormal Fourier basis:

ℓ n (θ, ξ) = ℓ {1 + a ln (ξ) cos 2(θ -θ σ )}, ℓ t (θ, ξ) = ℓ a lt (ξ) sin 2(θ -θ σ ), (12)
where a ln (ξ) and a lt (ξ) are the normal and tangential branch anisotropies in the ξ-networks. Note that, by construction we have a lt = 0 for disks (η = 0). The analytical form of ℓ t (θ, ξ) results from the orthonormal nature of the Fourier basis and the fact that the mean value of ℓ t vanishes due to disorder:

π 0 ℓ t (θ, ξ) P θ (θ, ξ) dθ = 0. (13)
Fig. 19 shows that this functional form provides a good approximation of the data.

For the calculation of a ln (ξ) and a lt (ξ), we introduce the following branch tensors [22]:

       χ ln αβ (ξ) = 1 ℓ π 0 ℓ n (θ, ξ)n α (ξ)n β (ξ)P θ (θ, ξ)dθ, χ lt αβ (ξ) = 1 ℓ π 0 ℓ t (θ, ξ)n α (ξ)t β (ξ)P θ (θ, ξ)dθ, (14)
The following relations are then easily obtained:

a ln (ξ) = 2[χ ln 1 (ξ) -χ ln 2 (ξ)]/tr[χ ln (∞)] -a c (ξ), a lt (ξ) = 2[χ l 1 (ξ) -χ l 2 (ξ)]/tr[χ l (∞)] -a c (ξ) -a ln (ξ), (15)
where χ l = χ ln + χ lt , and the subscripts 1 and 2 refer to the principal values of each tensor. By construction, we have trχ l = (χ l 1 + χ l 2 ) = ℓ . Note also that the two partial branch vector anisotropies a ln and a lt may take positive or negative values depending on the orientations θ ln and θ lt of the two tensors with respect to θ σ .

Figure 20 shows the branch-vector anisotropies a ln (ξ) and a lt (ξ) as a function of ξ in the critical state for all values of η. a ln (ξ) is positive for η = 0 and η = 0.1 and increases slightly with ξ, but for more elongated particles it takes negative values, which means that the particles tend to form longer branch vectors with their neighbors in the direction of extension. As ξ increases, this anisotropy increases in absolute value and reaches a plateau after passing by a peak value at a point in the range ξ ∈ [1,2]. This behavior suggests that the particles touch preferentially along their minor axes when the contact orientation is close to the compression axis (in the strong network), and along their major axis when the contact orientation is close to the extension axis (in the weak network), in agreement with the fact that the longest branches are in the weak network ; see Sec. III. As to a lt (ξ), its value is always negative and increases monotonically with ξ in absolute value. Note also that, for all values of ξ, a lt (ξ) is much higher than a ln (ξ) while both remain weak compared to a c (ξ).

B. Force anisotropies

We now consider the angle-averaged normal and tangential forces, f n (θ, ξ) and f t (θ, ξ) in the ξ-network. A second order Fourier expansion provides an adequate representation of these distributions for all values of ξ as shown in Fig. 21:

f n (θ, ξ) = f {1 + a f n (ξ) cos 2(θ -θ σ )} f t (θ, ξ) = f a f t (ξ) sin 2(θ -θ σ ), (16)
where a f n (ξ) and a f t (ξ) are the amplitudes of normal and tangential force anisotropies in the ξ-networks. No- tice that we have f t = 0 as a consequence of weak correlation between the branch vectors and contact forces as shown in Fig. 15 and the balance of force moments. Morevover, the orthogonality between the normal and tangential forces implies that the peak value of f t (θ, ξ) occurs at an angle rotated by π/4 with respect to that of are rotated to those of f n (θ, ξ).

As for the branch length vectors, the calculation of the anisotropy parameters a f n (ξ) and a f t (ξ) can be done by means of the following force tensors [11,22]:

       χ f n αβ (ξ) = 1 f π 0 f n (θ, ξ)n α (ξ)n β (ξ)P θ (θ, ξ)dθ, χ f t αβ (ξ) = 1 f π 0 f t (θ, ξ)n α (ξ)t β (ξ)P θ (θ, ξ)dθ. (17)
With these definitions, the following relationships can easily be established:

a f n (ξ) = 2 χ f n 1 (ξ) -χ f n 2 (ξ) tr[χ f n (∞)] -a c (ξ), (18)
a f t (ξ) = 2 χ f 1 (ξ) -χ f 2 (ξ) tr[χ f (∞)]
-a c (ξ) -a f n (ξ), (19) where χ f = χ f n + χ f t and the indices 1 and 2 refer to the principal values of each tensor. By construction, we have tr(χ f ) = χ f 1 + χ f 2 = f . The two partial force anisotropies a f n and a f t may take positive or negative values depending on the orientations θ f n and θ f t of the two tensors with respect to θ σ .

The normal and tangential force anisotropies are plotted in Fig. 22 as a function of ξ for all values of η. A remarkable feature of a f n (ξ) is that its value is negative in the weak network (ξ < 1) for all elongated particles, i.e. for all values of η in exception to η = 0 where it remains positive for all ξ. Hence, the weak forces in a packing of elongated particles occur at contacts preferentially oriented orthogonally to the principal stress direction θ σ whereas in a disk packing they are parallel. As we saw before, an inverse behavior occurs for the contact anisotropies, i.e. the weak contacts in the packings of elongated particles are parallel to the principal stress direction and orthogonal for the disk packings. a f n (ξ) increases in absolute value as ξ increases and passes by a peak at exactly ξ = 1, then declines as more contacts from the strong network with a positive contribution to the anisotropy are included in the ξ-network. At larger values (beyond ξ ≃ 2 for nearly all values of η), a f n (ξ) becomes positive as the strong forces tend to be parallel to the principal stress direction. This unmonotonic behavior of the partial force anisotropies for the elongated particles and the partial contact anisotropies for the disk packings underlies the differentiation between the weak and strong networks according to the values of the normal contact forces with respect to the mean force (ξ = 1). The difference between the elongated particle packings and disk packings reflects the formation of sideside contacts oriented along the principal stress direction tending to capture the strong force chains.

The tangential force anisotropy a f t (ξ) is an increasing function of both ξ and η. Its value is generally below a f n (ξ), but becomes comparable for the most elongated particles for which the friction mobilization plays a key role as discussed previously. This is plausible as the tangential force anisotropy represents friction mobilization at contacts oriented at π/4 with respect to the major principal stress direction.

C. Stress tensor

The physical importance of geometrical and mechanical anisotropies becomes clear when it is considered in connection with the stress tensor. As shown by Eq. 3, the stress tensor is a function of discrete microscopic parameters attached to the contact network. It is also possible to attribute a stress tensor to each ξ-network by restricting the summation to the corresponding contacts:

σ(ξ) = 1 V c∈V f c α (ξ)ℓ c β (ξ). (20)
For sufficiently large systems, the dependence of volume averages on individual discrete parameters vanishes [22,[START_REF] Rothenburg | [END_REF] and the discrete sums can be replaced by integrals as follows:

σ αβ (ξ) = n c Ω f α (ξ)ℓ β (ξ) P ℓf (ξ)df dℓ, (21)
where P ℓf is the joint probability density of forces and branch vectors in the ξ-networks, n c is the number density of contacts for the whole system and Ω is the integration domain in the space (ℓ, f ). The integral appearing in Eq. ( 21) can be reduced by integrating first with respect to the forces and branch vector lengths. Considering the components of the forces and branch vectors in contact frames (n, t), and neglecting the branch/force correlations (see Fig. 16), we get [16,22,[START_REF] Rothenburg | [END_REF]:

σ αβ (ξ) = n c π 0 { ℓ n (θ, ξ) n α (θ, ξ) + ℓ t (θ, ξ) t β (θ, ξ)} × { f n (θ, ξ) n α (θ, ξ) + f t (θ, ξ) t β (θ, ξ)} × P (θ, ξ) dθ. ( 22
)
The expression of the stress tensor by this equation makes appear explicitly the average directional functions representing the fabric and force states. Using the harmonic approximation introduced before, Eq. ( 22) can be integrated with respect to space direction θ and we get the following simple relation:

q(ξ) p ≃ 1 2 {a c (ξ)+ a ln (ξ)+ a lt (ξ)+ a f n (ξ)+ a f t (ξ)}, ( 23
)
where the cross products among the anisotropy parameters have been neglected. This relation expresses the normalized shear stress as a half-sum of texture and force anisotropies. Fig. 23 displays the partial shear stress q(ξ)/p as a function of ξ together with the approximation given by Eq. 23. As we see, equation (23) provides an excellent fit to the data for all values of ξ and η. Interestingly, q(ξ < 1)/p is zero for disk packings, implying that strong forces carry the whole deviatoric load. The partial stress deviator q(ξ = 1)/p in the weak network increases slightly with η but remains in all cases weak (below 0.1). This transition reflects a qualitative change in the condition of local force balance in the presence of clusters as shown in Fig. 11. In other words, for these packings the weak network sustains also partially the deviatoric load applied to the system. The weak values of q/p in the weak network is a consequence of the large positive value a c (ξ = 1) = 0.3 which compensates the negative values of a f n (ξ = 1), a ln (ξ = 1) and a lt (ξ = 1).

V. SUMMARY

In summary, using contacts dynamics simulations, we analyzed the granular texture and topology of forces chains in various packings composed of elongated particles under biaxial compression. As compared to disk packings, the effect of particle elongation is to enhance the heterogeneity of the packings by the clustering of the particles according to their contact modes. In particular, the side/side contacts tend to capture strong force chains and be oriented orthogonally to the major principal stress direction. These features are reinforced as the particle elongation is increased. The probability densities of the normal forces become broader with stronger force chains characterized by an exponential distribution as in disks packings, and with higher number of weak forces decreasing as a power law with the force.

An interesting finding of this work concerns the differentiation between the strong and weak force networks for elongated particles. In contrast to disks packings, where the contacts in the weak network are on the average perpendicular to the contacts in the strong network, the contacts in a packing of elongated particles are, on the average, oriented along the major principal stress direction both in the weak and strong networks. But, the weak forces in the case of elongated particles show a negative anisotropy in the sense that the average normal force in the weak network has its maximum value in the contacts perpendicular to the strong network. In other words, while in the disk packings the strong forces chains are propped by many weak lateral contact, for elongated particles the strong force chains are laterally sustain by less contact but larger weak forces. A harmonic decomposition of the stress tensor shows, however, that for both disks and elongated particles, the compensating effects of force and contact anisotropies lead to small shear stress deviator carried by the weak network.

Our simulation data indicate that the larger global shear strength of a packing of elongated particles increases with elongation mainly due to the increase of friction mobilization and friction force anisotropy. The normal force anisotropy is large but nearly independent of elongation. On the other hand, the correlation between contact forces and branch vectors joining particle centers reveal a sub-network of weak contacts with hight friction mobilization and small branch vector length.

In conclusion, the packings of elongated particles in 2D reveal a nontrivial texture allying the geometry of the particles with the preferred orientations of the contacts induced by shearing and equilibrium of particles. Some features are reminiscent of disk packings but are strongly modulated by the particle shape. More work is underway to clarify the effect of particle shape by focusing on the local structures. On the other hand, many aspects of the packings analyzed in this paper are specific to two dimensions. The side/side contacts in 3D between particles of spherocylindrical shape do not give rise to nematic ordering and the particle rotations and forces moments play a major role in the equilibrium of such particles. This point can only be analyzed by performing 3D simulations of large packings of sphero-cylinders of varying elongation. However, since the class of side/side contacts controls to a large extent the specific behavior of elongated particles in 2D, we believe that similar features should occur in 3D for platy particles, which may give spontaneously rise to geometrical chains of face/face contacts. Such simulations require, however, much more computational effort.

R∆RFIG. 1 :

 1 FIG.1: Shape of a Rounded-Cap Rectangle (RCR).

η = 0. 1 η = 0. 3 η = 0. 6 FIG. 3 :

 1363 FIG.3: Examples of the generated packings at the initial state.

FIG. 5 :

 5 FIG. 5: Probability distribution function of normal forces fn normalized by the average normal force fn in log-linear (a) and log-log (b) scales for different values of η.

FIG. 6 :

 6 FIG. 6: Probability distribution function of tangential forces ft normalized by the average tangential force ft in log-linear (a) and log-log (b) scales for different values of η.

FIG. 7 :

 7 FIG. 7: Friction mobilization IM averaged in the steady state as function of η.

  FIG. 8: A snapshot of the force-bearing particles at η = 0.2(a) and η = 0.7(b) and normal forces represented by the thickness of the segments joining the particle centers to the application point of the force. The strong and weak forces are in back and red, respectively. The diameter of yellow circle is proportional to Im at the contact.

7 FIG. 10 :

 710 FIG. 10: Probability distribution function of the friction mobilization index Im.

FIG. 12 :

 12 FIG. 12: Principals modes of contacts : cap-cap (a), side-side (b) and cap-side (c).

FIG. 15 :

 15 FIG. 15: Correlation C ℓ R f and C ℓf as a function of η averaged in the critical state.

FIG. 16 :

 16 FIG.16: Linear correlation between contact force f and branch length ℓ as a function of η.

FIG. 17 :

 17 FIG. 17: Distributions of contact orientations (symbols) in polar coordinates for η = 0.5 and several values of the force cutoff ξ together with their Fourier fits (11) (full lines).

7 FIG. 18 :

 718 FIG. 18: Partial fabric anisotropy ac as a function of force cutoff ξ normalized by the mean force f for different values of η.

FIG. 19 :

 19 FIG. 19: Distributions of ℓn (θ, ξ)(a) and ℓt (θ, ξ)(b) (symbols) in polar coordinates for η = 0.5 and several values of the force cutoff ξ together with their Fourier fits (19) (full lines).

7 FIG. 20 :

 720 FIG.20: Partial normal and tangential branch vector length anisotropies a ln and a lt as a function of force cutoff ξ normalized by the mean force f for different values of η.

FIG. 21 :

 21 FIG. 21: Distributions of fn (θ, ξ)(a) and ft (θ, ξ)(b) (symbols) in polar coordinates for η = 0.5 and several values of the force cutoff ξ together with their Fourier fits (16) (full lines).

7 FIG. 22 :

 722 FIG.22: Partial normal and tangential force anisotropies a f n and a f t as a function of force cutoff ξ normalized by the mean force f for different values of η.

7 FIG. 23 :

 723 FIG.23: Partial shear stress q/p as a function of force cutoff ξ for different values of η (plain line) together with approximation given by Eq. 23 (points).