Yannick Descantes
email: yannick.descantes@univ-eiffel.fr

A new discrete element modelling approach to simulate the behaviour of dense assemblies of true polyhedra

Keywords: Discrete element method, Molecular dynamics, Convex polyhedra, Granular packing, Granular flows, Rockfall simulation

A new molecular dynamics-like modelling approach aimed at simulating the mechanical behaviour of true polyhedra dense assemblies is presented. Thanks to an improved version of the Gilbert-Johnson-Keerthi contact detection algorithm called GJK -T D, this approach involves no edge or corner rounding and accounts for multipoint f ace -f ace or edge -f ace contacts between any set of two particles. Furthermore, torques equations are simply and efficiently solved using an improved leap-frog Verlet non-iterative method suggested by Omelyan [I.P. Omelyan, Molecular Simulation, 22 : 3, 213 -236 (1999)], in which no periodic renormalization of the quaternions is necessary. The potential of this new discrete element approach is then highlighted by simulating the gravity packing of frictionless polyhedra and the gavity flow of frictional polyhedra down an incline, and comparing the results with those reported in the literature. Of particular interest in the stationary flow regime, a linear decrease of the bulk solid fraction and a power law decrease of the bulk coordination number with increasing inertial number are observed, thus generalizing to polyhedra these observations initially made with disks in plane shear flow by da Cruz and coworkers [F.da Cruz, S. Emam, M. Prochnow, J.N. Roux and F. Chevoir, Phys.

)]. Beyond this regime, in the collisional regime, the granular temperature was found to achieve its maximum a few particle diameter above the rough bottom, suggesting the localization of significant particle

Introduction and background

Discrete element methods (DEMs) refer to numerical methods focussed on simulating the behaviour of assemblies of rigid macroscopic particles at micro to meso scale, by solving the equations of motion of individual particles interacting through contact or collision laws. Since early developments performed four decades ago by Cundall and Strack [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF], DEMs have been extensively used to investigate the behaviour of granular media such as rock blocks, natural sand, mineral aggregates, pharmaceutical tablets or agricultural seeds when subjected to packing [START_REF] Makse | Packing of compressible granular materials[END_REF][START_REF] Silbert | Geometry in Frictionless and Friction Sphere Packings[END_REF], flowing [START_REF] Silbert | Granular flow down an inclined plane[END_REF][START_REF] Da Cruz | Rheophysics of dense granular flows: Discrete simulation of plane shear flows[END_REF], crushing [START_REF] Cheng | Discrete element simulation of crushable soil[END_REF] or heat exchange [START_REF] Zhou | Particle scale study of heat transfer in packed and bubbling fluidized beds[END_REF] processes to mention just a few. At present, DEMs cover a variety of methods often broken down into smooth-DEM s and non -smooth -DEM s, which differ essentially in terms of description scales and contact laws regularization [START_REF] Radjaï | Contact dynamics as a nonsmooth discrete element method[END_REF]. In smooth -DEM s [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF], also known as molecular dynamics-like (MD-like) methods, particles are assumed to undergo local deflection at contacts (the rigidity assumption holds everywhere else), allowing to explicitly calculate contact forces from particles relative location and velocity vectors using regularized (visco-elastic) contact laws. By contrast, non -smooth -DEM s such as the contact dynamics (CD) method [START_REF] Radjaï | Contact dynamics as a nonsmooth discrete element method[END_REF][START_REF] Radjaï | Discrete-element Modeling of Granular Materials[END_REF] disregard time and displacement scales smaller than particle rearrangements, and reformulate the equations of motion as well as the contact laws without regularization in terms of velocity jumps and impulses (the integral of a force over one timestep), which are calculated using an implicit integration scheme.

Note that these methods are generally applied to dense granular assemblies, while other methods exist to account for extremely dense or dilute granular assemblies [START_REF] Radjaï | Discrete-element Modeling of Granular Materials[END_REF].

Since DEM simulation has begun, granular media have mainly been modelled as a collection of interacting disks or spheres [START_REF] Torquato | Is Random Close Packing of Spheres Well Defined ?[END_REF][START_REF] Makse | Packing of compressible granular materials[END_REF][START_REF] Silbert | Granular flow down an inclined plane[END_REF][START_REF] Zhang | A Simulation Study of the Effects of Dynamic Variables on the Packing of Spheres[END_REF][START_REF] Silbert | Geometry in Frictionless and Friction Sphere Packings[END_REF][START_REF] Da Cruz | Rheophysics of dense granular flows: Discrete simulation of plane shear flows[END_REF][START_REF] Agnolin | Internal States of Model Isotropic Granular Packings. i. Assembling Process, Geometry and Contact Networks[END_REF][START_REF] Farr | Close packing density of polydisperse hard spheres[END_REF].

However in recent years, DEM simulation of granular assemblies made of nonspherical particles (non-spherical DEM) has gained substantial interest [START_REF] Lu | Discrete element models for nonspherical particle systems: From theoretical developments to applications[END_REF].

This should be no surprise if one realizes that at least 70% of the raw materials feeding modern industries consist of such particles [START_REF] Zhong | Dem/cfd-dem modelling of non-spherical particulate systems: Theoretical developments and applications[END_REF]. Interestingly, nonspherical DEM has revealed significant differences from known sphere collection behaviour, which could to some extent be more representative of the behaviour of granular media. In the static regime, packings of flat and/or elongated particles tend to exhibit orientation order, which affects the force network and increases the solid fraction compared to sphere packings [START_REF] Delaney | The packing properties of superellipsoids[END_REF][START_REF] Camenen | Geometrical properties of rigid frictionless granular packings as a function of particle size and shape[END_REF]. In the flowing regime, significantly more energy is needed to trigger and sustain avalanches of dense assemblies of angular particles compared to spheres [START_REF] Wachs | Grains3d, a flexible dem approach for particles of arbitrary convex shape part i: Numerical model and validations[END_REF][START_REF] Azéma | Discrete simulation of dense flows of polyhedral grains down a rough inclined plane[END_REF].

As such, non-spherical DEM is far more challenging than DEM simulation performed with spherical particles (spherical-DEM) for several reasons [START_REF] Lu | Discrete element models for nonspherical particle systems: From theoretical developments to applications[END_REF].

First, the occurence of a contact between non-spherical particles depends on their orientation, which fluctuates as a function of time and shall hence be systematically recalculated. This difficulty increases when particle orientation is determined by a rotation matrix parameterized by Euler's angles, since such a matrix may be ill-defined [START_REF] Džiugys | An approach to simulate the motion of spherical and non-spherical fuel particles in combustion chambers[END_REF]. Second, contact between two non-spherical particles may not be a single point, nor located on the segment connecting their inertia centres, thus requiring more complex contact detection algorithms to determine both the contact frame and location needed respectively for contact forces and torques calculation. A third challenging issue may be added upon observing that the inertia of a non-spherical particle is a frame-dependent matrix, which causes an infortunate coupling between particle orientation and angular velocity when solving the torques equations of motion [START_REF] Allen | Computer simulation of liquids[END_REF].

Several non-spherical DEMs have been suggested and tested to circumvent these issues [START_REF] Lu | Discrete element models for nonspherical particle systems: From theoretical developments to applications[END_REF]. The most frequently used, known as the composite spheres method, consists in modeling a flat and/or elongated particle as a sphere cluster. This way, interparticle contacts are simply contacts between one or more couples of spheres belonging to different particles, and such contacts may be fully characterized by spherical DEMs. Another widely used method interpolates the shape of any particle with various aspect ratio and blockiness using a super -quadric, which is defined by a cartesian equation and may be seen as a generalized ellipsoid. Contact detection algorithms take advantage of the cartesian equation from which particle surface normal may easily be calculated.

Though far less frequent, particle shape may also be interpolated using polyhedra, either with sharp (true polyhedra) [START_REF] Radjaï | Contact dynamics as a nonsmooth discrete element method[END_REF] or rounded [START_REF] Wachs | Grains3d, a flexible dem approach for particles of arbitrary convex shape part i: Numerical model and validations[END_REF] edges, for which contact detection algorithms accounting for the presence of edges and flat surfaces shall be used.

Yet, non-spherical DEMs still face major difficulties to accurately account for true particle shape, with sharp edges and/or flat surfaces. The composite spheres method requires a large number of spheres of various diameters to build flat surfaces, which impacts the calculation efficiency, and still such clusters have failed to mimic the flow behaviour (e.g. arching effect) of true polyhedra in hopper discharge simulations [START_REF] Höhner | A numerical study on the influence of particle shape on hopper discharge within the polyhedral and multisphere discrete element method[END_REF]. Similarly, approaches implementing super-quadrics are not well-suited to mimic non-symmetrical particles with sharp edges, and they are known to yield computationally expensive simulations when designed to mimic flat or elongated particles [START_REF] Lu | Discrete element models for nonspherical particle systems: From theoretical developments to applications[END_REF]. Approaches implementing polyhedra appear more adequate, however some do not account for multiple contacts between two particles (e.g. edge/f ace contact) [START_REF] Wachs | Grains3d, a flexible dem approach for particles of arbitrary convex shape part i: Numerical model and validations[END_REF], whereas others [START_REF] Radjaï | Contact dynamics as a nonsmooth discrete element method[END_REF][START_REF] Radjaï | Discrete-element Modeling of Granular Materials[END_REF] may be biased by the use of Euler's angles and the common plane contact detection algorithm [START_REF] Descantes | Classical contact detection algorithms for 3d dem simulations: Drawbacks and solutions[END_REF].

