A parallel-GPU code for asteroid aggregation problems with angular particles

12/09/2019 ∙ by F. Ferrari, et al. ∙ NASA Politecnico di Milano 0

The paper presents a numerical implementation of the gravitational N-body problem with contact interactions between non-spherically shaped bodies. The work builds up on a previous implementation of the code and extends its capabilities. The number of bodies handled is significantly increased through the use of a CUDA/GPU-parallel octree structure. The implementation of the code is discussed and its performance are compared against direct N^2 integration. The code features both smooth (force-based) and non-smooth (impulse-based) methods, as well as a visco-elastic non-smooth method, to handle contact interaction between bodies. The numerical problem of simulating "rubble-pile" asteroid gravitational aggregation processes is addressed. We discuss the features of the problem and derive criteria to set up the numerical simulation from the dynamical constraints of the combined gravitational-collisional problem. Examples of asteroid aggregation scenarios that could benefit from such implementation are finally presented.



There are no comments yet.


page 12

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

The numerical resolution of the N-body problem is one of the most popular problems in high performance computing. The intrinsic complexity of the problem has challenged many astronomers and computer scientists since the beginning of the computing era. N-body algorithms have evolved side by side with the increase of computing power and availability of new resources, such as Graphic Processing Units (GPUs). State-of-the-art implementations include parallel N-body tree codes (Morbidelli, 2002; Stadel, 2001; Richardson et al., 2002b), hybrid codes (Aarseth, 1999, 2001; Wang et al., 2015), adaptive algorithms of optimal orders (Pruett et al., 2003), systolic algorithms (Dorband et al., 2003), or more generally symplectic codes (Wisdom and Holman, 1991; Duncan et al., 1998; Chambers, 1999). In order to be able to solve granular dynamics problems, the capability of gravitational N-body codes needs to be extended to include contact interaction between non-point-like bodies. Two main classes of methods are commonly used to deal with contact interactions: hard- and soft-body methods. The first consider impulsive contacts between nondeformable particles (J. and E., 1959; Jean and Moreau, 1987). This method is numerically very stable and allows to simulate effectively the dynamics of hundreds of thousands of bodies. In more recent years, the alternative class of soft-body methods have been developed (Cundall and Strack, 1979). Based on a force-driven approach, soft-body methods are more suitable than hard-body ones for the simulation of smooth dynamics. These could be used to complement some deficiencies of hard-body models. For example, due to the impulsive nature of contact exchanged, hard-body methods are not adequate to reproduce phenomena such as wave propagation in granular media (Gilardi and Sharf, 2002)

. However, the numerical solution of soft-body methods is more unstable and often requires extremely small time steps to adequately reproduce elastic forces at contact. Both hard- and soft-body models have their own advantages and shortcomings and their use shall be carefully weighted upon the application scenario. As a general rule, hard-body methods qualify to simulate non-smooth dynamics, and vice versa: soft-body is suitable for smooth problems. However, whether a problem should be classified in either one category is a controversial field of debate: both models rely on sets of non-directly measurable parameters, which often offer poor physical interpretation when compared to real-world scenarios 

(Dubois et al., 2018). To date, both hard- and soft-body methods have been successfully used to simulate a wide variety of planetary science scenarios: Richardson et al. (2002a) used their hard-body model to simulate gravitational aggregation of rubble-pile asteroids and planetary ring dynamics (Porco et al., 2008); Wada et al. (2006) used a soft-body approach to study cratering impacts under constant gravity field and more recently, soft-sphere methods with mutual gravity between bodies have been used to simulate rubble-pile scenarios (Sánchez and Scheeres, 2011; Tancredi et al., 2012; Schwartz et al., 2012). Another key aspect for a realistic evaluation of contact interactions concerns the shape of interacting bodies. All aforementioned codes are limited to spherically-shaped bodies. The use of spherical shapes is now considered to be a very relevant simplification for problems where contact between bodies happens. The use of spheres allows to drastically reduce the computational burden, but can heavily affect the realism of the model when solving for the contact interaction between bodies (Michel et al., 2004). As such, the determination of the final shape of the aggregates and fragments could not be addressed (Richardson et al., 2009). Few attempts have been made using hard-body polyhedral grains together with video game engines to solve for contact interactions (Movshovitz et al., 2012). Although polyhedral bodies provide a better modeling, compared to spheres, of some physical phenomena (e.g. related to interlocking between bodies) the computational accuracy provided by game engines is very low (single precision) and not able to satisfy basic energy or angular momentum conservation requirements. Recently Ferrari et al. (2017) implemented the N-body aggregation problem with non-spherical shapes and non-smooth contacts, using modules of Chrono::Engine (C::E) (Tasora, 2019; Tasora et al., 2016), an open software optimized for granular and multi-body engineering problems able to handle contacts and collisions of a large number of complex-shaped bodies. The results obtained were very promising and satisfactorily accurate in terms of energy and angular momentum conservation. In their work, Ferrari et al. (2017) made use of direct N-body simulations to solve for gravitational dynamics and only few thousands of bodies could be handled by the simulator.

The work presented here builds upon the work by Ferrari et al. (2017) and discusses the implementation of a parallel CUDA-GPU octree structure to extend the capabilities of code to a higher number of complex-shaped bodies. As its parent-code, the implementation is based on a re-work of the C::E multi-physics simulation engine to simulate collisions between non-spherical bodies and integrate the dynamics of the problem. The code features both non-smooth (impulse-based) and smooth (force-based) methods to handle contact interactions. Also, a compliant non-smooth method is featured and proposed for the first time for planetary science application, in the attempt to joint the advantages of numerical stability provided by its impulse-based formulation, together with its ability of simulating compliance at contact level. Section 2 discusses the general architecture of the code and presents methods available to solve gravitational and collisional dynamics. The performance of the code are discussed in Section 3, where the results of validation/test problems are reported. The applicability to asteroid aggregation problems is eventually discussed through simulation examples. We address the problem of rubble pile asteroid aggregation processes in section 4. The physics of gravitational and collisional processes involved are reviewed and their numerical modeling is discussed. We derive a set of qualitative and quantitative requirements from dynamical and numerical constraints of the gravitational-collisional problem.

2 Numerical implementation

The code is built using a modular approach, taking advantage of the object-oriented C++ architecture of C::E library. Modules represent interchangeable units of software, easily accessible by the user. Each module contains methods and routines to execute specific tasks. The main modules developed and/or directly retrieved by C::E libraries are:

  1. gravity

  2. contact

  3. rigid-body dynamics

  4. body creation

  5. numerical solvers

  6. data input/output

  7. visual interface

  8. post processing

Modules 257 are directly adapted from C::E libraries, modules 34 are based on C::E libraries and further extended with additional methods and routines, modules 168 are fully developed by the authors. Modules 6 and 7 provide routines that allow the user to interface with the code at a higher level than software implementation level. The visual interface is based on Irrlicht library (Gebhardt, 2016) and can be set to highlight desired geometrical or simulation features (forces, contacts), as well as to display real-time simulation outputs. The data input/output interface is straightforward: every user-tunable parameter of the model can be assigned via an input text file that is parsed and acquired by the code, while output data is typically saved into dedicated text files. This section presents a general overview of the numerical methods available in each module, focusing on those relevant to the gravitational aggregation problem.

2.1 Gravitational dynamics

Module 1 encloses methods and routines used to compute the gravity forces acting on each body. These are based on the evaluation of gravity field produced by point-mass or shape-based (e.g. polyhedron, ellipsoid) gravity sources. To simulate the problem of gravitational aggregation, the classical N-body problem with point-mass sources is considered. This is a common and reasonable assumption: for a high number of commensurate interacting bodies, the use of shape-based mutual gravity would impact dramatically on the computational effort and would not provide relevant benefits to the accuracy of results, as discussed in section 4.1. As detailed below, module 1 includes both direct N integration and a parallel-GPU octree method.

The N-body problem

The equations of motion of the N-body problem are typically written by using Newton’s law to compute the gravitational interactions between masses:



represents the position vector of the

-th body in an inertial frame, its mass, the universal gravitational constant and . A method implementing all N-to-N interactions is commonly referred as direct N or particle-particle method (Hockney and Eastwood, 1988; Hut and Makino, 2010). Such method is very simple to implement and provides exact computation of the motion of the N-bodies without any physical approximation. The accuracy of the result depends solely on the numerical errors of the solver. The simple dynamics of direct N method make it very easy to implement and parallelize. However, they are extremely inefficient from a computational point of view: interactions must be computed for each body involved in the simulation and this leads to a time complexity of . A direct consequence is that such codes are limited in terms of maximum number of bodies to be handled. The bottlenecks are simulation time that becomes prohibitively long, and hardware capabilities with computer running out of memory, both effects occurring very rapidly as N increases. The direct N method, even with advanced parallelism on processing units, is not suitable for efficient asteroid aggregation simulations that can involve more than few thousands of bodies.

