Numerical simulations of granular dynamics. I. Hardsphere discrete element method and tests
Derek C. Richardson
Kevin J. Walsh
Naomi Murdoch
Patrick Michel
Department of Astronomy
University of Maryland
College Park MD 207402421
University of Nice Sophia Antipolis, CNRS
Observatoire de la Côte d’Azur, B.P. 4229
06304 Nice Cedex 4, France
The Open University
PSSRI, Walton Hall
Milton Keynes, MK7 6AA, UK
Printed July 14, 2019
Submitted to Icarus
54 manuscript pages
8 figures including 4 in color (online version only)
Proposed running page head: Granular dynamics method and tests
Please address all editorial correspondence and proofs to:
Derek C. Richardson
Department of Astronomy
Computer and Space Sciences Building
Stadium Drive
University of Maryland
College Park MD 207422421
Tel: 3014058786
Fax: 3013149067
Email: dcr@astro.umd.edu
Abstract
We present a new particlebased (discrete element) numerical method for the simulation of granular dynamics, with application to motions of particles on small solar system body and planetary surfaces. The method employs the parallel body tree code pkdgrav to search for collisions and compute particle trajectories. Collisions are treated as instantaneous pointcontact events between rigid spheres. Particle confinement is achieved by combining arbitrary combinations of four provided wall primitives, namely infinite plane, finite disk, infinite cylinder, and finite cylinder, and degenerate cases of these. Various wall movements, including translation, oscillation, and rotation, are supported. We provide full derivations of collision prediction and resolution equations for all geometries and motions. Several tests of the method are described, including a model granular “atmosphere” that achieves correct energy equipartition, and a series of tumbler simulations that show the expected transition from tumbling to centrifuging as a function of rotation rate.

