E Azéma

F Radjai
email: radjai@lmgc.univ-montp2.fr

G Saussine
email: gilles.saussine@sncf.fr

Quasistatic rheology, force transmission and fabric properties of a packing of irregular polyhedral particles

Keywords: granular materials, polyhedral particles, contact dynamics method, shear strength, granular fabric, force chain PACS: 61.43.Bn, 81.05.Rm, 83.80.Fg, 45.70.Cc

By means of contact dynamics simulations, we investigate a dense packing composed of polyhedral particles under quasistatic shearing. The effect of particle shape is analyzed by comparing the polyhedra packing with a packing of similar characteristics except for the spherical shape of the particles. The polyhedra packing shows higher shear stress and dilatancy but similar stress-dilatancy relation compared to the sphere packing. A harmonic approximation of granular fabric is presented in terms of branch vectors (connecting particle centers) and contact force components along and perpendicular to the branch vectors. It is found that the origin of enhanced shear strength of the polyhedra packing lies in its higher force anisotropy with respect to the sphere packing which has a higher fabric anisotropy. Various contact types (face-vertex, face-face, etc) contribute differently to force transmission and fabric anisotropy. In particular, most face-face contacts belong to strong force chains along the major principal stress direction whereas vertex-face contacts are correlated with weak forces and oriented on average along the minor principal stress direction in steady shearing.

Introduction

Many recent numerical studies of granular media deal with model systems composed of spherical particles. The use of simplified particle shapes and contact interactions is needed in order to focus on the collective behavior of particles which is at the origin of many specific properties of granular materials. On the other hand, the numerical treatment of complex particle shapes by discrete element methods was until very recently out of reach due to demanding computational resources. There is presently, however, considerable scope for the numerical investigation of complex granular packings. This is not only due to available computer power and memory but also because during more than two decades of intense research in this field, many fundamental aspects of granular media have already been established for simplified particle shapes. In particular, various microscopic features such as fabric anisotropy [START_REF] Kruyt | Micromechanical definition of strain tensor for granular materials[END_REF], [START_REF] Bathurst | Micromechanical aspects of isotropic granular assemblies with linear contact interactions[END_REF], [START_REF] Rothenburg | Analytical study of induced anisotropy in idealized granular materials[END_REF], [START_REF] Radjai | Bimodal character of stress transmission in granular packings[END_REF], [START_REF] Kruyt | Kinematic and static assumptions for homogenization in micromechanics of granular materials[END_REF]), force transmission [START_REF] Liu | Force fluctuations in bead packs[END_REF], [START_REF] Radjai | Force distributions in dense two dimensional granular systems[END_REF], [START_REF] Coppersmith | Model for force fluctuations in bead packs[END_REF], [START_REF] Mueth | Force distribution in a granular medium[END_REF], [START_REF] Lovol | Force measurments on static granular materials[END_REF], [START_REF] Bardenhagen | Numerical study of stress distribution in sheared granular material in two dimensions[END_REF], [START_REF] Antony | Evolution of force distribution in three-dimensional granular media[END_REF], [START_REF] Silbert | Statistics of the contact network in frictional and frictionless granular packings[END_REF], [START_REF] Philip | Granular contact force density of states and entropy in a modified edwards ensemble[END_REF], [START_REF] Majmudar | Contact force measurements and stresseinduced anisotropy in granular materials[END_REF]) and friction mobilization [START_REF] Radjai | Bimodal character of stress transmission in granular packings[END_REF], [START_REF] Staron | Friction versus texture at the approach of a granular avalanche[END_REF]) have been analyzed for circular particles (in 2D) and spheres (in 3D). Hence, a recurrent issue today is how robust these findings are with respect to particle shape [START_REF] Ouadfel | stress-force-fabric' relationship for assemblies of ellipsoids[END_REF], [START_REF] Antony | Influence of particle shape on granular contact signatures and shear strength: new insights from simulations[END_REF], [START_REF] Cambou | Anisotropy in granular materials at different scales[END_REF], [START_REF] Nouguier-Lehon | Influence of particle shape and angularity on the behavior of granular materials: a numerical analysis[END_REF], Alonso-Marroquin and Herrmann [2002], [START_REF] Pena | Investigation of the asymptotic states of granular materials using a discrete model of anisotropic particles[END_REF]Pena et al. [ , 2006a,b],b], [START_REF] Azéma | Force transmission in a packing of pentagonal particles[END_REF]).

The issue of shape effect opens actually the door to a vast and substantial scientific domain given a multitude of potential particle morphologies. Several well-known examples are elongated and platy shapes (occurring in biomaterials and pharmaceutical applications), angular and facetted shapes (occurring in geomaterials) and nonconvex shapes (occurring in sintered powders). The macroscopic shear behavior is considerably influenced by particle shape. Rounded particles enhance flowability whereas angular shape is susceptible to improve shear strength, a factor of vital importance to civil-engineering applications [START_REF] Nouguier-Lehon | Influence of particle shape and angularity on the behavior of granular materials: a numerical analysis[END_REF]). In many engineering applications the particle shapes need to be optimized in order to increase performance [START_REF] Morgan | The effect of vibration on ballast beds[END_REF], [START_REF] Wu | The vibration behavior of railway track at high frequencies under multiple preloads and wheel interactions[END_REF], [START_REF] Lim | Discrete element modelling of railway ballast discrete element modelling of railway ballast[END_REF], [START_REF] Saussine | Modelling ballast behaviour under dynamic loading. part1 : A 2d polygonal discrete element method approach[END_REF], [START_REF] Lobo-Guerrero | Discrete element method analysis of railtrack ballast degradation during cyclic loading[END_REF], [START_REF] Lu | The importance of modelling ballast particle shape in the discrete element method[END_REF]).

In this paper, we employ the contact dynamics method to investigate the slow shear behavior of granular media composed of polyhedral particles. The facetted shapes give rise to a rich microstructure where the particles touch at their faces, edges and vertices. Face-face contacts are expected to play a major role in force transmission and statics of polyhedra by accommodating long force chains that are basically unstable in a packing composed of spheres. In order to isolate the effects arising from particle shape, the data from the polyhedra packing will be compared with a packing of spherical particles that, apart from the particle shape, is identical in all respects (preparation, friction coefficients, particle size distribution) to the polyhedra packing. Both packings are subjected to monotonous triaxial compression.