Hierarchical treecode algorithm

Algorithms based on tree data structures rely on more dynamic and adaptive computations that allow for a significant reduction of time complexity up to . The Barnes-Hut (also referenced as BH in the followings) algorithm (Barnes and Hut, 1986) groups particles using a hierarchy of cube structures. A node in the algorithm corresponds to a cube in physical space. Because of the use of octrees, each node has eight child nodes obtained by a simple homogeneous spatial subdivision performed along the three principal axis of the system. The tree is therefore built by recursive sub-division until each node of the tree contains zero or one particle. The structure is adaptive, implying that the size of the tree is not fixed but comes as a result of the repartition of the particles in the 3D space. The data structure can grow naturally to more levels in regions where the particle density is high. The octree is re-built at each time step 111See section 4.1 for more details on the choice of simulation time step(s) of the dynamic simulation and consists of bodies stored at terminal nodes, called leaf nodes, and intermediary internal nodes that behave as clusters of particles. The algorithm uses the newly obtained adaptive octree data structure to compute the center of mass of the cubes involved in the simulation and the forces applied to each body: body-to-body and body-to-node interactions are considered.

Figure 1: Schematics of Barnes-Hut clustering criteria


The algorithm relies on the idea that the force generated by a cluster of bodies can be approximated by treating the cluster as a single body. The accuracy of the approximations depends on the distance of the cluster from the body and the radius of the cluster of particles, as shown in Figure 1. We can therefore define the accuracy through , a user-tunable parameter used to set the error of the method. The position of center of mass and total mass of the node is computed for each leaf and internal node of the Barnes-Hut tree. For a given internal node :


where () are the coordinates of the position vector along the three principal axes of the system, while and (, , ) are respectively the mass and the coordinates of the position vector along the principal axes of the -th child node of .

Force computation

After the octree is built and all nodes contain the required information, parameter is computed for each body-node pair. (, ):


This is compared with the value

chosen by the user to set the accuracy and performance of the force estimation algorithm. A more detailed discussion on the effects of the user-tunable parameter

on both computational time and accuracy of the algorithm is provided in section 4.1. For each body , the tree is traversed from its root downwards along multiple branches and checks are performed at each encounter with a node , in order to compute force contributions:

  • If is a leaf node, the force contribution is evaluated as a classical body-to-body interaction:

  • If is an internal node and , the traversal is interrupted at this node and the force contribution is evaluated:


    where r is the distance from the particle to the center of mass of the cube softened according to a parameter :


    is called softening length and, as its name suggests, ensures the “softening” of gravitational forces when , to avoid very high values of speed and acceleration that the integrator could not handle. Studies show that, without opting for a time-dependent adaptive softening parameter, must be a factor of at least twice the minimum distance between a body and the attractive cluster involved (Rodionov and Sotnikova, 2005). In the case of finite-size body interaction, with contact and collisions, the softening parameter is not needed, since no overlap between center of mass of a body and a node could occur.

  • If is an internal node and , the traversal continues along all eight child nodes of .

CUDA implementation of GPU-parallel Barnes-Hut

The implementation of the Barnes-Hut algorithm on a GPU using CUDA language is inspired by the work of Burtscher and Pingali (Burtscher and Pingali, 2011). The physical domain is divided into sub-domains and the bodies are grouped following an octree structure. Similar to what done by Burtscher and Pingali (2011), the numerical tasks to follow the Barnes-Hut algorithm have been divided among five kernels, which are written to minimize memory accesses. The tasks are executed sequentially on the GPU as follows:

  • compute bounding box around all bodies (root of the octree)

  • build the octree

  • compute mass properties of each node

  • sort bodies by distance (to speed up force evaluation)

  • compute forces

Each body is assigned to a thread of the GPU, which traverses the octree to compute all force contributions acting on the body. Thanks to sorting, neighbor bodies require the traversal of a similar sub-part of the octree, whereas distant bodies may require completely different traversals.

2.2 Rigid-body dynamics and contact interactions

Modules 3 and 4 deal with the creation of bodies and the definition of their inertial and surface properties. The latter depend on the contact method the user choose to use, among those available in 2. The following paragraphs recall briefly the main features of the code concerning body handling and contact interactions, focusing on the parameters relevant to each method. References are provided for a more complete discussion of each method.

Creation of bodies and their properties

The shape of each body can be selected by the user among common geometries (sphere, box, cone,…) or directly input by the user as a triangulated mesh. The shape can also be generated as the convex envelope of a randomly created cloud of points. When dealing with a large number of bodies, methods exist to select their spatial distribution or pattern (regular grid, randomly uniform, Poisson-disk sampling,…) and bounding domain (common geometry or fitted in a triangulated mesh). Each body possesses 6 degrees of freedom, including both translational and rotational 3D motion. Random distribution routines can be used to set physical properties (size, density,…) and dynamical state (position, velocity, angular position, angular velocity) of bodies as well. Finally, the surface properties of each body can be defined. These are directly correlated to the contact method in use. Surface/contact parameters are discussed below both for smooth and for non-smooth methods. More details can be found in the documentation of 

Tasora (2019).

Collision detection

Collision detection is performed into two steps: a broad and a narrow phase. During the broad phase, pairs of bodies whose geometrical boundaries are close enough are identified. The bounding volume of a body is estimated based on its velocity and integration time step. If the bounding volume of two bodies overlap, the narrow phase is performed and contact points are precisely found using a GJK algorithm (Tasora and Anitescu, 2010; Ferrari et al., 2017).

Contact method: smooth contacts

The Smooth Contact Method (SMC) implemented in C::E (Fleischmann et al., 2015)

is a force-based soft-body DEM method. The equations of motion are formulated as a system of Differential Algebraic Equations (DAE): namely a system of Ordinary Differential Equations (ODE) to reproduce the dynamics and Algebraic Equations (AE) for the kinematic constraints. Forces are exchanged between bodies at contact using a two-way normal-tangent spring-dashpot system. The constitutive model of the spring-dashpot system can be selected between the simple Hooke and the more realistic Hertzian model. Relevant parameters of the SMC are coefficients of friction (static, dynamic, spinning), cohesion (as an attractive force at contact points), stiffness and damping. As any soft-body DEM, it is suitable for smooth problems with no discontinuities but requires very small time steps for cases with high stiffness to avoid numerical instability.

Contact method: non-smooth contacts

The Non-Smooth Contact (NSC) method implemented in C::E (Tasora and Anitescu, 2010; Anitescu and Tasora, 2010; Tasora and Anitescu, 2011) is an impulse-based hard-body method. Unlike the SMC, the equations of motion are formulated as Differential Variational Inequalities (DVI) and require the solution of a Cone Complementarity Problem (CCP) at each time step. The NSC shares with the SMC parameters defining friction (static, dynamic, spinning) and adhesion. The handling of contact interaction is however much simpler and solely relies on the restitution coefficient, defined as the ratio between velocity after and before the collision. Due to its impulse-momentum formulation, the NSC is best suited for problems with discontinuities or with nearly rigid contacts (high stiffness).

Contact method: non-smooth contacts with compliance

In the context of non-smooth dynamics, a method is available to simulate contacts with compliance and damping (Tasora et al., 2013). As for the NSC method, the equations of motion are formulated as DVI and are based on a impulse-momentum formulation. However, compliance and damping are enforced at the constraint level and the method is suitable not only for non-smooth problems but also to simulate the elastic and smooth behavior typical of soft-body DEM models. Also, since it does not rely on the solution of a DAE, but its solution is found after solving a CCP-based DVI, this method does not require the very small time step required by soft-body DEM. Its parameters are analogous to that of SMC: friction (static, dynamic, spinning), cohesion, stiffness and damping.

2.3 Numerical solvers

As discussed in the previous section, when using the SMC approach the dynamics are written as a system of DAE, and when using the NSC approach the dynamics are written as a DVI problem. From the numerical point of view, DAE require to numerically integrate ODEs, while DVI require the solution of an optimization problem at each time step. Despite the radical difference in approaching the solution, both approaches require a time stepper and a nonlinear solver. A large variety of time steppers and solvers is available in C::E (Tasora, 2017): the most interesting for gravitational-granular dynamics applications include symplectic methods (semi-implicit Euler, leapfrog) and Runge-Kutta methods (RK45, explicit Euler, implicit Euler, trapezoidal, Heun), to be complemented with either iterative or direct solvers. In particular, symplectic methods are very suitable for gravitational dynamics problem and semi-implicit Euler method provide the best performance among them, when paired with an iterative solver (Anitescu and Tasora, 2010; Mangoni et al., 2018). This is suitable and effective to solve both DAE and DVI problems.