The objective of the present paper is to introduce a new MD-like approach to simulate the behaviour of dense assemblies of true polyhedra. Its main benefits are the absence of edge or corner rounding and the management of multipoint f ace -f ace or edge -f ace contacts between two convex polyhedra; the use of an improved version of the Gilbert-Johnson-Keerthi algorithm for contact detection between convex polyhedra; and a simple yet efficient implementation of an improved leap-frog Verlet non-iterative method suggested by Omelyan [START_REF] Omelyan | A new leapfrog integrator of rotational motion. the revised angular-momentum approach[END_REF] to solve the torques equations of motion, in which no periodic renormalization of the quaternions is necessary. The paper is organized as follows. Section 2 outlines the main features of the new MD-like approach, which is assessed both in the static and flowing regimes in section 3. Assessment results are then discussed in section 4 and finally perspectives are drawn in section 5.

Main features of the new approach

Equations of motion

In the absence of physico-chemical reactions, a granular medium may be seen as a collection of N rigid particles of various sizes and shapes interacting through field and contact forces. As a consequence, the medium deformations derive from particles motion, each of which obeys Newton's second law for particle translation and rotation. When field forces reduce to gravity, Newton's second law yields the following equations:

                     m i d - → v i dt = - → F g i + ci j=1 αij α=1 - → F α ji i=1,. . . , N d - → J i dt = d(I i - → ω i ) dt = ci j=1 αij α=1 --→ G i P α ji ∧ - → F α ji (1)
where m i , I i , G i , -→ v i , -→ ω i and c i stand for the mass, inertia tensor, centre of mass and its velocity vector, angular velocity vector, and number of contacting neighbours of particle i respectively; -→ F g i and -→ J i stand for the gravity force and angular momenta vectors respectively; eventually P α ji and -→ F α ji stand for each of the α ij contacting points and their corresponding force vectors respectively (exerted by neighbouring particle j on particle i). For each particle i, the unknowns in equations system 1 are the centre of mass location and velocity vectors, the particle orientation and its rotational velocity vector. Note that unlike spheres, two contacting polyhedra may share more than one contact point.

Contact forces model

Modeling contact forces is all the more difficult that their measurement in three-dimensional static, opaque, dense assemblies of rigid frictional particles is still being investigated namely through X-ray or Neutron tomography/diffraction [START_REF] Hurley | Quantifying interparticle forces and heterogeneity in 3d granular materials[END_REF]. As a consequence, using a simple contact forces model consistent with basic physical principles seems preferable to more sophisticated models. Despite the rigid nature of true particles, micro-asperities present at their surface are likely to undergo elastic or plastic deformations. Hence, it is convenient that contact force models in granular media account for more or less elastic and dissipative collisions as well as frictional contacts between particles.

In classical models such as Hertz or Kelvin-Voigt, this is achieved by allowing rigid particles to slightly overlap, so that elastic-repulsive contact forces are simply proportional to a power function of the overlap depth. The Kelvin-Voigt model further combines such forces with parallel viscous forces opposing the relative motion between colliding particles, so that collisions dissipate energy in accordance with Newton's second law of partial energy restitution. Besides, the magnitude of contact forces in the plane perpendicular to the direction of maximum overlap (tangential plane) is limited by Coulomb friction. Hence, a tentative visco-elastic model to determine the contact forces between particles i and j could be:

           R n ji = min(0, k n ij δ ij -ν n ij δij ) along n ij R t ji = -min(-µR n ji , ν t ij || U ij -δij n ij ||) along t ij R s ji = 0 along s ij (2)
where (R n ji ,R t ji ,R s ji ) respectively stand for the normal and tangential components of the contact force R ji exerted by particle j on particle i, expressed in the local contact frame defined by the direction of maximum overlap n ij (normal direction, oriented from particle i to j), the direction t ij of the projection of the relative velocity at contact U ij onto the tangential plane ( U ij is the velocity of particle i relative to particle j), and s ij = n ij ∧ t ij ; δ ij and δij = U ij . n ij stand for the maximum overlap (in terms of absolute value since δ ij < 0) and its time derivative respectively, while k n ij , µ, ν n ij and ν t ij respectively stand for the contact normal elastic stiffness, Coulomb friction coefficient, and normal and tangential viscosities.

Note however that equations system 2 is only valid for couples of particles sharing a unique contact point, where the maximum overlap depth and relative velocity vector needed to determine t ij are calculated. As shown on Figure 1, this is the case for vertex/face and edge/edge types of contact between polyhedra, whereas edge/face and face/face types of contact require at least two and three non-aligned contact points respectively to be fully determined.

To sort this issue out, an arbitrary unit vector t ij of the tangential plane is considered to build the contact frame ( t ij , s ij , n ij ) and the following visco-elastic model is ultimately adopted to calculate the contact forces between particles i and j at each contact point α:

                   R αn ji = 1 α ij min(0, k n ij δ ij -ν n ij δα ij ) along n ij R αt ji = -min( ν t ij α ij , -µR αn ji || U α ij -δα ij n ij || ) U α ij . t ij along t ij R αs ji = -min( ν t ij α ij , -µR αn ji || U α ij -δα ij n ij || ) U α ij . s ij along s ij (3)
Observe that equation systems 2 and 3 are equivalent when α ij = 1 and

t ij = Uij || Uij || with || U ij || > 0.
Furthermore, note that each contact is assigned the maximum overlap depth between particles i and j to avoid unnecessary though timeconsuming overlap depth calculation at each contact. Besides, would the relative velocity tangential component at contact α be zero, then the term

-µR αn ji || U α ij -δα ij nij || would tend to +∞, yielding R αt ji = - ν t ij αij U α ij . t ij = R αs ji = - ν t ij αij U α ij . s ij = 0.
Last, contact force vectors in the laboratory and inertia frames relate through

F α ji = M ij R α ji
, where M ij is a 3 × 3 matrix whose columns are the coordinates of vectors ( t ij , s ij , n ij ) expressed in the laboratory frame.

Contact detection strategy and contacts location 2.3.1. Contact detection strategy

As evidenced by equations system 3, the maximum overlap between each pair of contacting particles shall first be determined to allow contact forces calculation. This implies that any contact between two polyhedra shall be identified and geometrically characterized to determine the overlap and contact frame.

Since checking the existence of a contact between any set of two particles in an assembly of N particles has a prohibitive calculation cost (proportional to N 2 ), an alternative strategy was adopted. This strategy combines as follows the twosteps gross detection algorithm suggested by Ogarko and Luding [START_REF] Ogarko | A fast multilevel algorithm for contact detection of arbitrarily polydisperse objects[END_REF] with a fine contact detection method adapted to polyhedra: 1) map each of the N particles into a cell of an optimized set of nested hierarchical grids based on its size and location in the granular assembly; 2) determine the list of neighbours of each particle by identifying in the cell to which it is mapped and in its surrounding cells the polyhedra whose bounding spheres intersect that of the particle; 3) use the modified version of the GJK algorithm suggested by Descantes et al [START_REF] Descantes | Classical contact detection algorithms for 3d dem simulations: Drawbacks and solutions[END_REF] to check the occurence of a contact between each polyhedron and its neighbours and, in this case, determine the contact geometrical characteristics.

