SIMULATING ASTEROID RUBBLE PILES WITH A SELF-GRAVITATING SOFT-SPHERE DISTINCT ELEMENT METHOD MODEL

In a given system of coordinates, the motion of the particles of a multiparticle system will be governed by Newton's second law of motion:

$$\begin{equation} m_i {\ddot{\vec{\bf r}}_i}=\vec{\bf f}_i, \end{equation} \tag{ 1 }$$

where m_i is the mass of particle i and

$$\begin{equation} {\vec{\bf f}_i}=-\nabla _{r_i}\mathcal V. \end{equation} \tag{ 2 }$$

$\vec{\bf f}_i$ is the force on that particle or the net force applied to the center of mass of a multiparticle body ($\mathcal V$ is the potential energy). Energy will be conserved in the system when neither the potential energy, $\mathcal V$, nor the kinetic energy ($\mathcal K$) depend explicitly on time. The energy conservation law applies whether or not an external potential exists; the essential condition is that no explicit time or velocity-dependent forces act on the system. If there is a smooth variation of the interaction potentials, a finite-difference method can be applied to solve Equation (1) (SSDEM); in contrast, if interacting potentials do not act smoothly, individual collisions must be identified leading to the class of HSDEM (also called Event-Driven Dynamics). Our current code makes use of an SSDEM method, and for our simulations, we use a "Mid-Step Velocity Verlet" algorithm, a variation of the original velocity Verlet algorithm. The former is mathematically equivalent, but numerically superior to the latter (Swope et al. 1982). The equations for this algorithm are as follows (Allen & Tildesley 1987).

• 1.  
  The new linear and angular positions and mid-step velocities are calculated from the previous values:

  $$\begin{eqnarray} {\vec{\bf r}} (t+\delta t)=& {\vec{\bf r}}(t)+{\vec{\bf v}}(t)\delta t+\frac{1}{2}{\vec{\bf a}}(t)\delta t^2 \end{eqnarray} \tag{ 3 }$$

  $$\begin{eqnarray} {\vec{\bf \theta }} (t+\delta t)=&{\vec{\bf \theta }} (t)+{\vec{\bf \omega }}(t)\delta t+\frac{1}{2}{\vec{\bf \alpha }}(t)\delta t^2 \end{eqnarray} \tag{ 4 }$$

  $$\begin{eqnarray} {\vec{\bf v}} \left(t+\frac{\delta t}{2}\right)= &{\vec{\bf v}}(t)+\frac{1}{2}{\vec{\bf a}}(t)\delta t \end{eqnarray} \tag{ 5 }$$

  $$\begin{eqnarray} {\vec{\bf \omega }} \left(t+\frac{\delta t}{2}\right)= &{\vec{\bf \omega }}(t)+\frac{1}{2}{\vec{\bf \alpha }}(t)\delta t. \end{eqnarray} \tag{ 6 }$$
• 2.  
  The forces, torques, and accelerations in the system are calculated in accordance with the collision and friction models:

  $$\begin{eqnarray} &&{\vec{\bf a}} (t+\delta t)= {\vec{\bf a}}\left({\vec{\bf r}} (t+\delta t),{\vec{\bf v}} \left(t+\frac{\delta t}{2}\right)\right) \end{eqnarray} \tag{ 7 }$$

  $$\begin{eqnarray} &&{\vec{\bf \alpha }} (t+\delta t)= {\vec{\bf \alpha }}\left({\vec{\bf \theta }} (t+\delta t),{\vec{\bf \omega }} \left(t+\frac{\delta t}{2}\right)\right). \end{eqnarray} \tag{ 8 }$$
• 3.  
  The final calculation of the linear and angular velocities is performed to complete the iteration:

  $$\begin{eqnarray} {\vec{\bf v}} (t+\delta t)=&{\vec{\bf v}} \left(t+\frac{\delta t}{2}\right)+\frac{1}{2}{\vec{\bf a}}(t)\delta t \end{eqnarray} \tag{ 9 }$$

  $$\begin{eqnarray} {\vec{\bf \omega }} (t+\delta t)=&{\vec{\bf \omega }} \left(t+\frac{\delta t}{2}\right)+\frac{1}{2}{\vec{\bf \alpha }}(t)\delta t. \end{eqnarray} \tag{ 10 }$$

  After this, the calculation is just repeated as many times as required. At this point, the kinetic energy at time t + δt is available and the potential energy at this time will have been evaluated in the force loop. One particularity of this algorithm is that after every iteration the positions, velocities, and accelerations on every particle are known at the same moment in time. The choice of time step depends on the nature of the forces, their models within the simulation, and their parameters. Our choice will be explained later.

When two particles collide in our simulations, normal and tangential forces are calculated. The calculation of the normal forces between colliding particles is modeled by a linear spring dashpot. The elastic force is modeled as

$$\begin{equation} {\vec{\bf f}}_e=k_n\xi {\bf \hat{n}}, \end{equation} \tag{ 11 }$$

and the damping force as

$$\begin{equation} {\vec{\bf f}}_d=-\gamma \dot{\xi }{\bf \hat{n}}, \end{equation} \tag{ 12 }$$

Then the normal force is calculated as ${\vec{\bf f}}_n={\vec{\bf f}}_e+{\vec{\bf f}}_d$. In these equations, k_n is the elastic constant, ξ is the deformation of the particles, γ is the damping constant (related to the dashpot), $\dot{\xi }$ is the rate of deformation, ${\bf \hat{n}}$ is the vector joining the centers of the colliding particles (see Figure 1, top), and ${\bf \hat{t}}$ is obtained as the unit vector of the relative approaching velocity minus its normal component. This unit vector is (by construction) parallel to the circular plane contained by the intersection of the spherical surfaces of the two particles. To show this, let particles k and m, of radii R_k and R_m, be in a collision, then

$$\begin{equation} {\bf \hat{n}}=\frac{{\vec{\bf r}}_k-{\vec{\bf r}}_m}{|{\vec{\bf r}}_k-{\vec{\bf r}}_m|}. \end{equation} \tag{ 13 }$$

To calculate ${\bf \hat{t}}$, we first obtain the total relative velocity of the points of contact of both particles. These points, p_k and p_m, respectively, are defined by the intersection of ${\bf \hat{n}}$ with the surface of each particle. Let us call ${\vec{\bf p}_k}$ and ${\vec{\bf p}_m}$ to their position vectors with respect to the center of each respective particle. Then,

$$\begin{equation} \begin{array}{l@{\,}l@{\,}l}\vec{\bf p}_k&=&-R_k{\bf \hat{n}}\\ \\ {\vec{\bf p}_m}&=&R_m{\bf \hat{n}.} \end{array} \end{equation} \tag{ 14 }$$

Then, the total collision velocity of the point ${\vec{\bf p}_k}$ with respect to ${\vec{\bf p}_m}$ will be

$$\begin{equation} {\vec{\bf v}}_{km}={\vec{\bf v}}_k-{\vec{\bf v}}_m+{\vec{\bf \omega }_k} {\times }{\vec{\bf p}_k}-{\vec{\bf \omega }_m} {\times }{\vec{\bf p}_m}. \end{equation} \tag{ 15 }$$