2.4 Post-processing: aggregate identification

Module 8 contains routines for the post-simulation analysis of results. We provide here one example that applies to the case of rubble-pile aggregation simulations. One of the tasks to be performed after such simulations is to identify the aggregate(s) of bodies and determine their properties such as shape, mass, inertia and void fraction (porosity).

Aggregate identification

The single or multiple aggregates obtained are identified through a method based on graph theory (Wilson, 1996). First, the graph is initialized to represent the whole system. Bodies are vertices of and the size of the vertex set is (number of bodies). Edges of represents physical contacts between objects: two bodies and , represented by vertices and of , are in contact if and only if there is an edge linking to . At the end of the simulation, a contact container is filled with pairs corresponding to pairs of bodies in contact with each other. The graph is filled with all existing contacts, drawn as edges in the graph. A schematic example is shown in Figure 2.

Figure 2: Example of global contact graph construction

Aggregates obtained correspond to the connected subgraphs of . Single bodies, that do not belong to any aggregate (e.g because they have escaped from the re-aggregation processes) are represented by isolated vertices. The number of aggregates is defined as the number of connected components of with a size greater than one, that is the number of the subset made of the connected subgraphs of with more than one vertex. A Depth First Search (DFS) algorithm (Wilson, 1996) is used to provide an effective way of traversing simple graphs. The output is a classification of the edges and a spanning tree that can be used for finding the connected components of a graph. Using such algorithm results in a complexity of (where and are respectively the number of vertices and edges of ) for the differentiation process.

Shape of the aggregate

The problem of defining the shape of a body made up of several discrete elements is not trivial and has a non-unique solution: infinite surfaces exist to envelope a given set of points. After identifying the aggregate and its children bodies, we address the problem of defining its overall shape by means of an alpha-shape algorithm (Edelsbrunner and Mücke, 1994). This computes the surface by enveloping a given set of points: in our case all vertices of children bodies. Unlike the simpler convex hull envelope, the alpha-shape surface is, in general, non-convex. The working principle of the alpha-shape algorithm can be visually represented as follows. We start from a cloud of points (vertices of children bodies). Let us imagine to roll a sphere of radius over the cloud of points. The alpha-shape envelope is the surface created by the sphere as it rolls over the cloud of points. In the limiting case where , the surface enveloping the cloud would be a convex hull. For a finite value of , non-convex envelopes can be obtained. Typically, should be on the same order of magnitude of the mean distance between vertices. As mentioned, the problem has a non-unique solution and the choice of makes the process arbitrary. This is affecting in a relevant manner the computation of global properties of the aggregate such as volume, void fraction (porosity) and shape-related quantities (elongation, axis ratios, …). To solve for this arbitrariness, we define an admissible range of enveloping surfaces, within a range of limiting cases. In general, two limiting surfaces can be defined: (i) the minimum volume surface, which is found using the minimum value of such to have no surfaces delimiting cavities inside the aggregate, and (ii) the maximum volume surface, which is the convex hull of the aggregate, obtained with . Global properties are then computed considering the uncertainty on their values due to the arbitrary definition of the surface, and range between values computed using minimum and maximum volume surfaces.

3 Validation scenarios

Original modules of C::E dealing with the modeling of granular dynamics have been widely tested and validated against benchmark and experimental tests (see “Validation Studies section in Tasora (2019)Heyn (2013) and others). The results of further test simulations are shown here to validate the gravity modules and their integration into the overall implementation. In particular, the accuracy of direct N and octree gravity models are assessed and their performance is compared in terms of computational time. Conservation of energy, linear and angular momentum are checked for the purely gravitational case and for the coupled gravity-collision problem.

3.1 Direct N-body integration

The first set of tests concerns the accuracy of force computation for the direct N method. Unlike the octree method, direct integration considers all gravitational interactions with no physical approximation. We performed several simulations of direct integration between mutually interacting bodies and in particular:

  • conic solutions of the two-body problem are reproduced as function of orbital energy;

  • orbits of planets in the Solar System are propagated for two-thousand years;

  • Halo and Lyapunov orbits around collinear libration points of the Sun-Earth and Earth-Moon circular restricted three-body problems (CR3BP) are reproduced.

When compared to the expected analytical (two-body problem) or numerical solution (N-body problem and CR3BP), the code is able to reproduce exactly the result, where the term exactly (here and from now on) implies to within round-off accuracy or numerical integration truncation errors. As mentioned in section 2.3, the semi-implicit Euler method is a very effective choice to solve gravitational dynamics, thanks to its symplectic nature that guarantees conservation of energy within the system. Energy, linear and angular momentum are conserved exactly for all cases, consistently with the nature of the direct N formulation that introduces no approximation errors in the model. The successful outcome of above tests guarantee that gravitational forces are computed correctly for each body. Figure 3 shows an example related to one of the aforementioned scenarios. The expected in-plane wobbling of the Sun around the barycenter of the Solar System (data from JPL’s Horizons ephemerides, solid line) is compared and matches well with the outcome of a direct N integration of the motion of Sun and planets during a 100-year timespan (cross markers).

Figure 3: Validation scenario: in-plane (x and y components in the ecliptic J2000 plane) wobbling of the Sun around the barycenter of the Solar System during a 100-year timespan. Comparison between JPL’s Horizons ephemerides333https://ssd.jpl.nasa.gov/(solid line) and numerical integration (cross markers)

3.2 GPU-parallel octree

We evaluate here the performance of the GPU-parallel implementation and compare it to direct N integration in terms of computational time and accuracy. Results presented here were obtained using a medium-range laptop with

  • Intel(R) Core(TM) i7-6500U CPU @3.1Ghz

  • Nvidia Geforce 940M (384 CUDA cores).

Although certainly dependent on hardware architecture, the results here shown identify trends and behaviors that are not dependent on hardware, but only on the algorithms used. To assess the performance of algorithms, two different sets of simulations are performed. For each simulation, bodies are created at time zero in random positions within a cubic domain and with a uniform distribution. Collision detection and interaction is disabled during the entire duration of the simulation, such that bodies only experience mutual gravity. The results of validation scenarios are in agreement with previous assessments concerning tree codes 

(Barnes and Hut, 1986; Hernquist, 1987).

Computational time

A first set of simulations is arranged to compare the computational time required by each algorithm. Simulations are performed for several scenarios, including ranging from 10 up to 10, and using both direct N integration and Barnes-Hut GPU algorithm with different level of accuracy (). We compare the time required by the code to compute one integration time step. In principle, results can be compared after computing only one time step. However, to avoid any bias and for the sake of a more accurate profiling, we integrate the dynamics forward for several time steps and take the average time required to compute one of them. Also, we want the dynamics to be steady such to not interfere with the evaluation of performance. To this goal we use a very short time step (10 s) and a total duration of the simulation of 10 s. Results are shown in Figure 4 using a logarithmic scale. As expected, the cost of direct integration increases with N. As discussed in section 2.1, the BH-GPU algorithm executes a sequence of kernels and the largest amount of computational time is shared between the creation of the octree and the evaluation of the gravitational body-body or body-cluster interactions. Both octree and gravity evaluation depend on the number of bodies in the system, but to very different extents. Figure 4 shows clearly that for a low number of bodies, the computational time required is dominated by the octree creation phase, which is independent from the choice of the accuracy parameter and increases very slowly with N. Conversely, for high number of bodies the time expense is dominated by the evaluation of the gravitational interactions, which instead depends on . When comparing direct N with BH-GPU, it clearly appears that for a low number of bodies (approximately N<700, the same for which is not affecting the time expense of BH-GPU), the process of building the octree makes the BH-GPU algorithm much more expensive than direct N integration. Conversely, for a high number of bodies, the BH-GPU is much faster than direct N. For N=1,000 the BH-GPU turns to be nearly two times faster and for N=10 it is faster of more than one order of magnitude. For a further increased number of bodies, direct N simulations become practically unfeasible and BH-GPU becomes the only feasible option. As discussed, BH-GPU algorithm can be tuned in accuracy: this affects the computational time as shown for different values used in Figure 4. In this high-N regime, the computational time of BH-GPU ranges from (with , i.e. no clustering of bodies, equivalent to direct integration) up to as increases.

Figure 4: Time to compute one integration time step as function of number of bodies in the simulation. Results are shown for direct N integration vs Barnes-Hut GPU-parallel algorithm ()
Figure 5: Maximum relative position error in time. The error is computed for simulations of 1,000 bodies with Barnes-Hut GPU-parallel algorithm (), and it is relative to a simulation with direct N integration