Steps 1 and 2 of this strategy differ from DEM classically used Verlet neighbour list and linked-cell methods [START_REF] Allen | Computer simulation of liquids[END_REF]. The Verlet list is built by periodically searching the list of possible neighbours of each particle among all others, which requires performing only periodically N 2 binary tests but nevertheless tends to be prohibitive when N is larger than a few thousand particles. The linked-cell method maps each particle to a unique grid whose cell size is at least larger that the largest particle, and then it searches the list of neighbours of each particle in the cell in which it is mapped and the surrounding cells. However in the case of highly polydisperse dense granular assemblies, using a unique grid yields prohibitive calculation time since hundreds of small particles may occupy the same cell in which tens of thousand binary tests will hence have to be performed. In contrast, the nested hierarchical cell space method suggested by Ogarko and Luding [START_REF] Ogarko | A fast multilevel algorithm for contact detection of arbitrarily polydisperse objects[END_REF] may be seen as a multi-grid linked-cell method, in which the number of grids and their cell sizes self-customize to the granular polydispersity so that the method performs better (CPU time for gross contact detection scales linearly with number of particles) than both the Verlet list and linked-cell methods. For this purpose, small particles are mapped to a grid with a small cell size, large particles are mapped to a grid with a larger cell size, and a first neighbour search is performed between particles of approximately the same size using the appropriate grid. Then a cross-level search is performed between large and small particles using the small cell size grid. Note that the nested hierarchical cell space method performs at least as well as if not better than other methods, such as computer graphics originated octrees whose CPU time for gross contact detection scales as O(N ) to O(N logN ) [START_REF] Vemuri | Fast collision detection algorithms with applications to particle flow[END_REF] Step 3 of this strategy is dedicated to fine contact detection between convex polyhedra, since any (non-convex) polyhedron may be built by assembling convex polyhedra. Though simple between two spheres -contact occurs when their centres are separated by a distance smaller than the sum of their radii and the maximum overlap occurs along the line joining their centres -maximum overlap and contact frame determination is much more complex and timeconsuming in the case of polyhedra. Several authors have investigated this issue and suggested original algorithms in particular in the frameworks of rock engineering [START_REF] Cundall | Formulation of a three-dimensional 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 algorithm for 3-d discrete element method[END_REF][START_REF] Boon | A new algorithm for contact detection between convex polygonal and polyhedral particles in the discrete element method[END_REF], robotics [START_REF] Lin | A fast algorithm for incremental distance cal-culation[END_REF] and computer graphics [START_REF] Gilbert | A fast procedure for computing the distance between complex objects in three-dimensional space[END_REF][START_REF] Mirtich | V-clip: Fast and robust polyhedral collision detection[END_REF]. Upon investigating these algorithms and their drawbacks, the GJK algorithm initially designed by Gilbert et al [START_REF] Gilbert | A fast procedure for computing the distance between complex objects in three-dimensional space[END_REF] and further modified by Descantes et al [START_REF] Descantes | Classical contact detection algorithms for 3d dem simulations: Drawbacks and solutions[END_REF] was implemented.

GJK algorithm looks for the closest feature (point, edge or face) between neighbouring polyhedra A and B (see Figure 2a), which is equivalent to finding the feature of the convex hull of their Minkowski difference A -B (see Figure 2b) that minimizes the distance to the coordinates origin O. Interestingly, observe on Figure 2b that when the polyhedra overlap, the coordinates origin is located inside the convex hull of the Minkowski difference A -B, so that point P of the convex hull located closest to the coordinates origin defines both the maximum overlap depth || --→ OP || and contact normal n = --→ OP of the overlapping polyhedra (pointing from polyhedron A to polyhedron B). Practically, point P is sought by building step by step a sequence of simplices whose vertices belong to the convex hull of the Minkowski difference, which converges to the simplex closest to the coordinates origin. When the polyhedra overlap, this simplex is a tetrahedron which encloses the coordinates origin and whose face point closest to the coordinates origin is taken as point P (see Figure 2c).

Finally, Descantes et al [START_REF] Descantes | Classical contact detection algorithms for 3d dem simulations: Drawbacks and solutions[END_REF] have investigated numerical unstabilities of the GJK algorithm reported by several authors [START_REF] Van Den Bergen | A fast and robust gjk implementation for collision detection of convex objects[END_REF][START_REF] Montanari | Improving the gjk algorithm for faster and more reliable distance queries between convex objects[END_REF], and they have observed that these unstabilities were likely caused by point P not belonging to the convex hull of the Minkowsky difference. Hence they have suggested a revised algorithm called GJK -T D, which rigorously relocates point P on the convex hull of the Minkowsky difference A -B of the polyhedra, so that overlap depth || --→ OP || and contact normal --→ OP of the overlapping polyhedra are correct.

Contacts location

Once the contact frame of each pair of overlapping polyhedra is known, determining the set of contact points P α ji is necessary in order to solve the second equation of system 1 for particles angular velocities ω i . Since the early 1990s, this contact point location issue has been discussed by some authors, mostly with the aim of achieving high computational efficiency using a single contact point while respecting the physics of granular media. Hogue and Newland [START_REF] Hogue | Efficient computer simulation of moving granular particles[END_REF] have investigated contact location between overlaping polygons and they reported physically correct results using the most burried vertex of one polygon into the other. Džiugys and Peters [START_REF] Džiugys | An approach to simulate the motion of spherical and non-spherical fuel particles in combustion chambers[END_REF] mention two methods, one in which the contact point is located in the middle of the overlap area and the other in which it is defined as the intersection between a line and a surface inside the overlap area. They notice that the second method is biaised in case one particle collides with a deformable wall. Wachs et al [START_REF] Wachs | Grains3d, a flexible dem approach for particles of arbitrary convex shape part i: Numerical model and validations[END_REF] round the corners and edges of their polyhedra with a curvature radius of magnitude r in order to obtain a unique contact point located in the middle of the overlap area, and they report physically correct behaviours in a rotating drum as well as good stability in time of a vertical column made of three cubes in face/face contact under their own weight.

More recently,Gay-Neto and Wriggers [START_REF] Neto | Discrete element model for general polyhedra[END_REF] observed that face/face and edge/face contact situations between polyhedra were improbable, hence giving them special treatment would be physically relevant without significantly extending the simulation duration. These authors suggested to approximate continuous edge/face and face/face contact between a pair of polyhedra by a set of point-wise contact interactions matching the singularities of the edges or faces (vertices located in the overlaping area and edges intersections). For this purpose, each particle face is split into triangular sub-regions, which are tracked for contact with other triangular sub-regions of neighboring particles. The same idea was implemented in the new approach as follows.

Consider two overlapping convex polyhedra, say A and B, whose faces are split into triangular sub-regions (figure 3a). Upon application of the GJK -T D algorithm [START_REF] Descantes | Classical contact detection algorithms for 3d dem simulations: Drawbacks and solutions[END_REF], point P is determined as the barycentre of the triangular face located closest to the coordinates origin on the convex hull of the Minkowsky difference A -B (see Fig. 2c). Note that applying the barycentric coordinates of point P = P A -P B to the vertices of the closest features between polyhedra A and B determines the closest points P A and P B between these polyhedra.

As a consequence, an affine contact plane is defined by its normal --→ OP and the midpoint of segment [P A , P B ]. In case of a face/face contact, the triangular sub-regions of both contacting faces are then projected on this affine plane (figure 3b), and singular points resulting from the intersection of two triangular sub-regions belonging to each of the polyhedra are identified and stored as pointwise contact interactions. The same applies to an edge/face contact, for which point-wise contact interactions consist of singular interaction points between the projection of each triangular sub-region of the face and the projection of the edge on the contact plane. Note that vertex/face and edge/edge types of contact (see Fig. 1) consist of a couple of contact points (P A ;P B ) determined by algorithm GJK -T D.

Particle-fixed frame to space-fixed frame rotation matrix

In dense granular assemblies, most particles share contacts with several neighbouring particles, characterized by different contact frames. Obviously, equations system 1 needs to be solved in a frame common to all contacts between particle i and its neighbours, namely a particle-fixed or space-fixed frame.

The particle inertia frame, which is a particle-fixed frame centered at the polyhedron center of inertia with its axes aligned with the axes of inertia of the polyhedron, is particularly convenient since in this frame the inertia tensor is diagonal and time-independent. Yet, this frame needs to be used in conjunction with a space-fixed frame common to all particles (laboratory frame), since calculating contact forces from equations system 3 requires determining the relative velocity of two contacting particles. A corollary is the need to determine the transformation between particle inertia and laboratory frames.

This transformation takes the form of a rotation matrix which is skewsymmetric and orthonormal, hence fully determined by three independent parameters. Among others, Džiugys and Peters [START_REF] Džiugys | An approach to simulate the motion of spherical and non-spherical fuel particles in combustion chambers[END_REF] have evidenced with details that classically used Euler's angles were unsuitable to parameterize this rotation matrix, since they cause particle orientation indeterminacy at specific Euler's angles values. In contrast, quaternions [START_REF] Hamilton | On quaternions[END_REF] are known as a reliable alternative to parameterize a singularity-free rotation matrix. Basically, a quaternion q is a set of four scalar variables which may define a rotation of a given angle around a fixed axis in three-dimensional space. According to Euler's rotation theorem, such a rotation is equivalent to any combination of rotations of a rigid body around a fixed point (e.g. its centre of inertia) in three-dimensional space. The quaternion-parameterized expression of the rotation matrix A i from particle i inertia frame to laboratory frame, whose columns stand for the components of the laboratory frame expressed in the inertia frame, is summarized e.g. by Allen

and Tildesley [START_REF] Allen | Computer simulation of liquids[END_REF]:

A i =      q 2 i0 + q 2 i1 -q 2 i2 -q 2 i3 2(q i1 q i2 + q i0 q i3 ) 2(q i1 q i3 -q i0 q i2 ) 2(q i1 q i2 -q i0 q i3 ) q 2 i0 -q 2 i1 + q 2 i2 -q 2 i3 2(q i2 q i3 + q i0 q i1 ) 2(q i1 q i3 + q i0 q i2 ) 2(q i2 q i3 -q i0 q i1 ) q 2 i0 -q 2 i1 -q 2 i2 + q 2 i3      (4)
where q i = (q i0 , q i1 , q i2 , q i3 ) T with q 2 i0 + q 2 i1 + q 2 i2 + q 2 i3 = 1 for a rotation quaternion. Note that a rotation quaternion also writes q i = (cos ζi 2 , a i sin ζi 2 , b i sin ζi 2 , c i sin ζi 2 ) T with ζ i standing for the angle of rotation around the axis defined by its unit vector (a i b i c i ) T . Furthermore, observe that a rotation quaternion relates to Euler's angles (θ i , φ i , ψ i ) according to the following set of equations:

                     q i0 = cos θ i 2 cos φ i + ψ i 2 q i1 = sin θ i 2 cos φ i -ψ i 2 q i2 = sin θ i 2 sin φ i -ψ i 2 q i3 = cos θ i 2 sin φ i + ψ i 2 (5)

Time integration scheme

Solving equations system 1 for particles location, orientation and velocity vectors requires an efficient time integration scheme. Several authors [START_REF] Allen | Computer simulation of liquids[END_REF][START_REF] Džiugys | An approach to simulate the motion of spherical and non-spherical fuel particles in combustion chambers[END_REF] have reviewed the advantages and drawbacks of the most commonly used schemes in DEM. Among these, the leap-frog Verlet algorithm was implemented for both equations due to its simplicity, stability and second order accuracy. Using this algorithm, equations system 1 rewrites as the following system of two pairs of equations at time t:

                             - → v i (t + ∆t 2 ) = - → v i (t - ∆t 2 ) + ∆t m i [ - → F g i (t) + αij α=1 - → F α ji (t)] + O(∆t 2 ) - → r i (t + ∆t) = - → r i (t) + ∆t. - → v i (t + ∆t 2 ) + O(∆t 2 ) i=1,. . . , N - → J i (t + ∆t 2 ) = - → J i (t - ∆t 2 ) + ∆t. αij α=1 ---→ G i P α ji (t) ∧ - → F α ji (t) + O(∆t 2 ) - → q i (t + ∆t) = - → q i (t) + ∆t. - → q i (t + ∆t 2 ) + O(∆t 2 ) (6)
where the first pair of equations determines the velocity vector v i at time t + ∆t 2 and location vector r i at time t + ∆t of the centre of mass of each polyhedron i from kinematic variables known at previous time step together with gravity and contact forces calculated at time t. The second pair of equations determines the angular momenta J i at time t + ∆t 2 and orientation quaternion q i at time t + ∆t of each polyhedron from kinematic variables known at previous time step and contact torques calculated at time t. Note that -→ q i stands for the time derivative of quaternion q i , which relates to angular velocity -→ ω p i = (ω p iX ω p iY ω p iZ ) T expressed in the particle inertia frame (denoted by superscript p) according to the following [START_REF] Allen | Computer simulation of liquids[END_REF]:

- → q i (t) = 1 2         q i0 (t) -q i1 (t) -q i2 (t) -q i3 (t) q i1 (t) q i0 (t) q i3 (t) -q i2 (t) q i2 (t) -q i3 (t) q i0 (t) q i1 (t) q i3 (t) q i2 (t) -q i1 (t) q i0 (t)                 0 ω p iX (t) ω p iY (t) ω p iZ (t)         (7)
Calculating the centre of gravity kinematic variables from system 6 first pair of equations associated with boundary and initial conditions raises no particular issue. Unfortunately, determining particles orientation and angular velocity vectors from the second pair of equations raises a major issue, since these two variables appear coupled in the expression of the quaternion time derivative (equation 7). To sort this issue out, Allen and Tildesley [START_REF] Allen | Computer simulation of liquids[END_REF] suggest using a modified leap-frog Verlet algorithm designed by Potter [START_REF] Potter | Computational physics[END_REF], but this algorithm involves guesses of angular momenta and time derivative of quaternions which cause bias. In contrast, Omelyan [START_REF] Omelyan | A new leapfrog integrator of rotational motion. the revised angular-momentum approach[END_REF] suggests an unbiased revised leapfrom algorithm for rotational motion which was implemented in the new approach.

This revised leapfrog algorithm is obtained by substituting rotation matrix A i for quaternion -→ q i in the second pair of equations of system 6 which rewrites:

                 - → J i (t + ∆t 2 ) = - → J i (t - ∆t 2 ) + ∆t. αij α=1 ---→ G i P α ji (t) ∧ - → F α ji (t) + O(∆t 2 ) i=1,. . . , N A i (t + ∆t) = A i (t) + ∆t. dA i (t + ∆t 2 ) dt + O(∆t 2 ) (8)
where:

dA i (t) dt = W [ - → ω p i (t)]A i (t) =      0 ω p iZ (t) -ω p iY (t) -ω p iZ (t) 0 ω p iX (t) ω p iY (t) -ω p iX (t) 0      A i (t) (9)
Then, upon observing that rotation matrix A i relates the angular momenta expressions in the inertia and laboratory frames as

I p i - → ω p i (t) = A i (t) - → J i (t)
i (t + ∆t 2 ) = 1 2 [A i (t) + A i (t + ∆t)], system 8 rewrites:                          - → J i (t + ∆t 2 ) = - → J i (t - ∆t 2 ) + ∆t. αij α=1 ---→ G i P α ji (t) ∧ - → F α ji (t) + O(∆t 2 ) i=1,. . . , N - → ω p i (t + ∆t 2 ) = (I p i ) -1 A i (t) - → J i (t + ∆t 2 ) + ∆t 2 (I p i ) -1 W [ - → ω p i (t + ∆t 2 )]
.

I p i - → ω p i (t + ∆t 2 ) + O(∆t 2 ) ( 10
)
Observe that the second equation of system 10 is implicit in -→ ω p i (t + ∆t 2 ), its unique unknown, -→ J i (t + ∆t 2 ) being calculated by the first equation. As a consequence, the second equation may easily be solved iteratively for -→ ω p i (t+ ∆t 2 ) upon initializing this vector to (

I p i ) -1 A i (t) - → J i (t + ∆t 2 )
in the right-hand side expression and calculating k successive vectors -→ ω

p(k) i (t + ∆t 2 ) until || - → ω p(k+1) i (t + ∆t 2 ) - - → ω p(k ) i (t + ∆t 2 )|| < .
In fact, Omelyan also suggests a non-iterative method to solve this equation. Both methods have been implemented in the new code and readers interested by the calculation details are invited to read ref. [START_REF] Omelyan | A new leapfrog integrator of rotational motion. the revised angular-momentum approach[END_REF].

An ultimate feature of Omelyan's revised leapfrog algorithm shall be underlined.

When interpolating

A i (t+ ∆t 2 ) by 1 2 [A i (t)+A i (t+∆t)
] and introducing equation 9 into the second equation of system 8, it is easy to come up with the following expression of A i (t + ∆t) as a function of A i (t):

A i (t + ∆t) = {I - ∆t 2 W [ - → ω p i (t + ∆t 2 )]} -1 {I + ∆t 2 W [ - → ω p i (t + ∆t 2 )]}A i (t) = (I -Θ i ) -1 (I + Θ i )A i (t) (11)
where I stands for the identity matrix. As stated by Omelyan, since W is skewsymmetric, Θ i is skewsymmetric and it is easy to see that (I -Θ i ) -1 (I +Θ i ) is orthonormal. As a consequence, if A i (t) is initially orthonormal, then it will remain so by construction at any further timestep. This is a valuable advantage over the use of quaternions, since these have to be frequently renormalized to avoid cumulative orientation errors as underlined by Allen and Tildesley [START_REF] Allen | Computer simulation of liquids[END_REF].

Validation of the new approach

In order to assess the new discrete element approach potential, simulations have been carried out. The following two classical problems have been simulated, for which numerical results with true polyhedra exist in the literature: [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF] the compaction of a biperiodic assembly of frictionless convex polyhedra under their own weight for comparison with the results reported by Camenen and

Descantes [START_REF] Camenen | Geometrical properties of rigid frictionless granular packings as a function of particle size and shape[END_REF] in particular; (2) the flow behaviour of frictional convex polyhedra down an inclined plane for comparison with the results reported by Azéma et al [START_REF] Azéma | Discrete simulation of dense flows of polyhedral grains down a rough inclined plane[END_REF] among others.

Compaction of frictionless polyhedra under their own weight