Consequently, the relative tangential velocity and tangential vector will be, respectively,

$$\begin{equation} {\vec{\bf v}}_t={\vec{\bf v}}_{km}-{\vec{\bf v}}_{km}{\times }{\bf \hat{n}} \end{equation} \tag{ 16 }$$

$$\begin{equation} {\bf \hat{t}}=\frac{{\vec{\bf v}}_t}{|{\vec{\bf v}}_t|}. \end{equation} \tag{ 17 }$$

We can now define ξ and $\dot{\xi }$ as

$$\begin{equation} \begin{array}{l@{\,}l@{\,}l}\xi &=&(R_k+R_m)-|{\vec{\bf r}}_k-{\vec{\bf r}}_m|\\ \dot{\xi }&=&{\vec{\bf v}}_{km}{\times }{\bf \hat{n}.} \end{array} \end{equation} \tag{ 18 }$$

Zoom In Zoom Out Reset image size

The tangential component of the contact force models friction, static and dynamic (see Figure 1, bottom). This is done by placing a linear spring at the contact point, attached to both particles, at the beginning of the collision (Herrmann & Luding 1998; Silbert et al. 2001), producing a restoring frictional force ${\vec{\bf f}}_{t}$. Let D be the total tangential elongation of this spring during the collision, then while there is a static contact, it is defined as

$$\begin{equation} D=\int _0^t|{\vec{\bf v}}_t|dt. \end{equation} \tag{ 19 }$$

In the code and during the static contact,

$$\begin{equation} D=\sum {|\vec{\bf v}}_t|\delta t. \end{equation} \tag{ 20 }$$

The restoring static frictional force is

$$\begin{equation} {\vec{\bf f}}_{t}=-k_t D\hat{t}. \end{equation} \tag{ 21 }$$

The magnitude of the elongation of this tangential spring is truncated in order to satisfy the local Coulomb yield criterion $|{\vec{\bf f}}_t|\le \mu |{\vec{\bf f}}_n|$. Thus,

$$\begin{equation} {\vec{\bf f}}_{fr}=\rm min(-k_t|D\hat{t}|,-\mu |{\vec{\bf f}}_n|)\hat{\bf t.} \end{equation} \tag{ 22 }$$

This provides an effective stick-slip friction force between the colliding particles. The value of k_t is taken as (2/7)k_n as suggested by Silbert et al. (2001).

Different authors have used values for k_n ranging from 1000 Nm⁻¹ (Milburn et al. 2005; Malone & Xu 2008; Biswas et al. 2003) up to 18,000Nm⁻¹ (Richard et al. 2008) for similar materials, the difference being the type of system that was simulated. The former dealt with a system in which collision speeds were damped due to an interstitial fluid and the fluid was the driving mechanism behind the phenomenon that was observed. The later dealt with collisional and frictional flows where collision speeds were much higher and drove the system. In fact, the choice of parameters seems to relate to the phenomena being analyzed much more than to the nature of the materials being simulated if a linear spring–dashpot model is used. The works of Yuu (Yuu et al. 1995) and Malone (Malone & Xu 2008) point to the differences and advantages of different choices in parameters. When a Hertzian spring is used the values for the elastic constants are much higher and can be calculated from the real Young's moduli and Poison's ratios of the materials (Teufelsbauer et al. 2009).

The linear spring–dashpot model for two colliding spherical particles can be solved analytically (Herrmann & Luding 1998; Luding 1998), leading to the following expressions for the contact duration and the restitution coefficient:

$$\begin{equation} \begin{array}{l@{\,}l@{\,}l}t_c&=&\pi /\omega \\ r&=&\exp {(-\pi \eta /\omega),} \end{array} \end{equation} \tag{ 23 }$$

where $\omega =\sqrt{\vphantom{A^A}\smash{\hbox{$\omega _0^2-\eta ^2$}}}$, ω²₀ = k_n/m_km, η = γ/m_km, and m_km = m_km_m/(m_k + m_m).

It is generally accepted that the time step should be less than a certain critical value, which is usually expressed as a fraction of the natural frequency of the equivalent mass–spring system; i.e., when γ = 0. However, the exact fraction of the natural frequency and indeed the relation used to obtain the natural frequency can differ significantly between authors (see Malone & Xu 2008 for a review). It has since become common practice in DEM codes to use a relation of the form

$$\begin{equation} \delta t\le C\sqrt{m_{km}/k_n,} \end{equation} \tag{ 24 }$$

where m_km corresponds to the reduced mass of the two least massive particles and the highest k_n (when dealing with mixtures, different k_n values and densities are sometimes used (Biswas et al. 2003)). Values between 0.2π and 2 have been used by many authors for C. This relation is used and maintained in this paper to calculate adequate time steps for the simulations.

Just like the hard-sphere model, however, this model has limitations, and unless an infinitely small δt is used, grazing collisions with an infinitely small overlap or collisions with very high speeds could go undetected. Therefore, there is a trade between accuracy and performance. To show this, let us consider two spherical particles of radii R_m and R_k, respectively, whose centers are at a distance R_m + R_k and located in points c_m and c_k at time t (see Figure 2). In the worst case scenario, for this grazing collision to be missed, particle k has to travel from point c_k to c'_k in a time lesser or equal to δt. It can be shown that for this to happen,

$$\begin{equation} \delta t \ge \frac{2\sqrt{2\Delta \left(R_m+R_k\right)-\Delta ^2}}{v}, \end{equation} \tag{ 25 }$$

where v is magnitude of the relative velocity of particle k with respect to particle m, Δ = [(R_m + R_k) − s] so that 0 < Δ ⩽ R_m + R_k and s is the minimum distance between these two particles, had the trajectory of particle k been found exactly through its equations of motion and not through a discrete numerical integration. Taking as reference the parameters and mean overlap reported by Milburn et al. (2005), grazing collisions could have been missed for $\delta t\ge 10^{-3}\,\rm s$. In their work, the particle speeds were ${\approx }0.1\,\rm m\,s^{-1}$ and they used a time step of $5\times 10^{-6}\; \rm s$. In fact, grazing collisions even 100 times faster could have been detected.

Zoom In Zoom Out Reset image size

We can re-write the last expression as

$$\begin{equation} \Delta ^2 -2\Delta \left(R_m+R_k\right)+\left(\frac{v\delta t}{2}\right)^2\ge 0. \end{equation} \tag{ 26 }$$

The roots of this expression are

$$\begin{equation} \Delta =(R_m+R_k)\pm \sqrt{\left(R_m+R_k\right)^2-\left(\frac{v\delta t}{2}\right)^2} \end{equation} \tag{ 27 }$$

and the critical point is found at Δ = (R_m + R_k). Notice that vδt is the distance between c_k and c'_k. If the discriminant (s² by construction) in Equation (27) is less or equal to zero, the distance that particle k can cover during one time step is large enough to completely miss even a head-on collision with particle m. If the discriminant is greater than zero, it is possible to conclude that collisions will be missed for