A second set of simulations is used to assess the accuracy of the BH-GPU algorithm. In this case, we wait for the bodies to complete a full orbit around the barycenter of the system and come back to their initial state at rest. Since collisions are disabled, all bodies are undisturbed in their gravity-driven path and the system exhibits a nearly perfect energy-conserved pulsating contracting-dilating behavior. This simulation time frame gives us a commensurate estimate of typical time duration of aggregation phenomena and provides us with every possible gravity-related condition, i.e. both, when bodies are far away and when they are close to each other. Compared to the first set of simulations, we use a larger time step and we integrate forward the dynamics for a much longer time. Results are here shown for the case of N=1,000 and for . A simulation using direct N integration is also performed and used as a real-world result to assess the accuracy of the BH-GPU algorithm. At each time step, the position of each body in the system is compared to its corresponding real-world one and the relative error is evaluated as percentage of the latter


where and are the positions of body at time when using, respectively, BH-GPU and N algorithms. Figure 5 shows the maximum position error among bodies in the system at each time step of the simulation. As expected, the error increases monotonically with , which is a direct mean and user-tunable parameter to select the desired level of accuracy. Few more features are worth of discussion. The error is shown to be time-dependent and exhibit a peak at half simulation time, i.e. at the end of contraction phase. This is expected, as the relative position error depends on the magnitude of accelerations acting on the bodies. The higher the acceleration, the higher the relative error when evaluating body-body vs body-cluster interactions. Accordingly, the error is shown to be higher in the middle part of the simulation, where bodies are closer to each other and it increases during contraction and decreases as they get further to each other. It is interesting to look at the worst case depicted in Figure 5. At the end of contraction phase, the error between the two extremes and is of about two orders of magnitude. However, the maximum relative error obtained using is still quite small, in the order of 0.3%. Also, it is worth noting that this simulation provides us with an extremely conservative estimate of the accuracy. This is because real-world scenarios include the dissipative effects introduced by collisions and often entail only a part of the contraction phase simulated here. For example, if collisions were enabled in the scenario considered here, they would become very relevant starting at approximately 2.510 s (transient phase) and the bodies would never get as close to each other as they are in the [3.5, 6.5]10 s range due to their physical extent (aggregate phase). They will therefore experience lower gravitational acceleration actions and position errors due to BH-GPU clustering: surely lower than the worst-case 0.3% value found here.

3.3 Collisions

A set of simulations is performed to test conservation laws during collisions and contact interactions between bodies. In this case, we simulated the settling dynamics of an already formed self-gravitating aggregate with 1,000 convex-hull bodies. The dynamics are propagated forward for one time step and contacts are solved using both smooth and non-smooth methods. The difference of physical and dynamical quantities are monitored. Several simulations are performed using the SMC method and different time steps, ranging between 10 and 10 seconds. The same is done using the NSC method which is able to handle higher time steps, between 10 to 10 seconds. A more detailed discussion on time step selection criteria is reported in Section 4.1. Unlike the purely gravitational case, energy, angular and linear momentum are not conserved exactly and collisions introduce approximation errors in the models. Since the semi-implicit Euler method is a first order integrator, the error is expected to depend linearly on the time step. This is confirmed by our simulations, as shown in Figure 6, which refer to the SMC case. It shows the angular momentum error in percentage, normalized to its initial value, for different time steps. As expected, a clear linear trend is identified for the error, as it scales with the time step.

Figure 6: Error on total angular momentum of the system after one simulation time step, with both gravity and contact interactions. Error is shown in percentage, normalized to the initial value of angular momentum. The linear trend of the error as function of time step duration is highlighted.

4 The asteroid aggregation problem

This section discusses how rubble-pile scenarios are simulated using the code described above (section 2). We discuss the appropriate choice of methods to use, how to tune their parameters and how to set up the numerical simulation.

4.1 Addressing the problem

To properly address the problem of simulating rubble-pile dynamics, we study the fundamental interactions that occur between bodies during gravitational aggregation phenomena. These can be associated to phases of the gravitational aggregation process:

  • gravitational phase: the dynamics of the bodies is only driven by mutual gravity and they do not experience contact interactions.

  • transient phase: this phase begins as bodies start to collide with each other and form small aggregates of few bodies each. This phase is characterized by an extremely chaotic motion of the bodies that are simultaneously subjected to both gravity and random collisions. From a macroscopic point of view, the overall dynamics of the system is dominated by mutual gravity between bodies and newly forming aggregates. However, at particle level, the dynamics is dominated by collisions, which creates impulsive (non-smooth contact) or quasi-impulsive (smooth) forces on a much faster scale than gravitational one. The combined effects at system and body level, and the dissipative nature of collisions leads eventually the system to a steady state condition.

  • aggregate phase: after transient phase, the system reaches a steady state equilibrium where bodies are either clustered in one or more aggregates, or dispersed.

The most challenging phase to be simulated is certainly the transient phase, where both gravitational and collision dynamics play a major role. Since gravity and collisions act at different scales, we explore the capability and suitability of the code by studying such dynamics separately. This leads to meaningful and general conclusions on how to numerically simulate such interactions, also applicable to the fully coupled problem.

Gravitational dynamics

In this early stage of the aggregation process, bodies are not in contact with each other and only interact gravitationally. As discussed in Section 3, the choice of the most suitable method depends on the number of bodies involved: direct N-body for few bodies, octree for a large number of bodies with a hardware-dependent threshold typically between 1,000 and 10,000 bodies. In order to correctly simulate the N-body dynamics, the integration time step must be consistent to the characteristic time of the problem and must sample the dynamics at least twice per . In particular, the following must be satisfied (Ferrari et al., 2017):


where is the universal gravitational constant and is the material density of the bodies. For typical values of asteroid material density ranging between 1 and 4 g/cm (Richardson et al., 2002a), maximum time steps are in the order of thousands of seconds.

Collisional dynamics

As done for gravitational dynamics, we study here the fundamental features related to collisional interactions between bodies. We analyze the different phases and numerical tasks involved when solving for collisions and derive numerical constraints. From the numerical and algorithmic point of view, collisional dynamics requires the solution of two different problems: collision detection and collision output determination.

The collision detection algorithm operates at each time step to identify contacts between body pairs. To ensure proper collision detection, the time step must be adequate: if the time step is too high, a collision could be missed. To quantify the dynamical constraint on the time step, we consider here the case of a collision between two spheres of radius . In particular, we consider the limiting case when a collision is not detected. It is easy to assess that, for the case of a direct collision occurring along the line connecting the centers of the spheres, the distance traveled during the time step must be lower than 4, and then


where is the relative velocity between bodies. Eq. (9) can be generalized for the case of grazing collisions (Sánchez and Scheeres, 2011):


where is the maximum overlap allowed between bodies: the smaller the overlap allowed, the smaller the time step required and the more precise the detection algorithm. Note that the case of direct collision in (9) can be retrieved from (10) when (complete overlapping). When considering a system of bodies, is the maximum relative velocity in the system and is the sum of the radii of the two smallest particles. Also, although (10) refers to a simplified case with spheres, the same relation can be conservatively used with convex hulls of any shape, where represents their minimum characteristic size.

The second task is to solve for the actual collisional dynamics when bodies are in contact, in order to compute the collision output. This task is performed by the contact method selected, as discussed in section 2.2. Hard- and soft-body methods deal with a very different modeling of the physics at contact. Hard-body methods model instantaneous collisions and deal with non-smooth dynamics based on impulse-momentum formulation. Since the collision is instantaneous, it does not make sense to speak of characteristic time involved in the dynamical process and then to derive a requirement on the time step from that. However, the time step still plays a major role and has to be selected properly to ensure the stability of the numerical solver. Hence, unlike gravity, collision detection and, as discussed below, soft-body methods, the choice of the time steps for hard-body methods is not driven by dynamical constraints, but by the properties of the numerical solver. We consider here the case of the semi-implicit Euler method. Looking at its stability region (Niiranen, 1999), we can derive a constraint on the time step to be used, as function of the roots of the characteristic equation. However, our problem involves many bodies and requires the solution of a CCP at each time step: the roots of the characteristic equation of such problem are not of easy determination. Alternatively, an empirical but very effective method is to determine a proper time step by comparison against known benchmark problems or known behavior of the system. For example, when benchmark problems are not available, a very convenient way is to tune the time step according to a desired level of accuracy in terms of conservation of total angular momentum of the system . As discussed in Section 3, the error depends linearly on the time step. This linear dependency provides us with a very effective means to estimate and select properly the time step. Further examples on how the time step can be set are provided in C::E technical documentation and validation studies (Tasora, 2019), with reference to non-smooth (NSC) and non-smooth visco-elastic methods.