The numerical procedures will be presented with a brief technical introduction to the detection and treatment of contacts between polyhedra in the framework of the contact dynamics method. We will consider the stress-strain and volume-change behavior. The harmonic approximation of the fabric and additive decomposition of the stress tensor into fabric and force anisotropies will be presented in detail. This will allow us to assess in clear terms the origins of shear strength in the polyhedra packing from fabric and force anisotropies in comparison to the sphere packing. The probability density functions of normal forces will be studied and compared between the two assemblies. Finally, we will focus on the contact networks of polyhedral particles and the role played by different contact categories with respect to force transmission.

Numerical method

In this section we briefly introduce the contact dynamics (CD) method with polyhedral particles and the numerical procedures used for sample preparation.

Contact dynamics method with polyhedra

The CD method is based on implicit time integration and nonsmooth formulation of mutual exclusion and dry friction between particles [START_REF] Jean | Unilaterality and dry friction in the dynamics of rigid body collections[END_REF], [START_REF] Moreau | Some numerical methods in multibody dynamics : application to granular[END_REF], [START_REF] Radjai | Etats internes des milieux granulaires denses[END_REF], [START_REF] Dubois | Lmgc90 une plateforme de développement dédiée à la modélisation des problèmes d'intéraction[END_REF]). The equations of motion are formulated as differential inclusions in which velocity jumps replace accelerations [START_REF] Moreau | Some numerical methods in multibody dynamics : application to granular[END_REF]). The unilateral contact interactions and Coulomb friction law are represented as set-valued force laws. The implementation of the time-stepping scheme requires the geometrical description of each potential contact in terms of contact position and its normal unit vector.

At each time step, all kinematic constraints implied by enduring contacts are simultaneously taken into account together with the equations of motion 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 that consists of solving a single contact problem, with other contact forces being treated as known, and iteratively updating the forces until a given convergence criterion is achieved. The method is thus able to deal properly with the nonlocal character of the momentum transfers resulting from the impenetrability of the particles and friction law.

The CD method is unconditionally stable due to its inherent implicit time integration method. The uniqueness of the solution at each time step is not guaranteed for perfectly rigid particles. However, by initializing each step with the forces calculated in the preceding step, the variability of admissible solutions shrinks to the numerical resolution. In the discrete element methods based on molecular dynamics, this "force history" is, by construction, included in the particle positions.

The treatment of a contact interaction between two particles requires the identification of the contact zone and a "common plane". For rigid particles it is possible to define this contact zone by a finite set of points. Before applying the contact detection algorithm between a pair of particles of irregular shapes, a "bounding box" method is used to compute a list of particle pairs potentially in contact. Then, for each pair, the first step is to determine if an overlap exists through a 3D extension of the "shadow overlap method" [START_REF] Saussine | Contribution à la modélisation de granulats tridimensionnels : application au ballast[END_REF], [START_REF] Dubois | Lmgc90 une plateforme de développement dédiée à la modélisation des problèmes d'intéraction[END_REF]). Several algorithms exist for overlap determination between convex polyhedra [START_REF] Cundall | Discrete numerical model for granular assemblies[END_REF], [START_REF] Cundall | Formulation of a three-dimensionnal distinct element modelpart i: a scheme to detect and represent contacts in a system composed of many polyhedral blocks[END_REF], [START_REF] Nezami | A fast contact detection for 3-d discrete element method[END_REF][START_REF] Nezami | Shortest link method for contact detection in discrete element method[END_REF], [START_REF] Dubois | Lmgc90 une plateforme de développement dédiée à la modélisation des problèmes d'intéraction[END_REF], [START_REF] Saussine | Contribution à la modélisation de granulats tridimensionnels : application au ballast[END_REF], [START_REF] Saussine | Modelling ballast behaviour under dynamic loading. part1 : A 2d polygonal discrete element method approach[END_REF], [START_REF] Pérales | Contribution à la modélisation des structures maconnées par approche discrete[END_REF]). When an overlap occurs, the contact plane is determined by computing the intersection between the two particles.

The contacts between polyhedral particles belong to different categories, namely face-face, edge-face, vertex-face, edge-edge, vertex-vertex, vertex-edge; see Fig. 1. The vertex-vertex and vertex-edge contacts are practically absent. In all cases, we determine one, two or three contact points which provide a good description of the contact zone. In this paper, the vertex-edge and edge-edge contacts are referred to as "simple" contacts whereas the edge-face and faceface contacts are treated as "double" and "triple" contacts since their representation involves 2 and 3 distinct points on the common plane, respectively.

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 by different resolution algorithms [START_REF] Dubois | Lmgc90 une plateforme de développement dédiée à la modélisation des problèmes d'intéraction[END_REF]).

Sample preparation

We generate two numerical samples. The first sample (S1) is composed of 36933 polyhedra. The particle shape are taken from a library of 1000 digitalised ballast grains provided by the French Railway Company SNCF. Each particle has at most 70 faces and 37 vertices and at least 12 faces and 8 vertices. Fig. 2 shows several examples of the polyhedral particles used in the simulations. The size of a particle is defined as two times the largest distance between the barycenter and the vertices of the particle, to which we will refer as "diameter" below. We used the following size distribution: 50% of diameter d min = 2.5 cm, 34% of diameter 3.75 cm, 16% of diameter d max = 5 cm. This distribution represents an approximation of that of railway ballast grains. The sample contains 7.1 10 5 vertices and more than 10 6 faces, the average numbers being 20 and 35, respectively. The second sample (S2) is composed of 19998 spheres with exactly the same size distribution as in S1. Fig. 3 shows a snapshot of the two samples in equilibrium state after deposition and isotropic compression under a constant stress of σ 0 = 10 4 Pa in a rectangular box at zero gravity.