$$\begin{equation} 0<\Delta \le (R_m+R_k)-\sqrt{\left(R_m+R_k\right)^2-\left(\frac{v\delta t}{2}\right)^2} \end{equation} \tag{ 28 }$$

In the case of Milburn et al. (2005), grazing collisions with an overlap of up to ${\approx }3.1\times 10^{-11}\,\rm m$ could have been missed. This is well within the accuracy of the method. From Equation (28), it is clear that v is the only variable that depends solely on the dynamics of the system, and should, therefore, serve as a guide to simulate it with the most appropriate method and parameters.

If the system to be studied has to handle particles with high velocities, algorithms with multiple time steps could be used (Tuckerman et al. 1990; Swindoll & Haile 1984; Teleman & Jonsson 1986). Still, if even with these methods the simulation is not practical or does not portray the system accurately, it might be recommended to use an HSDEM program instead of an SSDEM. However, that is beyond the scope of this paper.

In order to increase the speed of the program, a link-list method was implemented (Allen & Tildesley 1987). By means of this method, the space where the particles are confined is subdivided using a grid. Each cell (cell or cells for notation) in the grid contains only a small number of particles. In general, the size of the cells is chosen to be at least of the same length as twice the potential cutoff. The potential cutoff being the distance at which the interaction of any two particles is considered nonexistent or extremely small compared to other forces in the system. For contact forces, this potential cutoff is the sum of the radii of the two largest particles. The use of a spring–dashpot model, being a short range interaction, implies that each cell will be at least this size. One more point to consider is that within one time step two particles should not be able to reach the center of any cell. This is linked to the analysis made in the previous section referring to grazing collisions.

The number of the cell is a label for each particle and this allows the program to know if two particles are either in the same cell or in a neighboring cell. Only collisions between these particles are checked; this avoids the need to perform extra calculations when they are not close enough to interact.

The following is an explanation for a bi-dimensional system as the one shown in Figure 3. This algorithm first bins the particles in the different cells. Then it searches for collisions among the particles within each cell (the green one for this explanation) and forces are calculated when needed. In principle, collisions should be checked among the particles contained in all eight neighboring cells, but this would imply that for each pair of cells, collisions among their particles will be checked twice. To avoid this, only the four cells painted in light blue are checked; collisions among particles in the green cell and the orange cells will be checked when the later are the center of attention.

Zoom In Zoom Out Reset image size

For a three-dimensional (3D) simulation, a half-cube is considered, and the neighboring cells are the four shown in Figure 3 plus the nine immediately on them and out of the page, covering the green cell and the eight cells around it.

The initial conditions of the system refer to the choice of values of their geometrical and dynamical parameters and variables. The simulations presented here need to have values for the radii of the particles and a material density, from which the moment of inertia, volume, and mass will be later derived. Initial linear and angular positions and velocities must also be specified. The values for initial angular and linear accelerations will be derived from the initial values of force and torque. They all will, at some extent and in conjunction with the models implemented in the simulation code, define the kind of system that is simulated, i.e., if it is in one, two, or three dimensions, if the particles are spheres or disks, and if they themselves are solid or hollow. In our case, we have chosen to use solid spheres that move in three dimensions in a tri-dimensional space for all simulations.

The values for the radii of the particles are distributed randomly using a random number generator (Press et al. 1992) and allowed to vary by 20% from their mean value. Their initial positions are calculated as if they were in a 3D hexagonal close-packed (HCP) lattice, with a magnitude for the translation vector equal to 2.2 times the maximum of the radii of the particles. This gives the particles some room to move and avoids any initial overlapping that could render the simulation unstable from a numerical point of view. This adds a measure of randomness to the system and helps to avoid the artificial crystallization of the simulated grains. To further avoid crystallization due to the original HCP configuration, all the particles are given a random initially velocity and allowed to settle. The implementation of these two randomizing features made it possible to obtain filling fractions of ≈0.62, which is characteristic of a random packing of spheres.

For the simulations in the next section, a settling routine also had to be implemented. Usually, under 1g_e (g_e = 9.8 m s⁻², the average value for gravity on the surface of the Earth) a few seconds of simulation are needed to have the entire aggregate settled. However, given that we want to work under much lower and also higher levels of gravity, a settling routine cannot simply let them fall under their own weight until they rest motionless. If the simulations were to use a gravitational field higher than 1g_e, i.e., αg_e where α ⩾ 1, this routine would use this value to settle the particles and scale the value of the initial random velocities by $\sqrt{\alpha }$ to increase the kinetic energy accordingly.

If the value of gravity for the simulations was to be smaller than 1g_e, the particles would initially settle under 1g_e and initial velocities would not be scaled. When the settling under 1g_e has finished, the kinetic energy of the system is again calculated and the value of gravity is lowered by one order of magnitude. This unloading causes the particles to transform some of the energy stored in the springs into kinetic and potential energy. Then the settling process is started again with the just calculated value of kinetic energy. To avoid the excessive motion of the particles and to speed up the settling process, the particles are subjected to a Stokes-like drag. This process is repeated until the desired value of gravity is achieved and then the simulation starts.

During the settling process, the kinetic energy of the system is measured every 0.05 s. During the unload the aggregate behaves similarly to an overly damped spring. From the simulations it was seen that this time interval was sufficiently smaller than the oscillations of its kinetic energy.

The routine considers that the particles have settled when the variation in kinetic energy is less than 5 × 10⁻⁴ of its initial value. A smaller value did not change the results substantially, it only elongated the settling time.

As it was explained in the introduction, many of the observed asteroids are rubble piles; therefore, the understanding of the physical properties of these systems is essential. One of the most basic properties is the velocity of sound wave transmission (P-wave and S-wave velocities). Liu & Nagel (1992) have argued that the deformation of a bead under its own weight and thermal expansion could affect the transmission of sound; if this is so, an extremely low gravity environment could cause a change in this speed.

These simulations were carried out with 6000 and 24,000 glass-like particles with a density of 2500 kg m⁻³, 0.5 mm in mean diameter. In order to avoid any possible distortion of the results due to the interaction between the particles in the vertical walls of their container, they were placed in a rectangular box with periodic boundary conditions and solid bottom and top. The dimensions of the box were 6.25 × 50 × 6.25 mm for the simulations with 6000 particles and 12.5 × 50 × 12.5 mm for those with 24,000 particles; the height of the aggregate was approximately 17 mm in both cases. Measurements were made for gravity values of 1μg_e, 1mg_e 1g_e, 2g_e, and 3g_e. When the settling routine was finished, the first layer of particles at the bottom of the box was given an upward velocity of 0.5 m s⁻¹, while the velocities of the rest of the particles in the aggregate remained unchanged. This method approximated an initial plane wave, mimicking the experiments carried out by Teramoto & Yano (2005) and following the same procedures as Mouraille & Luding (2008). Moreover, the same set of simulations were carried out with a Hertzian spring and a dashpot for the normal part of the collisions (Mouraille & Luding 2008). This is,

$$\begin{eqnarray} &&{\bf f}_e=-k_n\xi ^{1.5}{\bf \hat{n}} \ \ \end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray} &&{\bf f}_d=-\gamma \xi ^{0.5}\dot{\xi }{\bf \hat{n}.} \end{eqnarray} \tag{ 30 }$$