Keywords: Asteroids, surfaces; Collisional physics; Geological processes; Regoliths
1 Introduction
The understanding of granular dynamics is taking an increasingly prominent place in the field of planetary science. Indeed, we now know that granular material, in the form of regolith, is covering the uppermost layer of most solid bodies in our solar system, from planets and their satellites to asteroids and comets. This presence of a regolith layer plays a particularly important role in the surface geology of asteroids. The same can be stated, although to a lesser extent, for bodies like Mars and the Moon, whose surface gravities are also smaller than that of Earth. Thus, flows of granular materials driven by different gravitational conditions are particularly important in the understanding of the geology of small bodies and planets. The exploration of the Moon and Mars in the next two decades will require deployment of landing vehicles on surfaces of loose granular material. Space missions to small bodies also involve measurements by landers (e.g., the Rosetta space mission of the European Space Agency) and sampling devices capable of coping with a wide range of surface properties. Therefore, understanding how granular materials, as a function of their properties (angle of friction, size distribution of their components, etc.), react to different kinds of stresses is of great interest for the design of landers and sampling devices of space missions. The same holds true in the framework of mitigation strategies against a potential impactor, which involve an interaction with the body’s surface.
The presence or relative absence of gravitational acceleration on granular flow is of importance for understanding the geology of small bodies and planets, and to clarify the environments that may be encountered during planetary exploration. Bodies with low surface gravity can be very sensitive to processes that appear irrelevant in the case of larger planetary bodies. For instance, seismic vibration induced by small impacts can occur throughout a small body and can be at the origin of motion of its granular surface. Such a mechanism has been proposed to explain the lack of very small craters on both asteroid 433 Eros (Richardson et al. 2004) and asteroid 25143 Itokawa (Michel et al. 2009). Shaking can also drive size sorting/segregation in granular media. This has been observed on Earth (Rosato et al. 1987), as well as on Eros (Robinson et al. 2002) and Itokawa (Miyamoto et al. 2007). In particular, on asteroid Itokawa it has been observed that the neck is formed mainly by centimeter or smallersized regolith but that both the head and body are covered by meterscale boulders. This finding has led some to suggest that subsequent impacts on an asteroid could provide the energy and driving mechanism for segregation to occur (Asphaug et al. 2001, Miwa et al. 2008), leading to a phenomenon commonly called the “Brazilnut effect.” However, it was then found for Itokawa that an erosion mechanism was not sufficient to explain the selective location of coarse material in the potential lows and highs, and that other factors, such as cohesion, rotation, or avalanches may be involved (Sanchez et al. 2010). In fact, rotational spinup of a small body as a result of the thermal YORP effect can lead to regolith motion on the surface (Scheeres et al. 2007), and even to the escape of material that can eventually result in the formation of a satellite (Walsh et al. 2008).
All these aspects demonstrate the great importance of understanding the dynamics of granular matter for planetary science applications. The study of granular dynamics requires dealing with complex mechanical processes and currently constitutes an entire field of research by itself. As examples of the complexity, it has been found that granular matter may strengthen under a variety of conditions, for instance as a result of a slow shift in the particle arrangement under shear stress or because of humidity (Losert et al. 2000b). Moreover, while granular materials are a discrete medium, their flow in response to stress can, under certain conditions, be described by continuum models (Savage 1998; Goddard 1990; Losert et al. 2000). However, many processes cannot be captured using classical continuum approaches. For example, highly stressed granular materials exhibit shear localization during failure (Mueth et al. 2000) but this localization cannot be described accurately by existing continuum models (Kamrin et al. 2007). In addition, situations have been identified in which material under shear stress weakens when the shear direction is changed/reversed (Toiya et al. 2004; Falk et al. 2008). Finally, at low gravity, other physical forces relevant to regolith in the asteroid environment, such as van der Waals forces, may exceed the particle weights and require consideration (Scheeres et al. 2010).
Improving our understanding of the dynamics of granular materials under a wide variety of conditions requires that both experimental and numerical work be performed and compared. Once numerical approaches have been validated by successful comparison with experiments, then they can cover a parameter space that is too wide for or unreachable by laboratory experiments.
There are several different approaches that have been used to perform modeling of granular materials (Mehta 2007). One commonly used technique is the discrete element method (DEM). DEM is a numerical method for computing the motion of a large number of particles of micronscale size and above. Though DEM is very closely related to molecular dynamics (in which atoms and molecules are allowed to interact for a period of time by approximations of known physics), the method is generally distinguished by its inclusion of rotational degrees of freedom, inelastic collisions, and oftencomplicated geometries (including polyhedra; e.g., Cleary and Sawley 2002; Fraige et al. 2008; Latham et al. 2008; Szarf et al. 2010).
However, the DEM method remains relatively computationally intensive, which limits either the length of a simulation or the number of particles. Several DEM codes, and related molecular dynamics codes, take advantage of parallel processing capabilities to scale up the number of particles or length of the simulation (e.g., Cleary and Sawley 2002; Kacianauskas et al. 2008).
Hardsphere particle dynamics have been used successfully in many granular physics applications. Hong and McLennan (1992) used hardsphere molecular dynamics to study particles flowing through a hole in a twodimensional box under the influence of gravity. Huilin et al. (2007) used an EulerianLagrangian approach coupled with a discrete hardsphere model to obtain details of particle collision information in a fluidized bed of granular material. Also Kosinski and Hoffman (2009) compared the standard hardsphere method including walls to a hardsphere model with walls that also accounts for particle adhesion. The van der Waals type interaction is presented as a demonstration case.
An alternative to treating all particles separately is to average the physics across many particles and thereby treat the material as a continuum. In the case of solidlike granular behavior, the continuum approach usually treats the material as elastic or elastoplastic and models it with the finite element method or a meshfree method (e.g., Elaskar et al. 2000; also see Holsapple 2001 and 2004, Holsapple and Michel 2006 and 2008, and Sharma et al. 2006 and 2009 for use of analytical and continuum approaches in modeling asteroid shapes). In the case of liquidlike or gaslike granular flow, the continuum approach may treat the material as a fluid and use computational fluid dynamics. However, as explained previously, the homogenization of the granularscale physics is not necessarily appropriate, and the discreteness of the particles and the forces between particles (and walls) need to be taken into account (Wada et al. 2006). Therefore limits of such homogenization must be considered carefully before attempting to use a continuum approach.
Recently, numerical codes have been developed to address specifically the dynamics of granular materials in the framework of planetary science. For instance, Sanchez et al. (2010) simulated the particles forming an asteroid by means of softsphere molecular dynamics. In this approach, particles forming the aggregate have short and longrange interactions, for contact and gravitational forces respectively, which are taken into account using two types of potentials (see Sanchez et al. 2010 for details).
In this paper, we present our method to simulate the dynamics of granular materials, along with a few basic tests that show its ability to address different problems involving these dynamics. A more complete validation by comparison with a welldocumented laboratory experiment of seismic shaking is the subject of a followup paper (Murdoch et al. 2011, in preparation). We use the body code pkdgrav (Stadel 2001), adapted for hardbody collisions (Richardson et al. 2000; Richardson et al. 2009). The granular material is therefore represented by hard spheres that interact via impulsive, pointcontact forces. Advantages of pkdgrav over other many discrete element approaches include full support for parallel computation, the use of hierarchical tree methods to rapidly compute longrange interparticle forces (namely gravity, when included) and to locate nearest neighbors (for shortrange Hooke’slaw type forces) and potential colliders, and options for particle bonding to make irregular shapes that are subjet to Euler’s laws of rigidbody rotation with noncentral impacts (cf. Richardson et al. 2009). In addition, collisions are determined prior to advancing particle positions, insuring that no collisions are missed and that collision circumstances are computed exactly (in general, to within the accuracy of the integration), which is a particular advantage when particles are moving rapidly.
Here we focus on the implementation of a wide range of boundary conditions (“walls”) that can be used to represent the different geometries involved in experimental setups, but more generally provide the needed particle confinement and possible external forcings, such as induced vibrations, for smallscale investigations of regolith dynamics in varying gravitational environments (different surface slopes, etc.). Our approach is designed to be general and flexible: any number of walls can be combined in arbitrary ways to match the desired configuration without changing any code, whereas many existing methods are tailored for a specific geometry. The longterm goal is to understand how scaling laws, different flow regimes, segregation, and so on change with gravity, and to apply this understanding to asteroid surfaces, without the need to simulate the surfaces in their entirety. We present the full detail of our numerical method in Section 2, describe basic tests of the method in Section 3, and offer discussion and conclusions in Section 4 and Section 5, respectively.
2 Numerical Method
This section details the numerical method employed in our granular dynamics simulations. We use pkdgrav, a parallel body tree code (Stadel 2001) that has been adapted for hardbody collisions (Richardson et al. 2000; Richardson et al. 2009). In what follows we present the general strategy of searching for and resolving collisions among particles, and between particles and the “walls” of the simulated apparatus. This is followed by detailed derivations of the collision equations for the four wall “primitives” that we have developed so far, namely the infinite plane, finite disk, infinite cylinder, and finite cylinder, plus degenerate cases of these.
2.1 General Strategy
In our approach, we predict when collisions will occur based on trajectory extrapolation at the beginning of each integration step. We use a secondorder leapfrog integrator in pkdgrav: each step consists of “kicking” particle velocities by a half step (keeping particle positions fixed), “drifting” particle positions at constant velocity for a full step, recomputing accelerations due to gravity, then performing a final halfstep velocity kick (see Richardson et al. 2009 for details). Collision searches are performed during the drift step by examining the trajectories of enough neighbors of each particle to ensure no collisions are missed. Since there is no interparticle gravity for the granular dynamics experiments discussed here (only a configurable uniform gravity field), the equations of motion are particularly simple, but we nonetheless use the full machinery of the body code for these simulations—the tree is still used to speed up neighbor searches, and parallelism reduces computation time for large simulations. Also, more complex gravity fields for which there is no analytical solution for particle motion will be used in the future.
The condition for two particles with initial vector positions and and velocities and to collide after a time interval (measured from the start of the drift step in this case) is
(1) 
where and are the radii of the particles (treated as perfect hard spheres), or
(2) 
where we have introduced the relative position and velocity (Fig. Figure Captions). To collide, the particles must be approaching one another, so . We assume the particles are not initially touching or overlapping, i.e., we require , where . Equation (2) can be solved for by squaring both sides and applying the quadratic formula:^{1}^{1}1In practice, a version of the quadratic formula optimized to reduce roundoff error is used in pkdgrav—cf. Press et al. 2007.
(3) 
where and the sign ambiguity is resolved by taking to be the smallest positive value of any real roots. If is between zero and the drift step time interval (i.e., the simulation time step), the collision takes place during this integration step. Collisions are treated as instantaneous events with no flexing and a configurable amount of energy loss due to restitution (parameterized by the normal coefficient of restitution , where 0 means perfect sticking and 1 means perfect bouncing) and surface coupling (parameterized by the transverse coefficient of restitution , where means reversal of transverse motion on contact, 0 means complete damping of transverse motion, and 1 means no surface coupling). For completeness, the collision resolution equations for perfect spheres are (Richardson et al. 2009, Appendix B; also see Richardson 1994, Fig. 1):
where the primes denote postcollision quantities, is the sum of the particle masses, is the total relative velocity at the contact point (with , , where is the spin vector of particle , , and ), , , is the reduced mass , and is the moment of inertia of particle ; the factors of 2/7 result from using spheres. These equations come from conservation of linear and angular momentum combined with the following statement of energy loss:
(4) 
[FIGURE Figure Captions GOES HERE]
The approach for handling wall collisions is similar: first the time to collision is determined, then the collision is resolved. The collision condition is that the distance of the particle center from the point of contact on the wall must equal the particle radius, i.e.,
(5) 
where is the vector position of the point of contact and we have dropped the subscript . (For the derivations that follow, we have adopted the convention that the wall corresponds to and the particle to ; this minimizes the number of minus signs in the equations.) Because depends on the particular wall geometry, it is best to consider each geometry in turn and exploit any symmetries to determine more conveniently when (if at all) this contact condition is satisfied for a given particle and wall pair. Note that often a wall consists of a combination of more than one geometry (for example, a finite disk is a round portion of a plane plus a surrounding ring, i.e., a cylinder of zero length)—the collision prediction routines for walls in pkdgrav return a list of possible collision times, one for each “face” (for the finite disk example, there are 3 faces: one flat side, the opposite flat side, and the perimeter ring), with the smallest positive time corresponding to the face that will be struck first.
In pkdgrav, after all neighbors of a given particle have been checked for possible collisions, every wall is checked also. Whatever body (whether particle or wall) that gives the shortest time to collision is taken to be the next collision event. That collision is resolved, collision times are updated depending on whether the collision that just happened changes future collision circumstances, and the next collision is carried out until all collisions during the current drift interval have been handled. Walls are treated as having infinite mass, so the restitution equations reduce to:
(6)  
(7) 
where
(8) 
(for a static wall), and are specific to the wall (see below), and , with being the unit vector in direction , i.e., from the point of contact (which depends on the wall geometry) to the particle center. To account for any wall motion, in Eq. (8) is adjusted as needed (see below).
The parameters of each wall type are specified by the user at run time via a simple text file. Information common to each wall that must be specified includes the origin (a vector reference point), orientation (a unit vector), and for the wall (which override the particle values during a particlewall collision), and color and transparency (for drawing purposes). In addition, each geometry has a unique set of parameters particular to that geometry, e.g., radius in the case of the finite disk and cylinders, and length for the finite cylinder. Also, simple motions are supported for each geometry; e.g., a plane or disk can have translational and/or sinusoidal oscillatory motion, while a cylinder can rotate around its symmetry (orientation) axis. Oscillations are treated by updating the wall position in stepwise fashion, holding the velocity constant between timesteps. In pkdgrav, a single function is used to get the position and velocity of a particle relative to any (possibly moving) wall at a given time. Finally, a wall can be “sticky,” meaning any particles that come into contact with it become rigidly held, or “absorbing,” meaning any particles that come into contact with it are removed from the simulation. These latter options are encoded as special cases of the wall normal coefficient of restitution (0 indicating sticky, indicating absorbing). The sticky wall is particularly useful for creating a rough surface (cf. Section 3.2); the absorbing wall is helpful when particles are no longer needed, such as when they have moved beyond the flow regime of interest.
For a particle hitting another particle that is rigidly stuck to a wall, Eqs. (6) and (7) also apply when resolving the collision (the stuck particle is treated as having infinite mass), with pointing from the center of the stuck particle to the center of the free particle, and and being the usual particle values (i.e., not the values specific to the wall). Detection of collisions with stuck particles is handled by setting the stuck particles’ velocities equal to the instantaneous velocities of their host walls at the start of the step. In the case of particles stuck on rotating cylinders, the component is added to the particle velocities, where is the cylinder’s rotation vector and is its radius, and in this case is the perpendicular from the cylinder rotation axis to the particle center at the start of the step (the sign of depends on whether the particle is on the outer [positive] or inner [negative] surface). The positions of particles stuck on moving/rotating walls are updated at the end of the drift step. Collisions with particles stuck on rotating cylinders can optionally be predicted to higher order using a quartic expression that accounts for the curvilinear motion (cf. Richardson et al. 2009, Eq. (A.4)).
For resolving collisions with moving walls, and particles stuck to moving walls, each wall’s instantaneous velocity at the start of the step is substracted from in Eq. (8) before applying Eqs. (6) and (7). For cylinder rotation, with the rotation axis parallel to the symmetry axis of the cylinder, is subtracted, where is as defined for those equations. For particles stuck to rotating walls, the stuck particle itself is given an extra speed component (see above), which is taken into account when resolving the collision outcome.
2.2 Specific Geometries
The following sections detail the implementation of the supported wall geometries in pkdgrav.
2.2.1 Infinite plane
The parameters for the infinite plane are the origin and normal , plus optional velocity , oscillation amplitude , and oscillation angular frequency (so the relative vector displacement after time due to oscillation, measured from the start of the simulation and evaluated at the start of the step, is ). The origin can be any point in the plane (the choice is arbitrary).
To simplify the equations in this and subsequent derivations, we define the relative position vector and separate it into perpendicular and parallel components, and , respectively, where , , and (so , which is only defined if and means generally has a time dependence; for completeness, we also define ). It is important to note that and are signed quantities (vector components), unlike defined previously, which is an unsigned magnitude. We similarly define the relative velocity , with corresponding perpendicular and parallel components.
Returning to the specific case of the infinite plane, the collision condition is , i.e., that the perpendicular distance from the particle center to the surface (the height above or below the plane) at the time of impact is equal to the particle radius (Fig. Figure Captions). For this geometry, and measured from the start of the drift step, the time to impact is
(9) 
where the first condition corresponds to impact with the “upper” face (the surface out of which the normal vector points) and the second corresponds to the opposite “lower” face. (We take this twocase approach to avoid having to solve a quadratic; cf. Section 2.2.3.) Note is a requirement for collision (otherwise the particle is moving away from or parallel to the plane). For completeness, the wall position at the start of the step is , where is the origin at the start of the simulation, and the wall velocity is .
[FIGURE Figure Captions GOES HERE]
2.2.2 Finite disk
This geometry uses the same parameters as the infinite plane (Fig. Figure Captions) but now includes the finite radius of the object. The origin is the geometric center of the disk. Collision prediction proceeds in 3 stages: first, Eq. (9) is used to determine the impact time (if any) with the infinite plane that contains the disk (i.e., taking ); second, the position the particle would have relative to the (possibly moving) disk at that impact time is checked to determine whether , i.e., that the separation of the particle center from the disk origin, projected onto the plane, would be less than the disk radius, indicating actual contact with the flat surface of the finite disk; third, the impact time (if any) with the disk perimeter, i.e., a ring, is checked (see Section 2.2.5)—the smallest positive time for impact with both faces and the perimeter is taken to be the next impact time with this object.
Collision resolution is the same as for the infinite plane, unless the particle has made contact with the perimeter, in which case the collision is treated as a ring bounce (Section 2.2.5).
2.2.3 Infinite cylinder
The parameters for the infinite cylinder are the origin , symmetryaxis orientation , radius , and optional angular frequency (with as the spin axis). The origin is any point along the symmetry axis. In the current implementation, cylinders are fixed in space, and since the rotation doesn’t change the orientation (it only affects the tangential motion and spin of particles that come into contact with the cylinder, assuming ), is constant (so ). The collision condition is , i.e., the perpendicular distance from the particle center to the cylinder wall is equal to the particle radius. Note that particles can either be inside (assuming ) or outside the cylinder; each face is considered in turn. Figure Figure Captions illustrates the geometry.
[FIGURE Figure Captions GOES HERE]
To find the collision time (if any), note , where recall , in this case, and is the impact time measured from the start of the drift step. Substituting into the collision condition equation,
(10) 
Taking the square,
(11) 
This is a quadratic equation for and can be solved in the usual way. For completeness,
(12) 
Note there are four solutions, two each for impacts with the inner face (corresponding to on the righthand side of Eq. (10)) and outer face (). Again, for each face, the smallest positive (real) value for is chosen. If , there will be no impact (the particle is either stationary with respect to the cylinder, or moving parallel to the cylinder’s symmetry axis). If the particle is outside the cylinder (so ), a collision can only occur if ; for an infinite cylinder, the collision can only be with the outer face. If the particle is inside the cylinder (so ), all solutions are real (i.e., there must be a collision, assuming ); furthermore, only collisions with the inner face are possible (even in the case of a finite cylinder). Together these facts can be used to reduce the computation time involved with this geometry. Note however, unlike other collision circumstances handled by pkdgrav, in the case of an interior bounce involving a cylinder (or a ring), it is possible for the particle to hit the same object two or more times in one integration step without an intervening collision with another object.
To resolve the collision, Eqs. (6) and (7) apply, with being either , if the outer face was struck (requiring ), or , if the inner face was struck (requiring ). Recall a cylinder can be spinning, in which case is subtracted from the total relative velocity in Eq. (8) before applying the collision equations.
2.2.4 Finite cylinder
This geometry uses the same parameters as the infinite cylinder, but now includes the finite length of the object. The origin is the geometric center of the cylinder. In a manner similar to the finite disk case (Section 2.2.2), collision prediction proceeds in 3 stages: first, Eq. (12) is used to determine the impact time (if any) with the infinite counterpart to the cylinder (i.e., taking ); second, the position the particle would have relative to the cylinder at that impact time is checked to determine whether , i.e., that the separation of the particle center from the cylinder origin, projected onto the cylinder’s orientation axis, would be less than half the cylinder length, indicating actual contact with the surface of the finite object (not the ends); third, the impact times (if any) with the cylinder ends are checked (see Section 2.2.5)—the smallest positive time for impact with both faces and the cylinder ends is taken to be the next impact time with this object. Three special cases for this geometry are supported: is a ring (Section 2.2.5); is a line (so there is no inner surface); and is a point, which is treated as a special case of an infinitesimal ring.
Collision resolution is the same as for the infinite cylinder, unless the particle has made contact with either end, in which case the collision is treated as a ring bounce (Section 2.2.5).
2.2.5 Ring
This is a special degenerate case of a zerolength cylinder, but it is also needed for both the perimeter of a finite disk and the ends of a finite cylinder. The parameters are inherited from whatever base type needs to consider ring collisions, i.e., in addition to the origin and orientation , the parameters are the radius of either the finite disk or the finite cylinder, and either , , and for the disk or for the cylinder—see the corresponding sections above.
Solving for the impact time of a sphere of radius with a thin ring is equivalent to finding the first intersection of a ray with a torus of tube radius , a common problem in computer ray tracing. We use the approach developed by Wagner (2004), for which the possible impact times are roots of the quartic equation
(13) 
where
(14)  
(15)  
(16)  
(17)  
(18) 
and
(19)  
(20)  
(21) 
for which and have the usual definitions (so the ray is given by ). However, this derivation only applies to the special case and , i.e., a ring centered at the origin of a Cartesian coordinate system and lying in the plane (so the normal is in the direction). For the general case, and must first be transformed (rotated) into the torus frame according to
(22)  
(23) 
where the primes indicate postrotation quantities,
(24) 
and and are constructed from using GramSchmidt orthonormalization (note rotation around the torus symmetry axis is unimportant, so we are free to choose any and so long as , , and are mutually orthogonal). Note , as required.
We solve the quartic using Laguerre’s method (Press et al. 2007), a rootfinding scheme optimized for polynomials that does not require prior bracketing of the roots.
There are several special cases that are checked for before committing to solving the quartic. In the case of (a point), Eq. (3) applies immediately, with and . In the case of (particle center on the ring symmetry axis) and (particle moving parallel to the symmetry axis), an arbitrary point is chosen on the ring and Eq. (3) is used to find the impact time. Also, as usual, if (), there can be no collision.
To resolve the collision, Eqs. (6) and (7) apply, with
(25) 
(So in the degenerate case of , i.e., contact with a point, .) For the special case of the particle center lying on the orientation axis, with (indicating the particle has struck an entire end perimeter perfectly and thus is undefined), , i.e., the same as for a plane bounce (Section 2.2.1).
2.3 Overlap Handling
The derivations above make the tacit assumption that particles are neither touching nor intersecting any geometric shapes (including other particles) at the start of the drift interval, i.e., legitimate impact times should always be positive. Unfortunately, due to finite computer precision, sometimes a particle ends up in contact with or interpenetrating another object at the start of the step, e.g., a collision during the previous step may have been missed, or roundoff error led to an overlap. The simplest fix is to allow to be negative (or zero) to correct such situations. It does, however, greatly complicate the decision logic when computing for the various wall geometries discussed above.
As an example, in the case of the infinite plane, a particle is touching or overlapping if , in which case the smallest nonpositive for impact is chosen. In practice, a modified version of Eq. (9) is used when computing : both faces are considered regardless of which side the particle center is on, allowing for the possibility of initial touching/overlap (which is revealed by one—but not both—of the ’s being negative, or one being negative and the other being zero). Note that it is still the case that collisions are ignored if or , regardless of whether the particle is touching or overlapping the surface—generally overlaps are a tiny fraction of the particle radius, whereas the condition with implies the particle has penetrated more than a full radius.
In the interest of brevity, we do not provide exhaustive casebycase touching/overlap remedies for all geometries here. It is sufficient to note that overlaps do happen, particularly in closepacked situations where a great many collisions can occur during a single timestep, but the errors are typically very small (many orders of magnitude smaller than the particle radius) and can be corrected by allowing to be negative (or zero). The test suite described in Section 3 was designed to stress the code and shows that any overlap artifacts encountered are negligible. Pkdgrav incorporates other overlap strategies for situations involving interparticle gravity and/or more complex aggregated particle shapes, but discussion of these strategies is beyond the scope of this paper.
3 Basic Tests
A suite of simple tests was developed during the coding stage to exercise each new feature/geometry as it was added. Examples include dropping particles into a cylinder bisected by an inclined plane, “rolling” a perfectly balanced sphere on the inside edge of a vertical ring, and bouncing hundreds of spheres on an oscillating disk inside a cylinder. Here we report on two tests designed to demonstrate correct dynamic behavior of the granular assembly. In the companion paper (Murdoch et al. 2011, in preparation), we present a detailed study of grain motions on a vibrating plate.
3.1 Model Atmosphere
For this test, a ball of approximately 1,000 closepacked particles is dropped from rest into the top half of an infinitely long vertical cylinder bisected by an infinite horizontal plane. A uniform gravity field points vertically downward (there is no interparticle gravity) and collisions among the particles and between particles and the walls are elastic (so there is no dissipation). The particles are roughly equally divided among three different masses: 1, 3, and 10 (arbitrary units; for this simulation, the cylinder radius is 1.0, the particle radius is 0.022 for all masses, the drop height of the center of the particle ball is 1.0, and the acceleration in the vertical direction is ; internally, the universal constant of gravitation ). The expected behavior, after the transient splash, is that the particles achieve energy equipartition, with the smallestmass particles reaching a scale height 3 times higher than that of the mediummass particles and about 10 times higher than that of the most massive particles.
Figure Figure Captions shows raytraced snapshots of the system, from the initial condition to the final equilibrium state. Figure Figure Captions shows the evolution of the mean heights of the three particle mass populations. The “height” of a particle is the distance separating its bottom surface from the plane, i.e., the coordinate of its center minus its radius, with the horizontal plane located at . Our hypothesis is that the height of a particle with mass is drawn from a probability distribution , where is the scale height associated with particles of mass . The best estimate of for our model (for a given ) is found by maximizing the likelihood
(26) 
Substituting for our model , we find is maximized when , where is the number of particles of mass ; in other words, the best estimate of is the average height of all particles with that mass. This is what is shown in Fig. Figure Captions (solid lines). The measured scale heights show the expected inverse proportionality with mass, but exhibit some fluctuations. In fact, a perfect match to the model is not expected due to finitesize effects: in the limit of large particle sizes compared to the cylinder diameter, the bottom layers may become fully occupied, which is not predicted by the functional form of used here. To test this, a few runs varying the cylindertoparticle diameter ratio were carried out and showed larger discrepancies from the ideal scale height ratios for smaller diameter ratios.
[FIGURE Figure Captions GOES HERE]
[FIGURE Figure Captions GOES HERE]
We also performed KolmogorovSmirnov (KS) tests to compare the height data with the expected distributions . Figure Figure Captions shows the KS curves (solid lines) for the final snapshot in Fig. Figure Captions compared to the normalized cumulative probability distributions (dotted lines), where and is the largest value of (among particles of mass ) for this snapshot. The normalization ensures the cumulative probability distribution has a value of unity at (each curve has a different , but the KS statistic is insensitive to scaling of the independent variable, which is why there are no tick marks on the horizontal axis of Fig. Figure Captions). For this particular snapshot, which represents an arbitrary instant after equilibrium has been reached, we find KS probabilities of 0.89, 0.08, and 0.70, corresponding to masses 1, 3, and 10, meaning we cannot reject the null hypothesis that the data is drawn from at those confidence levels. For the middle mass, this is not very convincing, but we find there is significant variation between snapshots. Taking the average probabilities from 50%, 60%, 70%, 80%, 90%, and 100% of the way through the simulation, we find mean KS probabilities of , , and (the errorbars denote 1 standard deviation from the mean). Moreover, if we repeat the analysis for a cylinder that is half as wide, the probabilities drop to , , and . We must caution that because the mean and maximum values measured from the data were used in the model for the KS test, the distribution of the KS statistic may not be fully consistent with the formula used to compute the null hypothesis probability (cf. Press et al. 2007). Repeating the analysis using the expected values of (see below), the mean KS probabilities for the radius 1.0 cylinder do not change significantly. Regardless, there is no way to avoid using to insure the cumulative probability is normalized, so the KS test probabilities must be considered formal values only.
[FIGURE Figure Captions GOES HERE]
Finally, we can apply the principles of hydrostatic equilibrium and the energy equipartition theorem to derive the expected scale heights and compare these with the measured values for each mass population. Treating the particles as an ideal gas at constant temperature (an ideal isothermal atmosphere), the number density of particles is proportional to the pressure . In hydrostatic equilibrium,
(27) 
where and is the (assumed constant) gravitational acceleration. Substituting from the ideal gas law, where is the Boltzmann constant, and integrating, we find
(28) 
which, since is proportional to , gives the desired scale height
(29) 
Now, the total energy of the system, which remained constant to better than 1 part in for the duration of all the model atmosphere tests we performed (the starting value was in system units), should be roughly equally divided among the particles once equilibrium is achieved. The energy of a given particle in our system is
(30) 
where is the particle speed. Applying the equipartition theorem, which states in part that for each degree of freedom that contributes to the energy, the average energy corresponding to that degree of freedom at thermal equilibrium is , the average energy per particle is
(31) 
(recall incorporates three degrees of freedom, one for each component of the velocity vector in three dimensions). Equating this to the total energy divided by the number of particles (exactly 991 in our simulation), solving for and substituting into Eq. (29), we get expected values of 1.88, 0.626, and 0.188 for our particle masses of 1, 3, and 10, respectively (notice these are in the correct 10:3:1 ratio; the values are shown as the dotted lines in Fig. Figure Captions). Averaging as we did previously for the KS probabilities, we find measured values of , , and (cf. Fig. Figure Captions), in quite good agreement. For the simulation where the radius is halved, the values are , , and , in worse agreement, as might be expected due to the greater severity of the finitesize effect.
In summary, we find that this test of the integrator with particle and wall bouncing gives results consistent with what would be expected in a real system.
3.2 Tumbler
Brucks et al. (2007) carried out a series of laboratory experiments to measure the dynamic angle of repose of glass beads in a rotating drum (tumbler) at various effective gravitational accelerations. The beads had diameter (not radius) mm and occupied about 50% of the volume of a tumbler of radius = 30 mm and length = 5 mm. A tumbler with = 45 mm was also used, but we restrict our comparison to the smaller one for this test. The inner wall of the cylinder was lined with sandpaper to provide a rough surface. A centrifuge was used to provide effective gravitational accelerations between 1 and 25 (where = 9.81 m s is Earthnorm gravity). The angular rotation speed of the tumbler was varied to achieve a range of Froude numbers (“Fr”) for the system, where Fr is defined as the ratio of centripetal to effective gravitational acceleration at the cylinder periphery:^{2}^{2}2Occasionally the square root of the quantity on the righthandside of Eq. (32) is found in the literature, i.e., ; we use the expression in Eq. (32) to compare our results more easily with Brucks et al. (2007).
(32) 
Figure 3 of Brucks et al. (2007) shows as a function of Fr for various obtained in their experiment. Our aim for this test was to compare simulated results for and . Smaller Fr numbers are computationally challenging due to the limited motion of the particles (but the curve is essentially flat for Fr 0.01 anyway). Different , particularly (a regime Brucks et al. 2007 could not explore), will be the subject of future work.
We constructed the simulated experimental apparatus using an infinite horizontal cylinder of radius (30 mm) cut into a short segment by two infinite vertical planes separated by a distance (5 mm). For the purpose of this test, the planes were “frictionless” (). To keep the tests tractable, a fixed particle radius of mm was used, i.e., twice the size of the glass beads used in the experiment (thereby providing nearly a factor of 10 savings in particle number, so each run would take of order a day, instead of more than a week). The initial conditions were prepared in stages, alternately filling the “bottom” of the cylinder (i.e., by removing one plane and having gravity point down the length of the cylinder), closing and rotating the cylinder to allow particles to stick uniformly to the inner cylinder wall (to mimic the roughness provided by the sandpaper in the real experiment), refilling to get the desired volume fraction (about 50%), rotating some more, and then allowing the particles to settle. There were 6513 particles total, of which 568 ended up stuck to the cylinder wall.^{3}^{3}3In order to allow for the finite diameter of the stuck particles, the cylinder was actually = 31.06 mm in radius. This starting condition was used for all the runs. For particle collisions, the normal coefficient of restitution (an estimate for glass beads derived by Heißelmann et al. 2009) and the tangential coefficient of restitution (to provide some roughness; this value is poorly constrained).
As an aside, initially this experiment was conducted with surface coupling (i.e., for the cylinder) instead of sticky walls. Although, for sufficiently fast cylinder spin, the particles also start to tumble in this setup, the particle spins quickly align with the cylinder spin (so that , within a factor of a few) and the particles settle to the bottom of the cylinder, sloshing around slightly. A rough surface instead provides a constant, somewhat random perturbation that ensures the tumbling behavior persists among the loose particles.
Figure Figure Captions shows raytraced snapshots after achieving equilibrium for 8 different Froude numbers (0.001, 0.01, 0.05, 0.1, 0.3, 0.5, 1.0, and 1.5). For the subcritical cases (), the total simulated run time was about 2.5 s, with a fixed timestep of about 5 s. For Fr = 1.0 and 1.5, the run time was 0.5 s with a timestep of 1 s (these latter cases achieve equilibrium faster, but the more rapid cylinder rotation dictated a smaller timestep). Also shown for the subcritical cases is the estimated slope, which was measured by fitting straight lines to sample surface points over several snapshots and taking the average. Only particles to the left of center were considered for these measurements, as particles tend to pile up on the right. As expected, for Fr = 1.0, it can be seen that centrifuging has begun; by Fr = 1.5, a thin layer of free particles persists along the inner surface of the upperright quadrant of the cylinder. For , depends on the rotation rate, with roughly the expected dependence from the experiments. A study by Taberlet et al. (2006) shows that the shape of the granular pile in these kinds of experiments depends on interactions with the end caps, with friction leading to “S” shapes and lack of friction leading to straighter slopes (as in our case). They also found that larger (i.e., using a longer cylinder, for a fixed particle size) leads to straighter slopes. Testing such dependencies in our approach will be the subject of future work.
[FIGURE Figure Captions GOES HERE]
Figure Figure Captions shows the measured values of from the simulations (filled squares), with 1 uncertainty errorbars, along with data taken from Fig. 3 of Brucks et al. (2007) for the = 30 mm tumbler. The experimental data show systematically larger dynamic slopes compared with the simulated data, but both follow a similar trend with Froude number. Some of the discrepancy in slope may arise from the different ratio for the data sets, namely about 110 for the experiment while only 57 for the simulations. Increasing the ratio for the simulations to match the experiments would require a factor of 8 more free particles (as well 4 times more stuck particles), which was beyond the scope of these tests. We note that the = 45 mm tumbler used in the experiments, with , exhibits equilibrium dynamic slopes about 10 deg higher than for = 30 mm (Brucks et al. 2007, Fig. 3).
[FIGURE Figure Captions GOES HERE]
There are other caveats to consider as well. The simulated particles are identical in size and are perfectly spherical; deviations from uniformity of the glass beads used in the experiments may increase the dynamic slopes (this could also give rise to a higher slope at small Froude numbers, where the experimental results did not appear to be as sensitive to ). Furthermore, the adopted and values for the simulations are mostly educated guesses; more experimentation is needed. The simulation used frictionless confining planes, whereas the flat walls in the experiment were comprised of a metal plate and a glass plate, both constrained to rotate with the cylinder, and neither of which were frictionless. In the simulations, the cylinder rotation rate was imposed instantaneously at the start, which may have caused spurious bulking as the assemblage reacted to the sudden shear stress. A better approach (not implemented here) would be to accelerate the tumbler gradually to the desired rotation rate. A related issue is that the dynamic angle of repose may be affected by the initial packing fraction of the particles, which depends on the dominant interparticle forces. In the absence of interparticle friction, random close packing is expected, whereas with friction, random loose packing is a more likely starting condition (Makse et al. 2000), and may lead to a different dynamic angle of repose. This difference was not explored here but will be part of a future investigation. Finally, certain minor run parameters required to keep the simulations tractable in the particlebased, collisionprediction approach may result in excessive particle energy, leading to shallower simulated slopes; experimentation with these thresholds is needed.
4 Discussion
The model atmosphere scenario presented in Section 3.1 is an example of a dilute granular flow in which particleparticle and particlewall interactions are dominated by nearinstantaneous, inelastic, twobody collisions. Our numerical approach, in the context of granular dynamics literature, can be characterized as hardsphere DEM (see Section 1), which is particularly wellsuited to simulating dilute granular flows (Mehta 2007). In the dilute regime, soft and hardparticle methods should give equivalent results, but socalled eventdriven simulations, in which no computations are needed between collision events, are typically much faster (Luding 2004). However, within dense granular flows in which manybody interactions dominate and particles have longlived contacts with many neighbors, it is possible that the assumptions of the eventdriven and hardsphere DEM models fail (Delannay et al. 2007).
In a rotating tumbler (cf. Section 3.2), there exists a region of longlived contacts but also a faster flowing layer. By varying the Froude number, it is possible to investigate many different regimes, from the fastflowing regime with a large Froude number to the regime when the Froude number is low and most contacts are longerlived. Thus the tumbler is an intermediate regime between dilute and dense flows, providing a convenient test of the range of applicability of our hardsphere method. Indeed, the tumbler has been studied in detail with many different numerical techniques; we provide a brief overview in the following.
McCarthy et al. (1996, 2000) and McCarthy and Ottino (1998) used a softsphere model to compare simulations of noncohesive granular materials in a slowly rotating drum to experimental data. In order to limit the number of particles necessary in the simulations they proposed a “hybrid” simulation technique. This method involves using a geometrical model to determine the bulk motion of the particles and then performing particle dynamics simulations on only the particles contained within and bordering the avalanching wedge. They found a favorable match between experiment and simulation, at low computational cost, but with the tradeoff that the method is tailored to a very specific geometry.
Khakhar et al. (1997) used a similar technique to separate the granular material in the rotating drum into a “rapid flow region” and a “fixed bed.” A continuum model in which averages are taken across the layer was used to analyze behavior of the flowing layer. The motion of grains on the free surface of a granular mixture in a rotating drum was also investigated by Monetti et al. (2001) using Monte Carlo simulations of a twodimensional lattice gas model. The model takes into account rotational frequency, frictional forces, and the gravitational field.
Softsphere methods have also repeatedly and successfully been used to describe granular material in a rotating drum (Dury et al. 1998; Dury and Ristow 1999; Rapaport 2002; Pohlman et al. 2006; Taberlet et al. 2006). This has also extended been to twodimensional simulations of irregularshaped particles by Poschel and Buchholtz (1995).
Hardsphere methods have not been commonly used for rotating drums, however Gui and Fan (2009) did perform numerical simulations of motion of rigid spherical particles within a 2D tumbler with an inner wavelike surface. The rotation of the tumbler was simulated as a traveling sine wave around a circle.
The observed behavior in our tumbler simulations shows that the code correctly models the transitions from global regimes with increasing Froude number. The differences between experiment and simulation are larger at lower Froude number and start to decrease at higher Froude number. This suggests that our hardsphere model is more suited to the dilute flow regimes of high Froude number rather than dense flow regimes of low Froude number. Softsphere DEM is typically used for dense granular flows (Zhu et al. 2007), as the detailed characterization of interparticle contacts is more suitable for this regime. However, hardsphere methods can approach experimental and softsphere simulation results for very low friction materials. For example, Pohlman et al. (2006) found an angle of repose of only 16.9, for low friction () chrome steel beads at a Froude number of . This is within a few degrees of our expected results at this Froude number.
In summary, our approach currently favors dilute or fluid flows over dense granular flows, but we are in the process of implementing a softsphere method. Whether regolith in lowgravity environments, such as the surface of an asteroid, when reacting to external stresses, is more appropriately modeled as a dilute or dense granular flow is the subject of ongoing investigation. Regardless, numerical methods will play an important role in the study of these environments.
5 Conclusions
We have added new features to an existing, welltested body code that allow convenient modeling of granular dynamics in a variety of conditions. Details of collision detection and resolution involving four principal supported wall geometries (infinite plane, finite disk, infinite cylinder, and finite cylinder) and associated degenerate cases (point, line, and ring) were provided. Two dynamical test suites were presented, one for a model “atmosphere” showing the correct equilibrium distribution of scale heights among a polydisperse population of dissipationless balls, and the other for a tumbler that showed correct qualititative behavior as a function of surface wall properties and cylinder rotation rate. Application of these new capabilities to numerical simulations designed to investigate collective grain behavior on vibrating plates will be presented in Murdoch et al. (2011, in preparation). Planned future code development includes adding support for rigid and semirigid particle aggregates (cf. Richardson et al. 2009) that would allow study of more complex particle shapes in simulations of granular dynamics, and implementing a softsphere model to better match the properties of dense granular flows.
Acknowledgments
This material is based upon work supported by the National Aeronautics and Space Administration under Grant No. NNX08AM39G issued through the Office of Space Science and by the National Science Foundation under Grant No. AST0524875. KJW was supported by a Poincaré Fellowship at Observatoire de la Côte d’Azur (OCA). PM and NM acknowledge financial support from the French Programme National de Planétologie. NM acknowledges financial support from Thales Alenia Space and The Open University. DCR thanks M. C. Miller for helpful discussions. We thank Dan Scheeres and an anonymous reviewer for their comments and suggested improvements. Some simulations in this work were performed on the CRIMSON Beowulf cluster at OCA. Raytracing for Figs. Figure Captions and Figure Captions was carried out using the Persistence of Vision Raytracer.^{4}^{4}4http://www.povray.org/ Data points from Fig. 3 of Brucks et al. (2007) were extracted using Plot Digitizer.^{5}^{5}5http://plotdigitizer.sourceforge.net/
References

