LMGC90: a Contact Dynamics open source code for the simulation of granular asteroid with realistic regolith shapes. Application to the accretion process

Granular asteroids are naturally occurring gravitational aggregates (rubble piles) bound together by gravitational forces. For this reason, it is reasonable to use the theoretical concepts and numerical tools developed for granular media to study them. In this paper, we extend the field of applicability of the Contact Dynamic (CD) method, a class of non smooth discrete element approach, for the simulation of three dimensional granular asteroids. The CD method is particularly relevant to address the study of dense granular assemblies of a large number of particles of complex shape and broad particles size distribution, since it does not introduces numerical artefacts due to contact sti↵ness. We describe how the open source software LMGC90, interfaced with an external library for the calculation of self-gravity, is used to model the accretion process of spherical and irregular polyhedral particles.


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" [1] 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 [2] and the two that are taking place at the time of this writing are the JAXA Hayabusa 2 mission to asteroid (162173) Ryugu [3] and the NASA OSIRIS-REx mission to asteroid (101955) Bennu [4]. 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 [2,5].
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 an important part of the research [6][7][8].
The common denominator of DEM is to consider the degrees of freedom associated with the 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: Smooth approaches and Non-smooth approaches Smooth approaches, also called Soft-Sphere ⇤ e-mail: diego.sanchez-lana@colorado.edu DEM (SSDEM) introduced in the 70s by Cundall [9] contact laws are modeled as continuous force linking contact forces to the geometric configuration of touching particles. In Non Smooth approaches, which main method is the Contact Dynamics (CD) introduced in the 90s by J.-J. Moreau [10,11], the interactions between particles are described by non-smooth contact laws instead of force laws.
In this paper, we use the CD method to model the accretion of spherical and polyhedral particles under selfgravity. The general ingredients of the CD method together with specific ingredients needed for modeling granular asteroids (polyhedral shape and self-gravity implementation) are first introduced. Then, the main results of accretion simulations are discussed and some perspectives are given regarding forthcoming research.

Main ingredients
The Contact Dynamics (CD) method integrates the equations of motion of a collection of rigid objects considering no elastic repulsive potential, no smoothing of the Coulomb friction law and dissipative collisions. The frictional contact interactions are reformulated as complementarity relations between the relative velocities between particles and the corresponding momenta at the contact points. The condition of geometrical contact between two particles is expressed by the following mutually exclusive equations, known as "Signorini's conditions": See [12,20] where f n is the contact normal force and v n , counted positive when the particles move away from each other, the contact relative velocity. In the same way, the Coulomb friction law can be reduced to the following set of inequalities: where v t is the sliding velocity at the contact, µ is the friction coe cient and f t is the friction force. Note that the Signorini Condition corresponds to a pure inelastic shock. The contact law is then completed by two coecients of restitution (normal e n and tangential e t ) controlling the amount of energy dissipated during collisions. The use of non-smoothed contact laws is an advantage of the method, since it allows the utilization of bigger time steps and greater numbers of grains in comparison with other discrete element methods. The CD method involves an e cient iterative algorithm based on a nonlinear Gauss-Seidel scheme for solving the equations of dynamics of the particles, satisfying at the same time the complementarity relations for contact forces and particle velocities at the end of each time step [12,13]. We use the LMGC90 open-source platform initially developed in Montpellier by Jean and Dubois [14] that implement the CD method. The software benefits today from a number of contributions such as contact detection between polyhedra [15], solver parallelization with OpenMP [16], thermal and electrical coupling e↵ects [17], coupling with fluids [18] or various type of contact interactions.

Regolith as polyhedral particles
For the contact detection between two polyhedral particles we use the common plane method introduced by Cundall [19] which consists to determine 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.
As shown in Fig.1, di↵erent situations may arise: contact point, contact line or contact surface. Point contacts include face-vertex, edge-edge, vertex-vertex or vertexedge contacts. Note that vertex-vertex and vertex-edge contacts are very rare. But when they occur, the common plane method can give the normal direction. Without any modification of the contact law, a face-edge contact (i.e. contact line) can be represented by two points whereas a face-face (i.e. contact surface) can be replaced by three points since they involve an equivalent number of geometrical unilateral constraints. The contact problem is then solved for the di↵erent points generated using the same strategy as mentioned above.

Gravitational forces
We used the Open Source Python library, written by Mike Grudić, called pykdgrav 1 . 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. In this routine, a kd-tree is implemented through the Barnes-Hut method [21] for computing the combined gravitational field and/or potential of n s particles. The package implements also an OpenMP multithreading.