The coefficient of friction is 0.5 between the particles in both samples and 0 with the walls. The normal and tangential coefficients of restitution are 0. The zero restitution simplifies the deposition and compaction process by The isotropic samples are subjected to vertical compression by imposing a constant downward velocity of 10 cm/s on the upper wall and a constant confining stress σ 2 = σ 3 = σ 0 on the lateral walls. Each simulation is stopped for a vertical deformation of 30%. The time step was 2.10 -4 s. The CPU time was 2 10 -3 s for S1 and 1 10 -3 s for S2, per particle and per time step on an Apple G5 computer. The deformation process can be considered to be quasistatic in view of the weak kinetic energy injected into the samples compared to the static pressure. This can be expressed more generally through the inertial number defined as [START_REF] Gdr-Midi | On dense granular flows[END_REF]):

I = ε m dp , (1)
where ε = Ḣ/H is the vertical strain rate, m is the total mass, p is the mean pressure and d is the mean particle diameter. In our simulations, we have I ≃ 10 -3 , corresponding to the quasistatic limit.

Stress-strain behavior

In this section, we compare the stress-strain and volume-change behavior between the packings of polyhedra (packing S1) and spheres (packing S2). The stress and strain variables are defined from numerical data. For the estimation of the stress tensor, we use the "tensorial moment" M i of each particle i defined by [START_REF] Moreau | Numerical investigation of shear zones in granular materials[END_REF], [START_REF] Staron | Friction versus texture at the approach of a granular avalanche[END_REF]):

M i αβ = c∈i f c α r c β , (2)
where f c α is the α component of the force exerted on particle i at the contact c, r c β is the β component of the position vector of the same contact c, and the summation runs over all contact neighbors of particle i (noted briefly by c ∈ i).

It can be shown that the tensorial moment of a collection of rigid particles is the sum of the tensorial moments of individual particles [START_REF] Moreau | Numerical investigation of shear zones in granular materials[END_REF]). The stress tensor σ for a packing of volume V is simply given by [START_REF] Moreau | Numerical investigation of shear zones in granular materials[END_REF], [START_REF] Staron | Friction versus texture at the approach of a granular avalanche[END_REF]):

σ = 1 V i∈V M i = 1 V c∈V f c α ℓ c β , (3)
where ℓ c is the branch vector joining the centers of the two touching particles at the contact c. Remark that the first summation runs over all particles whereas the second summation involves the contacts, each contact appearing only once.

Under triaxial conditions with vertical compression, we have σ 1 ≥ σ 2 = σ 3 , where the σ α are the stress principal values. Using the Cambridge representa-tion, we define the mean stress p and stress deviator q by [START_REF] Airey | The Cambridge true triaxial apparatus", Advanced Triaxial Testing of Soil and Rock[END_REF]) :

p = 1 3 (σ 1 + σ 2 + σ 3 ), (4)
q = 1 3 (σ 1 -σ 3 ).
(5)

For our system of perfectly rigid particles, the stress state is characterized by the mean stress p and the normalized shear stress q/p.

The cumulative strain components ε α are defined by

ε 1 = H H 0 dH ′ H ′ = ln 1 + ∆H H 0 , (6)
ε 2 = L L 0 dL ′ L ′ = ln 1 + ∆L L 0 , (7)
ε 3 = l l 0 dl ′ l ′ = ln 1 + ∆l l 0 , (8)
where H 0 , l 0 and L 0 are the initial height, width and length of the simulation box, respectively and ∆H = H 0 -H, ∆l = l 0 -l and ∆L = L 0 -L are the corresponding cumulative displacements. The volumetric strain is given by

ε p = V V 0 dV ′ V ′ = ln 1 + ∆V V 0 , (9)
where V 0 is the initial volume and ∆V = V -V 0 is the total volume change.

The cumulative shear strain is defined by

ε q ≡ ε 1 -ε 2 . ( 10
)
Figure 4 displays the evolution of q/p for the packings S1 and S2 as a function of ε q . For both packings, we observe a classical behavior characterized by a hardening behavior followed by (slight) softening and a stress plateau corresponding to the critical state of soil mechanics [START_REF] Mitchell | Fundamentals of Soil Behavior[END_REF]).

The critical-state strength in the case of polyhedra (≃ 0.46) is twice as high as that of spheres (≃ 0.23). This implies that the polyhedra packing has a

0 0.1 0.2 0.3 0.4 ε q 0.1 0.2 0.3 0.4 0.5 q/p S1 S2
Fig. 4. The strength parameter q/p as a function of shear strain ε q for the polyhedra packing S1 and sphere packing S2. Fig. 5. The volume change ε p as a function of shear strain ε q for the packings S1 et S2.

higher angle of internal friction ϕ defined by

sin ϕ = 3q 2p + q . ( 11
)
At the critical state, we have ϕ = ϕ 0 = 34 • for S1 and ϕ 0 = 18 • for S2.

Figure 5 shows the volumetric strain ε p as a function of shear strain ε q in S1 and S2. In both packings, we observe an early compaction slightly larger in S2 than in S1. The subsequent dilation is lower in S2 and the critical state with isochoric deformation is reached at ε q = 0.3. Dilation in S1 continues with a decreasing rate of volume change but the isochoric plateau is not fully reached. The dilatancy can be expressed in terms of the dilation angle ψ defined by

sin ψ = ε p ε q . ( 12
)
We have ψ ≃ 5 • for S1 and ψ ≃ 2.5 • for S2 at the stress peak state. The variation of ψ versus ϕ, a sort of stress-dilatancy diagram [START_REF] Wood | Soil behaviour and critical state soil mechanics[END_REF]), is displayed in Fig. 6 for polyhedra and spheres. For both packings, we have

ϕ ≃ kψ + ϕ 0 , ( 13
)
where k is a constant slightly smaller than 1 in both packings. This correlation between dilatancy and shear stress during stress-strain transients is a consequence of energy balance. The mechanical work performed on the system is partially dissipated in contact interactions and partially used in volume change [START_REF] Radjai | Contact dynamics study of 2d granular media : Critical states and relevant internal variables[END_REF]). Several stress-dilatancy relations have been proposed as flow rules for plastic deformations of granular media [START_REF] Wood | Soil behaviour and critical state soil mechanics[END_REF]). The relation ( 13) associates the peak state to the largest positive value of dilatancy and the critical state to zero dilatancy. It shows the "non associated" character of the flow rule in granular media (an associated flow rule implying ϕ = ψ).