Asphaug, E., King, P.J., Swift, M.R., Merrifield, M.R., 2001. Brazil nuts on Eros: Sizesorting of asteroid regolith. Proc. Lunar Sci. Conf. 32, 1708.

Brucks, A., Arndt, T., Ottino, J., Lueptow, R., 2007. Behavior of flowing granular materials under variable . Phys. Rev. E 75, 0323011–0323014.

Cleary, P.W., Sawley, M.L., 2002. DEM modelling of industrial granular flows: 3D case studies and the effect of particle shape on hopper discharge. Applied Mathematical Modelling 26, 89–111.

Delannay, R., Louge, M., Richard, P., Taberlet, N., Valance, A., 2007. Towards a theoretical picture of dense granular flows down inclines. Nature Materials 6, 99–108.

Dury, C.M., Ristow, G.H., 1999. Competition of mixing and segregation in rotating cylinders. Phys. Fluids 11, 1387–1394.

Dury, C.M., Ristow, G.H., Moss, J.L., Nakagawa, M., 1998. Boundary effects on the angle of repose in rotating cylinders. Phys. Rev. E 57, 4491–4497.

Elaskar, S.A., Godoy, L.A., Gray, D.D., Stiles, J.M., 2000. A viscoplastic approach to model the flow of granular solids. International Journal of Solids and Structures 37, 2185–2214.

