Paul Sánchez

Mathieu Renouf

Emilien Azéma

Rémy Mozul

Frédéric Dubois

A Contact Dynamics Code Implementation for the Simulation of Asteroid Evolution and Regolith in the Asteroid Environment

Keywords: celestial mechanics -granular -methods, numerical -minor planets, asteroids, general

Over the last decades, simulations by discrete elements methods (DEM) have proven to be a reliable analysis tool in various domains of science and engineering. By providing access to the local physical mechanisms, DEM allows the exploration of microscopic based phenomena related to particles properties and interactions in various conditions and to revisit constitutive equations consequently. The growing computer power and memory now allow us to handle large collections of grains of various shapes and sizes by DEM simulations and in particular, it offers new perspectives in the exploration of the behavior of asteroids seen as self-gravitating and cohesive granular aggregates. In this paper we describe the Contact Dynamics (CD) method, a class of DEM based on non-smooth mechanics, and its implementation in the open-source software LMGC90. In contrast to more classical approach, Hard-and Soft-Sphere DEM, the CD method is based on an implicit time integration of the equations of motion and on a non-regularized formulation of mutual exclusion between particles. This numerical strategy is particularly relevant to the study of dense granular assemblies (even of large size) because it does not introduce numerical artefacts due to contact stiffness. So that it can be used for Small Body research, we implement a parallelised kd-tree and monitor the performance of the code as it simulates a number of granular systems. We provide examples of the simulation of the accretion of self-gravitating aggregates as well as their rotational disruption. We use the routines at our disposal in the code to monitor their behaviour through the entire evolution and find agreement with previous research.

Introduction

During the last two decades, research about asteroids in general and small asteroids in particular has increased dramatically. The first near-Earth asteroid (NEA), (433) Eros was discovered 1898, but it would be a long time before its first pictures were taken by the NEAR Shoemaker mission almost 20 years ago. In the mean time, asteroids passed from being cursed as "the vermin of the sky" [START_REF] Friedman | Vermin of the sky[END_REF] by astronomers trying to observe distant stars and galaxies, to be the target of a number of space missions. The most recent finished one is the JAXA Hayabusa mission to asteroid (25143) Itokawa [START_REF] Fujiwara | The Rubble-Pile Asteroid Itokawa as Observed by Hayabusa[END_REF] and the two that are taking place at the time of this writing are the JAXA Hayabusa2 mission to asteroid (162173) Ryugu [START_REF] S.-I. Watanabe | Hayabusa2 mission overview[END_REF] and the NASA OSIRIS-REx mission to asteroid (101955) Bennu [START_REF] Lauretta | The osirisrex mission-sample acquisitions strategy and evidence for the nature of regolith on asteroid (101955)[END_REF]. One common objective in all these missions is the return of a sample of the particles on their surfaces for analysis. For the last two, the use of numerical methods to simulate the interaction of the spacecrafts with the asteroid surfaces has been a very important part of the research efforts [START_REF] Tardivel | Contact motion on surface of asteroid[END_REF][START_REF] Zacny | Chapter 8 -geotechnical properties of asteroids affecting surface operations, mining, and in situ resource utilization activities[END_REF][START_REF] Van Wal | Prearrival deployment analysis of rovers on hayabusa2 asteroid explorer[END_REF].

The pictures obtained of asteroid Itokawa revealed a small body covered in dust, pebbles, rocks and boulders going from micrometers to tens of meters in size [START_REF] Nakamura | Itokawa Dust Particles: A Direct Link Between S-Type Asteroids and Ordinary Chondrites[END_REF][START_REF] Tsuchiyama | Three-dimensional structure of hayabusa samples: Origin and evolution of itokawa regolith[END_REF]. The sample that was returned to Earth in 2010 revealed its chemical composition as well as the varied shapes of the few dust grains that were trapped in the sample canister after the failure of the sampling mechanism [START_REF] Yano | Touchdown of the Hayabusa Spacecraft at the Muses Sea on Itokawa[END_REF][START_REF] Nakamura | Itokawa Dust Particles: A Direct Link Between S-Type Asteroids and Ordinary Chondrites[END_REF]. These and previous findings confirmed the idea that asteroids had a granular structure and that their bodies were held together by gravitational attraction. This being so, the use of the numerical tools and theoretical approaches to study granular matter became a necessity, especially codes that implemented different Discrete Element Methods (DEM).

By definition, DEM aim to model the behaviour of a collection of distinct interacting particles. As any numerical method, it constitutes a significant support to experimentation in the sense that they give us access to information that is difficult to obtain experimentally (i.e., packing disorder, contacts orientation, force distributions). DEM make it possible to multiply numerical experiments in order to study the influence of particles characteristics (particle size, morphology), interactions (friction, cohesion) and loading, on the local and collective behaviour of the system. The common denominator of DEM is to consider the degrees of freedom associated with the elements (grains), considered as rigid objects, and to integrate the equations of motion for these degrees of freedom. DEM can be classified into two main families of approaches: Smooth and Non Smooth.

In Smooth approaches [START_REF] Cundall | A computer model for simulating progressive large scale movements of blocky rock systems[END_REF][START_REF] Allen | Computer simulation of liquids[END_REF], interaction forces between particles can be written as a function linking contact forces to contact kinematics. Typically, the normal reaction is taken as proportional to particle penetration (Hertzian contact model [START_REF] Hertz | On the contact of elastic solids[END_REF]). The particle motion is smooth (at least as twice differentiable) and, therefore, the equations of dynamics are written as ordinary differential equations that can be integrated by conventional methods such as the Gear method for Molecular Dynamics (MD) [START_REF] Allen | Computer simulation of liquids[END_REF], the Newmark [START_REF] Cundall | A computer model for simulating progressive large scale movements of blocky rock systems[END_REF], or Verlet scheme [START_REF] Iordanoff | Solid third body analysis using a discrete approach: Influence of adhesion and particle size on the macroscopic behavior of the contact[END_REF].

In Non Smooth approaches [START_REF] Hahn | Realistic animation of rigid bodies[END_REF][START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF][START_REF] Jean | Unilaterality and dry friction in the dynamics of rigid bodies collection[END_REF][START_REF] Baraff | Issues in computing contact forces for non penetrating rigid bodies[END_REF][START_REF] Dubois | The Contact Dynamics method: A nonsmooth story[END_REF], there is no explicit relation between contact forces and contact kinematics. Contact conditions are described by complementarity relations between the contact forces and displace- 1 Corresponding Author: diego.sanchez-lana@colorado.edu 2 ments or velocities (as the well know Signorini condition [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF]). No regularization is required. The equations of motion are rewritten as non-differentiable relations involving velocity jumps and impulse [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF]. Equations of motion can be driven using time stepping or Event Driven (ED) [START_REF] Radjai | Discrete-element modeling of granular materials[END_REF] and discretized by the θ-method [START_REF] Jean | Unilaterality and dry friction in the dynamics of rigid bodies collection[END_REF] or Leap Frog (LF) method [START_REF] Richardson | Direct large-scale n-body simulations of planetesimal dynamics[END_REF]. The ED method, associated with irregular time discretization, is well suited for collections of rigid bodies in which the time to cover the mean free path is much greater than the contact time of a collision between two bodies. The method then assumes that the collision time is effectively zero and so, by construction, only binary, and not multi-body, collisions can be simulated.

At the end fo the 90s, a first ED (Hard-Sphere DEM) model is proposed to the study of self-gravitating particles of meters in size [START_REF] Richardson | Direct large-scale n-body simulations of planetesimal dynamics[END_REF][START_REF] Richardson | Numerical experiments with rubble piles: equilibrium shapes and spins[END_REF]. The authors developed a code named PKDgrav [START_REF] Richardson | Direct large-scale n-body simulations of planetesimal dynamics[END_REF], originally used to study star clusters, which has been used also to study planetary rings, planetesimal formation, binary asteroid formation, rotational evolution of asteroids and asteroid collisions to name a few topics [START_REF] Walsh | Binary near-earth asteroid formation: Rubble pile model of tidal disruptions[END_REF][START_REF] Walsh | Rotational breakup as the origin of small binary asteroids[END_REF][START_REF] Cotto-Figueroa | Coupled spin and shape evolution of small rubble-pile asteroids: Self-limitation of the yorp effect[END_REF][START_REF] Walsh | Spin-up of rubble-pile asteroids: Disruption, satellite formation, and equilibrium shapes[END_REF][START_REF] Ballouz | Numerical simulations of collisional disruption of rotating gravitational aggregates: Dependence on material properties[END_REF][START_REF] Bagatin | Internal structure of asteroid gravitational aggregates[END_REF].