Harmonic representation of the fabric

The expression of stress tensor in Eq. ( 3) is an arithmetic mean involving the branch vectors and contact forces. Hence, in order to analyze the shear strength properties of the polyhedra packing compared to the sphere packing, we need a statistical description of the granular microstructure (texture or fabric) and force transmission.

In the presence of steric exclusions, the granular microstructure is highly disordered at the particle scale [START_REF] Troadec | Texture locale et plasticité des matériaux granulaires[END_REF], [START_REF] Troadec | Model for granular texture with steric exclusions[END_REF]). Since mechanical interactions are governed by contact and friction, the relevant descriptors of the microstructure are related to the contact network. At the lowest order, the contact network is characterized by the coordination number z which describes the compactness of a packing. This is a crude scalar information in view of the complex arrangement of the particles, but it is well-known that the compactness, generally expressed in terms of the solid fraction, controls the stress-strain behavior under monotonous shearing. Let us remark here that double and triple contact types (see section 2) are counted as single contacts for the coordination number although they are represented by two and three contact points, respectively, in the numerical treatment of interactions between polyhedra.

The evolution of z for polyhedra and spheres is shown in Fig. 7 as a function of ε q . It is remarkable that z is nearly constant in spite of the overall dilation in both packings. We have z ≃ 5.5 for polyhedra and z ≃ 4 for spheres. The connectivity of the contact network can be characterized in more detail by the fraction P (c) of particles with exactly c contact neighbors. The coordination number is the mean value of c : z = c cP (c). Fig. 8 shows P (c) for S1 and S2 in the critical state. The distribution is broader in S1 than in S2. This shows the wider range of potential equilibrium states in the polyhedra packing. For both packings, we observe a peak centered on c = 4 with a higher probability for S2.

Since the shear stress corresponds to the deviation of stress components from the mean stress p along different space directions, the coordination number z as a scalar quantity cannot account for the shear stress and its evolution with strain. Indeed, the expression of the stress tensor suggests that the useful information for the analysis of shear stress is the density and average force as a function of contact orientation. Such functions can be expanded in spherical harmonics in 3D [START_REF] Ouadfel | stress-force-fabric' relationship for assemblies of ellipsoids[END_REF]).

Let n be the unit vector along the branch vector ℓ ; Fig. 9. We set

ℓ = ℓn, ( 14
)
where ℓ is the length of the branch vector. We remark that the unit vector n does not coincide with the contact normal except in the case of spheres. We consider the components of the contact force in a local frame defined by n and an orthoradial unit vector t:

f = f n n + f t t, (15)
where f n and f t are the radial and orthoradial components of the contact force, respectively. The writing of Eq. ( 15) assumes that t is oriented along the orthoradial force.

We now define the angular averages associated with the branch vectors ℓ and contact force vectors f . Let A(Ω) be the set of branch vectors pointing in the direction Ω ≡ (θ, φ) up to a solid angle dΩ and N c (Ω) its cardinal. The angles θ and φ are shown in Fig. 10. The angular averages are defined as follows:

P Ω (Ω) = N c (Ω) N c , (16)
ℓ (Ω) = 1 N c (Ω) c∈A(Ω) ℓ c , ( 17
)
f n (Ω) = 1 N c (Ω) c∈A(Ω) f c n , (18)
f t (Ω) = 1 N c (Ω) c∈A(Ω) f c t , (19)
where N c = N c (Ω)dΩ is the total number of contacts, and ℓ c , f c n , and f c t are the actual values of branch vector length, radial force and orthoradial force for contact c, respectively.

Under the axisymmetric conditions of our simulations, the four functions defined in Eq. ( 19) are independent of φ. Fig. 11 displays a polar representation of these functions in the θ-plane for polyhedra (S1) and spheres (S2) at ε q = 0.3. We observe an anisotropic behavior in all cases except in ℓ (θ) for S2. A weak anisotropy can be seen for S1 in the latter case. The peak values occur along the compression axis except for f t (θ) in which the peaks are inclined at π/4 with respect to the vertical. The magnitude of anisotropy is larger for polyhedra compared to spheres except for P Ω (θ) which is weakly anisotropic for polyhedra.

The simple shapes of the above functions suggest that harmonic approximation based on spherical harmonics at leading terms captures their anisotropies. There are 9 second-order basis functions Y l m (θ, φ). But only the functions compatible with the symmetries of the problem, namely independent with respect to φ and π-periodic as a function of θ, are admissible. For P Ω (θ) as a scalar, and ℓ (θ) and f n (θ) as radial components of the vectors, the only admissible functions are Y 0 0 = 1 and Y 0 2 = 3 cos 2 θ -1. For f t (θ) as orthoradial component of a vector, the only function independent of φ and perpendicular to Y 0 0 = 1 and Y 0 2 = 3 cos 2 θ -1 is sin 2θ. Hence, within the harmonic model of fabric and force, we have

P Ω (θ) = 1 4π { 1 + a [3 cos 2 θ -1] }, (20)
ℓ (θ) = ℓ 0 { 1 + a l [3 cos 2 θ -1] } (21) f n (θ) = f 0 { 1 + a n [3 cos 2 θ -1] }, ( 22
) f t (θ) = f 0 a t sin 2[θ], (23)
where a, a l , a n and a t are the anisotropy parameters, ℓ 0 is the mean branch vector length, and f 0 the mean force. The probability density function P Ω (θ) is normalized to 1 ( S P Ω (Ω)dΩ = 1, where S is a sphere of unit radius). The values of the anisotropies a, a l , a n and a t can be calculated from generalized Fig. 12. Evolution of anisotropies a, a l , a n and a t with ε q for packings S1 and S2.

fabric tensors introduced in Appendix A.