Falk, M.L., Toiya, M., Losert, W., 2008. Shear transformation zone analysis of shear reversal during granular flow. ArXiv eprints 0802.1752.

Fraige, F.Y., Langston, P.A., Chen, G.Z., 2008. Distinct element modelling of cubic particle packing and flow. Powder Technology 186, 224–240.

Goddard, J.D., 1990. Nonlinear elasticity and pressuredependent wave speeds in granular media. Proc. Royal Soc. London A430, 105–131.

Gui, N., Fan, J., 2009. Numerical simulation of motion of rigid spherical particles in a rotating tumbler with an inner wavelike surface. Powder Technology 192, 234–241.

Heißelmann, D., Blum, J., Fraser, H., Wolling, K., 2009. Microgravity experiments on the collisional behavior of saturnian ring particles. Icarus 206, 424–430.

Holsapple, K.A., 2001. Equilibrium configurations of solid cohesionless bodies. Icarus 154, 432–448.

Holsapple, K.A., 2004. Equilibrium figures of spinning bodies with selfgravity. Icarus 172, 272–303.

Holsapple, K.A., Michel, P., 2006. Tidal disruptions: A continuum theory for solid bodies. Icarus 183, 331–348.

Holsapple, K.A., Michel. P., 2008. Tidal disruptions. II. A continuum theory for solid bodies with strength, with applications to the solar system. Icarus 193, 283–301.

