A contact dynamics code implementation for the simulation of asteroid evolution and regolith in the asteroid environment

Highlights

• We use an open-source, Non-Smooth Contact Dynamics code (LMGC90) for the simulation of granular asteroids.

• The pykdgrav module is implemented in the LMGC90 code to handle self-gravity.

• The gravitational accretion of spherical and polyhedral particles is simulated and analyzed.

• Spin-up evolution simulations of self-gravitating aggregates are run in order to validate the numerical approach.

Abstract

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 artifacts due to contact stiffness. So that it can be used for Small Body research, we implement a parallelized 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 behavior 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” (Friedman, 2013) 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 (Fujiwara et al., 2006) and the two that are taking place at the time of this writing are the JAXA Hayabusa2 mission to asteroid (162173) Ryugu (Watanabe et al., 2017) and the NASA OSIRIS-REx mission to asteroid (101955) Bennu (Lauretta et al., 2012). 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 (Tardivel et al., 2014, Zacny et al., 2018, Van wal et al., 2018).

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 (Nakamura et al., 2011, Tsuchiyama et al., 2011). 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 (Yano et al., 2006, Nakamura et al., 2011). 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 behavior 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 behavior 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 (Cundall, 1971, Allen and Tildesley, 1987), 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 (Hertz, 1882)). 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) (Allen and Tildesley, 1987), the Newmark (Cundall, 1971), or Verlet scheme (Iordanoff et al., 2002).

In Non Smooth approaches (Hahn, 1988, Moreau, 1988, Jean and Moreau, 1992, Baraff, 1993, Dubois et al., 2018), there is no explicit relation between contact forces and contact kinematics. Contact conditions are described by complementarity relations between the contact forces and displacements or velocities (as the well know Signorini condition (Moreau, 1988)). No regularization is required. The equations of motion are rewritten as non-differentiable relations involving velocity jumps and impulse (Moreau, 1988). Equations of motion can be driven using time stepping or Event Driven (ED) (Radjai and Dubois, 2011) and discretized by the $\theta$-method (Jean and Moreau, 1992) or Leap Frog (LF) method (Richardson et al., 2000). 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 of the 90 s, a first ED (Hard-Sphere DEM) model is proposed to the study of self-gravitating particles of meters in size (Richardson et al., 2000, Richardson et al., 2005). The authors developed a code named PKDgrav (Richardson et al., 2000), 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 (Walsh and Richardson, 2006, Walsh et al., 2008, Cotto-Figueroa et al., 2015, Walsh et al., 2012, Ballouz et al., 2015, Bagatin et al., 2018).

More recently (2008), an MD (Soft-Sphere DEM) has been introduced to model self-gravitating granular systems (Sánchez and Scheeres, 2009, Sánchez and Scheeres, 2011). 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 (Sánchez and Scheeres, 2012, Sánchez and Scheeres, 2014, Sánchez, 2015, Hirabayashi et al., 2015, Sánchez and Scheeres, 2016, Sánchez and Scheeres, 2018, Tardivel et al., 2018). This last point was largely applied to the simulation of the Touch-and-Go Sample Acquisition Mechanism (TAGSAM) for the OSIRIS-REx mission (Zacny et al., 2018). Note that a SSDEM version has been available in the PKDgrav code mentioned previously since 2012 (Schwartz et al., 2012).

Later on, in 2014, a first use of a Non Smooth approach was proposed to study regolith processes (Hartzell and Hunt, 2014). More recently, in 2018, the Contact Dynamics (CD) method, was applied to analyze the strength properties of self-gravitating aggregates of spheres (Azéma et al., 2018). A year later, Ferrari et al. (2019) and Ferrari and Tanga (2020) 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 (Dubois and Jean, 2006)¹ together with a direct application of self-gravity. In this paper, to underline the numerical efficiency of the CD method for modeling granular asteroids, the CD method is applied to model the accretion of spherical and polyhedral particles. Then, some specificities for the modeling 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.