On the other hand, unlike hard-body ones, soft-body methods directly rely on the physical modeling of collisional dynamics. The collision occurs in a finite time and foresee the visco-elastic behavior of the bodies at contact. Collisions are modeled by means of a spring-dashpot system at the contact point. We can therefore identify a characteristic time related to the fundamental frequency of the spring-dashpot system  (Tsuji et al., 1993):


When considering a system of bodies, is the highest frequency in the system, i.e. made of the combination of highest stiffness and lowest reduced mass 444 is the reduced mass of the two least massive bodies in the system. Accordingly, the time step must be at least two times smaller than the characteristic time


Practical applications typically make use of smaller time steps, ranging from to  (Sánchez and Scheeres, 2011), or even smaller than  (Herrmann and Luding, 1998).

To summarize, two different constraints on the simulation time step can be derived after the analysis of the collisional process. The first one is expressed by Eq. (10). It concerns the capability of the code to properly detect collisions and does not depend on the contact method in use. The second one depends on the contact method in use. In case of impulsive contacts, it depends on the stability of the numerical solver and can be determined empirically. In case of smooth dynamics, it depends on the dynamics of contact as expressed in Eq. (12). Both requirements shall be satisfied to properly simulate collisional processes. This duality is easily resolved by enforcing the most stringent requirement, i.e. by choosing the lowest time step required. As a general rule, collision detection will be the bottleneck for high velocity impacts and vice versa, an accurate collision output will be the driving criteria in case of low velocity collisions. In the attempt to quantify in a more rigorous way this general criterion, we investigate possible links between Eq. (10) and (12). Indeed, while the two effects appear not to be directly coupled for hard-body methods, some considerations can be made for the case of soft-body methods, where both (10) and (12) are built upon the dynamics of the problem. In particular, it is reasonable to assume that the maximum overlap allowed during collision detection must be consistently smaller than the maximum contraction of the spring-dashpot system. If this is not satisfied, the visco-elastic behavior at contact will not be accurately resolved. The maximum contraction of a spring-mass system can be computed based on energy conservation consideration, by equating the total energy right before contraction (only kinetic energy, no contraction) to the total energy at maximum contraction (only potential energy, zero velocity):


from which is easy to find that


In the case of a spring-dashpot system, the total energy is no longer conserved and a fraction of the initial kinetic energy is dissipated. The maximum contraction will be lower in this case, reduced by a factor <1


Now we define the maximum overlap allowed during collision detection as a small fraction of :


with , or equivalently


with . We can now substitute (17) into (10) and get


This relation correlates the two time steps involved in the numerical simulation of collision processes. In particular, we compare (18) with (12). The collision detection time step and the stiffness time step are equal if


which can be rewritten as


and finally, after trivial algebraic manipulation, as


This provides us with the quantification of the general criterion stated above. In particular, when selecting the time step, collision detection will be the driving criterion when , i.e. when


and vice versa. Eq. (22) gives us a useful condition to evaluate the tightest constraint through a simple comparison between maximum velocity and highest frequency in the system. Also, it correlates them with the accuracy of collision detection through . As expected, the problem of missing grazing collisions is an issue for high velocities and vice versa. More in detail, the break even point of the criterion is driven by the behavior , where is very small and is typically very high. Note that the dependency on shown in (22) might be misleading, since appears at the denominator of as well. In particular, the mass is proportional to and then the right hand side of Eq. (22) is proportional to .

The last point to discuss is about the choice of the contact method (SMC vs NSC vs NSC with compliance) to be used. As mentioned in section 2.2, SMC are naturally suited for smooth problems (low velocity contacts), NSC for non-smooth problems (high-velocity impacts) while in principle, the NSC method with compliance can be used for both. A more detailed comparison between the aforementioned methods, including simulation of different scenarios, evaluation of performance and criteria for selection of parameters is reported in C::E technical documentation (Tasora, 2019). Studies to assess the accuracy of such methods are also reported. Unlike gravity, where the exact dynamics of the system are known and can be computed, contact methods do not model in a comprehensive fashion the physics involved in collisional processes. While the results of gravity algorithms can be validated using an analytical exact formulation, the only way to assess the accuracy of contact methods is through experimental test. As mentioned in Section 3, the documentation of C::E (Tasora, 2019) reports the results of several benchmark scenarios, including cone penetration test, standard triaxial test and direct shear test. Further analytical and experimental validation scenarios are reported in Heyn (2013).

The gravity-collision problem

The full gravitational aggregation problem features both gravitational and collisional dynamics and, in principle, both of them constraint the choice of the simulation time step. However, gravity and collisions act on very different time scales and the characteristic time of gravity is typically several orders of magnitude higher than that of collisions. Also, their dynamics are not coupled. For these reasons, it is not required to perform the simulation using only one time step (which would be the smallest one, i.e. the one related to collisions ). Instead, as done in Sánchez and Scheeres (2012), the aggregation scenario can be simulated using two time steps: the gravity time step , determined using (8) and the collision time step , determined using either (10) or (12). In particular, we advance the simulation of to properly solve collisional dynamics, but we evaluate gravity (through interactions or by building the octree) every . This allows for consistent savings in terms of computational time, without jeopardizing the accuracy of results.

In this context, a more accurate modeling of the local gravity interactions between fragments would not be beneficial towards a more realistic representation of the granular N-body problem. Past astrodynamics studies have clearly shown that the mutual gravitational dynamics between objects of irregular shape, represented using shape-based (e.g. polyhedron) models, are very different from the case of point-mass sources interaction (e.g. Fahnestock and Scheeres, 2006; Ferrari and Lavagna, 2018). However, the granular N-body problem is very different from that of astrodynamics applications, which consider purely gravitational motion with no contact and collision interactions. In our case, due to the presence of both gravity and contact interactions, acting at different time scales, a noticeable effect of deviations due to enhanced gravity modeling would appear after a time when many contact/collision interactions have already occurred. In this case, the chaotic nature of the contact/collision interactions represents a much higher source of uncertainty compared to the accuracy of local gravity models. Because of this, all authors identify the better modeling of contact/collision interactions as the way to enhance the realism of simulations (e.g. Michel et al., 2004; Richardson et al., 2009; Movshovitz et al., 2012, and others) and agree to using octrees, at the cost of a lower accuracy of local gravity computations. Unlike astrodynamics applications, the goal here is not to accurately reproduce the coupled and orbital-attitude dynamics of each fragments, but rather to reproduce the global behavior of the system. Relevant parameters concerning gravity are indeed global ones, such as bulk density and total mass of the system, and lower-fidelity local models of gravity (as in octrees) are widely used.

4.2 Features of the problem

As discussed, the study of the dynamics of self-gravitating aggregates is a very complex problem and deals with challenges arising from the coupled and chaotic interactions occurring within the granular media. The study of asteroid shapes and rubble-pile dynamics has its fundamentals in the continuum theory, extended by Holsapple (2001, 2004, 2007, 2010) to the study of rubble-pile objects. However, self-gravitating granular systems are not strengthless fluid bodies and, for example, can spin faster than a perfect fluid before shedding mass (Richardson et al., 2005). Even for cohesionless aggregates, many effects can contribute to providing strength to the aggregate. These include surface interaction phenomena, such as friction, and other effects related to the geometry of contact interactions between particles such as interlocking (Sánchez and Scheeres, 2012). To summarize all these effects, the strength of a rubble pile can be quantified based on its effective angle of friction, which results from particle shapes, size distributions, packing and surface friction (Sánchez and Scheeres, 2012). Holsapple (2010) showed that the angle of friction is very important to establish the domain of admissible equilibrium shapes and that spins at limiting conditions become larger for a higher angle of friction, implying a larger strength of the aggregate. Accordingly, Sánchez and Scheeres (2016) observe that the angle of friction inhibits deformation. In Zhang et al. (2017) geometric interlocking is the main source of shear strength for crystal structures and higher values can be achieved for random packing and increased size heterogeneity of fragments. Among geometrical effects, extreme relevance is given by the angularity of non-spherical particles. As mentioned, spheres behaves very differently from angular particles, for which the contact dynamic problem is very different. Due to their non-smooth shape, angular particles can foster interlocking and between particles. In particular, spheres have only one contact point per pair, while non-spherical objects can have several. This has a great impact on all aspects related to contact dynamics and to the exchange of interactions/forces at contact points. For example, this affects how the value of parameters such as surface friction and restitution coefficients act on a global scale, since they both act at contact points. This means that their effect can be greatly amplified depending on relative geometry. In practice friction and restitution act on a larger scale and contact interactions results more dissipative with respect to the case of spheres. This is confirmed by experiments showing that restitution appears to be much lower for irregular objects (Hartmann, 1978) due to the loss of center-of-mass kinetic energy to rotational energy post-collision (Korycansky and Asphaug, 2006). In practice, angular shapes can increase the strength (Richardson et al., 2005; Jiang et al., 2015; Zhang et al., 2018) of the aggregate and, through a lower restitution coefficient, increase the likelihood of the formation of a stable aggregate during a gravitational re-accumulation process (Walsh et al., 2008). The resolution of the model, in terms of number of bodies in the rubble pile aggregate, is also very relevant. Richardson et al. (2005) show that coarse configurations consisting of a small number of larger particles are more resistant to tidal disruption or reshaping than fine configurations with many smaller particles, observing that the strength of the aggregate depends on the number of particles as well. A major role in the granular dynamics processes is played by friction. Sánchez and Scheeres (2012) found that friction increases substantially the limiting spin rate. For the case of angular bodies, where the overall effect of friction is higher (many contact points) the expected limiting spin rate would be even higher. Sánchez and Scheeres (2012) also observe that the addition of surface friction makes the aggregate more stable. They conclude that friction might modify the square root dependency between critical spin of a granular aggregate and bulk density. The same is found by Zhang et al. (2017), who conclude that higher interparticle friction can keep a spinning rubble pile stable at a higher spin rate. Furthermore, Pravec and Harris (2007) show how maximum spin rate increases (substantially) with friction. To summarize, the chaotic dependency between the dynamical history, initial shape/configuration of the particles, contact/surface parameters and properties of the final aggregate (Zhang et al., 2017) gives a clear picture of the complexity of the problem and the unpredictability of its dynamics.