Hong, D.C., McLennan, J.A., 1992. Molecular dynamics simulations of hard sphere granular particles. Physica A: Statistical Mechanics and its Applications 187, 159–171.

Huilin, L., Yunhua, Z., Ding, J., Gidaspow, D., Wei, L., 2007. Investigation of mixing/segregation of mixture particles in gassolid fluidized beds. Chem. Engin. Sci. 62, 301–317.

Kacianauskas, R., Maknickas, A., Kaceniauskas, A., Markauskas, D., Balevicius, R., 2010. Parallel discrete element simulation of polydispersed granular material. Advances in Engineering Software 41, 52–63.

Kamrin, K., Bazant, M.Z., 2007. Stochastic flow rule for granular materials. Phys. Rev. E 75, 041301.

Kharhar, D.V., McCarthy, J.J., Shinbrot, T., Ottino, J.M., 1997. Transverse flow and mixing of granular material in a rotating cylinder. Phys. Fluids 9, 31–43.

Kosinski, P., Hoffmann, A.C., 2009. Extension of the hardsphere particlewall collision model to account for particle deposition. Phys. Rev. E 79, 061302.

Latham, JP., Munjiza, A., Garcia, X., Xiang, J., Guises, R., 2007. Threedimensional particle shape acquisition and use of shape library for DEM and FEM/DEM simulation. Minerals Engineering 21, 797–805.