The evolution of the anisotropies with ε q are displayed in Fig. 12 for our packings S1 and S2. The fabric orientation anisotropy a increases with ε q and relaxes to a plateau after passing by a pronounced peak. Its value is systematically larger for spheres than for polyhedra (by a factor 3 in the critical state). The branch vector anisotropy a l is quite low compared to other anisotropies and its value all along shearing is negligible for spheres. It is remarkable that a l for polyhedra declines (as ε p , see Fig 5) at the beginning of shearing. The radial force anisotropy a n increases as the fabric anisotropy and tends to a plateau. But, in contrast to fabric anisotropy, its value is higher for polyhedra than for spheres. In other words, the aptitude of the polyhedra packing to develop large force anisotropy is correlated with particle shape rather than with fabric anisotropy (see section 7). The orthoradial force anisotropy a t has a similar behavior except that it takes considerably higher values in the case of polyhedra compared to spheres. In the following section, we study the relationship between the fabric and force anisotropies.

Origins of shear stress

In this section, we analyze the stresses in the framework of the harmonic approximation of granular microstructure introduced in the last section. Since this representation involves continuous functions of contact orientations, we need to express the stress tensor in integral form. The stress tensor as defined in Eq. ( 3) is an average:

σ αβ = n c ℓ k α f k β k , (24)
where n c = N c /V is the number density of contacts, ℓ k α is the α component of the branch vector at contact k and ℓ k β is the β component of the force vector at contact k. The average is taken over all contacts k in the control volume V . To express this mean as an integral, we introduce the joint probability density P Ωf ℓ (n, f , ℓ) of the force and branch vectors [START_REF] Bathurst | Micromechanical aspects of isotropic granular assemblies with linear contact interactions[END_REF], [START_REF] Rothenburg | Analytical study of induced anisotropy in idealized granular materials[END_REF], [START_REF] Ouadfel | stress-force-fabric' relationship for assemblies of ellipsoids[END_REF]). Then, from Eq. ( 24), we have

σ αβ = n c P Ωf ℓ (n, f , ℓ) ℓ(n) f β (n, ℓ) n α dΩ df dℓ, ( 25
)
where dΩ = sin θdθdφ.

Equation ( 25) can be simplified by integrating out the contribution of ℓ. Assuming that f is independent of ℓ (an assumption which is verified with a good approximation), we get

σ αβ = n c P Ωf (n, f ) ℓ (n) f β (n) n α dΩ df , (26)
where ℓ P Ωf = P Ωf ℓ dℓ.

Finally, integration of (26) over force vector yields the following expression for the stress tensor:

σ αβ = n c P Ω (n) ℓ (n) n α f β (n) dΩ, (27)
where f P Ω = P Ωf df . By introducing the average force components f n and f t in this equation, we get

σ αβ = n c P Ω (n) ℓ (n) { f n (n) n β + f t (n) t β }dΩ. ( 28
)
This writing of the stress tensor involves the functions previously introduced with the harmonic representation of the fabric (Eqs. ( 21), ( 22), ( 23) and ( 23)). Inserting these functions in the integral expression (Eq. 28) and given the definitions of mean stress p and stress deviator q in Eq. ( 5), one gets q/p, S1 0.4 (a+a n +a t +a l ) q/p, S2 0.4 (a+a n +a t )

p ≃ n c ℓ 0 f 0 , (29) q p ≃ 2 5 (a + a l + a n + a t ), ( 30
Fig. 13. The normalized shear stress q/p as a function of shear strain ε q for the packings S1 and S2 both from direct simulation data and theoretical prediction of Eq. ( 30).

where the cross products (aa l , aa n and aa t ) among the anisotropies have been neglected. Our simulation data are in quantitative agreement with this "stressforce-fabric" relation (a term coined by Rothenburg and Bathurst in [START_REF] Bathurst | Micromechanical aspects of isotropic granular assemblies with linear contact interactions[END_REF], [START_REF] Rothenburg | Analytical study of induced anisotropy in idealized granular materials[END_REF]) ) both for spheres and polyhedra, as shown in Fig. 13, all along the shear. We note that the theoretical fit would have been less satisfactory for polyhedra if the branch vector length anisotropy a l were omitted from the description.

Equation ( 30) is interesting as it exhibits the two origins of shear stress in a granular system: 1) the fabric anisotropies a and a l , related to the branch vector and 2) the force anisotropies a n and a t , related to the contact force.

Figure 12 shows that the values of these anisotropy parameters underlying the shear stress depend on the particle shape. In particular, the total force anisotropy a n + a t compared to the total fabric anisotropy a + a l is much higher in the case of polyhedra. In the critical state, we have a n + a t ≃ 0.88 and a+a l ≃ 0.2 for polyhedra, a n +a t ≃ 0.26 and a+a l ≃ 0.24 for spheres. The high value of the force anisotropy in the case of polyhedra comes from both radial and orthoradial components whereas in the sphere packing a t ≃ 0.05 is much less important than a n ≃ 0.21. This suggests that friction is more directly involved in force transmission in the polyhedral packing than in the sphere packing. The strong contribution of force anisotropy to the polyhedra packing is a particle shape effect related to the face-face contacts which carry most strong forces. This point will be analyzed in more detail below.

Force distributions