More recently (2008), an MD (Soft-Sphere DEM) has been introduced to model self-gravitating granular systems [START_REF] Sánchez | Granular Mechanics in Asteroid Regolith: Simulating and Scaling the Brazil Nut Effect[END_REF][START_REF] Sánchez | Simulating asteroid rubble piles with a self-gravitating soft-sphere distinct element method model[END_REF]. Based on a Smooth DEM formalism, the code developed by the authors has been used to study the rotational evolution of asteroids, internal structure and strength of small asteroids, binary asteroid and asteroid pair formation, and penetrometry in the asteroid environment [START_REF] Sánchez | Dem simulation of rotation-induced reshaping and disruption of rubble-pile asteroids[END_REF][START_REF] Sánchez | The strength of regolith and rubble pile asteroids[END_REF][START_REF] Sánchez | Asteroid evolution: Role of geotechnical properties[END_REF][START_REF] Hirabayashi | Internal structure of asteroids having surface shedding due to rotational instability[END_REF][START_REF] Sánchez | Disruption patterns of rotating self-gravitating aggregates: A survey on angle of friction and tensile strength[END_REF][START_REF] Sánchez | Rotational evolution of self-gravitating aggregates with cores of variable strength[END_REF][START_REF] Tardivel | Equatorial cavities on asteroids, an evidence of fission events[END_REF]. This last point was largely applied to the simulation of the Touch-and-Go Sample Acquisition Mechanism (TAGSAM) for the OSIRIS-REx mission [START_REF] Zacny | Chapter 8 -geotechnical properties of asteroids affecting surface operations, mining, and in situ resource utilization activities[END_REF]. Note that a SSDEM version has been available in the PKDgrav code mentioned previously since 2012 [START_REF] Schwartz | An implementation of the soft-sphere discrete element method in a high-performance parallel gravity tree-code[END_REF].

Later on, in 2014, a first use of a Non Smooth approach was proposed to study regolith processes [START_REF] Hartzell | Contact dynamics methods to study regolith processes[END_REF]. More recently, in 2018, the Contact Dynamics (CD) method, was applied to analyze the strength properties of self-gravitating aggregates of spheres [START_REF] Azéma | Scaling behavior of cohesive self-gravitating aggregates[END_REF]. A year later, Ferrrari, et al [START_REF] Ferrari | A parallel-GPU code for asteroid aggregation problems with angular particles[END_REF][START_REF] Ferrari | The role of fragment shapes in the simulations of asteroids as gravitational aggregates[END_REF] used the same method to study asteroid aggregation problems with angular particles.

The implicit formulation of the method and the introduction of nonsmooth laws in iterative or direct algorithms makes the CD method less accessible for computer implementation than other DEM methods based on explicit formulation. Thus, the aim of this paper is to present the spirit of the CD method, some details of its implementation in the LMGC90 open-source platform ( [START_REF] Dubois | The non smooth contact dynamic method: recent lmgc90 software developments and application[END_REF] 2 ) together with a direct application of self-gravity. In this paper, to underline the numerical efficiency of the CD method for modelling granular asteroids, the CD method is applied to model the accretion of spherical and polyhedral particles. Then, some specificities for the modelling of granular asteroid are introduced, including cohesion, particle shape and self-gravity. The main results of accretion simulations are discussed and some perspectives are given regarding forthcoming research avenues.

The Contact Dynamics Method

The Non-Smooth Philosophy

The objective here is to represent an asteroid as a collection of rigid solids subject to conditions of mutual noninterpenetrability and friction in case of contact. This situation can be described as non-smooth for at least three reasons: The geometrical conditions of non-interpenetration of the different objects (non-smoothness in space), the contact forces governed by non-regularized laws (non-smoothness in law) and velocity jumps due to collisions between bodies (non-smoothness in time).

One way to get rid of this non-smoothness is the use of regularization technics: This is what the smooth-DEM methods do when the non-interpenetrability of bodies is replaced by sufficiently stiff repulsion laws that apply when two of them approach each other. We are thus reduced to differential equations that can use classical numerical techniques. But, in each case, a compromise must be accepted between the requirement of precision and the stiffness of the approximate equations. This stiffness imposes very small discretization steps and often artificial inertias or viscosities are introduced to ensure numerical stability [START_REF] Iordanoff | Solid third body analysis using a discrete approach: Influence of adhesion and particle size on the macroscopic behavior of the contact[END_REF][START_REF] Nackenhorst | Formulation of a three-dimensional distinct element model -part I: a scheme to detect and represent contacts in a system composed of many polyhedral blocks[END_REF].

To overcome this kind of process, Moreau relies on the framework of convex analysis [START_REF] Rockafellar | Convex Analysis[END_REF] to develop what he will call the Contact Dynamics (CD) approach [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF]. It is based on two main ingredients: non-smoothed contact and non smooth dynamics. It implies:

1. An implicit time integration scheme to describe the dynamic evolution of the particles and to manage in the same time contacts and shocks; 2. A contact resolution algorithm ensuring the respect of constraints related to the choice of non-smooth interaction laws.

To use CD within a DEM philosophy, both previous ingredients are combined in a contact detection algorithm to deal with potential interactions between particles.

Dynamical Equations

In the body frame

The motion of a set of N p rigid particles is governed by the Newton-Euler's equations:

         M q = F ext + R, J ω + ω × Jω = M ext + M, (1)
involving the second time derivative of the configuration parameter q, the spin ω expressed in a frame attached to the body, the external force and moment resultants (for example gravitational forces), denoted F ext and M ext respectively, and the contact reaction forces (exerted by neighboring particles) called R. M and J represent respectively the mass and the inertia matrices.

Let us note that, in practice, the Euler equation is written in the main inertia frame of each body which allows to have a diagonal representation of J. In the specific case of spheres, ω becomes collinear to Jω thus allowing to remove the non-linear terms related to the cross product. On the other hand, for complex shapes, the Euler equation remains non-linear and depends implicitly on the orientation of the objects. We discuss in Sec 2.5.1 in more details the special treatment used to integrate the rotations when particles are non-spherical.

Thus, the velocity of a particle is the combination of a translation velocity of its center of gravity and an angular velocity vector expressed in the principal inertial frame attached to the particle. By excessive notation, the system [START_REF] Friedman | Vermin of the sky[END_REF] will be reduced to its first equation and the different variables used will contain the kinematics and forces related to translations and rotations:

M q = F ext + R, (2)
To deal with the non-smoothness tryptic mentioned previously, Equation ( 2) is reformulated in terms of equality of measure (in the mathematical sense) as follows:

Md q = F ext dt + dR, (3)
where d q is the differential measure of velocity, dt is the Lebesgue measure on the real space R, and dR is a measure of the impulse. In the following, for a sake of simplicity, R denotes the resultant of an impulse.

To compute the evolution of a system of particles, one should proceed to the time discretization. The multi-contact problem is solved over the interval ]t i , t i+1 ] of length h in the terms previously defined. Then successive approximations of both Equation (3) and first time derivative of the configuration parameter lead to the following system:

         qi+1 = q f ree i + M -1 R i+1 q i+1 = q i + hθ qi+1 + h(1 -θ) qi (4)
with

q f ree i = qi + M -1 h(θF ext i+1 + (1 -θ)F ext i ) (5)
where q f ree denotes the free velocity (velocity computed without contact forces). Quantities indexed by i (resp. i + 1) refer to time t i (resp. t i+1 ). The system is updated according to the second Equation of system (4). Note that, the update of the orientation of each particle when dealing with non-spherical particle will be discussed in details in Sec.2.5.1.