Losert, W., Bocquet, L., Lubensky, T. C., Gollub, J.P., 2000a. Particle dynamics in sheared granular matter. Phys. Rev. Lett. 85, 1428–1431.

Losert, W., Géminard, J.C., Nasuno, S., Gollub, J.P., 2000b. Mechanisms for slow strengthening in granular materials. Phys. Rev. E 61, 4060–4068.

Luding, S., 2004. Molecular dynamics simulations of granular materials. In: Hinrichsen, H., Wolf, D.E. (Eds.), The Physics of Granular Media. WileyVCH, Weinheim, pp. 299–323.REVISE CHECK END PAGE

Makse, H.A., Johnson, D.L., Schwartz, L.M., 2000. Packing of compressible granular materials. Phys. Rev. Let. 84, 4160–4163.

McCarthy, J.J., Ottino, J.M., 1998. Particle dynamics simulation: A hybrid technique applied to granular mixing. Powder Technology 97, 91–99.

McCarthy, J.J., Shinbrot, T., Metcalfe, G., Wolf, E.J., Ottino, J.M., 1996. Mixing of granular materials in slowly rotated containers. AIChE Journal 42, 3351–3363.

McCarthy, J.J., Khakhar, D.V., Ottino, J.M., 2000. Computational studies of granular mixing. Powder Technology 109, 72–82.