In this section, we study the probability density functions (pdf's) P (f n ) for sphere and polyhedra packings. Fig. 14 shows typical maps of normal forces in a portion of both packings in the critical state. The 3D force chains can be observed in both packings, but they seem more tortuous in the case of polyhedra.

The normal force pdf's are shown in Fig. 15 on log-linear and log-log scales at ε q = 0, 3. In both pdf's, the strong forces, i.e. forces above the mean normal force f n , fall off exponentially: P (f n ) ∝ e -βfn/ fn , with β ≃ 0.9 for S1 and β ≃ 1.1 for S2. In contrast, the shapes of the pdf's in the range of weak forces (f n < f n ) are radically different. In the sphere packing, the pdf slightly bends down as f n → 0 but does not tend to zero. We observe also a small peak close to the mean force. This is consistent with several other numerical and experimental observations for isotropic packings [START_REF] Lovol | Force measurments on static granular materials[END_REF], [START_REF] Bardenhagen | Numerical study of stress distribution in sheared granular material in two dimensions[END_REF], [START_REF] Antony | Evolution of force distribution in three-dimensional granular media[END_REF], [START_REF] Silbert | Statistics of the contact network in frictional and frictionless granular packings[END_REF], [START_REF] Philip | Granular contact force density of states and entropy in a modified edwards ensemble[END_REF], [START_REF] Majmudar | Contact force measurements and stresseinduced anisotropy in granular materials[END_REF]). In the case of polyhedra, the number of weak forces bends up as the force tends to zero. For both packings, the range of weak forces is well approximated by a power-law distribution :

P (f n ) ∝ [ fn fn ] α , (31)
with α = -0.24 for S1 and α = 0.05 for S2. The divergence of the number of weak forces in S1 should be attributed to the polyhedral shape of the particles favoring the arching effect an hence a higher fraction of weak forces. The coefficient of friction has a similar effect though to a lesser extent. We find, however, that in both systems the fraction of weak forces (f n < f n ) is about 60%.

Contact networks of polyhedral particles

In the case of the polyhedra packing, it is interesting to investigate the organization of the contact network in terms of simple, double and triple contacts.

The respective fractions of these contact types and their contributions to the structural anisotropy and force transmission are the key quantities for understanding the effect of particle shape on the shear strength properties of granular media. In fact, one expects that the triple (face-to-face) contacts play an essential role in force transmission. This feature was observed in the case of polygon packings for side-to-side contacts [START_REF] Azéma | Force transmission in a packing of pentagonal particles[END_REF]).

Considering the discrete expression of the stress tensor in Eq. ( 3) and restricting the summation to each contact type allows us to perform an additive decomposition:

σ = σ s + σ d + σ t , (32)
where the subscripts s, d and t design the respective contributions of simple, double and triple contacts. The corresponding stress deviators q s , q d and q t are then calculated and normalized by the mean stress p. Fig. 16 shows the evolution of partial shear stresses q s /p, q d /p and q t /p as a function of shear strain ε q . The contribution of simple contacts is larger for double and triple contacts. However, the double and triple contacts support together the largest portion of the overall shear stress, i.e. q d + q t > q s . The partial shear stress supported by each contact type depends on both the number of its contacts and their mean force. Fig. 17 shows the proportions k s , k d and k t of simple, double and triple contacts as a function of shear strain. k s declines during shear from 0.75 to 0.71 whereas k d and k t increase from 0.14 to 0.15 and from 0.11 to 0.14, respectively. Hence, the critical state is characterized by k s ≃ 0.7 et k t ≃ k d ≃ 0.15. Fig. 17 also shows the relative mean forces f s , f d and f t defined by

f s = k s f n s / f n , (33)
f d = k d f n d / f n , ( 34
) f t = k t f n t / f n , (35)
where f n s , f n d and f n t correspond to the mean normal forces of simple, double and triple contacts. We see that f s declines slightly with strain but is nearly two times larger than f t and 2.3 times larger than f d in the critical state. We have f s ≃ f d + f t . Hence, the lower contribution of triple contacts with respect to shear stress can be attributed to both the low level of the mean force (f t < 0.3) sustained by this class and to their weak number (< 15%). Following the same procedure as for the stress tensor, we now perform a similar decomposition of the fabric tensor F , defined by Eq. (A.4), into three terms:

F = F s + F d + F t , (36) (37)
where F s , F d and F t are the contributions of simple, double and triple contacts.

The corresponding anisotropies a s , a d and a t can be extracted, but since the principal directions of these partial fabric tensors are not necessarily identical to those of the overall fabric tensor, we define the "signed" anisotropies by multiplying each partial anisotropy a i by a phase factor cos 2(θ F -θ F i ):

a ′ i = a i cos 2(θ F -θ F i ). ( 38
)
Figure 18 shows the evolution of signed anisotropies of the three contact classes. We see that a ′ d and a ′ t increase with shear strain and tend to the limit value ≃ 0.04. As to a ′ s , we observe an initial increase followed by rapid decrease and change of sign at ε q ≃ 0.2. This evolution means that during shear the branch vectors of simple contacts tend to become increasingly perpendicular to the major principal direction (the direction of compression). A map of contact forces projected along the branch vectors is displayed in Fig. 19 in different colors according to the type of contact. The triple contacts, despite their lower proportion, appear clearly to be correlated in the form of long chains across the packing. These are mostly parallel to the direction of compression. We also observe a large number of weak forces mainly at simple contacts.

The pdf's of normal forces are shown in Fig. 20 separately for simple, double and triple contacts. The three contact types are involved in strong and weak networks. The strong forces have in all cases an exponential behavior but a major difference is observed in the range of weak forces where the proportion of simple contacts prevails. This correlation between simple and weak contacts is interesting as it clearly reveals the contrast between simple contacts, on one hand, and double and triple contacts, on the other hand, in the organization of the force network. In order to situate the simple, double and triple contacts with respect to the force network, we have plotted in Fig. 21

Conclusion

In this paper, granular materials composed of irregular polyhedral particles were numerically investigated. Macroscopic and microstructural properties were analyzed by (1) direct comparison with a similar packing composed of spherical particles and (2) characterization of contact networks and force transmission. A novel finding of this work is that the origin of enhanced shear strength in a polyhedra packing compared to a sphere packing lies in force anisotropy induced by particle shape. The fabric anisotropy associated with the network of branch vectors is lower in the polyhedra packing. This finding extends the results of a previous study of pentagonal particles in two dimensions to three dimensions [START_REF] Azéma | Force transmission in a packing of pentagonal particles[END_REF]. In other words, the force anisotropy, partially underlying shear strength, is mainly controlled by the fabric anisotropy in a sphere packing. This mechanism breaks down to some extent in a packing of polyhedra where force anisotropy results mainly from the "facetted" particle shape.