In order to be unconditionally stable, θ ranges from 0.5 to 1. When θ = 0, the contact detection is performed using the configuration parameter at time t i . On the contrary, when θ 1, the contact detection is performed using a prediction of the configuration parameter (at time t i + h(1 -θ)).

In the contact frame

Let us consider two touching particles, i and j, at the contact point α (cf. Fig. 1). The velocity of the point α attached to each particle is simply given by u α,i = qi + ω i × O i α, where O i α is the vector joining the center of particle i to contact point α (same notation fo particle j). The relative velocity at the contact point α is given by v α = u α,iu α, j . At the contact point α, it is possible to define a local frame (n α , t α , s α ) attached to the contact α where n α is oriented from particle j to particle i and where t α and s α define the orthogonal plane to n α ; The triplet (n α , t α , s α ) is an orthonormal basis (cf. Fig. 1).

Using such a frame, it is possible to build a linear operator H (and its transpose H T ) to connect the set of local variables, the relative contact velocity v and the contact impulse r, to the set of global variables, R and q, by the couple of relations:

         R = Hr v = H T q. (6)
A more detailed description of the construction of these matrices can be found in [START_REF] Radjai | Discrete-element modeling of granular materials[END_REF]. These two equations can be injected in the first equation of system (4), leading to the formulation of the equation of dynamics in the contact frames as:

Wr i+1 -v i+1 = -v f ree . (7)
In Equation [START_REF] Van Wal | Prearrival deployment analysis of rovers on hayabusa2 asteroid explorer[END_REF], W (= H T M -1 H) is called the Delassus operator while the right-hand-side represents the free relative velocity only accounting for the external forces given by Equation [START_REF] Tardivel | Contact motion on surface of asteroid[END_REF].

Equation [START_REF] Van Wal | Prearrival deployment analysis of rovers on hayabusa2 asteroid explorer[END_REF] presents two unknowns: the local contact impulse r and the contact relative velocity v (the time index is omitted to make pleasant reading). To find a solution to the contact problem, a relation between these two unknowns should be added through the definition of a contact law, detailed in the next section. Assuming such a law defined, Equation ( 7) can be solved contact by contact using the so-called Block Non Smooth Gauss-Seidel algorithm (NSGS) a splitting method exposed in [START_REF] Cottle | The linear complementarity problem[END_REF]. Using the global splitting scheme on Equation [START_REF] Van Wal | Prearrival deployment analysis of rovers on hayabusa2 asteroid explorer[END_REF], one obtains:

W αα r k+1 α -v k+1 α = -v α, f ree - n c β<α W αβ r k+1 β - n c β>α W αβ r k β ( 8
)
where the index k refers to the splitting method iterations and n c the number of contacts. The time index is omitted to make pleasant reading. This solver has proved to be very robust and efficient on a large collection of heterogeneous problems [START_REF] Saussine | Modelling ballast behaviour under dynamic loading. part 1: A 2d polygonal discrete element method approach[END_REF][START_REF] Azéma | Internal Structure of Inertial Granular Flows[END_REF][START_REF] Zeng | Numerical investigation on the electrical transmission ability of a shearing powder layer[END_REF] and benefits of a parallel version [START_REF] Renouf | A parallel version of the Non Smooth Contact Dynamics algorithm applied to the simulation of granular media[END_REF] to ensure reduced simulation time.

It should be noted that the problem of frictional contact between rigid bodies can have a multiplicity of solutions (inherent to the rigid formulation). Each solution will depend on the reading direction of the contact list. Nevertheless, whatever the solution obtained, it respects all the constraints imposed and all the solutions always present the same physical properties [START_REF] Renouf | A parallel version of the Non Smooth Contact Dynamics algorithm applied to the simulation of granular media[END_REF][START_REF] Alart | On inconsistency in frictional granular systems[END_REF].

Contact laws

When considering a system composed of rigid bodies, the physical behaviour of the system depends entirely on the laws of interaction between particles. Therefore, depending on the physics one wishes to represent, the contact law may be more or less complex, from a simple unilateral contact law [START_REF] Jean | The non smooth contact dynamics method[END_REF] to cohesive zone models [START_REF] Raous | A consistent model coupling adhesion, friction, and unilateral contact[END_REF]. The case of a unilateral frictional contact law with Coulomb friction is presented therefore.

The contact unilaterality is described by the well-known Signorini Condition [START_REF] Signorini | Questioni di elasticità non linearizzata e semilinearizzata[END_REF] relating the normal component of the contact relative velocity v n and the normal component of the contact impulse r n . As long as v n remains positive, no force is acting and r n is equal to 0. When v n = 0, r n becomes positive. Such conditions are summarized by the following set of inequalities:

r n ≥ 0 v n ≥ 0 r n .v n = 0. (9)
Note that the Signorini condition corresponds to a pure inelastic shock. According to the set of inequalities [START_REF] Tsuchiyama | Three-dimensional structure of hayabusa samples: Origin and evolution of itokawa regolith[END_REF], when a collision occurs between two bodies, the relative velocity vanishes. In order to allow bounces in the simulation, one can introduce, for example, a restitution coefficient, linking the velocity before and after shock (v i+1 n = -e n v i n ). In the formalism presented above, the contact problem with restitution becomes:

r n ≥ 0 (1 + e n )v n ≥ 0 r n .(1 + e n )v n = 0. ( 10
)
where v = v i+1 n + e n v i n . When e n = 1 the shock is fully elastic, while when e n = 0 the shock is plastic and equivalent to the Equation ( 9).

The Amontons-Coulomb friction law is one of the simplest way to consider a resistance in the tangential direction.

When the norm of the sliding velocity v t is equal to 0, the norm of the tangential impulse r t belongs to the interval [-µr n , µr n ] where µ is called friction coefficient. When the norm of v t is no longer equal to 0, the direction of r t is opposed to the sliding motion and equal to the friction threshold. This friction law can be summarized by the following set of inequalities:

         v t = 0 ⇒ -µr n r t µr n v t 0 ⇒ r t = -µr n * v t /||v t || . ( 11
)

Resolution of the local frictional contact problem

Before trying to solve the frictional contact problem, it is necessary to rewrite Equation 8 differently. If we make the difference of two successive iterations appear to the left hand side, then the equation can be written as follows:

W αα (r k+1 α -r k α ) = v k+1 α -v α, f ree - n c β<α W αβ r k+1 β - n c β≥α W αβ r k β ( 12
)
The convergence of the Gauss-Seidel algorithm makes the right-hand side member tend towards zero. It is therefore possible to replace the Delassus operator W by any matrix with better properties, such as a diagonal matrix. In practice, the matrix W is replaced by its diagonal W [START_REF] Radjai | Discrete-element modeling of granular materials[END_REF]. Thus, the modified Equation ( 8) can be explicitly written as:

Wαα r α,k+1v α,k+1 = -bα [START_REF] Hertz | On the contact of elastic solids[END_REF] where bα = b α + ( Wαα -W αα )r α,k and b α is the right hand side of Equation ( 8). As a result, the three components of the local impulse are thus decoupled and can be explicitly expressed through a case study, illustrated by Fig. 2. If b n is positive then r vanishes (no contact). On the contrary, due to Equation ( 13), r * = -W -1 nn b n . This value allows us to determine the circle C t , intersection between the tangent plane P t and the Coulomb cone (cf. Fig. 2a). If the norm of the vector b t is less than µW tt r * , i.e. that the solution is within the friction cone, then r t = -W -1 tt b t (cf. Fig. 2b). Otherwise, the solution is outside the friction cone and should be projected on its boundary (cf. Fig. 2c). This discussion is summarized in the Table 1. This resolution method is independent of the particle geometry and, in the case of collections of rigid bodies, the approximation performed on the W matrix does not disturb the convergence of the Gauss-Seidel algorithm [START_REF] Radjai | Discrete-element modeling of granular materials[END_REF]. It can be noted that there are methods that do not approximate the matrix W [START_REF] Alart | Méthode de newton généralise en mécanique du contact[END_REF][START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems: Applications in Mechanics and Electronics[END_REF], but they are slower and can be put in default when using high coefficients of friction.

if bn > 0 then r = 0 (status = no contact) else tt bt (status = sliding contact)

r n = -W-1 nn bn if bt ≤ µ Wtt r n then r t = -W-1 tt bt (status = sticking contact) else β = µ Wtt
Table 1: Pseudo-code of the resolution of the 3D frictional contact problem.

Extensions of interest to model granular asteroids

The CD method could be directly used to model the granular behaviour of asteroids by considering the previous ingredients. However, to take full advantage of the potential of the method, different ingredients can be introduced, both in terms of geometry, interaction law and boundary conditions.

Indeed, granular asteroids are gravitational aggregates made of regolith, i.e. a competent rocks of various sizes (from microns to meters) and irregular shapes, with low, though non-negligible, internal cohesive strength. The shapes of the regolith can be angular or more rounded, usually elongated and generally convex following the detailed studies carried out on the Itokawa dust particles [START_REF] Fujiwara | The Rubble-Pile Asteroid Itokawa as Observed by Hayabusa[END_REF] brought by the Hayabusa mission.

This size dispersion induces interactions between particles that cannot be reduced to simple laws of frictional contact. More representative laws need to be introduced in order to take into account this complexity, without neglecting the presence of the gravitational forces. The next subsections will detailed these specific features in the framework of Contact Dynamics.

Modelling regolith as polyhedral particles

Grain angularity modeling is of major importance to model realistic granular asteroids in order to better predict their behavior. This importance is underlined by the pionneering work of Korycansky and Asphaug, who have studied the low-speed impact of rubble piles modeled by polyhedra [START_REF] Korycansky | Low-speed impacts between rubble piles modeled as collections of polyhedra[END_REF]. Recent reports show that particle shape strongly affects the strength and dilatancy properties of the granular media. A basic result is that angularity increases the shear strength [START_REF] Mitchell | Fundamentals of Soil Behavior[END_REF], but it is only very recently that systematic studies have been reported on the effect of particle shape angularity. In this case, a nonlinear dependence of the shear strength and packing fraction has been evidenced

numerically [START_REF] Azéma | Packings of irregular polyhedral particles: Strength, structure, and effects of angularity[END_REF].

The introduction of polyhedral particles requires the use of specific detection algorithms and a method to integrate the rotations.

First, concerning the contact detection. While the coarse detection (sorting method) can be the same as for spheres, the fine detection cannot be based only on a distance test between the particle centres. Coarse detection, also called proximity search, is a sorting method allowing to build a neighborhood list from the bounding spheres of the different objects constituting the study system. For polyhedral particles, fine detection consists of two parts. The first, called rough, allows the elimination of useless contacts by a technique of axis separator (shadow overlap) [START_REF] Saussine | Modelling ballast behaviour under dynamic loading. part 1: A 2d polygonal discrete element method approach[END_REF]. The second, called accurate, can use a method of intersection of triangular faces [START_REF] Saussine | Modelling ballast behaviour under dynamic loading. part 1: A 2d polygonal discrete element method approach[END_REF] or a method of plane separator [START_REF] Nackenhorst | Formulation of a three-dimensional distinct element model -part I: a scheme to detect and represent contacts in a system composed of many polyhedral blocks[END_REF]. Here, the choice is made on the second method which consists in determining the distance between each pair of polyhedra by computing the separating plane. This is an iterative method based on the perturbation of the orientation of the normal vector where the process is initialized by the vector joining the polyhedron centers. The main advantage of this method is its speed when optimized as proposed in [START_REF] Radjai | Discrete-element modeling of granular materials[END_REF]. Concerning the first point, and assuming that the rotation speeds of the grains are "low" (see test cases below) we use a semi-implicite integrator for the Euler equations as follows:

(a) (b) (c) (d)
ω i+1 = ω f ree + J -1 M i+1 , (14)
where J -1 is the inverse of the inertia matrix, M is the resultant of the contact impulse torque, and ω f ree = J -1 M ext i+1ω i × Jω i + ω i the free rotation speed. Concerning the second point, i.e. for the integration of the rotation speed ω, there are several approaches in the literature to deal with this problem [START_REF] Boyle | The Integration of Angular Velocity[END_REF], but there is no "perfect" one. The default 235 choice that we made in the LMGC90 platform consists in a sort of "linearized" method. Basically, each vector of the inertia frame e k i at the step i is transformed by a rotation increment hω as follows:

e * k+1 = e k i + (1 -θ)hω i × e k i ( 15
)
which are used to predict the orientation. After solving the contact forces, Equation ( 14) is used to update the velocities:

e k i+1 = e * k+1 + θhω i+1 × e * k . (16)
Equations ( 15) and ( 16) require a re-orthogonalization of the main inertial frame. Note that an alternative method

In order to assess the accuracy of the strategy used to integrate rotations, two test cases are considered: (1) the free rotation of a polyhedron combining different main axes of inertia and (2) a collision test between two polyhedra in the presence of gravitational force considering a normal restitution coefficient equal to 1. For both cases, we evaluate, at each moment, the total energy of the system, noted E T . We note E T,0 the total energy of the system at the initial instant (t=0). The relative error, noted ∆E T , expressed as a percentage, is given by:

∆E T = 100 E T -E T,0 E T,0 . (17)
Figure 4(a) shows the evolution of ∆E T for the first benchmark (reduced in this case to rotational energy) after 150 minutes of simulated time, as a function of the time step and the initial rotation speed. As expected, ∆E T increases almost linearly with the time step due to the linearization proposed by the equation [START_REF] Hahn | Realistic animation of rigid bodies[END_REF]. It can also be noted that ∆E T increases with the increase of ω 0 (two orders of magnitude at the same time step).

The figure 4(b) shows the maximum error ∆E max T when two rotating polyhedra collide (measured over a period of 150 minutes) depending on the time step and the initial rotation speed. Again, as the θ method is a first-order integrator, the total energy is not exactly conserved since, in addition to the time integration of the rotations, collisions introduce approximation errors into the models. Note that the vertical dashed line shown in both figures materialize the time step used in the simulations presented in Sec.3 (see corresponding discussion below).

For a more technical explanation of the implementation of the method, see the chapter 9 in [START_REF] Radjai | Discrete-element modeling of granular materials[END_REF]. The reliability and robustness of this method have been tested for several years in previous applications to granular materials, see for example [

Cohesive contact laws

It has been proposed that the smallest particles in an aggregate could constitute a sort of granular bridges between the larger particles [START_REF] Sánchez | The strength of regolith and rubble pile asteroids[END_REF]. As such, the finest material in an asteroid would form a weak cohesive matrix that binds the larger rocks and boulders. By implicitly taking into account this phenomenon, cohesion between large particles can be implicitly represented by a Maugis-Dugdale law, based on a threshold force at separation, denoted r 0 . In the first order, it can be considered constant but it could be dependent on the particle size, shape, cleanliness, rugosity or others characteristics of the systems such as temperature, quantities of water and the presence of volatiles among the most important [START_REF] Mitchell | Fundamentals of Soil Behavior[END_REF]. Such an attractive force r 0 is activated when particles are close enough, leading to the definition of an activation distance, denoted d 0 .

In the CD framework, cohesive forces are easily introduced by shifting upwards the Signorini's condition (see Fig. 5b) and by replacing r n by r * n = r n + r 0 both in Eq. ( 9) leading to a new set of inequalities:

(r n + r 0 ) ≥ 0 v n ≥ 0 (r n + r 0 ).v n = 0 , (18)
for δ n ≤ d 0 . When δ n > d 0 , r n vanishes. shifted upwards from a value of -r 0 (b).

Gravitational forces

The implementation of gravitational forces is of paramount importance when simulating granular asteroids. These forces will be distinctly important for asteroid-size systems, but different algorithms exist depending of the number of particles that form the system. The first form is a direct calculation of complexity O(N 2 ) which is best suited for small system of only up to a few thousands of particles. Next come tree codes algorithms, originally developed by Barnes

and Hut [START_REF] Barnes | A hierarchical O(N log N) force-calculation algorithm[END_REF] and later modified in order to improve the formation of the hierarchical tree (O(N log N)). A variation of a tree-code algorithm was implemented by Sánchez [START_REF] Sánchez | Simulating asteroid rubble piles with a self-gravitating soft-sphere distinct element method model[END_REF] which used a static, regular-grid to take advantage of the packing of a granular system in a condensed phase. This method has been called a "static domain decomposition technique." This method is of order O(N 2 ), but is at least an 10 times faster than the direct calculation.

Regardless of the algorithm used to compute gravitational forces on the particles, applying them within the structure of the CD-method it is very simple. In this context, gravitational forces are seen as external forces, which are taken into account directly through the explicit calculation of the free velocities v α, f ree as shown in Eq.7.

In this work the Open Source Python library written by Mike Grudić, called pykdgrav, is used3 [START_REF] Guszejnov | Evolution of giant molecular clouds across cosmic time[END_REF]. This routine is not to be confused with the PKDgrav code which has been used for many years in the Planetary and Space Science community to carry out research on asteroids, comets and planetary rings. Grudić affirms however, having taken the idea of implementing a kd-tree instead of an octree from PKDgrav. pykdgrav, which has been validated in [START_REF] Guszejnov | Evolution of giant molecular clouds across cosmic time[END_REF][START_REF] Gurvich | Pressure balance in the multiphase ism of cosmologically simulated disc galaxies[END_REF], is a package that implements the Barnes-Hut method [START_REF] Barnes | A hierarchical O(N log N) force-calculation algorithm[END_REF] for computing the combined gravitational field and/or potential of N particles. The kd-tree is implemented as a numba jitclass to achieve much higher performance than the equivalent pure Python implementation. The package implements OpenMP multithreading, but no support for higher parallelism is implemented at the moment.

LMGC90 Implementation

The CD framework has been implemented in the LMGC904 open-source platform under open-source license CECILL, initially developed by Jean and Dubois ([43]5 ) to address multi-contact problems. The software benefits today from a number of contributions such as contact detection between polyhedra by Saussine et al. [START_REF] Saussine | Modelling ballast behaviour under dynamic loading. part 1: A 2d polygonal discrete element method approach[END_REF], solver parallelization with OpenMP by Renouf et al. [START_REF] Renouf | A parallel version of the Non Smooth Contact Dynamics algorithm applied to the simulation of granular media[END_REF], thermal and electrical coupling effects by Renouf [START_REF] Zeng | Numerical investigation on the electrical transmission ability of a shearing powder layer[END_REF], coupling with fluids [START_REF] Constant | Implementation of an unresolved stabilised FEMDEM model to solve immersed granular flows[END_REF] or various type of contact interactions such cohesive zone models [START_REF] Blal | Criteria on the artificial compliance inherent to the intrinsic cohesive zone[END_REF] for simulating, for example, breakable particles.

In order to facilitate the evolution and favour the adaptability of the code, it has been designed using an object oriented approach and is managed via a python interface. Naturally, particles are objects, however, they are not the only entities that benefit from the pseudo-orientation object of the code. In fact, contacts between two particles and their behaviour are also considered as objects. Each of these objects has data, as well as methods for input, output, display, pre-detection and contact detection or explicit and implicit resolution. Practically, the code proposes a global pattern to model and solve the problem, with various functionalities, which the user may enrich with its own routines through plug points. Basically the user needs to define the modelling ingredients (the bulk behaviour of the components of the system, the kind of interactions between these components among others) and the numerical strategy to simulate the evolution of the system.

Table 2 shows how the CD method is used in the platform. The main steps exposed in the previous sections are presented as a pseudo algorithm.

Application to the Simulation of Granular Asteroids

In order to highlight the potential of the CD method, through the LMGC90 software, in the modelling of asteroids as self-gravitating and cohesive granular assembly, different types of simulations have been carried out: 1. the gravitational accretion of spherical and polyhedral particles and, 2. the spin-up processes of spherical self-gravitating aggregates. First, the accretion process is presented, detailing the different parameters and initial conditions used, both for aggregates made of spheres or polyhedra. In a second step, starting from the final configuration of the accretion process of aggregates of spheres, simulations reproducing the spin-up process are performed and their results discussed.

Packing preparation

Our numerical samples are composed of n s solid particles with density ρ = 5 164kg.m -3 and with diameters (i.e., circumscribed spheres in the case of polyhedral shape particle) d varying between d min = 7.5m and d max = 10m with a uniform distribution of particle volume fractions, where d min (resp. d max ) is the minimum (resp. maximum) particle

System initialization

                                                                                                       Time step loop: i = i + 1
Position prediction:

q * = q i + h(1 -θ) qi Rotation prediction: e * k+1 = e k i + h(1 -θ)ω i × e k i
Gravitational Force computation: F G Free Velocity computation: ( q f ree , ω f ree )

Contact detection in configuration q * , e * k+1

Contact problem initialization

                                 NSGS iterations: k = k + 1                  Contact loop: α = α + 1 RHS computation: b k α Local contact problem resolution: (v k+1 α , r k+1 α ) solution                 

Convergence test

Position correction: qi+1 = q f ree + M -1 Hr i+1

Rotation correction:

ω i+1 = ω f ree + J -1 MHr i+1
Table 2: Pseudo code of the CD approach diameter. Particles are randomly deposited in a cubic grid of mesh size equal to d max ; see Fig. 6. We consider two groups of particles assemblies: one made of spheres and a second made of polyhedra. Polyhedral particles are build by randomly generating n v vertex on a unit sphere randomly distributed between 10 to 70. The convex hull of these points is created by associating three vertices for each face. This condition implies that the number of faces is simply given by n f = 2n v -4 ( [START_REF] Azéma | Packings of irregular polyhedral particles: Strength, structure, and effects of angularity[END_REF]). Particle platyness and elongation are defined from the ratio of the principal axis length of the particles. In order to avoid undesirable effects induced by anisotropic shape of the particles, platyness and elongation parameters are set in the range [0.5, 1]. This means that, numerically, the sets of vertices are generated until this condition is satisfied. Figure 6(c) shows a few examples of the polyhedral particles generated for the simulations. In order to evaluate the correctness of the physical behavior and the scaling of the simulations, n s is varied in the set {125, 512, 1 000, 3 375, 8 000, 15 625, 27 000} for spheres assembly, and in the set {125, 512, 1 000, 3 375, 8 000} for polyhedra assembly. Thus, in total 12 packings were build, 7 composed of spheres and 5 composed of polyhedral particles. For all the simulations presented in this paper, the restitution coefficients are set to zero. This choice is motivated by experimental observation that suggests that, in a multicontact system, multiple shocks can dissipate large amounts of kinetic energy in very short times, even if the restitution coefficient is large [START_REF] Luding | Studies of columns of beads under external vibrations[END_REF][START_REF] Da Cruz | Rheophysics of dense granular materials: Discrete simulation of plane shear flows[END_REF]. So, since in CD the characteristic observation time (i.e. the time step) is large, compared to the characteristic time of successive collisions, the choice of a restitution coefficient equal to 0 is appropriate. The time step is fixed to 5 × 10 -2 s which is sufficient to obtain less than 1% error on the calculation of the rotations (see dashed line in Fig. 4a) since the maximum rotational velocities of the particles in the simulations presented in section 3.2 do not exceed 5 × 10 -2 rad/s over a few seconds while the mean rotation speed is of order 10 -3 .

Particle accretion process.

In this section we discuss the particle accretion process and the effect of parameters such as particle number, intergranular friction and grain shape on different physical quantities. Numerical efficiency is also discussed in terms of CPU time spent in each phase of the CD-algorithm. Note that, during the accretion simulation cohesion between the particles is set to 0.

Kinetic energy

Figure 7 shows different snapshot taken during the accretion for the 8000-particle polyhedron assembly. As classically observed, particles are approaching each other in a non-homogeneous manner. We observe some zones denser than others with, in particular, the formation of "branches of grain" aligned along the vector joining each corners of the englobing cube to the potential aggregate center. Indeed, near the box corners the grains are farthest from the center and therefore have to travel a longer distance in order to be agglomerated. At the end of the accretion process, in absence of inter-particle cohesion and anisotropic shape of the particles or anisotropic shape of the containing box, the shape of the aggregate is spherical as illustrated in Fig. 7(e) for aggregate composed of polyhedral particles. Figure 7(f) shows also the shape of aggregate made of spheres at the end of the accretion process for different values of n s .

(f)

Connectivity and packing fraction

Figure 9(a) shows the evolution of the coordination number Z (defined as the mean number of contact per particle) as a function of time in assembly of spheres. Note that for the sake of clarity the evolution of Z in aggregates of polyhedra is not shown as the curves are basically the same. First, Z quickly increases from 0 to 1.5 on the time interval [0, 20min], which coincide with the rapid increase of E c . In this interval we do not observe any effects of the number of particles. Binary shocks dominate the dynamics and energy is mainly dissipated through inelastic collisions (we recall that the coefficients of restitution are fixed to 0 in all the simulations). Then, after a very small decrease, we observe a non-linear increase of the coordination number in all simulations, which tends asymptotically to a maximum value Z max . Fluctuations around the mean evolution decline also as the number of particles is increased. asymptotically tends to values close to 5 as n s increases in aggregates of spheres, and from 5 to 7 in aggregates of polyhedra. Note that for polyhedra, face-face contacts are counted three times, face-edge contacts are counted twice while face-vertex and edge-edge contacts once since, as discussed in section 2.5.1, each type of contact involves a different number of geometric constraints. At the same time, we see that Z max declines from 6 to 4 as the inter-particle friction increases in aggregates of spheres. Basically, as described in several earlier studies, in the idealised case of a confined granular assembly without stress gradient, for the so-called "jammed state", Z jammed declines as µ increases tends to 6 as µ tends to 0 for aggregates of spheres (12 for non-spherical particles assembly), and to 4 for large friction values, both in spherical and non-spherical particle aggregates ( [START_REF] Agnolin | Internal states of model isotropic granular packings. I. Assemblies process, geometry and contact networks[END_REF][START_REF] Van Hecke | Jamming of soft particles: geometry, mechanics, scaling and isostaticity[END_REF][START_REF] Azéma | Packings of irregular polyhedral particles: Strength, structure, and effects of angularity[END_REF]). Obviously, in the presence of bulk forces (i.e. gravitational forces in our case), we expect Z max < Z jammed since less geometrical and mechanical constrains are necessary in order to reach acceptable equilibrium positions, in particular for particles surrounding the aggregates. But we see here, that for gravitational aggregates, Z max is largely dependent on n s ; therefore, we must carefully calibrate the number of particles in order to carry out the most realistic simulations.

(a) (b) with n s . For a same number of particles, the size of the aggregate made of spheres is larger than that of the aggregate made of polyhedra since the polyhedra volume is lower to that of spheres. Nevertheless, it is interesting to note that, by normalizing D max by d * , the mean diameter of equivalent spheres of the same volume to that of polyhedra, the data collapses on a same curve (see inset Fig. 10(a)). Finally, we show that D max first increases with the inter-particle friction and then saturates for µ > 0.4 in spheres aggregates (see inset Fig. 10(b) for n s = 3375).

Along the same line, the variations of the final state value of the packing fraction ν (defined from the total volume of the particles divided by the volume of the encompassing sphere and averaged over t > 125min), is shown in Fig. 10(b) as a function of the number of particles n s , both for frictionless spheres and polyhedra aggregates. In contrast to the evolution of the final-state coordination number, the packing fraction remains slightly dependent on the number of particles. Also, as often observed in the literature, and consistently with the variation of D max with µ, ν declines simulation data via the contacts position and forces. Let's start with the internal moment tensor M p of each particle p, defined by ( [START_REF] Moreau | Some numerical methods in multibody dynamics: application to granular materials[END_REF][START_REF] Andreotti | Granular media: between fluid and solid[END_REF]):

M p i j = c∈p f c i r c j , (19)
where f c i is the i component of the force exerted on particle p at contact c, r c j is the j component of the position vector of the same contact (vector joining the particle center to the contact point), and the summation runs over all contacts c of particle p. The average stress tensor σ in a given volume V is then defined by ( [START_REF] Moreau | Some numerical methods in multibody dynamics: application to granular materials[END_REF][START_REF] Andreotti | Granular media: between fluid and solid[END_REF]):

σ = 1 V p∈V M p , (20)
In the case of particle accretion, the definition of the volume in which particles evolve is neither trivial nor unique. On the contrary, once agglomerated, and given the previous discussion on the average shape and size of the agglomerates at the end of the accretion, it is quite natural to consider the volume of the encompassing sphere. Thus, and just for the sake of convenience, the volume we consider is that of the encompassing sphere at each instant t. The mean stress in 3D is then given by P = (σ 1 + σ 2 + σ 3 )/3, where σ 1 , σ 2 and σ 3 are the principal stresses.

Figure 11(a) shows the evolution of P as a function of time in an assembly of spheres at zero friction. Again, for the sake of clarity, the evolution of P in assembly of polyhedra is not shown as the curves remain the same as for assembly of spheres. We see that P quickly increases to a peak value which depend on n s . As discussed before, in the first instance the dynamics are mainly governed by rapid collisions within the assembly. The number of collisions increases as the number of particles increases too, explaining the increase of the value of the peak-pressure as n s is increased. Then, P quickly decreases around a mean constant value which fluctuations declines as the time is increased.

Figure 11(b) shows P , defined from the mean values of P beyond t > 125min, as a function of D max , both for spheres and polyhedra assembly. The inset shows a log-log representation. We observe that P increases linearly with D max for both, spheres and polyhedra. Interestingly, the data collapses along a single common straight line regardless of grain shape. Let's consider the ideal case of a "non-granular" asteroid, i.e. an asteroid made of a single monolithic rock with a homogeneous distribution of its density. Following the Gauss theorem, saying that the gravitational field g(r) at distance r from the center of a spherical mass of diameter D max is identical to that of a material point, where all the mass of the corresponding sphere would be concentrated in this point, we get:

g(r) = g 0 r D max , (21)
where g 0 = GM a /D 2 max with G the gravitational constant and M a = p∈V m p the total mass of the asteroid, with m p the mass of a particle p. Then, if we assume that the interior of the asteroid is in hydrostatic equilibrium like a liquid mass floating in a vacuum, the differential of the pressure can be related to the gravitational field by dP = -ρ 0 g(r)dr, which leads to the following equation for the stress gradient within the granular asteroid:

P(r) = 1 4 ρg 0 D max (1 -4 r 2 D 2 max ). ( 22
)
As a result, the mean pressure P is given by:

P = 2 D max D max /2 0 P(r)dr = 1 6 ρg 0 D max . (23)
Finally, assuming a small particle size distribution, the total mass M a of the asteroid can be rewritten as a product of characteristic density ρ * and the volume of the encompassing spheres as M a = ρ * πD 3 max /6, with ρ * = n s ρ 0 d 3 . By replacing this relation in the expression of g 0 we deduce that P ∝ D 2 max . As we can see in Figure 11, the prediction given by Eq.23 is in good agreement with our numerical data both for aggregates of spheres and polyhedra. The discrepancy we see for the largest diameters can be explained by the very simplifying assumptions made.

Numerical efficiency

Finally, in this section we discuss more numerical details about the time spent in each step of the Contact Dynamics algorithm in order to evaluate the numerical method (and its implementation). The average time spent in the main time-consuming parts of a time step is measured as a function of the number of particles. Figure 6 presents, in loglog scale, the average CPU time used during contact detection and resolution (Gauss-Seidel algorithm) as well as in the pykdgrav library as a function of the number of particles present in the simulation in spheres assembly (a) and polyhedra assembly (b).

In both systems, the CPU-time spent in contact detection and gravitational forces calculation is almost linear (in log-log scale) with the number of particles. In spheres assemblies the CPU-time spent for gravitational forces calculation is longer than the CPU-time spent in contact detection. We observe the contrary for polyhedra assemblies in which the CPU-time spent for contact detection is larger to the CPU-time spent in gravitational forces calculation.

Indeed, the time spent in the pykdgrav library is naturally independent of particle shape, while the procedure to detect contact between particles is much more complex to that between spheres as discussed in Sec.2.5.1.

In both cases, for spheres and polyhedra assemblies, most of the time is spent solving the contact problem. This

Spin-up process

The test conditions are close to those used in the literature [START_REF] Sánchez | Dem simulation of rotation-induced reshaping and disruption of rubble-pile asteroids[END_REF]. All the simulations consist of 3 375 spheres and have as initial state the final state obtained after the accretion process. An initial rotation speed of 2.7 10 -4 rad.s -1 is given to the sample around the Oz axis passing through the centre of gravity of the sample. Friction and normal restitution are set to zero. Then, every 50 minutes, an increment of 9 × 10 -5 rad.s -1 is given to the sample. The difference between each simulation is related to the value of the cohesive force used. The latter is varied in the set f 0 ∈ {1N, 10N, 100N, 1 000N}. The behaviour of the system over time is observed through different quantities.

Figure 13 shows four configurations of the most cohesive aggregate during its evolution over time. In order to better appreciate the deformation of the aggregate, it is presented in the (O; y, z) (right) and (O; y, z) (left) planes.

When the speed of rotation remains low, the aggregate keeps a more or less spherical shape (a). Above a certain threshold, the aggregate becomes more deformed. This results in a collapse in the direction of the axis of rotation and a relatively uniform expansion in the orthogonal plane (b and c). Then, high porosities appear in the sample, leading to a complete disruption (d).

The effect of local cohesion modifies the time at which the phenomena described above appear. Basically, the collapse of the cohesive self-gravitating aggregate occurs earlier as f 0 tends to 0. Figure 14 shows the evolution of As cohesion is increased, a slower decrease in the value of the semi-axis in the direction of the axis of rotation and a reduction in the gap between the two main semi-axes occurs. Finally, the simulation time from which the values of the semi-axes diverge by more than 3 mean diameters increases with the value of cohesion, from 290min for f 0 = 1N to 354min for f 0 = 1 000N. As expected, cohesion has the effect of limiting radial expansion and preserving symmetry for the aggregate as long as a critical spin velocity is not reached, in agreement with previous research [START_REF] Sánchez | Disruption patterns of rotating self-gravitating aggregates: A survey on angle of friction and tensile strength[END_REF][START_REF] Zhang | Rotational Failure of Rubble-pile Bodies: Influences of Shear and Cohesive Strengths[END_REF][START_REF] Leisner | An extended parameter space study of the effect of cohesion in gravitational aggregates through spin-up simulations[END_REF].

Our data reveals also that the critical spin rate increases with cohesion.

The loss of symmetry in the radial plane is shown in Figure 15 showing the velocity field within the aggregate for the four cohesion values at two different simulation times. After 400min of simulated time, the velocity field remains homogeneous whatever the cohesion value used. After 560min of simulated time, the different systems no longer present the same appearance. For the weakest cohesion values, some particles are ejected from the aggregate (see (a)

and (b) in Fig. 15). For the largest cohesions, strong deformations are observed (see (c) and (d) in Fig. 15), but no particles are yet in free flight. It can also be noted that the aggregate takes on a three-pointed star shape, with the length of the branches being more or less pronounced as a function of the cohesion intensity.

Figure 16 shows the evolution of the mean kinetic energy of the system (a) and the evolution of the mean pressure P within the assembly (b) both as a function of time and for the different values of f 0 . The mean pressure P is computed following Eq.20 in which the total volume V is computed as the volume of the encompassing ellipse. The small 500 jumps in kinetic energy observed in Fig. 16(a) all along the simulation are directly related to the rotation increments.

Basically, the evolution of the kinetic energy is almost identical for all cohesion values, which is : first a non linear increases, until a divergence is observed from t > 800min. Conversely to the mean kinetic energy, P declines as the time is increased and reach a value of 0 precisely from t > 800min. We note also that fluctuations around a mean tendency increases as f 0 increases too. Finally, Fig. 17 shows the evolution of (a) the solid fraction ν, and (b) the evolution of the coordination number Z, both as a function of time. The solid fraction remains constant during the first third of the simulation and then starts to decrease during the second third (see Fig. 17a). This decrease is characterized by small oscillations related to the pulses given to the aggregate to spin it up (see insert in Fig. 17a). This is followed by a rapid decrease, characterizing the disruption of the aggregate. A strong cohesion has the effect of erasing these oscillations. This results in greater variations in the evolution of the solid fraction, leading to a less compact system than low cohesive systems.

The evolution of Z is also affected by the value of the particle-particle cohesive forces within the system. For the two lowest values, there is an almost rapid decrease up to half of the simulation time, then this decrease is combined with a periodic fluctuation, where the period is related to the velocity increments given to the aggregate. For the cohesive forces of 100 N, the behaviour is similar but the coordination number remains greater. Finally, for the greatest cohesion, it remains constant over three quarters of the simulation and then decreases and increases over the last quarter. The velocity increments thus induce a destabilisation of the sample, leading to a radial expansion of the sample. When the cohesion is sufficiently large, it remobilises the contacts ensuring a re-compaction of the system.

Otherwise, the particles move more freely around the periphery, probably to the point of sample rupture.

Conclusion

In this paper, we have presented a numerical method and its implementation in an open-source code (LMGC90)

for the simulation of gravitational aggregates. The Contact Dynamics method is one of a family of Discrete Element Methods that have been developed during the last four to five decades for the simulation of granular materials, or divided media in general terms. As it has been shown, these methods have also been applied for the simulation of galaxies, gases and man-made structures among the most prominent. In the context of asteroids however, its use is still very rare with less than a handful of papers that have shown its potential as a research tool for Planetary Science research. We have provided a simplified description of the mathematical formalism that was introduced by the creators of the method, its implementation in the LMGC90 simulations code, the parallelisation method that has been implemented and a few examples of gravitational accretion and rotational disruption of asteroid size cohesive aggregates.

One of the main advantages of this code over others that have been in use for the last two decades is the relative ease to simulate non-spherical particles. As it can be readily observed, naturally occurring granular materials are formed by a large number of non-spherical particles with a very wide size distribution. Which means that, contrary to what happens with spherical particles, like the ones commonly used in simulation, any pair of particles can have multiple points of contact and they can all have a different geometry. This, together with the added complexity of the implementation of self-gravity turned the problem of simulating self-gravitating aggregates into a challenging and complex task.

In order to test the performance of the implementation of the numerical method in the code, now that self-gravity was implemented, we ran a number of simulations in which we changed the number of particles of the aggregates as well as their shape and cohesive forces. In these examples, though not fully representative of real asteroids, we wanted to show not only the potential of the numerical method, but also of the computational implementation. In order to do that, the examples were of aggregates formed by spherical and not spherical particles, we varied the friction and cohesive forces between pairs of particles in contact and monitored their evolution through the accretion and disruption processes. Quantities such as internal pressure, which comes from the evaluation of the stress tensor, coordination number, solid fraction, kinetic energy, angular velocity, semi axis of rotation, velocity fields were directly measured or evaluated. All of them were in agreement with previously published research.

This being so, and given the open source platform used by the developers at the University of Montpellier, we now feel confident that this simulation code is ready to be used as a research tool for Planetary Sciences, especially on topics related to small planetary bodies.

Acknowledgments

P.S. would like to acknowledge support from NASA grant 80NSSC18K0491, support from NASA's SSERVI program and support from the Observatory of Paris. We acknowledge the support of the High-Performance Computing Platform MESO@LR. We would like to thank the anonymous referees for their valuable and constructive comments.

Figure 1 :

 1 Figure1: Contact frame (n α , t α , s α ) defined at contact point α between two particles i and j.

Figure 2 :

 2 Figure 2: (a) Coulomb cone of the different possible scenarios for the solution of the 3D frictional contact problem. The circle C t is defined from the intersection between the tangent plane P t and the Coulomb cone. (b) the solution r * (at the iteration k+1) is inside the friction cone. (c) The solution r * is outside the Coulomb cone and should be projected on its boundary.

  r t = -β W-1

Figure 3 :

 3 Figure 3: Different configurations of contacts between two polyhedra leading to a different number of contact points: (a) vertex to face; (b) edge to face; (c) face to face; (d) edge to edge; See[START_REF] Radjai | Discrete-element modeling of granular materials[END_REF][START_REF] Azéma | Packings of irregular polyhedral particles: Strength, structure, and effects of angularity[END_REF]

Figure 4 :

 4 Figure 4: (a) Relative error ∆E T , in percent, on the total energy conservation for the free non-principal rotational motion of one polyhedron, at the end of a simulated time of 150 minutes, as a function of the time step and for increasing rotation speed ω 0 imposed on each axis of the inertia frame. (b) Maximum relative error ∆E max T during the collision of two rotating polyhedrons (for a duration of 150 minutes) as a function of the time step and for increasing rotation speed. Note, the vertical dashed line materialize the time step used in the simulations presented in Sec.3.

Figure 5 :

 5 Figure 5: Signorini condition re-written in terms of impulse and formal velocity for a cohesive contact (a) with a threshold value of -r 0 and (b)

Figure 6 :

 6 Figure 6: Initial position for a packing of a) 1 000 spheres and b) 1 000 polyhedrons. Examples of polyhedra used in the simulations (c)