Mehta, A.J., 2007. Granular Physics. Cambridge Univ. Press.New York

Michel, P., O’Brien, D., Abe, S., Hirata, N., 2009. Itokawa’s cratering record as observed by Hayabusa: Implications for its age and collisional history. Icarus 200, 503–513.

Miwa, Y., Yano, H., Morimoto, M., Mori, O., Kawaguchi, J., 2008. A study on constitution of asteroids with Brazilnut effect. Proc. 26th International Symposium on Space Technology and Science 2008, 2008k03.

Miyamoto, H. and 14 colleagues, 2007. Regolith migration and sorting on asteroid Itokawa. Science 316, 1011–1014.

Monetti, R., Hurd, A., Kenkre, V.M., 2001. Simulations for dynamics of granular mixtures in a rotating drum. Granular Matter 3, 113–116.

Mueth, D.M., Debregeas, G.F., Karczmar, G.S., Eng, P.J., Nagel, S.R., Jaeger, H.M., 2000. Signatures of granular microstructure in dense shear flows. Nature 406, 385–389.

Murdoch, N., Berardi, C., Michel, P., Richardson, D.C., Wolfgang, L., Green, S.F., 2011. Numerical simulations of granular material dynamics. II. Comparison with shaking experiments. Icarus, in preparation.

Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P., 2007. Numerical Recipes: The Art of Scientific Computing, third ed.. Cambridge Univ. Press, New York.

Pohlman, N.A., Severson, B.L., Ottino, J.M., Lueptow, R.M., 2006. Surface roughness effects in granular matter: Influence on angle of repose and the absence of segregation. Phys. Rev. E 73, 031304.