The face-face contacts were shown to belong mostly to the strong force network. The local equilibrium structures involving face-face and edge-face contacts accommodate force lines that are basically unstable with spherical par-ticles. Hence, the term "arching" seems to be more adapted to the description of force patterns in an assembly of polyhedra than in an assembly of spheres. This effect is crucial for the probability density function of normal forces in the range of weak forces that is well approximated by a decreasing power-law in the case of polyhedra.

In this investigation the polyhedra were irregular with a given number of faces, edges and vertices. These shape parameters can now be changed and the resulting packings can be analyzed along the same lines as in the present investigation. Since the face-face contacts seem to play a key role, it would be interesting to consider irregular polyhedra with less faces in number but with larger areas. From a mechanical point of view, there should be little difference between a small face and a vertex. The best shape from the shear strength viewpoint can be obtained with a large number of faces of large area, but these two conditions can not be realized at the same time. It seems thus that an optimal polyhedral shape should exist with a number of faces of not two low areas. The work is under way to elucidate this point and other aspects of the problem concerned with packing structure by systematically changing the particle shape parameters.

We acknowledge assistance by F. Dubois with the LMGC90 platform employed for the simulations, as well as the precious help of V. Richefeu with 3D visualization of forces. This work was funded by the French Railway Society, the SNCF, and the Région Languedoc-Roussillon of France.

A Fabric tensors

The anisotropies a, a n , a t and a l can be calculated from the tensors F , H (n) , H (t) and H (l) defined by [START_REF] Bathurst | Micromechanical aspects of isotropic granular assemblies with linear contact interactions[END_REF], [START_REF] Rothenburg | Analytical study of induced anisotropy in idealized granular materials[END_REF], [START_REF] Ouadfel | stress-force-fabric' relationship for assemblies of ellipsoids[END_REF]) : Using the equations ( 21), ( 22), ( 23) and ( 23), it is then easy to show that the corresponding anisotropies are :

Fig. 1 .

 1 Fig. 1. Different types of contacts between two polyhedra.

Fig. 2 .

 2 Fig. 2. Examples of polyhedra used in the simulations.

Fig. 3 .

 3 Fig. 3. Snapshots of the two packings S1 (polyhedra) and S2 (spheres). The walls are not shown enhancing dissipation during dynamics rearrangements. The initial value of the solid fraction is ρ ≃ 0.6 in both samples. Both samples have a nearly square bottom of side such that L ≈ l and an aspect ratio H/L ≃ 2, where H is the height. The initial configuration is defined by H 0 ≃ 30D M for S1 and for S2 with D M the mean diameter.

Fig. 6 .

 6 Fig. 6. The stress-dilatancy diagram representing the relation between the internal angle of friction and the dilation angle for polyhedra and spheres.

Fig. 7 .Fig. 8 .

 78 Fig. 7. Evolution of the coordination number z as a function of the cumulative shear strain ε q for polyhedra (S1) and spheres (S2)

Fig. 9 .

 9 Fig. 9. Geometry of a contact between two polyhedra.

  Fig. 10. Spherical coordinates.

Fig. 11 .

 11 Fig. 11. Polar representation of density probability function P Ω (θ), f n (θ), f t (θ) and ℓ (θ) for S1 et S2 in residual state.

Fig. 14 .Fig. 15 .

 1415 Fig. 14. Force maps in a portion of the packings S1 (right) and S2 (left). The segments are branch vectors with thickness proportional to the normal force, and gray level proportional to the depth of field.

Fig. 16 .Fig. 17 .

 1617 Fig. 16. Evolution of partial shear stresses as a function of shear strain for simple (s), double (d) and triple (t) contacts, as well as the total shear stress (s+d+t).

Fig. 18 .

 18 Fig. 18. Evolution of the signed anisotropies a ′ of simple (s), double (d) and triple (t) contacts as a function of shear strain ε q .

Fig. 19 .

 19 Fig.19. Map of contact forces projected along branch vectors at ε q = 0.4. Line thickness is proportional to the force. The simple, double and triple contacts are in red (dark gray), in green (light gray) and in blue (black).

Fig. 20 .

 20 Fig. 20. Probability distribution functions of radial forces at simple (s), double (d) and triple (t) contacts on log-linear (a) and log-log (b) scales.

  trH (n) = f , trF = 1 et trH (l) = ℓ 0 .


O Ravel

The Codalema Collaboration

30TH INTERNATIONAL COSMIC RAY

The principle of the CODALEMA experiment is based on an original approach of the detection of radio transients associated with extensive air showers induced by ultra high-energy cosmic rays. Since September 2006, CODALEMA has been under operation with a new setup at the Nançay Radio Observatory, France. It uses 16 broadband dipole antennas associated with 13 particle detectors generating the trigger and allowing the primary cosmic ray energy estimation. We will present evidence for the radio detection of cosmic rays above 10 17 eV, based on an event-by-event analysis and we will discuss the radio characteristics of these showers.

Introduction

Radio emission associated with the development of Extensive Air Showers (EAS) initiated by high-energy cosmic rays was investigated in the 1960's [START_REF] Weekes | Proc. of the first int. workshop on "Radiodetection of high Energy Particles[END_REF]. Since 2002, with the availability of fast electronics, several groups have undertaken the task of reinvestigating the phenomenology of radio pulses. The experiments CODALEMA, located at the Nançay radio Observatory in France and LOPES in Germany obtained firm evidences for a radio emission counterpart of cosmic ray air showers [2,3,4].

The new CODALEMA setup including arrays of active dipoles and ground particle detectors is described in details in an other contribution of this conference [5].

Event selection

The results thereafter presented have been deduced from a sample of events taken since the beginning of December 2006 and corresponding to 104 effective days of data acquisition. The trigger of the experiment is generated by coincidences of the scintillators stations. For the period of study we have a sample of 15029 internal events air showers (an event is classified as "internal event" if the 5 central particle detectors are flagged and the central station gets the maximum signal). Among these events, 101 events detected by the antennas meet the following criteria: a time coincidence and an angular difference of the arrival direction within respectively 200 ns and 20° between the information reconstructed from the particle and radio detectors. The direction of the particle showers is determined by triangulation independently made with the two arrays, using arrival times from the digitised signals. The shower energy of the "internal events" is determined by the Constant Intensity Cut (CIC) method [6] with a precision of around 30 %. Finally if we take into account the "internal events" criterion we have 31 radio-detected cosmic rays for the period.