Accretion of polyhedral particles
Our numerical samples are composed of n s solid particles with density ⇢ = 5,164 kg m 3 and with diameters (i.e., circumscribed spheres in the case of polyhedral shape particle) d varying between d min =7.5 m and d max =10 m with a uniform distribution of particle volume fractions, where d max is the maximum diameter. Particles are randomly deposited in a cubic grid of mesh size equal to d max ; see Fig.2(a). We consider two groups of particles assemblies: one made of spheres and a second made of polyhedra. Polyhedral particles are built by randomly generating n v vertex on a unit sphere randomly distributed between 10 to 70. The platyness and elongation parameters of the particles are set in the range [0.5, 1]. Figure  2

Kinetic energy
Figs.2(a) to 2(d) shows di↵erent snapshots taken during the accretion for the 8000-particle polyhedron assembly.  The inset shows the variations of the maximal value of < E c > as a function of n s in assembly of spheres and polyhedra. Different colors correspond to di↵erent particle numbers. (b) Z max as a function of 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 shows Z max as a function of µ for n s = 3375 in spherical particles.
At the end of the accretion process, in absence of interparticle 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.2(e) for aggregate composed of polyhedral particles. Fig.3 shows the evolution of the mean kinetic energy hE c i (in spherical assemblies), as a function of time for di↵erent values of n s and an inter-particle friction equal to 0. The evolution of hE c i is almost identical for all the simulations. We observe a quick jump of the kinetic energy to a peak at t ' 20 min followed by a rapid fall-o↵ over a more or less wide time interval depending on the number of particles. The initial rapid increases in hE c i is 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 > 125 min for all n s . The inset in figure 3(a) shows the variation of the kinetic energy peak both in assemblies of spheres and polyhedra at µ = 0. We see that hE c i increases exponentially with the number of particles n s regardless of the shape of the particles. The evolution of the kinetic energy is closely related to various dissipative mechanisms at the particles scale, which can be understood through the evolution of the particle connectivity within the assembly.  Di↵erent colors correspond to di↵erent particle numbers. Figure 3(b) shows Z max (where Z is the connectivity number as defined in [20]), averaged in the stabilized regime (i.e. for t > 125min), as a function of n s in aggregates of spheres ( ) and polyhedra (⌅) at µ = 0. Basically, Z max increases from 3.5 and asymptotically tends to values close to 5 as n s increases in aggregates of spheres, and from 5 to 7 in aggregates of polyhedra. At the same time, we see that Z max declines from 6 to 4 as the inter-particle friction increases in aggregates of spheres (inset). In fact, in the idealized case of a confined granular assembly without a stress gradient, Z jammed (i.e. the value of Z in the jammed state) declines with µ and 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. In the presence of gravitational forces, we expect Z max < Z jammed since less geometrical and mechanical constrains are necessary in order to reach equilibrium positions. Here, we see also that Z max is largely dependent on n s . This is because Z is not constant between the bulk and the free surface. Fig.4(a) shows the variation in aggregate diameter D max , normalized by the mean diameter hdi as a function of n s . Note that, D max increases as a power law with n s . For the same number of particles, the size of the aggregate of spheres is larger than that of the aggregate of polyhedra since the polyhedra volume is lower to that of spheres. We note also 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 (inset Fig.4(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.4(b) for n s = 3375). The variations in the final state value of the packing fraction ⌫ (defined from the total volume of the particles divided by the volume of the encompassing sphere), is shown in Fig.4(b) as a function of n s (for frictionless spheres and polyhedra aggregates) and as a function of ⌫ (for spheres in the inset). In contrast to the evolution of the final-state Z max , the packing fraction remains slightly dependent on the number of particles. Also, we see that ⌫ declines quickly as the local friction is increased and tends toward a plateau at larger values of friction (µ > 0.4).

Conclusions
In this paper, we have presented the general ingredients of the Contact Dynamics (CD) method for the simulation of gravitational aggregates, and its implementation in the open-source code LMGC90. The CD method is certainly the most popular discrete element approach based on non-smooth dynamics and have been mostly used for the simulation of granular materials, or divided media in general terms. 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. In particular, one of the main advantages of this code (over others) is the relative ease to simulate non-spherical particles.
In order to highlight the potential of the CD method, through the LMGC90 software, in the modeling of asteroids as self-gravitating granular assemblies, we have presented a parametric study on the aggregation process of spheres and polyhedral particles. We have discussed on the e↵ects of di↵erent parameters such as particle number, inter-granular friction and grain shape. The applicability of CD method in a Planetary Sciences context open new perspectives in the modeling of realistic self-gravitating systems by incorporating a wide range of particle sizes and shapes, and various sources of cohesion (which are generally coupled with particle size) for di↵erent dynamical scenarios (not discussed here).