Poschel, T., Buchholtz, V., 1995. Complex flow of granular material in a rotating cylinder. Chaos, Solitons & Fractals 5, 1901–1905, 19071912.

Rapaport, D.C., 2002. Simulational studies of axial granular segregation in a rotating cylinder. Phys. Rev. E 65, 061306.

Richardson, D.C., 1994. Tree code simulations of planetary rings. Mon. Not. R. Astron. Soc. 269, 493–511.

Richardson, D.C., Quinn, T., Stadel, J., Lake, G., 2000. Direct largescale body simulations of planetesimal dynamics. Icarus 143, 45–59.

Richardson, D.C., Michel, P., Walsh, K.J., Flynn, K.W., 2009. Numerical simulations of asteroids modeled as gravitational aggregates. Plan. & Space Sci. 57, 183–192.

Richardson, J.E., Melosh, H.J., Greenberg, R., 2004. Impactinduced seismic activity on asteroid 433 Eros: A surface modification process. Science 306, 1526–1529.

Robinson, M.S., Thomas, P.C., Veverka, J., Murchie, S.L., Wilcox, B.B., 2002. The geology of 433 Eros. Meteorit. Planet. Sci. 37, 1651–1684.

Rosato, A., Strandburg, K.J., Prinz, F., Swendsen, R.H., 1987. Why the Brazil nuts are on top: Size segregation of particulate matter by shaking. Phys. Rev. Lett. 58, 1038–1040.

Sánchez, P., Scheeres, D.J., Swift, M.R., 2010. Impact driven size sorting in selfgravitating granular aggregates. Lunar Planet. Sci. 41. Abstract 1533.

Savage, S.B., 1998. Analyses of slow highconcentration flows of granular materials. J. Fluid Mech. 377, 1–26.

Scheeres, D.J., Abe, M., Yoshikawa, M., Nakamura, R., Gaskell, R.W., Abell, P.A., 2007. The effect of YORP on Itokawa. Icarus 188, 425–429.

Scheeres, D.J., Hartzell, C.M., Sánchez, P., Swift, M., 2010. Scaling forces to asteroid surfaces: The role of cohesion. Icarus 210, 968–984.

Sharma, I., Jenkins, J.T., Burns, J.A., 2006. Tidal encounters of ellipsoidal granular asteroids with planets. Icarus 183, 312–330.

Sharma, I., Jenkins, J.T., Burns, J.A., 2009. Dynamical passage to approximate equilibrium shapes for spinning, gravitating rubble asteroids. Icarus 200, 304–322.

Stadel, J., 2001. Cosmological body simulations and their analysis. Thesis, University of Washington, Seattle. 126 pp.

Szarf, K., Combe, G., Villard, P., 2010. Polygons vs. clumps of discs: A numerical study of the influence of grain shape on the mechanical behaviour of granular materials. Powder Technology, in press.

Taberlet, N., Richard, P., Hinch, J.E., 2006. S shape of a granular pile in a rotating drum. Phys. Rev. E 73, 050301(R).

Toiya, M., Stambaugh, J., Losert, W., 2004. Transient and oscillatory granular shear flow. Phys. Rev. Lett. 93, 088001.

Wada, K., Senshu, H,, and Matsui, T., 2006. Numerical simulation of impact cratering on granular material. Icarus 180, 528–545.

Wagner, M. 2004. Ray/torus intersection. Independent online publication.^{6}^{6}6http://www.emeyex.com/site/projects/raytorus.pdf

Walsh, K.J., Richardson, D.C., Michel, P., 2008. Rotational breakup as the origin of small binary asteroids. Nature 454, 188–19.

Zhu, H.P., Zhou, Z.Y., Yang, R.Y., Yu, A.B., 2007. Discrete particle simulation of particulate systems: Theoretical developments. Chem. Engin. Sci. 62, 3378–3396.
Figure Captions

Figure Figure Captions: Diagram illustrating the quantities used for predicting the collision of two spheres (cf. Eqs. (1)–(3)). Although the diagram depicts a collision in two dimensions, the derivations apply to the full threedimensional case.

Figure Figure Captions: Geometry of sphereplane collision. Here the plane is depicted as a finite, thick disk for illustration only; in reality the plane has zero thickness and infinite extent. For simplicity, velocity vectors are not shown. The unit vector is in the plane, while is perpendicular to it (for clarity of presentation, the vector representing is given length , i.e., the distance from the plane origin to the contact point). Note and , where in this scenario must equal , the radius of the impacting sphere. The geometry is similar for a spheredisk collision (Section 2.2.2), where the disk in the figure would have radius but still have no thickness.

Figure Figure Captions: Geometry of spherecylinder collision. Shown is a finite cylinder of radius and length , with origin at the geometric center. For an infinite cylinder, and can be anywhere along the symmetry axis. For ease of illustration, only the case of exterior impact is shown; interior impact is similar, except has opposite sign. The unit vector is the perpendicular pointing from the cylinder axis to the point of contact, while points along the cylinder axis (for clarity of presentation, the vectors representing and are given length and , respectively, where is the radius of the impacting sphere). Note (for exterior impact) and , where in this scenario must equal ; is simply .

Figure Figure Captions: Snapshots from the model atmosphere test (Section 3.1) showing (a) the initial condition, (b) splash, and (c) final equilibrated state. The particle colors red, magenta, and blue correspond to particle masses 1, 3, and 10, respectively (particle radii are the same for all masses). The scale heights of the 3 populations are expected to be inversely proportional to their masses. The full vertical extent of the smallestmass particle distribution is not shown. Gravity is directed down in the figure.

Figure Figure Captions: Evolution of the mean particle heights in the model atmosphere test. The colors red (top line), magenta (middle line), and blue (bottom line) correspond to particle masses 1, 3, and 10, respectively. The dotted horizontal lines are the corresponding expected scale heights assuming hydrostatic equilibrium, the ideal gas law, and energy equipartition. The timestep for these simulations was in system units. The snapshots of Fig. Figure Captions correspond to timesteps 0, 30, and 1,000, respectively.

Figure Figure Captions: Cumulative probability distributions for the height data (solid lines) and model (dotted lines) corresponding to the final snapshot in Fig. Figure Captions. The colors red (middle lines at small ), magenta (bottom lines at small ), and blue (top lines) correspond to particle masses 1, 3, and 10, respectively. The KS statistic measures the maximum vertical deviation of the data from the model (see text).

Figure Figure Captions: Snapshots of the simulated tumbler experiments taken at the end of each run. The Froude number (Eq. (32)) is indicated. View size is about 7.5 cm across each frame. The view is down the cylinder axis, and rotation is in the clockwise rotation direction (at this viewing angle, the short cylinder itself cannot be seen). The blue particles along the inner surface of the cylinder are “glued” in place to provide roughness. Where shown, the yellow line indicates an estimate of the surface slope (also see Fig. Figure Captions).

Figure Figure Captions: Plot of dynamic repose angle as a function of Froude number for the simulated tumblers (filled squares) and from the experiments of Brucks et al. (2007; open squares connected with dashed line, from their Fig. 3). Errorbars for the simulated points are 1 uncertainties from averaging slope estimates over several snapshots. The simulated points are shifted down and to the right compared with the experiments, which may be due partly to the difference in particletotumbler radius ratio. See main text for discussion.
Figure 1
Figure 2
Figure 3
Figure 4
Figure 5
Figure 6
Figure 7
Figure 8