Radio detection efficiency

The energy distribution as above determined for internal events measured by the ground particle array and detected by the dipole array are displayed in Figure 1. Moreover we select only showers inside a radius of 200 m around the central scintillator station. The energy threshold of the ground particle array is around 10 15 eV; the two distributions converge when energy increases suggesting a rise in the radio-detection efficiency. The efficiency plot of the Figure 2 allows to estimate the radio-detection threshold at 10 % at 9.10 16 eV. Because of the low number of events available we do not establish the efficiency at higher energy, nevertheless we should expect a full efficiency around 10 18 eV.

Geomagnetic effect

Figure 3 shows the sky map of the 101 events reconstructed using the antenna information.

Being strongly correlated with their particle detector information counterpart, they are clearly associated to air showers. A noticeable deficit of events coming from the southern direction can be observed whereas the azimuth angle distribution measured by particle detectors is compatible with a uniform distribution. In Figure 4, the zenith angle distribution shows the classical zenithal dependence.

. The distribution of the angle α between shower axis and geomagnetic field (Figure 5) associated with each radio-detected event seems to confirm the radio-detection suppression for small values of α (below 15°). Moreover, we observe an obvious dependence of the radio event rate on the geomagnetic angle. Figure 6 displays the shower energy as a function of the α angle. We notice that we detect less showers of low energy (~10 17 ) when α is small (~20°) and on the contrary more events when α is close to 70°, i.e. when the arrival direction of the shower is close to the perpendicular direction to the geomagnetic field, where the shower charge deviation is maxim With the previous CODALEMA setup [2] using conic helix log-periodic antennas, this effect has not been highlighted. The helix antenna gain pattern is maximum for a zenith angle of 20° towards south, which allows a best detection of showers of low energy (10 17 eV) for the south arrival direction. The current CODALEMA setup features only active dipoles. The gain pattern, in contrary to conic helix antennas is symmetric for southern and northern arrival directions. That could explain why we observe this deficit in the south close to the geomagnetic orientation. This result seems to be consistent with the contribution of the geomagnetic emission process in the generation of radio signal from air shower induced by cosmic rays. However, this observation alone does not exclude other possible radio emission contributions.

Electric field lateral profile

As in the first configuration of CODALEMA using conic helix log-Periodic antennas we find again a satisfactory dependence of the electric field amplitude to the shower axis distance with the 2 parameters (E 0 , d 0 ) fit . Figure 7 displays one of the most energetic event measured by CODALEMA. This event is unfortunately not "internal"; despite that, we can estimate with confidence its energy because the position of the shower core (evaluated by radio) allows to constrain the fit. For this event, E ) / ( exp . However, we need at least 5 antennas tagged to be able to fit the E profile. We do not have enough "internal event" data to seek a correlation between the value of E 0 and the estim of the primary c

The first results of the CODALEMA experiment with the new set-up (using only active dipoles) confirm the field characteristics from our previous measurement [2,3]. The measurement of showers arrival time and direction coincidence is clearly seen. Moreover the addition of a dense network of particle detectors makes possible the shower measurement and the correlation of the measured energy with radio observables. We measure the CODALEMA radio-detection efficiency threshold at 9.10 16 eV (100% efficiency is expected at 10 18 eV). On the other hand, the new result is the apparent depletion of the number of radio-detected events for south directions, which seems consistent with the decreasing of the electric field geosynchrotron contribution at small α. This dependence in α doesn't call into questions the cosmic rays radio-detection capacities. Indeed we can suppose that the Coulombian emission processes, presently believed to be less important will be dominant close to the terrestrial magnetic field direction. I.e. for a given radio detector the energy threshold should be higher for the South arrival direction. Improved statistics for "internal events" and an experimental study of the signal polarization will allow confirming this assumption. Moreover the exponential dependence of the electric fields to the shower axis distance is also observed and should make possible to establish the correlation between the energy of the primary cosmic ray with the electric field amplitude. We will increase at the end of the year the size of the array at Nançay by addition of autonomous antenna stations. The CODALEMA collaboration also explores the radio detection beyond 10 18 eV in colla

Figure 1 :

 1 Figure 1: Distribution of the energy measured by the ground particle detector array (top, blue dots) and the distribution of the radiodetected showers for "internal events" ottom, red dots) (b .

Figure 2 :

 2 Figure 2: Efficiency of radio detection versus energy deduced from the scintillator analysis.

Figure 3 :

 3 Figure 3: Sky map of the arrival directions of 101 cosmic ray showers recorded by the CODALEMA dipoles in coincidences with particle detectors, after time and angular selections. The red triangle represents the geomagnetic field orientation. Dotted internal circles refer to indicated values of zenith angles up to the horizon; the external values stand for azimutal angles.Such north-south asymmetry was previously observed at the end of the 60' at Haverah Park[6]. A probable explanation of this radio-detection deficit in the south celestial hemisphere could issues from the geomagnetic field oriented at 27° from the zenith and south azimuth. Indeed, one can expect a suppression of the geosynchrotron

Figure 4 :

 4 Figure 4: Zenith (top) and azimut (bottom) angular distributions of the internal showers and E > 10 16 eV detected by the particle etector array used as a trigger. d

Figure 5 :

 5 Figure 5: Distribution of α, this distribution is corrected for the acceptance distribution of the trigger calculated from histograms of the Fig. 2.

Figure 6 :

 6 Figure 6: Primary cosmic ray energy (eV) versus the α angle. With the present counting, no radio signal is detected below 17°.

  /m and d 0 =135 m. Response of the 12 dipoles simultaneous tagged are corrected for the antenna and amplifier gains. The energy is evaluated at 2.10 18 eV.

Figure 7 :

 7 Figure 7: EAS electric field amplitude versus distance to the shower axis for the higher energetic event. The inset displays the topology of the ground field with the econstructed direction and core location of r the shower.