4.3 Simulation examples

A set of simulations is performed here to extend results presented by Ferrari et al. (2017) to an increased number of bodies and to verify the correlations between parameters of the simulation, initial conditions of the N bodies and aggregation outcome. To simulate realistic aggregation scenarios, it is important to carefully select the physical properties of the N bodies and their initial dynamics. Initial conditions play a crucial role to the formation of the aggregates and their properties. The initial dynamical state includes the position and velocity of the center of mass of each body, as well as the angular position and spin rate of each body. Typical asteroid aggregation scenarios consider these bodies as fragments created from a parent larger body after, for example, a collision event (Michel et al., 2001; Campo Bagatin et al., 2018). In such case, the fragments could share a common residual angular momentum due to spinning motion of the parent body. Parameter is introduced to simulate such effect, and it is defined as the initial residual angular velocity of the system of particles: when the fragments have a nonzero orbital velocity around the center of mass of the system. In order to compare and extend the results in Ferrari et al. (2017), similar scenarios are set up, using non-smooth dynamics and a time step of 10 s. In our simulations, bodies are convex hulls whose shape is created as the geometrical envelope of randomly generated points. In particular, we use 16 points and corresponding convex hulls have, on average, 10 vertices. Surface friction is based on a simple Coulomb model with coefficient of 0.6. To foster the formation of aggregates, restitution is set to zero (perfectly inelastic collisions). Results presented here refer to simulations with 10 bodies, with average characteristic size of each fragment of 3 km and maximum initial distance between bodies in the order of 100 km. The physical properties of fragments are chosen among typical values of objects belonging to the main asteroid belt or Near Earth Asteroids population. The material density is set to 3 g/cm, a common choice for asteroid aggregation simulations (Sánchez and Scheeres, 2011), suitable to produce values of bulk density between 1.2 and 2.5 g/cm, typical of the asteroid population (Richardson et al., 2002a). At initial time, the angular and center of mass position of bodies are randomly generated with uniform spatial distribution, as well as the directions of their linear and angular velocity vectors.

Figure 7: Aggregation sequence example with 1,000 bodies

An example of aggregation sequence is shown in Figure 7, for a simplified case with 10 bodies (to facilitate the visualization of images). While not in contact, the motion of the bodies is driven solely by their mutual attraction. When fragments start to interact at closer range, few small bodies are scattered away due to collision mechanisms. Depending on initial conditions, one, multiple or no stable aggregates are eventually formed. Aggregation is not observed when relative velocities between bodies are too high and in particular when parameters , and are above certain threshold values, above which no aggregation is possible. As for the case under study, these values are identified approximately with m/s, rad/s and rad/s. After reaching a steady state, aggregates are considered as single asteroids. Different kind of aggregates have resulted from the simulation performed.

Figure 8 shows the properties of the largest aggregate formed after the numerical simulation, as function of the initial dynamics of the system. In particular it shows (a) the inertial elongation

, defined as the ratio between maximum and minimum moment of inertia of the aggregate 

(Ferrari et al., 2017), (b) the number of fragments in the aggregate and (c) its rotation period , as function of initial conditions (red diamonds), (blue asterisks) and (green stars). The values of are shown on the abscissa normalized to their maximum values, so that each parameter ranges from 0 to 1. Polynomial (for and ) and exponential (for ) fitting functions are used to highlight the trend of data distribution.

Figure 8: Properties of the final aggregate as functions of normalized initial conditions: (a) inertial elongation, (b) number of fragments in the main aggregate and (c) rotation period of the main aggregate. Uncertainties due to definition of final aggregate’s surface are highlighted.

The measure of the global properties and geometry of the final aggregate is affected by the uncertainty on the computation of its enveloping surface. As discussed in section 2.4, quantities such as volume, porosity, bulk density and size of the aggregate are better defined within a range of values computed between a minimum and a maximum volume surface. As a consequence, the inertial elongation and the angular velocity (or equivalently the rotation period ) of the aggregate are affected by the uncertainties of this measurement and should also be defined in within a range of values. In particular, the angular velocity of the aggregate is measured after computing its total angular momentum and inertia matrix , as , as in Richardson et al. (2005)

. The computation of the inertia tensor of the aggregate (and its principal inertia axes) rely on the definition of the enveloping surface. The uncertainty of this measurement depends on the size of the single particle with respect to the total aggregate. Hence, the accuracy depends on the number of particles in the aggregate and it is higher for smaller aggregates (those that spin faster in our simulations).

Richardson et al. (2005) estimate that the error in axis measurement could be as large as on particle radius along each axis, with errors of about 10% for the case of 1,000 bodies. When propagated to the computation of angular velocity, they estimate errors larger than 20%. Similarly, we estimated the error in the computed values of and , based on the number of bodies in the final aggregate, and show their uncertainties around the computed mean value in Figures 8 and 9. In our case the smaller aggregates are the fastest spinning ones: they are made of tens or few hundreds bodies, meaning potentially very large measurement errors. In other words, the resolution of our smaller aggregates is too coarse and does not allow to estimate accurately their spin rate, elongation, porosity and bulk density. Also, due to their low number of bodies, it is rather inaccurate to use the term rubble pile when referring to such small aggregates. To a closer inspection, the structure of these aggregates is closer to a monolith rather than a rubble pile: they are made of one (or few) big boulders and tens of smaller particles resting on their surface. This structure is by no means a typical self-gravitating aggregate and it’s rather closer to the “test particle on a rigid sphere” case reported in Weidenschilling (1981) and Richardson et al. (2000), where critical spin may be estimated using the simple balance between gravitational and centrifugal force (Hestroffer et al., 2019). This is indeed the result we find for them. The spin rate for aggregates with many more particles are comparable to what has been found by other authors within the error prescribed by the uncertainties mentioned above. Among these, our fastest spinning aggregates have spin rates in the order of 4.20.9 h.

General trends of Figures 8 and 9 agree with expected behavior and are coherent with what discussed in section 4.2. As expected, higher elongations are obtained for higher values of , and . Maximum elongation is reached for values between 75%-85% of their limiting value. For higher values of initial conditions (75%-85%) the simulations result in the formation of multiple aggregates of smaller size, which are more typically regularly shaped than larger ones. It is indeed observed that maxima in Figure 8(a) represent thresholds after which the formation of more than one significant aggregate occurs. As confirmed in Figure 8(b), the higher the velocities, the lower number of bodies is found in the main aggregate, which is smaller for higher values of , and . Figure 8(c) confirms a well-known result, with larger aggregates having a higher rotation period (lower spin rate) with respect to smaller ones. Figure 9 shows the rotation period of the final aggregate as function of (a) inertial elongation and (b) mass of the aggregate. As expected, fast rotators are more elongated and less massive, while slow rotators have a more regular rounded shape and are more massive.

Figure 9: Rotation period of the final aggregate as functions of (a) inertial elongation, (b) mass of the aggregate. Uncertainties due to definition of final aggregate’s surface are highlighted.

5 Conclusion

