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

TL;DR: The Contact Dynamics method is described, a class of DEM based on non-smooth mechanics, and its implementation in the open-source software LMGC90, and a parallelized kd-tree is implemented to monitor the performance of the code as it simulates a number of granular systems.

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.

## Summary (4 min read)

### 1. Introduction

- During the last two decades, research about asteroids in general and small asteroids in particular has increased dramatically.
- 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 [8, 9].
- 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.
- In Smooth approaches [11, 12], interaction forces between particles can be written as a function linking contact forces to contact kinematics.
- Then, some specificities for the modelling of granular asteroid are in-60 troduced, including cohesion, particle shape and self-gravity.

### 2.1. The Non-Smooth Philosophy

- 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).
- This stiffness imposes very small discretization steps and often artificial inertias or viscosities are introduced to ensure numerical stability [14, 44].75.
- 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.

### 2.2.1. In the body frame

- On the other hand, for complex shapes, the Euler equation remains non-linear and depends implicitly on the orientation of the objects.
- The system is updated according to the second Equation of system (4).
- Note that, the update110 of the orientation of each particle when dealing with non-spherical particle will be discussed in details in Sec.2.5.1. 2.2.2.
- 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.
- The time index is omitted to135 make pleasant reading.

### 2.3. 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.
- The contact unilaterality is described by the well-known Signorini Condition [54] relating the normal component of the contact relative velocity vn and the normal component of the contact impulse rn.
- According to the set of inequalities (9), 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 (vi+1n = −envin).
- The Amontons-Coulomb friction law is one of the simplest way to consider a resistance in the tangential direction.

### 2.4. Resolution of the local frictional contact problem

- If the authors make the difference of two successive iterations appear to the left hand side, then the equation can be written as follows: Wαα(rk+1α − rkα) = vk+1α − vα, f ree − nc∑ β<α Wαβrk+1β − nc∑ β≥α Wαβrkβ (12) The convergence of the Gauss-Seidel algorithm makes the right-hand side member tend towards zero.
- 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 [20].
- 5. Extensions of interest to model granular asteroids180.
- The next subsections will detailed these specific features in the framework of190 Contact Dynamics.

### 2.5.1. Modelling regolith as polyhedral particles

- Grain angularity modeling is of major importance to model realistic granular asteroids in order to better predict their behavior.
- Recent reports show that particle shape strongly195 affects the strength and dilatancy properties of the granular media.
- The second, called accurate, can use a method of intersection of triangular faces [47] or a method of plane separator [44].
- Without any modification of the contact law, a face-edge contact (i.e. contact215 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.
- (17) Figure 4(a) shows the evolution of ∆ET 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.

### 2.5.3. Gravitational forces

- The implementation of gravitational forces is of paramount importance when simulating granular asteroids.
- This method has been called a “static domain decomposition technique.” pykdgrav, which has been validated in [80, 81], is a package that implements the Barnes-Hut method [78] 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 parallelism290 is implemented at the moment.

### 2.6. 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.[47], solver295 parallelization with OpenMP by Renouf et al.[50], thermal and electrical coupling effects by Renouf [49], coupling with fluids [82] or various type of contact interactions such cohesive zone models [83] for simulating, for example, breakable particles.
- In fact, contacts between two particles and their behaviour are also considered as objects.
- 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.
- First, the accretion process is presented, detailing the different parameters and initial conditions used, both for aggregates made of spheres or polyhedra.

### 3.2. Particle accretion process.

- Intergranular friction and grain shape on different physical quantities.the authors.
- 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.345.

### 3.2.1. Kinetic energy

- Figure 7 shows different snapshot taken during the accretion for the 8000-particle polyhedron assembly.
- 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.
- The rapid increases of 〈Ec〉 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 ns.
- 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.

### 3.2.2. Connectivity and packing fraction

- Figure 9(a) shows the evolution of the coordination number Z (defined as the mean number of contact per particle)370 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.
- But390 the authors see here, that for gravitational aggregates, Zmax is largely dependent on ns; therefore, they must carefully calibrate the number of particles in order to carry out the most realistic simulations.
- 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 ns, 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.

### 3.2.3. Stress transmission

- It is also interesting to evaluate the evolution of the mean pressure within the assembly during the accretion process.
- Thus, and just for the sake of convenience, the volume the authors consider is that of the encompassing sphere at each instant t.
- Interestingly, the data collapses along a single common straight line regardless430 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.
- As the authors can see in Figure 11, the prediction given by Eq.23 is in good agreement with their numerical data both for aggregates of spheres and polyhedra.

### 3.2.4. Numerical efficiency445

- Finally, in this section the authors 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).
- In spheres assemblies the CPU-time spent for gravitational forces calculation is longer than the CPU-time spent in contact detection.
- This point has been identified for a long time and has led to the implementation of parallel computing techniques to reduce the time consumed by this part of the method [50, 90].
- The evolution of the time spent in solving the contact problem is460 less linear (in log-log scale) than the other two curves.
- Indeed, the number of iterations of the Gauss-Seidel algorithm is not constant over time, strong variations of this number can disturb the mean value, which explains the non-linear behaviour observed on the Figs 12a and Fig. 12b.

### 3.3. Spin-up process

- The test conditions are close to those used in the literature [31].
- As cohesion is increased, a slower decrease in the value of the semi-axis in the direction of the axis of rotation and485 a reduction in the gap between the two main semi-axes occurs.
- This results in greater510 variations in the evolution of the solid fraction, leading to a less compact system than low cohesive systems.
- 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.
- The authors 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.

Did you find this useful? Give us your feedback

...read more

##### Citations

17 citations

##### References

20,549 citations

3,553 citations

3,283 citations

2,936 citations