This simulation was inspired by numerous authors who studied the compaction of particle assemblies under their own weight [START_REF] Macrae | Significance of The Properties of Materials in the Packing of Real Spherical Particles[END_REF][START_REF] Zhang | A Simulation Study of the Effects of Dynamic Variables on the Packing of Spheres[END_REF][START_REF] Silbert | Geometry in Frictionless and Friction Sphere Packings[END_REF], and particularly by Camenen and Descantes [START_REF] Camenen | Geometrical properties of rigid frictionless granular packings as a function of particle size and shape[END_REF] who investigated the ability of frictionless pinacoids -a variety of convex polyhedra comprised of eight vertices, fourteen edges, and eight faces as shown in Fig. 4 -to achieve their random close packed state (RCP) when randomly dropped onto a bottom wall and left to relax under the action of their own weight. The RCP state of granular particles, also defined as their maximally randomly jammed state [START_REF] Torquato | Is Random Close Packing of Spheres Well Defined ?[END_REF], corresponds to the stable mechanical equilibrium state of rigid, frictionless grains under an isotropic confining pressure with no traces of crystallization [START_REF] Agnolin | Internal States of Model Isotropic Granular Packings. i. Assembling Process, Geometry and Contact Networks[END_REF]. The objective here is to confirm that the new discrete element approach allows to achieve such a state, then to determine the packing solid fraction as well as microstructural characteristics, and finally to confront these results with those reported in particular by Camenen and Descantes [START_REF] Camenen | Geometrical properties of rigid frictionless granular packings as a function of particle size and shape[END_REF]. oriented isometric pinacoid was introduced into each sphere so that the sphere is circumscribed to the pinacoid. Finally, the spherical shells were removed, gravity was set downward, periodic boundary conditions were substituted for the four lateral walls and the loose pinacoid assembly was left to relax under its own weight in the absence of Coulomb friction until a steady-state is achieved.

Figure 5 shows the pinacoid assembly at initial (loose) and final (dense) stages.

Although it has been shown that, in the absence of Coulomb friction, a granular assembly will achieve its RCP state regardless of the densification process (e.g. regardless of the falling height of individual particles) [START_REF] Agnolin | Internal States of Model Isotropic Granular Packings. i. Assembling Process, Geometry and Contact Networks[END_REF], some damping is necessary to mitigate particles velocities which would otherwise cause excessive overlap. Hence, the local non-viscous damping suggested by Potyondy and Cundall [START_REF] Potyondy | A bonded-particle model for rock[END_REF] was applied to each grain up to roughly 6.10 6 timesteps. This damping takes the form of force and torque components, which are proportional to the magnitude of the unbalanced force/torque components applied to each particle and oppose its velocity components. In addition, various combinations of normal stiffness and viscosity at contacts were tested in order to examine their effect on particles overlap and kinetic energy, until a combination of these parameters was ultimately selected to finalize the densification process.

Table 1 gathers the simulated system parameters. Note that 1 has been assigned to both the density ρ and gravity g, so that all parameters are fully determined upon assigning the diameter d sph of the spherical shell circumscribed to the mean pinacoid. In addition, the mechanical parameters have been fixed upon calculating the inertia characteristic time T inertia from both the gravity and diameter d eq of the sphere of same volume as the mean pinacoid, and then separating the inertia, stiffness T stif f ness and simulation ∆t time scales.

Besides, note that only the final combination of mechanical parameters is summarized in the table (other combinations will be reported in section 3.1.2). The normal viscosity γ n was taken critical to mitigate non-realistic oscillations at contacts, whereas both the tangential viscosity γ t and Coulomb friction µ were set to zero. The non-viscous damping coefficient α P ot was fixed according to Potyondy and Cundall [START_REF] Potyondy | A bonded-particle model for rock[END_REF]. A unique simulation running over 14.10 6 timesteps was performed.

Compaction results

The results presentation illustrates first the quality of mechanical equilibrium achieved by the simulation, then focuses on macroscopic properties of the packing, and finally investigates its microstructure. timesteps, as a consequence of particles undergoing dissipative collisions while they accelerate from initially loose state to achieve a dense state (see kinetic energy variations on Fig. 6b). Note that the mean overlap never exceeds 10 -3 d sph and varies inversely with the mean contact stiffness k mean n as recalled by Wachs et al [START_REF] Wachs | Grains3d, a flexible dem approach for particles of arbitrary convex shape part i: Numerical model and validations[END_REF]. Besides, the mean overlap stabilizes at roughly 3 × 10 -5 d sph in zone 4 of the simulation, which is almost two orders of magnitude better than the results reported by Camenen and Descantes [START_REF] Camenen | Geometrical properties of rigid frictionless granular packings as a function of particle size and shape[END_REF]. Last, observe on Fig. 6b that at the end of the simulation, the total kinetic energy stabilizes in the range 10 -8 E p to 10 -7 E p , which is similar to the results reported by Camenen and Descantes [START_REF] Camenen | Geometrical properties of rigid frictionless granular packings as a function of particle size and shape[END_REF].

Macroscopic properties have been calculated once the mechanical equilibrium is reached, that is in zone 4 of the simulation when the total kinetic energy stabilizes to its minimum. These properties have been averaged over N configurations spread every 2000 timesteps.

Figure 7 shows the mean solid fraction profile along the z -axis. This profile was calculated from N = 5 configurations using the Monte-Carlo method with 10 5 shots randomly fired in the granular packing broken down into 0.1d eqthick horizontal slices. Away from the free surface and bottom wall, where particles layering is known to cause large solid fraction fluctuations visible up to z = 2d eq [START_REF] Suzuki | Study of the wall effect on particle packing structure using x-ray micro computed tomography[END_REF], the solid fraction is homogenous. The mean solid fraction calculated in the homogenous zone (between z = 3d eq and z = 6d eq inclusive) yields ν = 0.67, which compares very well with the results reported by Camenen and Descantes [START_REF] Camenen | Geometrical properties of rigid frictionless granular packings as a function of particle size and shape[END_REF].

Figure 8 depicts the profiles of four stress components along the z -axis, as well as that of the weight-to-packing surface area ratio. Each profile was calculated layer by layer using the symmetrized stress tensor expression suggested by Moreau [START_REF] Moreau | Numerical investigation of shear zones in granular materials[END_REF][START_REF] Cambou | Micromécanique des matériaux granulaires, Traité MIM, série géomatériaux[END_REF]:

σ(z) = 1 L x L y δz sym   z≤zn≤z+δz m∈T actn - → F mn ⊗ ---→ G n P mn   - 1 L x L y δz z≤zn≤z+δz [W ( - → ω n )W ( - → ω n ) + W ( - → ω n )]I n (12)
where z n stands for the z-coordinate of the centre of inertia of any particle n inside the [z; z + δz] layer, with G n , I n , -→ ω n and -→ ω n being the particle centre of inertia, inertia tensor, rotational velocity and acceleration vectors respectively;

T act n defines the set of particles m sharing a contact P mn with particle n, and W stands for the matrix expression of the rotational velocity and acceleration vectors respectively (see equation 9). In addition, each profile was averaged from N = 100 configurations. Figure 8 shows that, in the absence of Coulomb friction, the shear stresses equal zero (e.g. see σ yz ), the normal stresses σ xx , σ yy and σ zz are almost identical, and their mean trend superimposes with the weight profile. This isotropic stress state suggests that the granular packing behaves like a fluid at rest. Like the solid fraction profile, the normal stress profiles undergo large fluctuations up to z = 2d eq as a consequence of the particles layering close to the bottom wall.

Ultimately, the packing microstructure was investigated upon achievement of the mechanical equilibrium, in order to verify the absence of crystallization and check the contacts network.

The absence of translational arrangement was checked using the classical pair correlation function g(r) [START_REF] Allen | Computer simulation of liquids[END_REF], whereas the absence of orientational order of the pinacoids was verified upon calculating the nematic order parameter Q 2 00 [START_REF] Allen | Computer simulation of liquids[END_REF][START_REF] Camenen | Effect of confinement on dense packings of rigid frictionless spheres and polyhedra[END_REF]. These function and parameter were averaged from N = 5 configurations. Figure 9a displaying the pair correlation function fluctuations confirms the absence of translational order, since g(r) becomes homogenous and clause to 1 beyond 4R min . With Q 2 00 = 0.18, no orientational order has been observed either. These conclusions are consistent with those of Camenen and Descantes [START_REF] Camenen | Geometrical properties of rigid frictionless granular packings as a function of particle size and shape[END_REF].