For the linear spring model, two sets of parameters were used (L₁ and L₂; Richard et al. 2008; Mouraille & Luding 2008), while for the nonlinear (H) model only one set of parameters was used (Mouraille & Luding 2008). Table 1 summarizes the parameters used for the simulations.

In order to analyze the transmission of the sound wave through the aggregate, it was necessary to obtain a speed profile during the process. The container was divided into slices of height 0.5 mm parallel to the bottom of the container and the instantaneous speed of the particles inside each slice was averaged. This method provided a space averaged instantaneous speed profile.¹ Knowing that the experimental sound speed in glass was about 170 m s⁻¹ (Teramoto & Yano 2005) and that the height of the aggregate was ≈17 mm, this average was calculated every 5 × 10⁻⁶ s during the simulation. This allowed us to have ≈100 speed profiles during the transmission of the sound wave. In order to understand the transmission process and how the natural granularity of the system affected wave transmission, a 3D speed profile was also calculated following the procedure described above, but this time each slice was subdivided into cubes of dimensions 1.5 × 0.5 × 1.5 mm. The time interval that the sound wave took to go through the aggregate was calculated from the moment the first layer of particles was set in motion until when a given monitored layer started increasing its upward velocity steadily (P-waves).

Figure 4 shows an example of the speed of the P-waves obtained through simulations. From these plots it is evident that the L₁ and L₂ set of parameters produce a very similar sound speed behavior regardless of the value of gravity; L₁ being the closest to the trend of values observed by Teramoto (Teramoto & Yano 2005). More importantly the linear model using L₁ parameters and the simulations using a Hertzian spring show very similar results for gravity values of 1g_e and 2g_e; the black horizontal line is only a guide to the eye to show this agreement between simulations and experiments. Disagreements in the numerical values could be rightfully attributed to a number of different causes including the wrong choice of parameters or the correctness of either model. However, it is also evident that for extremely low or for higher values of gravity the behavior of the models diverges and the literature does not contain experimental evidence to validate either model or any set of parameters for the values of gravity here used. Given that both, linear and Hertzian, contact models are generally used in DEM simulations of granular materials, their disagreement, when it comes to sound wave transmission, should be investigated in greater depth. Furthermore, it has been experimentally observed that when 1 m spheres collide, the coefficient of restitution is independent of the impact velocity, a typical behavior of a linear spring contact law (Durda et al. 2010).

Zoom In Zoom Out Reset image size

Another interesting finding regarding sound wave transmission through granular materials is shown in Figure 5. It shows a 3D speed profile of the aggregate exactly when the peak of the wave is passing through this layer of particles. From it one can appreciate that the sound wave is not plane anymore. This indicates that particles in different regions of the aggregate absorb and transmit sound (pressure information) at different rates, a feature that in a cohesive system could produce internal cracks. This is a normal consequence of the random nature of rubble piles, and its implications will require a more detailed study (Scheeres et al. 2010).

Zoom In Zoom Out Reset image size

Let us suppose a spherical asteroid of radius R_A rotating at a speed ω; we study the dynamics of the particles on its surface when the asteroid increases its spin velocity beyond a critical value ω_c that should cause them to go into orbit. In order to find out what this critical value is, and for programming purposes, we will need to make some assumptions and simplifications. Figure 6 is a representation of the system we will be simulating.

Zoom In Zoom Out Reset image size

We will be looking at the aggregate from a frame of reference attached to it, the rotation velocity will be $\vec{\omega }=\omega \hat{k}$ and the y coordinate will be measured from the very surface and not from the bottom of the depicted container. In such a system the motion of a particle i would be governed by the following equations, the same that include a Coriolis term due to the change of frame of reference:

$$\begin{eqnarray} &&f_x^{(i)}=2\omega \dot{r}_y^{(i)}m^{(i)} \end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray} &&f_y^{(i)}=\left[-2\omega \dot{r}_x^{(i)}+\omega ^2\big(R_A+r_y^{(i)}\big)+g_y\right]m^{(i)} \end{eqnarray} \tag{ 32 }$$

$$\begin{equation} f_z^{(i)}=0, \end{equation} \tag{ 33 }$$

where

$$\begin{equation} g_y=\left\lbrace \begin{array}{@{}ll}-\frac{GM_A}{\big(R_A+{r_y^{(i)}}\big)^2} & \mbox{if }\; r_y^{(i)}> 0\\[11pt] -\frac{GM_A}{R_A^3}\big(R_A+r_y^{(i)}\big)& \mbox{if }\; r_y^{(i)}\le 0. \end{array}\right. \end{equation} \tag{ 34 }$$

M_A is the total mass of the asteroid, R_A is the radius of the asteroid, ω is the angular velocity around the z axis, G is the universal gravitational constant, and g_y is the magnitude of gravity. Let us suppose that the particles in the very surface of the aggregate feel no attractive force, that is, they are basically in orbit, but still at rest, similar to the particles in the equator of asteroid 1999-KW4. For this condition to be true the critical spinning velocity of the asteroid has to be

$$\begin{equation} \omega _c=\sqrt{\frac{GM_A}{R_A^3}}. \end{equation} \tag{ 35 }$$

Equation (35) shows a relation between the critical spinning velocity of an asteroid and R_A. If ω>ω_c⇒f⁽ⁱ⁾_y>0 ∀i, which means that the asteroid will desintegrate.

For these simulations used 7998 particles, 2 mm in diameter, with a density of ρ_p = 3200 kg m⁻³ plus two 4 mm intruders on the surface of the aggregate (same density). The physical parameters are as follows: k_n = 9000Nm⁻¹, γ = 0.01Nsm⁻¹, and k_t = 2.57 × 10³Nm⁻¹; only the linear spring–dashpot model was used, the time step δt = 4 × 10⁻⁶ s. The value for k_n was chosen in relation to the value of k_n of the previous section. The particles here are twice as big, so it seems reasonable to deform them twice as much. In addition, given the low gravity environment, the forces at which they are subjected should cause an even smaller deformation than that expected under 1g_e. The particles were contained in a rectangular box of dimensions 20 × 260 × 20 mm with periodic boundaries. This arrangement provides a bulk density of ≈2 kg m⁻³; the total mass of an asteroid with a radius of 8 km, formed with a material of such characteristics has an approximate mass of 4.26 × 10¹⁵ kg. The system was settled as specified in the section about the simulation details. This arrangement attained an average depth of h = 136 mm. The nominal radius of the simulated asteroid was 8000 m which has an surface gravity of 0.0044 m s⁻²; however, the critical angular velocity of the system was calculated for an asteroid of 7999.950 m, which meant that ω_c = 0.00074453 rad s⁻¹. The angular speed was incremented by 0.1ω_c every 10 s of simulation for a total of 110 s; this method would allow us to see how this minimal variation in the angular velocity initiates the disruption process. Calculating an ω_c for a radius just 5 cm smaller than the nominal radius of the asteroid will allow us to see if it is only the upper 5 cm of regolith that is lifted or the entire surface.

3.2.1. Calculation of the Stress Tensor

