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

Paul Sánchez, Mathieu Renouf, Emilien Azéma, Rémy Mozul

## ► To cite this version:

Paul Sánchez, Mathieu Renouf, Emilien Azéma, Rémy Mozul. LMGC90: a Contact Dynamics open source code for the simulation of granular asteroid with realistic regolith shapes. Application to the accretion process. Powders & Grains, Jul 2021, Buenos Aires, Argentina. pp.14007, 10.1051/epj-conf/202124914007. hal-03259852

HAL Id: hal-03259852

https://hal.science/hal-03259852

Submitted on 14 Jun 2021

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# LMGC90: a Contact Dynamics open source code for the simulation of granular asteroid with realistic regolith shapes. Application to the accretion process

Paul Sánchez<fnref n="1" />,*, Mathieu Renouf<fnref n="2" />, Emilien Azéma<fnref n="2" />,<fnref n="3" />, and Rémy Mozul<fnref n="2" />

**Abstract.** 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 stiffness. 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.

## 1 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–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

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 self-gravity. 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.

## 2 The Contact Dynamics method for granular asteroid

### 2.1 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”:

$$f_n \ge 0 \quad v_n \ge 0 \quad f_n \cdot v_n = 0, \quad (1)$$

A video is available at https://doi.org/10.48448/ymga-dr73

<footnote n="1">¹Colorado Center for Astrodynamics Research, University of Colorado Boulder, 3775 Discovery Dr, Boulder, CO 80303, USA</footnote>

<footnote n="2">²LMGC, Université de Montpellier, CNRS, Montpellier, France</footnote>

<footnote n="3">³Institut Universitaire de France (IUF), Paris, France</footnote>

<footnote>*e-mail: diego.sanchez-lana@colorado.edu</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Different types of contacts between two polyhedra: (a) vertex to face; (b) edge to face; (c) face to face; (d) edge to edge; See [12, 20]</figcaption></figure>

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:

$$ \left\{ \begin{array}{l} \|\mathbf{v}_t\| = 0 \quad \Rightarrow \quad -\mu f_n \le \|\mathbf{f}_t\| \le \mu f_n \\ \|\mathbf{v}_t\| \ne 0 \quad \Rightarrow \quad \|\mathbf{f}_t\| = \mu f_n \end{array} \right. \qquad (2) $$

where $v_t$ is the sliding velocity at the contact, $\mu$ is the friction coefficient 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 coefficients 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 efficient 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 effects [17], coupling with fluids [18] or various type of contact interactions.

## 2.2 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, different situations may arise: contact point, contact line or contact surface. Point contacts include face-vertex, edge-edge, vertex-vertex or vertex-edge 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

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Accretion evolution of 8 000 polyhedra (snapshots a, b, c and d). Shape of the final state aggregate made of polyhedra (e). Examples of polyhedra used in the simulations (bottom).</figcaption></figure>

solved for the different points generated using the same strategy as mentioned above.

## 2.3 Gravitational forces

We used the Open Source Python library, written by Mike Grudić, called *pykdgrav*<fnref n="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.

# 3 Accretion of polyhedral particles

Our numerical samples are composed of $n_s$ solid particles with density $\rho = 5,164 \text{ kg m}^{-3}$ and with diameters (i.e., circumscribed spheres in the case of polyhedral shape particle) $d$ varying between $d_{min}=7.5 \text{ m}$ and $d_{max}=10 \text{ 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(bottom) shows a few examples of the polyhedral particles generated for the simulations. The number of particles is varied in the set {125, 512, 1 000, 3 375, 8 000, 15 625, 27 000} for the spherical assembly, and in the set {125, 512, 1 000, 3 375, 8 000} for the polyhedral assemblies. For all the simulations presented in this paper, the restitution coefficients are set to zero. Particles velocities are also set to zero at the beginning of the simulations.

### 3.1 Kinetic energy

Figs.2(a) to 2(d) shows different snapshots taken during the accretion for the 8000-particle polyhedron assembly.

<footnote n="1">¹https://github.com/mikegrudic/pykdgrav</footnote>
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

**Figure 3.** (a) Evolution of $\langle E_c \rangle$ as a function of time at $\mu_s = 0$. The inset shows the variations of the maximal value of $\langle E_c \rangle$ as a function of $n_s$ in assembly of spheres and polyhedra. Different colors correspond to different particle numbers. (b) $Z_{max}$ as a function of $n_s$ averaged in the stabilized regime (i.e. for $t > 125\,\text{min}$) for frictionless assembly of polyhedra (square symbols) and spheres (disks symbols). The inset shows $Z_{max}$ as a function of $\mu$ for $n_s = 3375$ in spherical particles.

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.2(e) for aggregate composed of polyhedral particles.

Fig.3 shows the evolution of the mean kinetic energy $\langle E_c \rangle$ (in spherical assemblies), as a function of time for different values of $n_s$ and an inter-particle friction equal to 0. The evolution of $\langle E_c \rangle$ is almost identical for all the simulations. We observe a quick jump of the kinetic energy to a peak at $t \approx 20$ min followed by a rapid fall-off over a more or less wide time interval depending on the number of particles. The initial rapid increases in $\langle E_c \rangle$ 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 $\mu = 0$. We see that $\langle E_c \rangle$ 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.

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

**Figure 4.** (a) Log-log representation of $D_{max}/d$, as a function of $n_s$ in assembly of spheres (disks symbols) and polyhedra (square symbols) at $\mu = 0$. (b) $\nu$ as a function of $n_s$ in assembly of spheres (●) and polyhedra (■) at $\mu = 0$. Inset: variation of both, $\nu$ and $d_{max}$ as a function of $\mu$ in aggregates of spheres for $n_s = 3375$. Different colors correspond to different particle numbers.

## 3.2 Connectivity and packing fraction

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 > 125\,\text{min}$), as a function of $n_s$ in aggregates of spheres (●) and polyhedra (■) at $\mu = 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 $\mu$ and tends to 6 as $\mu$ 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 $\langle d \rangle$ 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 $\mu > 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 $v$ (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 $v$ (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 $v$ declines quickly as the local friction is increased and tends toward a plateau at larger values of friction ($\mu > 0.4$).

## 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 effects of different 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 different dynamical scenarios (not discussed here).

## References

[1] L.D. Friedman, Vermin of the sky, http://www.planetary.org/blogs/guest-blogs/lou-friedman/20130219-vermin-of-the-sky.html (2013)

[2] A. Fujiwara, J. Kawaguchi, D.K. Yeomans, M. Abe, T. Mukai, T. Okada, J. Saito, H. Yano, M. Yoshikawa, D.J. Scheeres et al., Science **312**, 1330 (2006)

[3] S.i. Watanabe, Y. Tsuda, M. Yoshikawa, S. Tanaka, T. Saiki, S. Nakazawa, Space Sci. Rev. **208**, 3 (2017)

[4] D. Lauretta, M. Barucci, E. Bierhaus, J. Brucato, H. Campins, P. Christensen, B. Clark, H. Connolly, E. Dotto, J. Dworkin et al., *The Osiris-Rex Mission-Sample Acquisitions Strategy and Evidence for the*

*Nature of Regolith on Asteroid (101955) 1999 RQ36* (ACM, Japan, 2012), Vol. 1667, p. 6291

[5] H. Miyamoto, H. Yano, D.J. Scheeres, S. Abe, O. Barnouin-Jha, A.F. Cheng, H. Demura, R.W. Gaskell, N. Hirata, M. Ishiguro et al., Science **316**, 1011 (2007)

[6] S. Tardivel, D.J. Scheeres, P. Michel, S. Van wal, P. Sánchez, J. Spacecraft Rockets (2014)

[7] K. Zacny, E.B. Bierhaus, D.T. Britt, B. Clark, C.M. Hartzell, L. Gertsch, A.V. Kulchitsky, J.B. Johnson, P. Metzger, D.M. Reeves et al., *Geotechnical Properties of Asteroids Affecting Surface Operations, Mining, and In Situ Resource Utilization Activities* in Primitive Meteorites and Asteroids, edited by N. Abreu (Elsevier, 2018), pp. 439–476, ISBN 978-0-12-813325-5

[8] S. Van wal, Y. Tsuda, K. Yoshikawa, A. Miura, S. Tanaka, D. Scheeres, J. Spacecraft Rockets **55**, 797 (2018)

[9] P.A. Cundall, *A computer model for simulating progressive large scale movements of blocky rock systems*, in Proceedings of the symposium of the international society of rock mechanics, **1**, 132–150 (1971)

[10] J.J. Moreau, *Unilateral contact and dry friction in finite freedom dynamics*, in Non Smooth Mechanics and Applications, CISM Courses and Lectures, edited by J. Moreau, e. P.-D. Panagiotopoulos, **302** (Springer-Verlag, Wien, New York), 1–82 (1988)

[11] M. Jean, J.J. Moreau, *Unilaterality and dry friction in the dynamics of rigid bodies collection*, in Contact Mechanics International Symposium, edited by A. Curnier (Presses Polytechniques et Universitaires Romanes, 1992), pp. 31–48

[12] F. Radjai, F. Dubois, *Discrete-element modeling of granular materials* (Wiley-Iste, 2011)

[13] F. Dubois, V. Acary, M. Jean, Comptes Rendus - Mecanique **346**, 247 (2018)

[14] F. Dubois, M. Jean, *The non smooth contact dynamic method: recent LMGC90 software developments and application* in Analysis and Simulation of Contact Problems, edited by P. Wriggers, U. Nackenhorst (Lec Not. App. and Comp. Mech., 2006), pp. 375–378, ISBN 978-0-12-813325-5

[15] G. Saussine, C. Cholet, F. Dubois, C. Bohatier, P. Gautier, J. Moreau, Comput. Methods Appl. Mech. Engrg. **195**, 2841 (2006)

[16] M. Renouf, F. Dubois, P. Alart, J. Comput. Appl. Math. **168**, 375 (2004)

[17] C. Zeng, M. Renouf, Y. Berthier, R. Hamdi, Granul. Matter **18**, 19 (2016)

[18] M. Constant, F. Dubois, J. Lambrechts, V. Legat, Comp. Part. Mech. **6**, 213 (2018)

[19] C. P.A., Int. J. Rock Mech. Mining Sci. Geomech. **25**, 107 (1988)

[20] E. Azéma, F. Radjai, F. Dubois, Phys. Rev. E **87**, 062203 (2013)

[21] J. Barnes, P. Hut, Nature **324**, 446 (1986)