The paper presents the numerical implementation of a code for coupled gravity-granular dynamics problems and reports numerical methods available to solve gravitational and collisional dynamics. We build upon the work by Ferrari et al. (2017) and extend the capabilities of the code by introducing a parallel CUDA-GPU octree structure, to be able to evaluate mutual gravity for a higher number of bodies. Unlike existing N-body codes, the implementation presented here has the capability to handle non-spherically shaped bodies to a full double precision accuracy. Compared to classical sphere-based codes, this allows for a more realistic simulation of contacts between bodies and open to new opportunities and scenarios to be simulated. After presenting the results of test scenarios to validate newly implemented modules, we discuss the features and performance of numerical methods and derive requirements arising from the dynamics of asteroid aggregation scenarios. The problem of numerical simulation of gravitational aggregation is addressed: we provide guidelines and quantify criteria to properly set up such numerical simulations. We show examples of possible applications by performing numerical simulations of gravitational aggregation. The parameter space explored includes initial conditions of bodies in terms of relative velocity between them, spinning rate and residual angular velocity of the system. The outcome of the analysis is discussed showing the properties of the final self-gravitating aggregate after its stable formation in terms of inertial elongation, period of rotation and number of bodies in the aggregate.


This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 800060. Part of the research work was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. The authors would like to thank the anonymous reviewers for their comments and suggestions that helped to increase the quality of the paper.


  • S.J. Aarseth (2001) NBODY2: a direct n-body integration code. New Astronomy 6, pp. 277–291. Cited by: §1.
  • S. J. Aarseth (1999) From NBODY1 to NBODY6: the growth of an industry. 111 (765), pp. 1333–1346. External Links: Document, Link Cited by: §1.
  • M. Anitescu and A. Tasora (2010) An iterative approach for cone complementarity problems for nonsmooth dynamics. 47 (2), pp. 207–235. Cited by: §2.2, §2.3.
  • J. Barnes and P. Hut (1986) A hierarchical o(n log n) force-calculation algorithm. Nature 324, pp. 446–449. External Links: Link Cited by: §2.1, §3.2.
  • M. Burtscher and K. Pingali (2011) GPU computing gems. pp. 75 – 92. Cited by: §2.1.
  • A. Campo Bagatin, R. A. Alemañ, P. G. Benavidez, and D. C. Richardson (2018) Internal structure of asteroid gravitational aggregates. 302, pp. 343 – 359. External Links: ISSN 0019-1035, Document, Link Cited by: §4.3.
  • J.E. Chambers (1999) A hybrid symplectic integrator that permits close encounters between massive bodies. Monthly Notices of the Royal Astronomical Society 304, pp. 793–799. Cited by: §1.
  • P. A. Cundall and O. D. L. Strack (1979) A discrete numerical model for granular assemblies. 29 (1), pp. 47–65. External Links: Document, https://doi.org/10.1680/geot.1979.29.1.47, Link Cited by: §1.
  • E. Dorband, M. Hemsendorf, and M. D. (2003) Systolic and hyper-systolic algorithms for the gravitational n-body problem, with an application to brownian motion. Journal of Computational Physics 185, pp. 484–511. Cited by: §1.
  • F. Dubois, V. Acary, and M. Jean (2018) The contact dynamics method: a nonsmooth story. Comptes Rendus Mécanique 346 (3), pp. 247 – 262. Note: The legacy of Jean-Jacques Moreau in mechanics / L’héritage de Jean-Jacques Moreau en mécanique External Links: ISSN 1631-0721, Document, Link Cited by: §1.
  • M.J. Duncan, H.F. Levison, and M.H. Lee (1998) A multiple time step symplectic algorithm for integrating close encounters. Astronomical Journal 116, pp. 2067–2077. Cited by: §1.
  • H. Edelsbrunner and E.P. Mücke (1994) Three-dimensional alpha shapes. Transactions on Graphics 13 (1), pp. 43–72. Cited by: §2.4.
  • E. G. Fahnestock and D. J. Scheeres (2006) Simulation of the full two rigid body problem using polyhedral mutual potential and potential derivatives approach. 96 (3), pp. 317–339. External Links: ISSN 1572-9478, Document, Link Cited by: §4.1.
  • F. Ferrari, A. Tasora, P. Masarati, and M. Lavagna (2017) N-body gravitational and contact dynamics for asteroid aggregation. Multibody System Dynamics 39 (1), pp. 3–20. External Links: ISSN 1573-272X, Document, Link Cited by: §1, §1, §2.2, §4.1, §4.3, §4.3, §5.
  • F. Ferrari and M. Lavagna (2018) Ballistic landing design on binary asteroids: the aim case study. 62 (8), pp. 2245 – 2260. Note: Past, Present and Future of Small Body Science and Exploration External Links: ISSN 0273-1177, Document, Link Cited by: §4.1.
  • J. Fleischmann, R. Serban, D. Negrut, and P. Jayakumar (2015) On the importance of displacement history in soft-body contact models. 11 (4), pp. 044502–044502–5. External Links: ISSN 1555-1415, Document, Link Cited by: §2.2.
  • N. Gebhardt (2016) IrrLicht release 1.8.4. Note: http://irrlicht.sourceforge.netLast accessed 2016-12-01 Cited by: §2.
  • G. Gilardi and I. Sharf (2002) Literature survey of contact dynamics modelling. Mechanism and Machine Theory 37 (10), pp. 1213 – 1239. External Links: ISSN 0094-114X, Document, Link Cited by: §1.
  • W. K. Hartmann (1978) Planet formation: mechanism of early growth. 33 (1), pp. 50 – 61. External Links: ISSN 0019-1035, Document, Link Cited by: §4.2.
  • L. Hernquist (1987) Performance characteristics of tree codes. 64, pp. 715–734. External Links: Document Cited by: §3.2.
  • H.J. Herrmann and S. Luding (1998) Modeling granular media on the computer. 10 (4), pp. 189–231. External Links: ISSN 1432-0959, Document, Link Cited by: §4.1.
  • D. Hestroffer, P. Sánchez, L. Staron, A. C. Bagatin, S. Eggl, W. Losert, N. Murdoch, E. Opsomer, F. Radjai, D. C. Richardson, M. Salazar, D. J. Scheeres, S. Schwartz, N. Taberlet, and H. Yano (2019) Small solar system bodies as granular media. 27 (1), pp. 6. External Links: ISSN 1432-0754, Document, Link Cited by: §4.3.
  • T. Heyn (2013) On the modeling, simulation, and visualization of many-body dynamics problems with friction and contact. Ph.D. dissertation, University of Wisconsin-Madison. Cited by: §3, §4.1.
  • H.W. Hockney and J.W. Eastwood (1988) Computer simulation using particles. CRC Press. Cited by: §2.1.
  • K.A. Holsapple (2001) Equilibrium configurations of solid cohesionless bodies. 154 (2), pp. 432 – 448. External Links: ISSN 0019-1035, Document, Link Cited by: §4.2.
  • K. A. Holsapple (2004) Equilibrium figures of spinning bodies with self-gravity. 172 (1), pp. 272 – 303. Note: Special Issue: Cassini-Huygens at Jupiter External Links: ISSN 0019-1035, Document, Link Cited by: §4.2.
  • K. A. Holsapple (2007) Spin limits of solar system bodies: from the small fast-rotators to 2003 el61. 187 (2), pp. 500 – 509. External Links: ISSN 0019-1035, Document, Link Cited by: §4.2.
  • K. A. Holsapple (2010) On yorp-induced spin deformations of asteroids. 205 (2), pp. 430 – 442. External Links: ISSN 0019-1035, Document, Link Cited by: §4.2.
  • P. Hut and J. Makino (2010) Moving stars around. Cited by: §2.1.
  • Alder,B. J. and Wainwright,T. E. (1959) Studies in molecular dynamics. i. general method. 31 (2), pp. 459–466. External Links: Document, https://doi.org/10.1063/1.1730376, Link Cited by: §1.
  • M. Jean and J. J. Moreau (1987) Dynamics in the presence of unilateral contacts and dry friction: a numerical approach. In Unilateral Problems in Structural Analysis — 2, G. Del Piero and F. Maceri (Eds.), Vienna, pp. 151–196. External Links: ISBN 978-3-7091-2967-8 Cited by: §1.
  • M. Jiang, Z. Shen, and J. Wang (2015) A novel three-dimensional contact model for granulates incorporating rolling and twisting resistances. 65, pp. 147 – 163. External Links: ISSN 0266-352X, Document, Link Cited by: §4.2.
  • D.G. Korycansky and E. Asphaug (2006) Low-speed impacts between rubble piles modeled as collections of polyhedra. 181 (2), pp. 605 – 617. External Links: ISSN 0019-1035, Document, Link Cited by: §4.2.
  • D. Mangoni, A. Tasora, and R. Garziera (2018) A primal–dual predictor–corrector interior point method for non-smooth contact dynamics. 330, pp. 351 – 367. External Links: ISSN 0045-7825, Document, Link Cited by: §2.3.
  • P. Michel, W. Benz, and D.C. Richardson (2004) Catastrophic disruption of asteroids and family formation: a review of numerical simulations including both fragmentation and gravitational reaccumulations. Planetary and Space Science 52, pp. 1109–1117. Cited by: §1, §4.1.
  • P. Michel, W. Benz, P. Tanga, and D.C. Richardson (2001) Collisions and gravitational reaccumulation: forming asteroid families and satellites. Science 294 (5547), pp. 1696–1700. Cited by: §4.3.
  • A. Morbidelli (2002) Modern integrations of solar system dynamics. Annual Reviews of Earth and Planetary Sciences 30 (1), pp. 89–112. Cited by: §1.
  • N. Movshovitz, E. Asphaug, and D. Korycansky (2012) Numerical modeling of the disruption of comet d/1993 f2 shoemaker-levy 9 representing the progenitor by a gravitationally bound assemblage of randomly shaped polyhedra. 759 (2), pp. 93. External Links: Link Cited by: §1, §4.1.
  • J. Niiranen (1999) Cited by: §4.1.
  • C. C. Porco, J. W. Weiss, D. C. Richardson, L. Dones, T. Quinn, and H. Throop (2008) Simulations of the dynamical and light-scattering behavior of saturn’s rings and the derivation of ring particle and disk properties. 136 (5), pp. 2172. External Links: Link Cited by: §1.
  • P. Pravec and A.W. Harris (2007) Binary asteroid population: 1. angular momentum content. 190 (1), pp. 250 – 259. External Links: ISSN 0019-1035, Document, Link Cited by: §4.2.
  • C.D. Pruett, J.W. Rudmin, and J.M. Lacy (2003) An adaptive n-body algorithm of optimal order. Journal of Computational Physics 187, pp. 298–317. Cited by: §1.
  • D. C. Richardson, P. Michel, K. J. Walsh, and K. W. Flynn (2009) Numerical simulations of asteroids modelled as gravitational aggregates. Planetary and Space Science 57, pp. 183–192. Cited by: §1, §4.1.
  • D.C. Richardson, Z.M. Leinhardt, H.J. Melosh, W.F. Bottke Jr., and E. Asphaug (2002a) Gravitational aggregates: evidence and evolution. In W.F. Bottke Jr., A. Cellino, P. Paolicchi, R.P Binzel (Ed) Asteroids III, Tucson (Ed.), Vol. 160, pp. 501–515. Cited by: §1, §4.1, §4.3.
  • D.C. Richardson, T. Quinn, J. Stadel, and G. Lake (2002b) Direct large-scale N-body simulations of planetesimal dynamics. Planetary and Space Science 30 (1), pp. 45 – 59. Cited by: §1.
  • D. C. Richardson, P. Elankumaran, and R. E. Sanderson (2005) Numerical experiments with rubble piles: equilibrium shapes and spins. 173 (2), pp. 349 – 361. External Links: ISSN 0019-1035, Document, Link Cited by: §4.2, §4.3.
  • D. C. Richardson, T. Quinn, J. Stadel, and G. Lake (2000) Direct large-scale n-body simulations of planetesimal dynamics. 143 (1), pp. 45 – 59. External Links: ISSN 0019-1035, Document, Link Cited by: §4.3.
  • S. A. Rodionov and N. Ya. Sotnikova (2005) Optimal choice of the softening length and time step in n-body simulations. Astronomy Reports 49 (6), pp. 470–476. External Links: ISSN 1562-6881, Document, Link Cited by: 2nd item.
  • D. P. Sánchez and D. J. Scheeres (2012) DEM simulation of rotation-induced reshaping and disruption of rubble-pile asteroids. 218 (2), pp. 876 – 894. External Links: ISSN 0019-1035, Document, Link Cited by: §4.1, §4.2.
  • P. Sánchez and D. J. Scheeres (2011) Simulating asteroid rubble piles with a self-gravitating soft-sphere distinct element method model. 727 (2), pp. 120. External Links: Link Cited by: §1, §4.1, §4.1, §4.3.
  • P. Sánchez and D. J. Scheeres (2016) Disruption patterns of rotating self-gravitating aggregates: a survey on angle of friction and tensile strength. 271, pp. 453 – 471. External Links: ISSN 0019-1035, Document, Link Cited by: §4.2.
  • S. R. Schwartz, D. C. Richardson, and P. Michel (2012) An implementation of the soft-sphere discrete element method in a high-performance parallel gravity tree-code. 14 (3), pp. 363–380. External Links: ISSN 1434-7636, Document, Link Cited by: §1.
  • J. Stadel (2001) Cosmological N-body simulations and their analysis. Ph.D. Thesis, University of Washington, Seattle. Cited by: §1.
  • G. Tancredi, A. Maciel, L. Heredia, P. Richeri, and S. Nesmachnow (2012) Granular physics in low-gravity environments using discrete element method. Monthly Notices of the Royal Astronomical Society 420 (4), pp. 3368–3380. External Links: Document, /oup/backfile/content_public/journal/mnras/420/4/10.1111_j.1365-2966.2011.20259.x/2/mnras0420-3368.pdf, Link Cited by: §1.
  • A. Tasora and M. Anitescu (2011) A matrix-free cone complementarity approach for solving large-scale, nonsmooth, rigid body dynamics. 200, pp. 439–453. Cited by: §2.2.
  • A. Tasora, M. Anitescu, S. Negrini, and D. Negrut (2013) A compliant visco-plastic particle contact model based on differential variational inequalities. 53, pp. 2 – 12. Note: Multibody System Dynamics: A Selective Perspective External Links: ISSN 0020-7462, Document, Link Cited by: §2.2.
  • A. Tasora (2019) Chrono::engine project. Note: http://projectchrono.org/ Cited by: §1, §2.2, §3, §4.1, §4.1.
  • A. Tasora (2017) External Links: Link Cited by: §2.3.
  • A. Tasora and M. Anitescu (2010) A convex complementarity approach for simulating large granular flows. 5 (3), pp. 031004–031004–10. External Links: ISSN 1555-1415, Document, Link Cited by: §2.2, §2.2.
  • A. Tasora, R. Serban, H. Mazhar, A. Pazouki, D. Melanz, J. Fleischmann, M. Taylor, H. Sugiyama, and D. Negrut (2016)

    Chrono: an open source multi-physics dynamics engine

    In High Performance Computing in Science and Engineering, T. Kozubek, R. Blaheta, J. Šístek, M. Rozložník, and M. Čermák (Eds.), Cham, pp. 19–49. External Links: ISBN 978-3-319-40361-8 Cited by: §1.
  • Y. Tsuji, T. Kawaguchi, and T. Tanaka (1993) Discrete particle simulation of two-dimensional fluidized bed. 77 (1), pp. 79 – 87. External Links: ISSN 0032-5910, Document, Link Cited by: §4.1.
  • K. Wada, H. Senshu, and T. Matsui (2006) Numerical simulation of impact cratering on granular material. Icarus 180 (2), pp. 528 – 545. External Links: ISSN 0019-1035, Document, Link Cited by: §1.
  • K. J. Walsh, D. C. Richardson, and P. Michel (2008) Rotational breakup as the origin of small binary asteroids. 454, pp. 188. External Links: Link Cited by: §4.2.
  • L. Wang, R. Spurzem, S. Aarseth, K. Nitadori, P. Berczik, M. B. N. Kouwenhoven, and T. Naab (2015) nbody6++gpu: ready for the gravitational million-body problem. 450 (4), pp. 4070–4080. External Links: ISSN 0035-8711, Document, http://oup.prod.sis.lan/mnras/article-pdf/450/4/4070/5782448/stv817.pdf, Link Cited by: §1.
  • S.J. Weidenschilling (1981) How fast can an asteroid spin?. 46 (1), pp. 124 – 126. External Links: ISSN 0019-1035, Document, Link Cited by: §4.3.
  • R.J. Wilson (1996) Introduction to graph theory. Prentiss Hall. Cited by: §2.4, §2.4.
  • J. Wisdom and M. Holman (1991) Symplectic maps for the n-body problem. Astronomical Journal 102, pp. 1528–1538. Cited by: §1.
  • Y. Zhang, D. C. Richardson, O. S. Barnouin, C. Maurel, P. Michel, S. R. Schwartz, R. Ballouz, L. A.M. Benner, S. P. Naidu, and J. Li (2017) Creep stability of the proposed aida mission target 65803 didymos: i. discrete cohesionless granular physics model. 294, pp. 98 – 123. External Links: ISSN 0019-1035, Document, Link Cited by: §4.2.
  • Y. Zhang, D. C. Richardson, O. S. Barnouin, P. Michel, S. R. Schwartz, and R. Ballouz (2018) Rotational failure of rubble-pile bodies: influences of shear and cohesive strengths. 857 (1), pp. 15. External Links: Document, Link Cited by: §4.2.