The stress tensor is not necessary for an SSDEM simulation; however, it was implemented here in order to gain greater insight into the possible role that interlocking of the particles could play in the disruption process. In order to calculate the stress tensor in our simulations, and given that in all the simulations the particles are contained in a box with or without periodic boundary conditions, the box was divided by planes perpendicular to the coordinate axes. The stress tensor was calculated over these planes, periodically and at a given instant. The stress tensor has nine components in three dimensions, where the trace represents the normal stresses along the directions of the axes. The other six components represent the shear stress on these planes. In a liquid at rest, these components are very close to zero; in a granular material, due to friction and the discrete nature of a granular system, this is usually not true.

To calculate the stress tensor on a given plane the code looks at the particles that are in contact and in different sides of this plane (see Figure 7). Then, the stress tensor σ^s_kl is

$$\begin{equation} \sigma _{kl}^s=\frac{1}{S_c}\sum _{i=1}^{N_1}\sum _{j=1}^{N_2}{\vec{\bf f}_{ji}}, \end{equation} \tag{ 36 }$$

where k = x, y, z represents the axis to which the plane over which the calculation is carried out is perpendicular; l = x, y, z represents the component of the total force on that plane; N₁ and N₂ are the number of particles whose centers are in either side of the plane. $\vec{\bf f}_{ji}$ is the force exerted from particle j to particle i given that they are in contact and in opposite sides of the surface S_c. Figure 7 shows a scheme of the forces that are considered for this calculation.

Zoom In Zoom Out Reset image size

The calculation of the stress tensor requires the binning of the particles in volumes; the walls of these volumes are the surfaces on which stress is going to be calculated. In principle the shape of the volume can be arbitrary; however, due to the geometry of the systems to be studied, the geometry of the imposed gravitational field, the coordinate system used (Cartesian) and ease of coding, the walls of these volumes are always parallel to the walls of the container.

Obviously, due to Newton's second law, every force ${\vec{\bf f}_{ji}}$ has its counter part ${\vec{\bf f}_{ij}}$. We have chosen to account for the forces that the particles on top exert over the particles beneath them (in agreement with the direction of gravity), the particles on the right over the particles on the left, and the particles on the front over the particles at the back. This should be changed if a self-gravitational field is to be implemented so that the decision of which set of particles is pushing the other is based solely on the direction of gravity. This is not necessarily in agreement with the usual convention in Continuum Mechanics (normal stress positive in traction, negative in compression).

As in any DEM simulation, the particles are never completely static unless it is so specified, which implies that contact forces between particles are never constant. However, the calculated stress tensor should vary around a stable configuration; in particular, in these simulations there is also an initial, transient variation due to the change in angular velocity that rapidly dies down to become a "solid line." The inset in Figure 8 shows an example such variation of the σ_yy component of the stress tensor for ω = 0.4ω_c with respect to the height of the aggregate during the 10 s that this angular velocity was maintained. In this, and the rest of figures in this section, the height (H) is expressed in terms of the mean size of the particles (d). The pressure goes down to zero when ω_c is reached. Simulations ran for longer periods of time, beyond ω_c, show that all the particles begin to rise and all contact is lost.

Zoom In Zoom Out Reset image size

Equations (32) and (34) can be used to find these values analytically. To calculate σ_yy we need to find the weight per unit area that the particles are exerting over the bottom of the container:

$$\begin{eqnarray} \sigma _{yy}&=&\frac{1}{A}\sum _{i=1}^N f_y^{(i)}\\ &=&\frac{1}{A}\sum _{i=1}^N\left[\omega ^2\big(R_A+r_y^{(i)}\big)-\frac{GM_A}{R_A^3}\big(R_A+r_y^{(i)}\big)\right]. \end{eqnarray} \tag{ 37 }$$

We will substitute the mass of the individual particles for the average particle mass given that the radii of the particles have a uniform distribution between 0.9 and 1.1 mm. In addition, we will substitute r⁽ⁱ⁾_y by the average depth of the particles over the first layer where the stress tensor was calculated (located 5 mm over the bottom of the container). This leads to

$$\begin{equation} \sigma _{yy}=\phi \rho _p h^{\prime }\left(\omega ^2-\frac{GM_A}{R_A^3}\right)\left(R_A-\frac{h^{\prime }}{2}\right), \end{equation} \tag{ 39 }$$

where ϕ ≈ 0.6 is the average filling fraction of the aggregate and h' = h − 0.005. This equation holds as long as there is a contact between the particles and the bottom of the container. For ω = 0, σ_yy = −1.13Pa, in agreement with the simulations. For ω = 0.00077445 rad s⁻¹, the maximum spin rate used in the simulations, σ_yy = 9.3 × 10⁻³Pa. As the particles are not cohesive, a positive value implies that the particles are no longer held by gravity; therefore, making the pressure at the bottom of the container equal to zero, as the simulations clearly show. Equation (39) shows the ω² dependency that the values of pressure at the bottom of the container show. Figure 8 shows a comparison between the values obtained through simulation and by applying Equation (39).

Figure 9 shows the velocity profiles of the aggregate when ω>ω_c at t = 100 s, t = 105 s, and t = 110 s. This figure clearly shows that at this point the particles acquire an upward velocity that depends on their height. In addition, all the particles have an upward velocity and have lost all contact with the bottom of the container, despite the fact that an asteroid with a radius only 5 cm smaller should be still stable at this spin rate. With these simulations we cannot see if in a real asteroid the particles at the center could still be under some pressure; however, it seems plausible that the outermost layer of particles will be disrupted quite rapidly in the absence of cohesive forces.

Zoom In Zoom Out Reset image size

The inset shows the instantaneous value of σ_yy for 100 s < t < 110 s at intervals of 0.1 s. During this last 10 s of simulation, there are almost no particle–particle contacts as there are only three instances where σ_yy ≠ 0 and with values up to six orders of magnitude inferior to those found at the previous spin rate. This shows that friction is not locking the particles together and therefore cannot stop the disruption process.

Given that the particles in a real asteroid are not spherical, as those here simulated, one would expect an amount of interlocking due to their shape; however, it must also be considered the fact that the outward motion of the particles will cause the dilation of the aggregate and a subsequent increase in porosity. This will in turn cause the particles to loose their contacts with other particles, possibly reflecting a behavior similar to that observed here.

As gravity is a long-range interaction, an exact algorithm should calculate it between all the particles in the system. Thus, if a granular system is formed by N particles, the complexity of the force calculation algorithm, which is where the program will spend most of its time, is O(N²). To support computationally tractable simulations when the number of particles is large, we wish to reduce the total number of computations. To do this, we must remember that at sufficiently long distances the mass of two bodies can be considered as condensed at its center of mass and their gravitational interaction can be calculated as for two point masses. Just like for a Tree code, we will take advantage of this approximation. Given the geometry of the calculation of the contact forces, we will divide the available space for the simulation again, using bigger cells (CELLS for notation). These CELLS must contain an integer number of cells to facilitate calculations. So the basic simplification behind this algorithm and a Tree code is the same, but the implementation is different as here we will use a regular and static grid to divide space and in a Tree code the subdivision of space is dynamic and depends on the spatial distribution of the particles.