Figure 7 :

 7 Figure 7: Accretion evolution of a packing made of 8 000 polyhedra (snapshots a, b, c and d). Shape of the final state aggregate made of polyhedra (e) and spheres for n s ∈ [125, 1 000, 8 000, 27 000] (f).

Figure 8 Figure 8 :

 88 Figure8shows the evolution of the mean kinetic energy E c (in spheres assemblies), defined as the total kinetic energy divided by the number of particles, as a function of time for (a) different values of n s and inter-particle friction equal to 0, and for (b) n s = 3375 and inter-particle friction varied from 0 to 0.8. The evolution of E c is almost identical for all the simulations. We observe a quick jump of the kinetic energy to a peak at t 20min followed by a rapid fall-off over a more or less wide time interval depending on the number of particles. Basically this fall-off occur for t ∈ [∼ 20min, ∼ 25min] for n s = 27000 and t ∈ [∼ 20min, ∼ 40min] for n s = 125. The inter-particle friction do not modify the maximum value of the kinetic energy can reach, but the increase in friction has the effect of shortening the time interval over which the first drop in kinetic energy occurs. The rapid increases of E c is mainly due to the very loose nature of the initial-state packing. Then, the kinetic energy relaxes as the time is increased and reach a constant plateau with very small values from t > 125min for all n s . The residual-state value of the kinetic energy is, in average, independent with n s . The inset in figure8(a) show the variation of the kinetic energy peak both in assembly of spheres

