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.

About: This article is published in Icarus.The article was published on 2021-07-15 and is currently open access. It has received 6 citations till now.

Summary

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.

Figures

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 [20, 59]

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: Deviation of the semi-axes of the dynamically equivalent equal-volume ellipsoid during the spin-up process: (a) f0 = 1N, (b) f0 = 10 N, (c) f0 = 100 N and (d) f0 = 1 000 N.

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 10: (a) Log-log representation of the evolution of the aggregate diameter, normalized by particle mean diameter d, as a function of ns in assembly of spheres (disks symbols) and polyhedra (square symbols) at µ = 0. (b) Packing fraction ν of the aggregate as a function of ns in assembly of spheres (disks symbols) and polyhedra (square symbols) at µ = 0. The inset in (b) shows the variation of both, ν and dmax as a function of the inter-particle friction µ in aggregates of spheres for ns = 3375. Different colors correspond to different particle numbers.

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

Table 1: Pseudo-code of the resolution of the 3D frictional contact problem.

Figure 5: Signorini condition re-written in terms of impulse and formal velocity for a cohesive contact (a) with a threshold value of −r0 and (b) shifted upwards from a value of −r0 (b).

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 ns ∈ [125, 1 000, 8 000, 27 000] (f).

Figure 11: (a) Evolution of the stress P as a function of time in assembly of spheres for different values of ns 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 Dmax 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 8: Evolution of the mean kinetic energy, noted < Ec >, during the accretion of spheres (a) for inter-particle friction fixed to 0 and (b) for ns = 3 375 and different values of the inter-particle friction. The inset in (a) shows the relation between the maximal value of < Ec > and ns in assembly of spheres and polyhedra. Different colors corresponde to different particle numbers.

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 4: (a) Relative error ∆ET , 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 ∆EmaxT 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 16: Evolution of (a) the mean kinetic energy, noted < Ec > and (b) the pressure at the centre of the aggregate, as a function of time for different values of f0.

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

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

Figure 2: (a) Coulomb cone of the different possible scenarios for the solution of the 3D frictional contact problem. The circle Ct is defined from the intersection between the tangent plane Pt 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.

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 Zmax as a function of the number of particles ns 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 Zmax as a function of µ for ns = 3375 spherical particles. Error bars represent the standard deviation in the stabilised regime. Different colors correspond to different particle numbers.

Table 2: Pseudo code of the CD approach

Frequently asked questions

Q1. What have the authors contributed in "A contact dynamics code implementation for the simulation of asteroid evolution and regolith in the asteroid environment" ?

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. In this paper the authors 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, Hardand 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, the authors implement a parallelised kd-tree and monitor the performance of the code as it simulates a number of granular systems. The authors provide examples of the simulation of the accretion of self-gravitating aggregates as well as their rotational disruption.