These two characteristics of the gravitational grid, namely, its regularity and staticity, can be exploited to speed up the calculation of gravity. The position and size of the CELLS will not change throughout the entire simulation. Therefore, it is possible to calculate and store in memory the relative distances ($|\vec{R}|$) between the centers of any two given CELLS and other related quantities that will be needed by the numerical approximation.

As our code must already verify the distance between close neighbors, we can take advantage of this and use it to calculate a "local" gravitational term. First, a gravitational force is calculated among the particles inside each individual cell. Then, we perform the same calculation among the particles in neighboring cells as long as they are in the same CELL. This avoids completely the double calculation of forces.

In order to calculate a long-range gravitational term, we assume that the total mass of the particles inside each CELL can be concentrated in one point, their center of mass (CM), and treat it as a virtual particle. Now, in general, the CM will not coincide exactly with the geometrical center of the CELL, but it will be very close. Close enough as to approximate its location and relative position with respect to the CMs of the other CELLS with some corrections to the relative position vectors that were previously calculated and stored. This numerical approximation is explained in the next section.

Figure 10 shows a 2D depiction of the procedure we have implemented. The red circles represent the virtual particles. The black arrows in Figure 10 represent the long-range interactions to be calculated.

Zoom In Zoom Out Reset image size

Obviously, as the algorithm does not calculate, but approximates, the exact gravitational forces some of it is either over or underestimated. Analytically this is the force that the N_B particles that are contained inside a CELL B are going to exert over a particle j that is an element of a CELL A (Figure 11 shows a 3D scheme of the arrangement and notation we are going to use.) The force over particle j has two components, contact forces and gravitational attraction; therefore,

$$\begin{equation} m\ddot{\vec{r}}_j=\sum \vec{\bf F}_{\rm cont}+\sum \vec{\bf F}_{\rm grav}^{(i)}. \end{equation} \tag{ 40 }$$

For the gravitational component, let $\vec{\bf r}_{ji}=\vec{\bf r}_i-\vec{\bf r}_j$, then

$$\begin{equation} \vec{\bf F}_{ji}=\frac{-Gm_im_j}{r_{ji}^3}\vec{\bf r}_{ji}, \end{equation} \tag{ 41 }$$

which is the gravitational force of m_i over m_j. If we now sum over all the particles contained in CELL B, we obtain

$$\begin{eqnarray} &&\vec{\bf F}_{jB}=\sum _{i=1}^{N_B}\left(\frac{-Gm_im_j}{r_{ji}^3}\right)\vec{\bf r}_{ji}\hphantom{\vec{\bf F}_{BVB}==}\, \end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray} &&=-Gm_j\left(\sum _{i=1}^{N_B}\frac{m_i(\vec{\bf r}_i-\vec{\bf r}_j)}{|\vec{\bf r}_i-\vec{\bf r}_j|^3}\right). \end{eqnarray} \tag{ 43 }$$

From Figure 11, it can be seen that $\vec{\bf r}_{ji}=\vec{\bf r}_i-\vec{\bf r}_j=\vec{\bf R}+\vec{\bf \rho }_i-\vec{\bf \rho }_j$. If we let ${\bf \Delta }{\vec{\bf \rho }}_{ji}=\vec{\rho }_i-\vec{\rho }_j$ we obtain

$$\begin{equation} \vec{\bf r}_i-\vec{\bf r}_j=R\left({\bf \hat{R}}+\frac{{\bf \Delta }\vec{\bf \rho }_{ji}}{R}\right). \end{equation} \tag{ 44 }$$

Then,

$$\begin{equation} \vec{\bf F}_{jB}=-Gm_j\left(\sum _{i=1}^{N_B}\frac{m_i(\vec{\bf R}+{{\bf \Delta }\vec{\bf \rho }}_{ji})}{R^3\left|{\bf \hat{R}}+\frac{{\bf \Delta }\vec{\bf \rho }_{ji}}{R}\right|^3}\right). \end{equation} \tag{ 45 }$$

Now we can work only on the second term of the denominator of the fraction, this is

$$\begin{eqnarray} &&\left|{\bf \hat{R}}+\frac{{\bf \Delta }\vec{\bf \rho }_{ji}}{R}\right|^{-3}=\left(1+\frac{2{{\bf \Delta }\vec{\bf \rho }}_{ji}\cdot {\bf \hat{R}}}{R}+\frac{|{{\bf \Delta }\vec{\bf \rho }_{ji}}|^2}{R^2}\right)^{-\frac{3}{2}} \hphantom{=sd=}\, \end{eqnarray} \tag{ 46 }$$

$$\begin{eqnarray} &&\sim \left(1-3\frac{{\bf \Delta }\vec{\bf \rho }_{ji}\cdot {\bf \hat{R}}}{R}+\cdots \right), \end{eqnarray} \tag{ 47 }$$

assuming $|{\bf \Delta }\vec{\bf \rho }_{ji}|\ll R$ the gravitational force over m_j becomes

$$\begin{eqnarray} \nonumber \vec{\bf F}_{jB}&=&{-}Gm_j\sum _{i=1}^{N_B}\frac{m_i}{R^2}\left({\bf \hat{R}}+\frac{{\bf \Delta }\vec{\bf \rho }_{ji}}{R}\right)\left(1-3\frac{{\bf \Delta }\vec{\bf \rho }_{ji}\cdot {\bf \hat{R}}}{R}+\cdots \right)\\ \nonumber &=&{-}Gm_j\sum _{i=1}^{N_B}\frac{m_i}{R^2}\left({\bf \hat{R}}+\frac{{\bf \Delta }\vec{\bf \rho }_{ji}}{R}-3\frac{{\bf \Delta }\vec{\bf \rho }_{ji}}{R}{\bf \hat{R}\hat{R}}+\cdots \right)\\ \nonumber &=&{-}Gm_j\sum _{i=1}^{N_B}\frac{m_i}{R^2}\left[{\bf \hat{R}}+({\bf I}-3{\bf \hat{R}\hat{R}})\cdot \frac{{\bf \Delta }\vec{\bf \rho }_{ji}}{R}+\cdots \right]\nonumber \\ &\approx& {-}\frac{Gm_jM_B}{R^2}{\bf \hat{R}}-\frac{Gm_j({\bf I}-3{\bf \hat{R}\hat{R}})}{R^3}\sum _{i=1}^{N_B}m_i(\vec{\bf \rho }_i-\vec{\bf \rho }_j). \nonumber \\ \end{eqnarray} \tag{ 48 }$$

Zoom In Zoom Out Reset image size

The summation over the particles in B can be developed as

$$\begin{eqnarray} \sum _{i=1}^{N_B}m_i(\vec{\bf \rho }_i-\vec{\bf \rho }_j)=&\sum _{i=1}^{N_B}m_i\vec{\bf \rho }_i-\sum _{i=1}^{N_B}m_i\vec{\bf \rho }_j \end{eqnarray} \tag{ 49 }$$

$$\begin{eqnarray} &&\hspace{47pt}=M_B\vec{\bf \rho }_B-M_B\vec{\bf \rho }_j, \end{eqnarray} \tag{ 50 }$$