Figure 9b depicts various coordination number profiles along the z -axis, namely the total mean number of contacts per particle z tot as well as its distribution by contact type. These profiles look homogenous at least two particle diameters away from both the free surface and bottom wall. The mean value of the total coordination number in this central zone (between z = 3d eq and z = 6d eq inclusive) is z tot = 8.3, which is close to z tot = 8.4 found by Camenen and Descantes [START_REF] Camenen | Geometrical properties of rigid frictionless granular packings as a function of particle size and shape[END_REF]. Although this value is well above the classical threshold value of 6 contacts per particle -twice its number of degrees of freedom -corresponding to isostatic sphere assemblies [START_REF] Alexander | Amorphous solids: their structure, lattice dynamics and elasticity[END_REF][START_REF] Silbert | Geometry in Frictionless and Friction Sphere Packings[END_REF][START_REF] Agnolin | Internal States of Model Isotropic Granular Packings. i. Assembling Process, Geometry and Contact Networks[END_REF], it is irrelevant for checking the isostaticity of non-spherical particles [START_REF] Jaoshvili | Experiments on The Random Packing of Tetrahedral Dice[END_REF][START_REF] Jiao | Maximally random jammed packings of platonic solids: Hyperuniform long-range correlations and isostaticity[END_REF][START_REF] Camenen | Effect of confinement on dense packings of rigid frictionless spheres and polyhedra[END_REF]. Indeed, upon considering a particle located away from the free surface and bottom wall, and assigning 1, 2 and 3 constraints to its z s single-point contacts, z d edge/f ace contacts and z t f ace/f ace contacts respectively, the total number of constraints on this particle writes z iso = z s + 2z d + 3z t [START_REF] Camenen | Geometrical properties of rigid frictionless granular packings as a function of particle size and shape[END_REF]. Figure 9b shows the profile of z iso , which falls about 11 though short of twice the 6 degrees of freedom of the pinacoid. As a consequence, and by contrast with Camenen and Descantes [START_REF] Camenen | Geometrical properties of rigid frictionless granular packings as a function of particle size and shape[END_REF], the present pinacoid packing seems to be slightly hypostatic. Last, it is worth pointing out that roughly 76% of the contacts consist of a single point (z s , either vertex/f ace or edge/edge), whereas 21% consist of a particle edge (z d ) and 3% are f ace/f ace contacts (z t ). Even if f ace/f ace contacts remain scarce as pointed out by Wachs et al [START_REF] Wachs | Grains3d, a flexible dem approach for particles of arbitrary convex shape part i: Numerical model and validations[END_REF], these results justify considering contact situations with more than a single point. Besides, these results seem slightly different from those reported by Camenen and Descantes [START_REF] Camenen | Geometrical properties of rigid frictionless granular packings as a function of particle size and shape[END_REF], who noted respectively 66% of single point contacts, 25% of edge/f ace contacts and 9% of f ace/f ace contacts.

Dense flow of frictional convex polyhedra down an inclined plane

Dense flow of granular assemblies down an inclined plane has received wide experimental and numerical attention over the last decades [START_REF] Savage | The motion of a finite mass of granular material down a rough incline[END_REF][START_REF] Pouliquen | Scaling laws in granular flows down a rough inclined plane[END_REF][START_REF] Silbert | Granular flow down an inclined plane[END_REF][START_REF] Jop | A constitutive law for dense granular flows[END_REF]. Here, emphasis is placed on the ability of the new discrete element approach to replicate the steady flow simulation of frictional pinacoids (see Fig. 4) down a rough inclined plane as reported by Azéma et al [START_REF] Azéma | Discrete simulation of dense flows of polyhedral grains down a rough inclined plane[END_REF]. In particular, the range of tilt angles for which a steady flow is observed will be investigated and benchmarked against Azéma et al's findings in terms of solid fraction, mechanical stress state and microstructural characteristics.

Flow simulation procedure and parameters

In order to build an approximately two-layer thick rough substrate of pinacoids, 2681 spheres of diameter d sph ± 0.1d sph were geometrically deposited into a parallelepiped container of length L y = 30d sph and width L x = 26d sph using the same procedure as in section 3.1.1. Then a ramdomly oriented isometric pinacoid was placed into each sphere, and the spherical shells were subsequently discarded. Finally the gravity was set downwards, walls perpendicular to x and y axes were replaced by biperiodic boundary conditions and the pinacoids were left to relax under their own weight for 10 5 timesteps (the timestep was set to five times the value of ∆t used in the flow simulation, see table 2). Figure 10a shows the resulting rough substrate.

Next, the same procedure was used to geometrically deposit into the same container an assembly of 20000 randomly oriented pinacoids inscribed into spheres of diameters d sph ± 0.1d sph . This pinacoid assembly was further translated upward along the z -axis, so that the lowest particle is located above the highest particle of the rough substrate. Figure 10b shows the resulting initially loose assembly of pinacoids on top of the rough substrate.

Finally, gravity was applied in the yz-plane with an initial tilt angle θ = 50 • to trigger the flow, and the tilt angle was gradually reduced to the target value, the idea being to ease the achievement of a steady flow by strongly accelerating the 20000 pinacoid assembly in the flow direction y before it had fully relaxed on the fixed rough substrate. First, 768000 timesteps were simulated with a tilt angle θ = 50 • before the tilt angle was reduced to θ = 47.5 • and θ = 45 • , then another 214000 timesteps were run before the tilt angle was reduced from 45 • to 42 • and 39 • , then another 126000 timesteps were run before the tilt angle was reduced from 39 • to 37 • , and finally another 255000 timesteps were run before the tilt angle was reduced from 37 • to 36 • and 30 • . The tangential viscosity γ t was assigned the same critical value as the normal viscosity γ n and the Coulomb friction was set to µ = 0.4 in accordance with Azéma et al [START_REF] Azéma | Discrete simulation of dense flows of polyhedral grains down a rough inclined plane[END_REF]. The various simulations were run over 4.10 6 to 5.10 6 timesteps.

Flow results

The results presentation first checks the absence of excessive overlap and then attempts to identify the steady flow regime by examining how the tilt angle θ influences the profiles of various quantities in the yz-plane. These quantities are the packing structure (solid fraction and coordination number), grains velocity and stresses. These profiles were averaged both in 0.1d eq -thick horizontal layers and in time over N configurations spread every 1000 timesteps.

Figure 11a displays the mean overlap between particles as a function of time.

It can be concluded that the mean overlap never exceeds 10 -4 d sph whatever the tilt angle, which is reasonable [START_REF] Cleary | Dem prediction of industrial and geophysical particle flows[END_REF]. For all the tilt angles tested, Figure 12a depicts the solid fraction fluctuation profiles in the yz-plane, each averaged from N = 3 configurations. A polynomial fit of these profiles is also represented using symbols to help distinguish between the curves. As reported by Azéma et al [START_REF] Azéma | Discrete simulation of dense flows of polyhedral grains down a rough inclined plane[END_REF], these profiles evidence a vertical expansion of the packing when the tilt angle increases as a consequence of the Reynolds dilatancy. This expansion obviously results in a solid fraction decrease. Note that in a narrow range of tilt angle values bounded by [36 •

; 39 • ],
the solid fraction is reasonably uniform in the bulk (viz. ν 0.45 for z > 5d eq ) though 10% smaller than the figure reported by Azéma et al for polyhedra (ν = 0.5), and it increases close to the rough bottom in full agreement with these authors. However for tilt angles θ = 42 • and higher, the solid fraction profile tends to be increasingly reversed S -shaped, with its minimum reached a few particle diameters above the rough bottom (say z = 4d eq to z = 5d eq ). For each tilt angle tested, figure 13a shows the velocity profile averaged from N = 100 configurations. For tilt angles θ = 39 • and below, these profiles evidence no sliding at the rough bottom interface as reported by Azéma et al, whereas such a sliding occurs for tilt angles θ = 42 • and above. Furthermore, these profiles look convex close to the rough bottom and then concave in the bulk, but the convex zone tends to reduce with increasing θ and ultimately disappears for θ = 50 • . With v max y / gd eq 7.5 compared to 12, the maximum velocity for θ = 37 • looks 38% smaller than that reported by Azéma et al. The inset confirms the reasonably Bagnoldian shape of the profiles at least around