Figure 9 (

 9 Figure 9(b) shows Z max , averaged in the stabilized regime (i.e. for t > 125min), as a function of the number of particles in aggregates of spheres (disks symbols) and polyhedra (square symbols) at µ = 0. Z max as a function of the inter-particle friction in spheres assembly is shown in the inset of Fig.9(b). Basically, Z max increases from 3.5 and

Figure 9 :

 9 Figure 9: (a) Evolution of the coordination number Z as a function of the time (i.e during the accretion process) in assembly of spheres at µ = 0. For the sake of clarity the variation of Z in assembly of polyhedra is not shown as the evolution remain basically the same. (b) Coordination number Z max as a function of the number of particles n s averaged in the stabilized regime (i.e. for t > 125min) for frictionless assembly of polyhedra (square symbols) and spheres (disks symbols). The inset in (b) shows averaged values of Z max as a function of µ for n s = 3375 spherical particles. Error bars represent the standard deviation in the stabilised regime. Different colors correspond to different particle numbers.

Figure 10 (

 10 Figure 10(a) displays the variation of aggregates diameter D max , normalized by the mean diameter d (diameter of the circumscribed spheres for polyhedra), as a function of the number of particles n s . Let us recall that, by construction, d is identical between spherical and polyhedral particles aggregates. Basically, D max increases as a power law

Figure 11 :

 11 Figure 11: (a) Evolution of the stress P as a function of time in assembly of spheres for different values of n s in frictionless assembly of spheres. (b) Variation of the mean stress P averaged in the steady state (i.e. for t > 125min) as a function of D max for both, assembly of spheres and polyhedra at µ = 0. The inset show the same data in log-log representation. The red lines are the approximations given by Eq.23. Different colors corresponde to different particle numbers

Figure 12 :

 12 Figure 12: Computation time vs. n s log(n s ) during accretion process for (a) spheres assembly and (b) polyhedra assembly. Different colors corresponde to different particle numbers.

Figure 13 :

 13 Figure 13: Different simulation snapshots for the aggregate with the highest cohesion value in the (O; x, y) and (O; y, z) planes at time equal to a) 160, b) 320, c) 480 and d) 640 min.