Where M_B is the total mass of all the particles contained in CELL_B and $\vec{\bf \rho }_B$ is the position of their center of mass. Putting all together we obtain

$$\begin{equation} \frac{\vec{\bf F}_{jB}}{m_j}\approx -\frac{GM_B}{R^3}[\vec{\bf R}+({\bf I}-3{\bf \hat{R}\hat{R}})\cdot (\vec{\bf \rho }_B-\vec{\bf \rho }_j)]. \end{equation} \tag{ 51 }$$

Let $\vec{\bf {\it F}}_{jB}$ be the approximated value of $\vec{\bf F}_{jB}$, then

$$\begin{equation} \vec{\bf {\it F}}_{jB}=-\frac{GM_Bm_j}{R^3}[\vec{\bf R}+({\bf I}-3{\bf \hat{R}\hat{R}}){\times }(\vec{\bf \rho }_B-\vec{\bf \rho }_j)]. \end{equation} \tag{ 52 }$$

This is the equation that will be used in simulations. The first term in the last equation, multiplied times $\vec{\bf R}$ is the usual expression to calculate the gravitational attraction between two point masses. The remaining terms give us a first-order correction due to the discrepancy between the geometric center of the CELL and the center of mass of the particles in it. This correction reflects the discrete nature of the particles and the inhomogeneity of the system. If we further develop the last equation, a summation over all m_j will lead to

$$\begin{equation} \vec{\bf {\it F}}_{AB}=-\frac{GM_BM_A}{R^3}[\vec{\bf R}+({\bf I}-3{\bf \hat{R}\hat{R}})\cdot (\vec{\bf \rho }_B-\vec{\bf \rho }_A)]. \end{equation} \tag{ 53 }$$

If the system were completely homogeneous, the difference between the position of the CMs with respect to the center of each respective CELL, $(\vec{\bf \rho }_B-\vec{\bf \rho }_A)$, would be equal to zero so recovering the original form of the gravitational attraction between two point masses.

It is feasible to continue the expansion in Equation (52) to attain higher acuracy, but that is left to future development.

From Equation (52), it can be seen that the values for $|\vec{\bf R}|$, $1/|\vec{\bf R}|$, $1/|\vec{\bf R}|^2$, and $1/|\vec{\bf R}|^3$ will be needed. Therefore, these are the quantities that we will need to store in memory for a quick calculation of gravitational forces. This is how the staticity of the grid is exploited, the grid will not change during the entire simulation; therefore, these quantities are not going to change either.

Let us have a virtual space of dimensions L_x × L_y × L_z where our simulation takes place. In order to conserve the isotropy and homogeneity of space, we should divide it into so many CELLS such that L_x/Gn_x = L_y/Gn_y = L_z/Gn_z, thus producing cubic CELLS, where Gn_x, Gn_y, and Gn_z $\in \mathbb {N}$.

If we were to calculate $\vec{\bf R}$ between the centers of any two given CELLS in the simulation, we would end up with Gn²_xGn²_yGn²_z possible permutations, many of which would result in the same relative position vector. This number can grow very rapidly and can easily surpass the machine memory limit, turning the algorithm useless for practical purposes. This problem can be overcome by noticing that the distance we need to calculate depends only on the absolute values of the differences of the components of the position vectors of each CELL.

Mathematically, let $\vec{\bf R}_A$ and $\vec{\bf R}_B$ be the positions of the centers of CELL_A and CELL_B, respectively, and $\vec{\bf R}=\vec{\bf R}_B-\vec{\bf R}_A$. If we express

$$\begin{eqnarray} \vec{\bf R}_A=&a\hat{i}+b\hat{j}+c\hat{k} \end{eqnarray} \tag{ 54 }$$

$$\begin{eqnarray} \vec{\bf R}_B=&d\hat{i}+e\hat{j}+f\hat{k}, \end{eqnarray} \tag{ 55 }$$

where

$$\begin{eqnarray} {-}Gn_x-1\le &a,d\le Gn_x-1 \end{eqnarray} \tag{ 56 }$$

$$\begin{eqnarray} {-}Gn_y-1\le &b,e\le Gn_y-1 \end{eqnarray} \tag{ 57 }$$

$$\begin{eqnarray} {-}Gn_z-1\le &c,f\le Gn_z-1{,} \end{eqnarray} \tag{ 58 }$$

so that {a, b, c, d, e, f} $\in \mathbb {Z}$, then the square of the distance between the two CELLS will be

$$\begin{equation} |\vec{\bf R}|^2=|d-a|^2+|e-b|^2+|f-c|^2. \end{equation} \tag{ 59 }$$

Therefore, $|\vec{\bf R}|^2\in \mathbb {Z}\ge 0$.

Given that any two CELLS in the simulation space cannot be further apart than the origin to the opposite corner:

$$\begin{equation} 0\le |\vec{\bf R}|^2\le (Gn_x-1)^2+(Gn_y-1)^2+(Gn_z-1)^2. \end{equation} \tag{ 60 }$$

It is the regularity of the grid we are using that allows us to do this. In memory we will only need to store 4[(Gn_x − 1)² + (Gn_y − 1)² + (Gn_z − 1)² + 1] quantities including zero.

At this stage of the development of the algorithm we have chosen to leave the unit vectors to be calculated whenever they are needed. We have done this as many values for $|\vec{\bf R}|^2$ may correspond to more than one unit vector. The determination of the corresponding unit vector from CELL_A to CELL_B would be more time consuming that the multiplication of $1/|\vec{\bf R}|^2$ times a given component of $\vec{\bf R}$.

In order to calculate the gravitational forces with this structure, the particles are binned into the CELLS. This allows the program to know which CELL contains which particles and to store this information in another 2D matrix. An added benefit is that this lets the program know which CELLS contain particles and which do not, so avoiding an extra check over empty space. After doing this, the total mass of the virtual particles and the center of mass of all the particles in each CELL with respect to its geometrical center is calculated.

The algorithm now goes through all the CELLS that contain particles and uses Equation (52) to calculate the force that the particles in CELL_B exert over each particle j in CELL_A. A similar equation is used to calculate the force that the particles in CELL_A exert over each particle i in CELL_B.

Here, we need to show that the procedure and approximations used are valid and the results do not deviate from the exact calculations in a meaningful way. To do this we have ran a test simulation using this algorithm while at the same time calculating the exact gravitational forces on each particle only for comparison. This will allow us to see how close we are and if there is a further refinement that needs to be implemented.