θ 39 • .
Figure 13b displays angular velocity profiles for various tilt angles. These profiles have been averaged from N = 100 configurations. As noticed by Azéma et al, these profiles tend to zero at the rough bottom interface and reach a maximum in absolute value a few layers away from this interface. Furthermore, each of these profiles agrees reasonably well with its corresponding -γ/2 profile [START_REF] Da Cruz | Rheophysics of dense granular flows: Discrete simulation of plane shear flows[END_REF][START_REF] Azéma | Discrete simulation of dense flows of polyhedral grains down a rough inclined plane[END_REF], in particular in the range of tilt angles [36 Figure 14b displays the friction coefficient profiles µ(z) calculated for all tested tilt angles from averaged σ yz (z) to σ zz (z) ratios (N = 100 configurations). A polynomial fit of these profiles is also represented using symbols to help distinguish between the curves. In the stationary flow regime, Savage and

Hutter [START_REF] Savage | The motion of a finite mass of granular material down a rough incline[END_REF] recall that µ must be uniform along z and equal to the tangent of the tilt angle. At least in the [36

Discussion

The ability of the new discrete element approach to simulate the quasistatic as well as the dynamic behaviour of dense polyhedra assemblies has been assessed.

In the quasistatic regime, this approach allowed to achieve the random close packing state described by Camenen and Descantes [START_REF] Camenen | Effect of confinement on dense packings of rigid frictionless spheres and polyhedra[END_REF], that is with a solid fraction ν 0.67 free of crystal nucleus, a mean coordination number in the bulk slightly in excess of 8, a hydrostatic stress distribution and even less steric hindrance violation than Camenen et al. In fact, apparent discrepancy was found solely on the proportions by contact types, with only 66% of single-point contacts reported by Camenen et al against 76% observed in the present paper.

In the dynamic regime, a stationary flow was identified for a range of tilt angle values consistent with the one reported by Azéma et al [START_REF] Azéma | Discrete simulation of dense flows of polyhedral grains down a rough inclined plane[END_REF], inside which the solid fraction profiles were found reasonably uniform in the bulk (viz. for z > 5d eq ), the translational velocity profiles clearly Bagnoldian, the angular velocity profiles similar to the -γ/2 profiles, and the stress profiles consistent with the mascroscopic expressions derived from the momentum balance equations describing the behaviour of the equivalent continuum. In addition, granular flows simulated by the new discrete element approach were found to be free of excessive particle overlap. A few discrepancies were also observed, in particular the coordination number values as well as the mean solid fraction values were found smaller than those reported by Azéma et al respectively by a factor of 1.5 to 2 and by 10%.

These discrepancies are briefly discussed in the next subsections.

Contact type discrepancies in the quasistatic regime

In order to investigate the contact type discrepancies from those reported by Camenen and Descantes [START_REF] Camenen | Effect of confinement on dense packings of rigid frictionless spheres and polyhedra[END_REF], it should be observed that both their results and those reported by Azéma et al [START_REF] Azéma | Discrete simulation of dense flows of polyhedral grains down a rough inclined plane[END_REF] were obtained from numerical simulations using the contact dynamics method [START_REF] Radjaï | Contact dynamics as a nonsmooth discrete element method[END_REF][START_REF] Radjaï | Discrete-element Modeling of Granular Materials[END_REF]. As already pointed out in section 1, this method disregards small sub-particle scales and focusses on the scale of particles rearrangements to calculate the unknowns, velocities and impulses (the integral of a force over one timestep), under the constraint of non-regularized contact laws reformulated to account for shocks and velocity jumps. Unlike this method, MD-like methods resolve contact situations at sub-particle scale, thus requiring a regularized mechanical model of the contact zone to explicitly calculate contact forces, here a regularized spring-dashpot model (see equation system 3).

With MD-like methods, attention shall be paid to avoid unphysical oscillations at interparticle contact, which may be done upon taking the dashpot viscosity critical. In the case of spherical particles, with a unique contact point between two spheres, the critical viscosity may easily be calculated (see Table 1). In the case of polyhedra for which edge -f ace and f ace -f ace contacts occur, adopting the same expression as for spheres for the critical viscosity may not totally eliminate unphysical oscillations in a f ace -f ace contact. A direct consequence of these oscillations is contact type variability as shown on Figure 15, where the contact type alternates between f ace -f ace and edge -f ace. This variability may explain the differences in terms of proportions by contact types between

Camenen et al and the present paper, as well as the fact that z iso was found slightly short of 12 (see Figure 9b). istics (e.g. energy stability for similar tilt angle values, see section 3.2.2). To the author's best knowledge, comparisons between contact dynamics and MDlike simulation results, which could shed some light on these discrepancies, have been performed in the static [START_REF] Richefeu | Force distribution in cohesive and non cohesive granular media[END_REF][START_REF] Camenen | Effect of confinement on dense packings of rigid frictionless spheres and polyhedra[END_REF] and flow regimes [START_REF] Azéma | Discrete simulation of dense flows of polyhedral grains down a rough inclined plane[END_REF], but unfortunately exclusively for sphere assemblies, which probably reflects the lack of MD-like codes accounting for true polyhedra. Although a detailed study of these discrepancies falls beyond the scope of the present paper, at least assumptions on their origin may be suggested.

Wachs et al [START_REF] Wachs | Grains3d, a flexible dem approach for particles of arbitrary convex shape part i: Numerical model and validations[END_REF] have used their in-house MD-like code to study the flow of cubes and tetrahedra assemblies in a rotating drum. They identified a rolling regime for drum rotational velocities in the range [20; 65] seems plausible for pinacoids. The corresponding solid fraction decrease is similar to that reported by da Cruz et al, but twice their value when expressed as a percentage of the maximum solid fraction (achieved for I = 0). A tentative explanation of this apparent discrepancy could be higher restitution coefficient values in the present study compared to Azéma et al, that would increase particle agitation and thus slightly decrease the solid fraction. As already underlined in section 4.1, the critical viscosity expressions adopted in the present study are valid for spheres, but could yield slightly subcritical damping when used for polyhedra. The level of particle agitation may be simply assessed upon calculating the particles translational velocity variance profile, also called granular temperature profile, defined as follows [START_REF] Richard | Influence of lateral confinement on granular flows: comparison between sheardriven and gravitydriven flows[END_REF]:

   T = T xx + T yy + T zz T jj = v j i , v j i δz,t -v j i 2 δz,t (13)
where ., . δz,t stands for the time-averaged mean value over particles i located in horizontal layer of thickness δ z and j = x, y, or z component. Figure 16d shows eight granular temperature profiles corresponding to the tested tilt angles, each averaged from N = 100 configurations. Note that these profiles differ significantly from those reported for gravity-driven flows of spheres by Richard 780 et al [START_REF] Richard | Influence of lateral confinement on granular flows: comparison between sheardriven and gravitydriven flows[END_REF]: here, the granular temperature first increases from bottom wall to roughly five particle diameters where it reaches its maximum. Observe that this maximum temperature strongly increases with the tilt angle value. Then, the granular temperature remains reasonably uniform up to the free surface in the stationary flow regime (viz. for tilt angle values in the [36 • ; 39 • ] range), whereas it decreases linearly up to the free surface in the collisional regime. Further investigation falling beyond the scope of the present study should hence focus on strategies to better control particles restitution coefficient when simulating true polyhedra flow using MD-like strategies, all the more that various contact types are possible.

Conclusion and perspectives

A new molecular dynamics-like approach aimed at simulating the mechanical behaviour of true polyhedra assemblies has been presented. This approach is innovative insofar as it combines an efficient two-steps gross contact detection algorithm suggested by Ogarko and Luding [START_REF] Ogarko | A fast multilevel algorithm for contact detection of arbitrarily polydisperse objects[END_REF] with the following three features. First, thanks to the use of an improved version of the GJK algorithm called GJK -T D [START_REF] Descantes | Classical contact detection algorithms for 3d dem simulations: Drawbacks and solutions[END_REF], this approach involves no edge or corner rounding and accounts for multipoint f ace -f ace or edge -f ace contacts between any set of two particles. Second, quaternion parameterization of the particle rotation matrices prevents the occurrence of conditioning issues commonly reported when using Euler angle parameterization. Third, torques equations are simply and efficiently solved using an improved leap-frog Verlet non-iterative method suggested by Omelyan [START_REF] Omelyan | A new leapfrog integrator of rotational motion. the revised angular-momentum approach[END_REF], in which no periodic renormalization of the quaternions is necessary.

The potential of this new discrete element approach has then been assessed by comparing the results of numerical simulations carried out both in the quasistatic and flowing regimes with those reported in the literature.

In the quasistatic regime, this approach allowed to achieve the random close packing state of polyhedra assemblies described by Camenen and Descantes [START_REF] Camenen | Effect of confinement on dense packings of rigid frictionless spheres and polyhedra[END_REF],

that is with a solid fraction ν 0.67 free of crystal nucleus, a mean coordination number in the bulk slightly in excess of 8 and a hydrostatic stress distribution.

In the dynamic regime, a polyhedra stationary flow down an incline was identified for a range of tilt angle values [36 • ; 39 • ] consistent with the one reported by Azéma et al [START_REF] Azéma | Discrete simulation of dense flows of polyhedral grains down a rough inclined plane[END_REF], inside which the solid fraction profiles were found reasonably uniform in the bulk (viz. for z > 5d eq ), the translational velocity profiles clearly Bagnoldian, the angular velocity profiles similar to the -γ/2 profiles, and the stress profiles consistent with the mascroscopic expressions derived from the momentum balance equations describing the behaviour of the equivalent continuum. Besides, upon averaging the coordination number, solid fraction and inertial number values in the bulk, a linear decrease of the bulk solid fraction and a power law decrease of the bulk coordination number with increasing inertial number were observed, thus generalizing to polyhedra these observations initially made with disks in plane shear flow by da Cruz et al [START_REF] Da Cruz | Rheophysics of dense granular flows: Discrete simulation of plane shear flows[END_REF].

In both regimes, particle overlap was found to be of order 10 -4 d eq which is two orders of magnitude better than in previous works.

A few discrepancies were also observed, regarding slightly differing proportions by contact types in the quasistatic regime, as well as slightly lower coordination number and mean solid fraction values in the flowing regime. These discrepancies could be caused by difficulties in controling particles restitution coefficient when simulating true polyhedra behaviour using molecular dynamics-like strategies, all the more that various contact situations are possible. This will be investigated further in a coming paper.

Several perspectives may also be drawn. In the short term, the new approach will account for spherical particles to simulate the behaviour of combined sphere and polyhedra assemblies. OpenMP parallelization to minimize calculation expenditures is also underway, as well as lattice-Boltzmann coupling in order to account for a fluid phase.

  115

3. 1 . 1 .

 11 Compaction simulation procedure and parameters 2500 spheres of diameter d sph ± 0.05d sph were successively geometrically deposited following guidelines from Camenen and Descantes [17] into a parallelepiped container of length L x = 10d sph and width L y = 10d sph , at random locations corresponding to local potential energy minima. Then, a randomly

Figure 6a displays the

  Figure 6a displays the mean and maximum overlap between two particles as a function of time, while Fig. 6b depicts corresponding kinetic energy fluctuations. Both figures evidence four zones corresponding to different combinations of contact mechanical parameters (their values appear in the figure caption). The overlap is expressed as a fraction of the mean diameter d sph of pinacoidcircumscribed spheres. Overlap fluctations are maximum during the first 2.10 6

  Fig. 6b depicts corresponding kinetic energy fluctuations and shows that the flow is about to stop for θ = 30 • , whereas it accelerates for θ ≥ 42 • . The [36 • ; 39 • ] range seems to frame the tilt angle value interval for which the energy curves are constant, suggesting the achievement of a stationnary flow. This range looks consistent with the [36 • ; 38 • [ range of tilt angle values for which Azéma et al [19] reported stationnary flows.

Figure

  Figure 12b displays eight coordination number profiles corresponding to the tested tilt angles, each averaged over N = 100 configurations. In these profiles, each contact is counted once regardless of its type. As reported by Azéma et al, coordination number decreases when tilt angle increases, and the former increases close to the rough bottom. In the range of tilt angle values [36 • ; 39 • [, deemed to frame the stationary regime, the coordination number decreases from 1.8 to 1.4 and looks more or less homogenous in the bulk, whereas it looks more homogenous in the tilt angle range [42 • ; 50 • ] where fluctations lie between 1

Figure 14a depicts the

  Figure14adepicts the normal and shear stress profiles σ zz (z) and σ yz (z) corresponding to the tested tilt angle values. These profiles were calculated using the stress tensor expression suggested by Moreau (equation 12) and averaged from N = 100 configurations. The trend of these profiles is clearly linear, and the normal stress profiles are hydrostatic as evidenced e.g. by the inset for a tilt angle θ = 39 • (the weight by unit surface superimposes very well with the trend of σ zz ). These results are consistent with those reported byAzéma et

4. 2 .

 2 Solid fraction and coordination number discrepancies in the stationary flow regime Whereas in the quasistatic regime both the density and coordination number (regardless of contact type) calculated from the new discrete element approach compare remarkably well with those determined by Camenen et al from contact dynamics simulations (see section 3.1.2), this does not seem to be the case in the stationary flow regime despite general agreement on other flow character-

  rounds per minute, and reported in this range a decreasing trend for the coordination number from Z = 3 to Z = 2 with increasing rotational velocity. These coordination number values seem somewhat more consistent with those reported by Azéma et al for pinacoids in stationary flow. It should however be observed that matching the rolling regime to the stationary flow on an incline is disputable, given that the angle of repose measured in the former are significantly higher than the tilt angle values reported by Azéma et al in the latter, and they fluctuate with the solid volume inside the drum. Besides, with the in-house code used, the edges of polyhedra have to be rounded for contact detection reasons, which may impact both the solid fraction and coordination number in the flowing regime compared to true polyhedra. Da Cruz et al[START_REF] Da Cruz | Rheophysics of dense granular flows: Discrete simulation of plane shear flows[END_REF] have reported that the coordination number of frictional disks in 2D plane shear flow was a strongly decreasing function of the inertial number, more precisely according to a power law of the form Z = Z max -cI γ with Z max the maximum coordination number for a given interparticle friction µ, c and γ being two fit parameters. They partly explained this upon showing that the solid fraction was a linearly decreasing function of the inertial number in the stationary flow range. They also observed for an interparticle friction µ = 0.4 that the power law c coefficient was affected by particles restitution coefficient e (ratio of particle velocities after to before shock), whereas γ 0.5 remained constant. Besides, accross their stationary flow regime defined by the [10 -3 ; 10 -1 ] range of inertial number values, they observed a coordination number decrease from Z = 4 to Z = 1.5 or Z = 0.5 for low (e = 0.1) to high (e = 0.9) restitution coefficient values respectively, as well as a 5% drop of the solid fraction (from 0.82 to 0.78). It is hence interesting to check whether such conclusions apply to assemblies of polyhedra in their stationary flow regime, since coordination number discrepancies could then be justified at least by higher restitution coefficient values in the present study compared to Azéma et al.

Figure 16a displays eight 1 2

 1 Figure16adisplays eight inertial number profiles corresponding to the tested tilt angles, each averaged from N = 125 configurations. A polynomial fit of these profiles is also represented using symbols to help distinguish between the curves. Although these profiles look slightly more scattered than those reported by Azéma et al, their shape is similar in the [36 • ; 39 • ] range of tilt angle values, with I reasonably uniform in the bulk (away from rough bottom and free surface). Note that the polynomial fit tends to be convex for θ = 30 • , whereas it tends to be more and more S-shaped for tilt angles θ = 42 • and higher, thus suggesting a quasistatic to stationary regime transition between θ = 30 • and θ = 36 • , and a stationary to collisional regime transition between θ = 39 • and θ = 42 • . Last, note that the stationary flow of pinacoids seems to take place in the [0.5; 0.7[ range of inertial number values, which is significantly higher than the range reported by da Cruz et al for disks, though 0.2 to 0.3 points lower than the [0.7; 1] range observed by Azéma et al.

Figure 1 :

 1 Figure 1: Contact types between polyhedra (a) Vertex/face (b) edge/edge, (c) edge/face, (d) face/face.

Figure 2 :

 2 Figure 2: GJK algorithm principle, (a) Two contacting particles A and B, (b) Minkowksi difference of the particles, (c) final simplexe (tetrahedra) containing origin O with closest point P = P A -P B and contact normal n.

Figure 3 :

 3 Figure 3: (a) Intersection between triangulated particle faces, (b) In-common plane zoom of intersections between triangles.

Figure 4 :

 4 Figure 4: Pinacoid, a four-parameter polyhedron (length L, width G, thickness E and angle α).

Figure 5 :

 5 Figure 5: Pinacoid packing at (a) initial and (b) final timestep.

Figure 9 :Figure 10 :

 910 Figure 9: (a) Pair correlation function of the 2500 pinacoids packing, and (b) Coordination number profiles broken down into vertex/face & edge/edge (zs), edge/face (z d ) and face/face (zt) contacts, with ztot their sum and z iso = zs + 2z d + 3zt.

Figure 11 :

 11 Figure 11: (a) Particles mean overlap δ (normalized by the mean diameter d sph of pinacoidcircumscribed spheres) as a function of time, and (b) Ratio of the total kinetic energy to the total potential energy of the granular assembly.

Figure 15 :

 15 Figure 15: Time oscillations between f ace -f ace (1) and edge -f ace (2) at contact between a polyhedron and the bottom wall in the 2500 pinacoids packing simulation; inset is a zoom over 10 timesteps.

Table 2 :

 2 Parameters used to simulate the flow of 20000 pinacoids down a rough inclined plane. cos θ, -sin θ) T θ {30 • ; 36• ; 37 • ; 39 • ; 42 • ; 45 • ; 47.5 • ; 50 • } T inertia d eq /g T stif f ness T inertia

Figure 16 :

 16 Figure 16: (a) Inertial number profiles and their polynomial fits, (b) Coordination number as a function of inertial number (symbols) and its power law fit, Ztot(I) 4.5 -3.62I 0.5 with R 2 = 0.96 (dotted line), (c) Solid fraction as a function of inertial number (symbols) and its linear fit, ν(I) 0.50 -0.07I with R 2 = 0.99 (dotted line), (d) Granular temperature profiles for various tilt angle values.

Table 2

 2 gathers the simulated system parameters. This table incorporates a few differences from Table1. Particle polydispersity has been doubled to 0.1d sph as in Azéma et al[START_REF] Azéma | Discrete simulation of dense flows of polyhedral grains down a rough inclined plane[END_REF] to avoid crystallisation in the flowing layer. The gravity has been gradually inclined to eight discrete values spanning a significantly larger range than the [35.7 • ; 38 • ] steady flow range reported by Azéma et al.

Acknowledgements

The author thanks his former master's students F. Tricoire and A. Haidar for contributions which have proved useful for developping the present molecular dynamics-like approach. He also acknowledges interesting discussions with P.

Richard at various stages of this development.

eq / α P ot = 0.7 in zone 4. Table 1: Parameters used to simulate the compaction of 2500 pinacoids under their own weight.

Parameters Value


Oleg Lookin

Anastasia Khokhlova

Tatiana Myachina

Xenia Butova

Olivier Cazorla

Pieter De Tombe

Contractile State Dependent Sarcomere Length Variability in Isolated Guinea-Pig Cardiomyocytes

come