Figure 14 :

 14 Figure 14: Deviation of the semi-axes of the dynamically equivalent equal-volume ellipsoid during the spin-up process: (a) f 0 = 1N, (b) f 0 = 10 N, (c) f 0 = 100 N and (d) f 0 = 1 000 N.

Figure 15 :Figure 16 :

 1516 Figure 15: Representation of the velocity field in the (O; x, y) plane at two different simulation times and for the different cohesion values: (a) 1 N, (b) 10 N, (c) 100 N and (d) 1 000 N

Figure 17 :

 17 Figure 17: Evolution of (a) the solid fraction ν and (b) coordination number Z, both as a function of time and for all values of f 0 .

  [START_REF] Azéma | Short-time dynamics of a packing of polyhedral grains under horizontal vibrations[END_REF][START_REF] Rafi | Modelling and analysis of the nimes arena and the arles aqueduct subjected to a seismic loading, using the non-smooth contact dynamics method[END_REF][START_REF] Azéma | Quasistatic rheology, force transmission and fabric properties of a packing of irregular polyhedral particles[END_REF][START_REF] Azéma | Discrete simulation of dense flows of polyhedral grains down a rough inclined plan[END_REF][START_REF] Azéma | Packings of irregular polyhedral particles: Strength, structure, and effects of angularity[END_REF][START_REF] Quezadas | Penetration test in coarse granular material using Contact Dynamics Method[END_REF][START_REF] Mollon | A numerical framework for discrete modelling of friction and wear using Voronoi polyhedrons[END_REF][START_REF] Taforel | Finite element-discrete element coupling strategies for the modelling of ballast-soil interaction[END_REF][START_REF] Lancionni | Seismic vulnerability of ancient stone arches by using a numerical model based on the non-smooth contact dynamics method[END_REF][START_REF] Ouhbi | Railway ballast: Grain shape characterization to study its influence on the mechanical behaviour[END_REF][START_REF] Ferrelec | Analysis of compaction of railway ballast by different maintenance methods using dem[END_REF][START_REF] Gazzani | Modal parameters identification with environmental tests and advanced numerical analyses for masonry bell towers: a meaningful case study[END_REF][START_REF] Artoni | Fracture of granular materials composed of arbitrary grain shapes: A new cohesive interaction model[END_REF][START_REF] Orozco | Discrete-element model for dynamic fracture of a single particle international[END_REF][START_REF] Cantor | Microstructural analysis of sheared polydisperse polyhedral grains[END_REF][START_REF] Clementi | Damage assessment of ancient masonry churches stroked by the central italy earthquakes of 2016 by the non-smooth contact dynamics method[END_REF].

https://git-xen.lmgc.univ-montp2.fr/lmgc90/lmgc90 user

proposed by Hugues and Winget[START_REF] Hugues | Finite rotation effects in numerical integration of rate constitutive equations arising in large-deformations analysis[END_REF], which preserves both the orthogonality and the metrics of the main frame of inertia is also implemented in the LMGC90 platform. Both methods give essentially the same results.

https://github.com/mikegrudic/pykdgrav

In French "Logiciel de Mécanique Gérant le Contact" written in Fortran90

https://git-xen.lmgc.univ-montp2.fr/lmgc90

Stress transmission

It is also interesting to evaluate the evolution of the mean pressure within the assembly during the accretion process. This can be done via the granular stress tensor σ at any stage of accretion, which is calculated from the