The simulations use a continuum size distribution for the particles with diameters ranging from 8 to $9\,\rm m$ and a density of 3200 kg m⁻³. The following parameters could be considered small compared to the real material parameters and to the parameters used in other SSDEM simulations (Mishra & Rajamani 1992); however, as explained before, the choice of parameters in general has more to do with the specific phenomenon to be analyzed than with the simulated materials. Here, we will only address the performance of the algorithm to calculate gravitational forces. The normal component of the force during collisions is simulated through a linear spring–dashpot model with k_n = 1.8 × 10⁵Nm⁻¹ and η = 3 × 10⁴Nsm⁻¹ (Herrmann & Luding 1998; Richard et al. 2008). Tangential forces, static and dynamic friction, are simulated through a tangential spring (Silbert et al. 2001) with k_t = 5.14 × 10⁴Nm⁻¹. This provides a coefficient of restitution of ≈0.84 and contact duration for the two smallest particles of ≈4.86 s. With these parameters we have chosen a time step for the integration of the equations of motion (δt) of 5 × 10⁻² s allowing about 100 iterations per collision. The particles are initially positioned in a HCP lattice, with a translation vector ($\vec{t}_v$) whose magnitude was 2.2 times the radius of the biggest particle in the distribution. For this and the following simulations the dimensions of the simulation box are 400 × 400 × 400 m and Gn_x = Gn_y = Gn_z = 20. The number of particles is 3680; this was chosen to stack the particles in two complete parallel layers. The particles have no initial velocity for these simulations.

The code so implemented simulated 10 hr of real time. Figure 12 shows the initial and final configuration of the aggregate. It is noticeable that the self-gravitating aggregate formed does not show a perfectly spherical, though stable, configuration. This clearly shows its capacity to sustain internal shear stress which is also a characteristic observed in real asteroids and a very desirable feature in simulations. It must be noted that despite the low value for k_n, the average overlap was 0.006% of the mean particle size, that is, 5 × 10⁻⁴ m, and a maximum overlap of ≈5 × 10⁻³ m.

Zoom In Zoom Out Reset image size

In order to measure the agreement between the exact direct calculation and the approximation of our algorithm, we can define a vector $\Delta \vec{\bf f}_g=\vec{\bf f}_d-\vec{\bf f}_a$, the difference between the exact gravitational force ($\vec{\bf f}_d$) and the force this algorithm has calculated ($\vec{\bf f}_a$). Then we will be able plot $|\Delta \vec{\bf f}_g|/|\vec{\bf f}_d|$ and the angle θ between $\vec{\bf f}_d$ and $\vec{\bf f}_a$. Figure 13 represents the relative and cumulative histograms of these two quantities. From them, it is observable that for about 70% of the particles the error in magnitude is under 5% and for 90% of the particles, θ is under 3^o. These two results show that, even though the forces are not exact, they are well within an acceptable error. It must also be noted that the center of mass of the system remains static through out the entire simulation.

Zoom In Zoom Out Reset image size

In order to test the performance of this algorithm, we have carried out a number of runs with different number of particles. Except for the changing number of particles in the simulation, all the other parameters are kept as before. The simulations ran for 10,000 s of real time to settle the aggregate and then 2000 s during which CPU time was monitored; the number of particles used were 100, 200, 400, 800, 1600, 3200, 6400, and 12,800. The results are shown in Table 2; their time dependence can be fit with an equation of the form

$$\begin{equation} t=4.07096\times 10^{-5}N^2+0.0707153N+17.4444. \end{equation} \tag{ 61 }$$

This shows that the complexity of the code as a whole is O(N²), albeit with a small dependency on N².

As this algorithm works with a static grid, then the less dispersed the particles are in space, the less occupied CELLS in which calculations will have to be performed. To verify this, we simulated 1000 particles inside a box with periodic boundary conditions, which $|\vec{t}_v|$ was increased so that the particles could occupy more than the minimum possible number of CELLS. The program simulates 400 s of real time. The results are tabulated in Table 3.

Table 3. Algorithm Dependence on the Number of CELLS in the Gravitational Grid

$|\vec{t}_v|/r_{\rm max}$	Number of CELLS	Time (s)
2.2	513–306	99
2.5	512–471	137.1
3.0	514–514	149.4
3.5	1026–996	500.6
4.0	1058–1058	420.1
4.5	1506–1506	946.1
5.0	2006–2006	1603.2

Note. The number of CELLS change as the particles aggregate over time.

Download table as:  ASCIITypeset image

Now we can take one more step, and to probe the previous finding, we can keep the size of the box constant and have two different initial conditions. In the first batch of simulations, the particles have a random initial velocity between 0 and 0.3 m s⁻¹ along any axis. In the second batch, this range is between 0 and 0.005 m s⁻¹, just enough as to move the particles so that they are not all in the same plane and the system is essentially bi-dimensional, but not sufficiently large as to move them to another CELL. We do this while increasing the number of particles for comparison. The box is a cube of 400 m per side with periodic boundary conditions. The program simulates 400 s of real time. The results are summarized in Table 4. These results show that the more clustered the particles are, the faster the algorithm becomes. This means that the algorithm here explained is better suited for self-gravitating granular aggregates in a dense phase than in a very dilute one, which is exactly what a rubble-pile asteroid is.

Apart from the previous examples, there are two numerical experiments that seem relevant to the asteroid environment: (1) asteroid formation from a system with no angular momentum ($\vec{L}=0$) and (2) asteroid formation from a system with a given initial angular momentum ($\vec{L}={L_0}\hat{k}$). The former is exactly the same as the one we have just analyzed in the previous section and the later, having exactly the same parameters, has also been given an initial angular velocity. Figure 14 shows the configuration of the rotating aggregate at t = 0 and t = 10 hr.

Zoom In Zoom Out Reset image size

Figure 15 shows the variation of the magnitude of the angular velocities of both experiments in rad s⁻¹. For both plots, a mean radius of 150 m was assumed; the total mass of the aggregate is ≈3.74 × 10⁹ kg and the moment of inertia was calculated as for a solid sphere.

Zoom In Zoom Out Reset image size

As expected for a DEM simulation, angular momentum is only conserved within a degree of accuracy, given that this is a DEM simulation and we are using approximations of the real forces that are involved. However, it is remarkable that the variations are in the order of 10⁻⁷ rad s⁻¹, orders of magnitude smaller than the typical rotation rate for near-Earth objects (NEOs) and main belt asteroids (MBAs) (Scheeres et al. 2004), even for slow rotators.

This algorithm, as it is presented here, is a first attempt, albeit not unrefined, to build a simple way to calculate self-gravitating forces in a granular system. Therefore, further improvements are always possible. Here we will only enumerate some that seem to be the most obvious.

• 1.  
  It should be possible to periodically calculate the exact gravitational forces and a vector $\Delta \vec{\bf f}_g$ to correct the approxi- mated force. This should minimize the error of the approxi- mation and would also allow us not to have to calculate gravitational forces at every time step, still maintaining the same level of exactitude.
• 2.  
  Here all forces are calculated at every time step. An avenue of further improvement may come from an adaptive updating method, dependent on the kinetic activity of the system, its temperature or some other magnitude that depends on the motion of the particles and the aggregate as a whole. This would reduce the complexity of the simulation to the complexity of the contact detection algorithm, the link list, which is O(Nlog N), when gravity is not calculated.
• 3.  
  Another way to improve the algorithm is by improving the numerical approximation of the gravity expansion without increasing the simulation time or the complexity of the code.
• 4.  
  At the moment the algorithm uses two grids to calculate gravity. It should be possible to repeat the process with a third grid of larger proportions that could be used for particles far away from the main aggregate. If done recursively, this would effectively transform this algorithm into a static version of an octree configuration.