John E. Lane

ASRC Aerospace, Kennedy Space Center, Florida

Philip T. Metzger

Kennedy Space Center, Kennedy Space Center, Florida

R. Allen Wilkinson

Glenn Research Center, Cleveland, Ohio

A Review of Discrete Element Method (DEM) 

Particle Shapes and Size Distributions for 

Lunar Soil

NASA/TM—2010-216257

December 2010
NASA STI Program . . . in Profi le

Since its founding, NASA has been dedicated to the 

advancement of aeronautics and space science. The 

NASA Scientifi c and Technical Information (STI) 

program plays a key part in helping NASA maintain 

this important role.

The NASA STI Program operates under the auspices 

of the Agency Chief Information Offi cer. It collects, 

organizes, provides for archiving, and disseminates 

NASA’s STI. The NASA STI program provides access 

to the NASA Aeronautics and Space Database and 

its public interface, the NASA Technical Reports 

Server, thus providing one of the largest collections 

of aeronautical and space science STI in the world. 

Results are published in both non-NASA channels 

and by NASA in the NASA STI Report Series, which 

includes the following report types:

• TECHNICAL PUBLICATION. Reports of 

completed research or a major signifi cant phase 

of research that present the results of NASA 

programs and include extensive data or theoretical 

analysis. Includes compilations of signifi cant 

scientifi c and technical data and information 

deemed to be of continuing reference value. 

NASA counterpart of peer-reviewed formal 

professional papers but has less stringent 

limitations on manuscript length and extent of 

graphic presentations.

• TECHNICAL MEMORANDUM. Scientifi c 

and technical fi ndings that are preliminary or 

of specialized interest, e.g., quick release 

reports, working papers, and bibliographies that 

contain minimal annotation. Does not contain 

extensive analysis.

• CONTRACTOR REPORT. Scientifi c and 

technical fi ndings by NASA-sponsored 

contractors and grantees.

• CONFERENCE PUBLICATION. Collected 

papers from scientifi c and technical 

conferences, symposia, seminars, or other 

meetings sponsored or cosponsored by NASA.

• SPECIAL PUBLICATION. Scientifi c, 

technical, or historical information from 

NASA programs, projects, and missions, often 

concerned with subjects having substantial 

public interest.

• TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientifi c and 

technical material pertinent to NASA’s mission.

Specialized services also include creating custom 

thesauri, building customized databases, organizing 

and publishing research results.

For more information about the NASA STI 

program, see the following:

• Access the NASA STI program home page at 

http://www.sti.nasa.gov

• E-mail your question via the Internet to help@

sti.nasa.gov

• Fax your question to the NASA STI Help Desk 

at 443–757–5803

• Telephone the NASA STI Help Desk at

 443–757–5802

• Write to:

 NASA Center for AeroSpace Information (CASI)

 7115 Standard Drive

 Hanover, MD 21076–1320
John E. Lane

ASRC Aerospace, Kennedy Space Center, Florida

Philip T. Metzger

Kennedy Space Center, Kennedy Space Center, Florida

R. Allen Wilkinson

Glenn Research Center, Cleveland, Ohio

A Review of Discrete Element Method (DEM) 

Particle Shapes and Size Distributions for 

Lunar Soil

NASA/TM—2010-216257

December 2010

National Aeronautics and

Space Administration

Glenn Research Center

Cleveland, Ohio 44135
Available from

NASA Center for Aerospace Information

7115 Standard Drive

Hanover, MD 21076–1320

National Technical Information Service

5301 Shawnee Road

Alexandria, VA 22312

Available electronically at http://gltrs.grc.nasa.gov

Trade names and trademarks are used in this report for identifi cation 

only. Their usage does not constitute an offi cial endorsement, 

either expressed or implied, by the National Aeronautics and 

Space Administration.

Level of Review: This material has been technically reviewed by technical management. 
NASA/TM—2010-216257 1 

A Review of Discrete Element Method (DEM) Particle Shapes and 

Size Distributions for Lunar Soil 

John E. Lane 

ASRC Aerospace 

Kennedy Space Center, Florida 32899 

Philip T. Metzger 

National Aeronautics and Space Administration 

Kennedy Space Center, Florida 32899 

R. Allen Wilkinson 

National Aeronautics and Space Administration 

Glenn Research Center 

Cleveland, Ohio 44135 

Summary 

This report reviews two topics that are important to discrete element method (DEM) modeling the 

behavior of soils (such as lunar soils): (1) methods of modeling particle shapes and (2) analytical 

representations of particle size distribution. These two topics are somewhat unrelated, but both are 

reviewed in this report as a part of ongoing efforts to develop models of lunar soil mechanics. The choice 

of particle shape complexity is driven primarily by opposing tradeoffs with total number of particles, 

computer memory, and total simulation computer processing time. Also, the choice is dependent on 

available DEM software capabilities. For example, PFC2D/PFC3D and EDEM support clustering of 

spheres; MIMES incorporates superquadric particle shapes; and BLOKS3D provides polyhedra shapes. 

Most commercial and custom DEM software supports some type of complex particle shape beyond the 

standard sphere. Convex polyhedra, clusters of spheres and single parametric particle shapes such as the 

ellipsoid, poly-ellipsoid, and superquadric, are all motivated by the desire to introduce asymmetry into the 

particle shape, as well as edges and corners, in order to better simulate actual granular particle shapes and 

behavior. 

An empirical particle size distribution (PSD) formula is shown to fit the example desert sand data 

from Bagnold, using a coarse-grade exponent c and small-grade exponent s. An additional empirical 

parameter γ has been introduced to connect the two (log-log) linear fits in the small- and coarse-grade 

particle size regions. Note that the PSD has been defined so that the integral from size 0 to ∞ is equal to 

100 in units of percent. Particle size data of JSC-1a obtained from a fine particle analyzer (FPA) at NASA 

Kennedy Space Center is also fitted to a similar empirical PSD function in order to extract the equivalent 

small- and coarse-grade exponents. 

1.0 Introduction 

The starting point with any discrete element method (DEM) software is to model grains as spheres of 

a specified diameter range. The spherical particle shape greatly simplifies simulations and accommodates 

the maximum number of particles for any given central processing unit (CPU) execution time budget. 

However, it is well known that simulations of most real granular assemblies, such as sand or lunar 

regolith, are not accurately modeled by spheres because of the rolling of contact surfaces versus the actual 

frictional sliding of real particle surfaces. Most researchers active in granular modeling agree that a 

particle shape more complex than the sphere is required to generate simulations that can match the 

observable properties of real granular systems (Refs. 1 to 4). Section 2.0 of this report addresses particle 

shape. 
NASA/TM—2010-216257 2 

The particle size distribution is one of the primary factors affecting how a soil behaves mechanically. 

In order to model soil mechanics, it is important to start out with an understanding of the size distribution 

of the soil that is being modeled. Section 3.0 of this report examines analytical models that can 

approximate a lunar soil particle size distribution. Appendix A shows magnified images and elemental 

spectra of JSC-1a lunar regolith stimulant, developed by the NASA Johnson Space Center. It should be 

noted that particle size distributions have been measured by several methods, which do not always agree 

with one another because each method is measuring a different aspect of particle “size.” Therefore, these 

analytical functions should not necessarily be considered definitive for lunar soil or lunar simulants until 

further work is accomplished. Instead, these should be considered as tentative “models” of the lunar soil 

particle size distribution that will be refined with further research. 

2.0 Methods of Modeling Particle Shapes 

Here, Subsections 2.1 to 2.7 cover different particle constructs. Subsection 2.8 describes a contact 

detection scheme, given nonspherical particles. Appendix B presents a summary of particle shape types 

associated with various DEM software packages, comparing some of the pros and cons of each of these 

idealized particle characteristics. 

2.1 Clusters of Spheres 

Popular DEM software, such as that developed and marketed by the companies Itasca and DEM 

Solutions Ltd., make use of sphere clusters (Fig. 1). Normally, the spheres making up a sphere cluster 

form a rigid body where the adhesion force between spheres is set to infinity. In some cases, researchers 

prefer to connect the spheres using finite adhesion forces in order to model crushing and breaking of large 

complex granular particles. 

Another variation of particles composed of sphere clusters uses overlapping spheres (Fig. 2) (Refs. 7 

to 9). Even though this is the most general case of sphere clusters, it is not necessarily the most common 

approach because of complications caused by its implementation in most of the commercial DEM 

software. DEM software typically models particles using a constant mass density. If particles of equal 

density overlap, the density doubles. Therefore, it is often necessary to implement a careful scaling of 

mass density of each of the component spheres so that the final cluster density uniformity is optimized. 

An algorithm to adjust mass density in two dimensions is described by Ashmawy (Ref. 8). 

Figure 1.—Particle shape defined by cluster of 

nonoverlapping spheres. 

Figure 2.—Particle shape defined by cluster of 

overlapping spheres. 
NASA/TM—2010-216257 3 

The equation of the surface a single sphere of the cluster contains in general, four parameters, the 

radius Rk and the vector offset rk = (xk , yk , zk ) relative to a local particle origin: 

1

2 2 2

 = 

 



 −  + 

 



 −  + 

 



 −

k

k

k

k

k

k

R

z z

R

y y

R

x x (1) 

For a cluster of N spheres, there will be 1 + 4(N – 1) parameters since one of the spheres is defined to be 

centered at the origin of the local particle coordinate system. If the specific DEM software requires a mass 

density scaling of the component spheres, an additional parameter per sphere is needed, for a total of 

5N – 3 parameters. If however, a density adjustment algorithm is utilized, such as that proposed by 

Ashmawy, only the final density of the cluster is required for a total number of parameters equal to 4N – 2. 

Equation (1) is valid for both overlapping and nonoverlapping clusters. A more useful representation 

of the surface equation is given in spherical coordinates as 

( , ) sin cos

( , ) sin sin

( , ) cos for 0 2 , 0

k k

k k

k

xR x

yR y

z R z 

θ φ = θ φ+

θ φ = θ φ+

θ φ = θ+ ≤φ< π ≤θ<π

 (2) 

2.2 Ellipsoid 

Other popular particle types are the two-dimensional ellipse and three-dimensional ellipsoid (Fig. 3) 

(Ref. 10). According to many researchers, the asymmetry provided by the ellipsoid yields the best 

tradeoff between particle complexity and total simulation execution time. The equation for the ellipsoid 

surface is similar to Equation (1): 

1

2 2 2

 = 

 



  +  



  +  





c

z

b

y

a

x (3) 

where a, b, and c are the half-widths along each of the principal axes. The parametric form of 

Equation (3) is 

( , ) sin cos

( , ) sin sin

( , ) cos for 0 2 , 0

x a

y b

z c 

θφ = θ φ

θφ = θ φ

θφ = θ ≤φ< π ≤θ<π

 (4) 

Figure 3.—Particle shape defined by an ellipsoid 

(Eq. (3)), with a = 1, b = 0.5, and c = 2. 
NASA/TM—2010-216257 4 

Neither the ellipse or ellipsoid lends itself to clustering strategies, since the attraction to clustering is 

based on a simple unit particle—the sphere. If the unit particle is an ellipse, the difficulties that arise from 

the ellipse tend to negate the advantages of clustering. 

The ellipsoid and sphere are subsets of the superquadric shapes (Ref. 11). 

2.3 Superquadric 

The equation for the superquadric surface is similar to Equation (3): 

+ + =1

K L M

c

z

b

y

a

x (5) 

where again a, b, and c are the half-widths along each of the principal axes (Fig. 4). The parametric form 

of Equation (3) is: 

( )

( )

( )

≤ φ < π ≤ θ < π

θ φ = θ θ

θ φ = θ φ θ φ

θ φ = θ φ θ φ

 for 0 2 , 0

( , ) cos sgn cos

( , ) sin sin sgn sin sin

( , ) sin cos sgn sin cos

2 /

2 /

2 /

M

L

K

z c

y b

x a

 (6) 

where sgn(x) is the sign of the argument x. The exponent parameters K, L, and M adjust the blockiness of 

the shape along the three principal axis. If the exponent is less than 2, the shape becomes concave along 

that axis; it becomes convex for an exponent greater than or equal to 2. The superquadrics provides an 

additional degree of asymmetry over the ellipsoid. The greatest degree of asymmetry is achieved by 

defining multiple superquadric parameters in different regions. This strategy leads to the asymmetric 

superquadric. 

Figure 4.—Superquadric particle shape defined 

by Equation (5) where a = 1, b = 1, c = 1, 

K = 0.6, L = 4, and M = 2. 
NASA/TM—2010-216257 5 

2.4 Asymmetric Superquadric 

When Equation (6) is used with P sets of parameters in P regions, an arbitrarily asymmetric shape can 

be generated (Fig. 5): 

( )

( )

( )

( )

( )

( ) 













θ θ φ ≤ φ < π θ ≤ θ < π

θ θ ≤ φ < φ ≤ θ < θ

θ φ =















θ φ θ φ φ ≤ φ < π θ ≤ θ < π

θ φ θ φ ≤ φ < φ ≤ θ < θ

θ φ =















θ φ θ φ φ ≤ φ < π θ ≤ θ < π

θ φ θ φ ≤ φ < φ ≤ θ < θ

θ φ =

− −

− −

− −

1 1

1 1 1

1 1

1 1 1

1 1

1 1 1

cos sgn cos for 2 , 

cos sgn cos for 0 , 0

( , )

sin sin sgn sin sin for 2 , 

sin sin sgn sin sin for 0 , 0

( , )

sin cos sgn sin cos for 2 , 

sin cos sgn sin cos for 0 , 0

( , )

2 /

2 / 1

2 /

2 / 1

2 /

2 / 1

P P P

P P P

P P P

M P

M

LP

L

KP

K

c

c

z

b

b

y

a

a

x

 

 

 

(7) 

For example, Figure 5 is generated by defining P = 2 regions: 

Region 1: { } ≤ φ < π ≤ θ < π2

1 0 2 , 0

Region 2: { } ≤ φ < π π ≤ θ < π 2

1 0 2 , 

In the general case, the slopes are not continuous and are characterized by a sharp surface ridge. 

Figure 5.—Asymmetric superquadric particle shape 

defined by Equation (7) where a1 = a2 = 0.7, 

b1 = b2 = 1, c1 = 0.75, c2 = 2, K1 = K2 = 4, 

L1 = L2 = 2, M1 = 4, and M2 = 0.7. 
NASA/TM—2010-216257 6 

2.5 Eight-Quadrant Superquadric 

Defining the eight quadrants of the asymmetric superquadric as three-dimensional regions generates a 

useful particle shape definition: 

( )

( )

( )









θ ≤ θ < π

θ ≤ θ <

θ φ =







≤ θ < π

θ φ π ≤ φ < π

θ φ ≤ φ < π θ φ =

≤ θ < π 







θ φ ≤ φ <

θ φ − ≤ φ <

θ φ =

π

π

π π

π π

2

2 /

2

2

2 /

1

2 /

2

2 /

1

2

3

2

2 /

2

2 2

2 /

1

cos for 

cos for 0

( , )

 0

sin sin for 2

sin sin for 0 ( , )

 0

sin cos for 

sin cos for 

( , )

2

1

2

1

2

1

M

M

L

L

K

K

c

c z

b

b

y

a

a

x

 (8) 

The description of the eight-quadrant superquadric form Equation (8) is implemented in Mathematica 

(Ref. 12) and is shown in Figure 6, as well as in Appendix C where the eight-quadrant superquadric shape 

parameters are randomly modified to demonstrate the corresponding range of particle shapes. 

Figure 6.—Eight-quadrant superquadric video sequence using Equation (8). 
NASA/TM—2010-216257 7 

Figure 7.—Poly-ellipsoid particle shape defined 

by Equation (9) where a1 = 2, a2 = 0.7, b1 = 0.5, 

b2 = 1, c1 = 0.75, c2 = 2, and K1 = K2 = L1 = L2

= M1 = M2 = 2. 

2.6 Poly-Ellipsoid 

Poly-ellipsoids, as defined by Peters (Ref. 13), can be viewed as a special case of the eight-quadrant 

superquadric discussed above. In this case (Fig. 7), there are again eight distinct regions corresponding to 

the eight spatial quadrants: 







φ ≤ θ < π

φ ≤ θ <

θ φ =







θ φ π ≤ φ < π

θ φ ≤ φ < π

θ φ =







θ φ ≤ φ <

θ φ − ≤ φ <

θ φ =

π

π

π π

π π

2 2

2 1

2

1

2

3

2 2

2 2 1

cos for 

cos for 0

( , )

sin sin for 2

sin sin for 0

( , )

sin cos for 

sin cos for 

( , )

c

c

z

b

b

y

a

a

x

 (9) 

where, from Equation (8), L1 = L2 = M1 = M2 = N1 = N2 are set equal to 2 only in the special case of the 

poly-ellipsoid. This special case defines the poly-ellipsoid and guarantees that the slopes at all quadrant 

boundaries are continuous. 

Care must be taken in defining the parameter sets so that the surface is continuous and does not 

contain voids. A logical method to ensure surface continuity is accomplished by use of the poly-ellipsoid 

particle shape. 

2.7 Polyhedra 

All particle shapes up to this point have been described by closed-form expressions, which is a great 

convenience for computing first order (mass and volume) and second-order (moments of inertia) particle 

parameters of interest for DEM simulations, as well as computing particle interaction (contact forces). All 

simulated particle shapes from compound sphere clusters to poly-ellipses have one goal in common: to 

include shape asymmetry and sharper edges and corners for the purpose of modeling more realistic 

particle shapes. 
NASA/TM—2010-216257 8 

A convex polyhedron can be defined by V vertices, E edges, and F faces, described conveniently by 

the Euler characteristic 

V − E + F = 2 (10) 

The most convenient and compact method to define a general convex polyhedron is by a linear 

system of inequalities, which can be expressed in matrix form as follows: 

A⋅r ≤ b (11) 

where r is an N × 1 position vector in N-dimensional space; A is an F × N matrix (F is the number of 

polyhedron faces); and b is an N × 1 vector. 

In three-dimensional space ( )T

r = x y z , where T denotes “transpose.” The inequality of 

Equation (11) discriminates between points in the polyhedron from those outside of the surface, which are 

not part of the polyhedron. Replacing the inequality with an equality for each of the ith rows of A, the 

plane of the ith face of the polyhedron is then expressed by 

ai1 x + ai2 y + ai3 z − bi = 0 (12) 

In this way, the region bounded by the intersection of all F infinite planes defined by Equation (11) 

defines the interior region and surface of the polyhedron. Some specific examples are shown below in 

Equations (13) to (16) and Figures 8 to 11. 

Figure 8.—Four-faced polyhedron (tetrahedron) 

particle shape. 





























=





























− −

− −

− − =

1

1

1

1

1 1 1

1 1 1

1 1 1

1 1 1

b

A

 (13)

 
NASA/TM—2010-216257 9 

Figure 9.—Five-faced polyhedron (pyramid) particle shape. 

































=

































−

−

−

=

2

1

2

1

2

1

2

1

2

1

1

1

1

1

0 0 1

0 1

0 1

1 0

1 0

b

A

 (14)

Figure 10.—Six-faced polyhedron (cube) particle shape. 









































=









































−

−

−

=

1

1

1

1

1

1

0 0 1

0 0 1

0 1 0

0 1 0

1 0 0

1 0 0

b

A

 (15)
NASA/TM—2010-216257 10 

Figure 11.—Eight-faced polyhedron (octahedron) 

particle shape. 

111

11 1

1 11

1 11

11 1

11 1

1 11

111

1

1

1

1

1

1

1

1

    −  

  −   − −   =   −

 

  − −

  − −  

  −−−

   

 

 

 

  =  

 

 

 

     

A

b

 (16)

Simulated particle geometries based on polyhedra are a radical departure from the previous shape 

methodologies. General polyhedra can easily be defined with numerous faces, edges, and vertices, as well 

as asymmetries to closely agree with actual granular particle shapes. Traditionally, polyhedra shapes have 

been difficult to implement in DEM software, since the tradeoff between particle shape complexity, 

number of particles, memory, and CPU speed has driven most commercial DEM software towards 

utilizing nonpolyhedra particle shapes. Advances in computer hardware in recent years, however, has 

made utilization of polyhedra particles a practical strategy. The DEM software BLOKS3D developed at 

the University of Illinois by Zhao (Ref. 14), is an example of this strategy. Simulations using 800 

particles composed of tetrahedrons (four vertices), pyramids (five vertices), cubes (eight vertices), and 

rhombic-dodecahedrons (14 vertices) have been reported by Zhao. 

2.8 Contact Detection and Force Model 

The most difficult parts of DEM from a CPU resources perspective are contact detection and contact 

force model implementation (Ref. 15). The complexity of a contact detection algorithm increases rapidly 

with particle shape complexity. Particle shapes that are simple spheres require the least complex detection 

strategy. Poly-ellipsoids, super-ellipsoids (also called superquadrics), asymmetric superquadrics, and 

general polyhedra require more sophisticated contact detection algorithms. 

On the two extremes of contact detection algorithms, there are (1) continuous math algorithms based 

on analytic models of particle shape and (2) discrete math algorithms, which use gridded representations 

of particles. Numerous combinations of analytical and gridded algorithms, along with ad hoc algorithms, 

led to many possible approaches to contact detection. The details of the approach used may be intimately 

dependent on other unrelated aspects of the DEM code. 

 
NASA/TM—2010-216257 11 

One method of determining the contact force direction in a continuum domain for arbitrarily shaped 

particles, which is loosely based on the method of Hogue (Ref. 16), makes use of the equations of the 

continuous particle surfaces. In this approach, two particle surfaces, described by f(x, y) = 0 and 

g(x, y) = 0 , are said to overlap when a set of coordinates can be found that simultaneously satisfy both 

equalities. For example, the two-dimensional case of Equation (3) defines an ellipse, so that 

2

1

2

1

f( , ) 1 



 



  −  



 = − b

v

a

u

u v

 

( )

2

2

0

2

2

0 g , 1 



 



 −  −  



 − = − b

v v

a

u u

u v

 (17) 

where u and v are local coordinates of the particle coordinate system. Since the particles are rotated and 

translated relative to the origin of the global coordinate system: 





 



 ⋅



 



  =



 





y

x

v

u

0 1

1 0 for particle 1 (18a) 





 



 ⋅



 





− φ φ

φ φ  =



 





y

x

v

u

sin cos

cos sin for particle 2 (18b) 

then 

2

1

2

1

f( , ) 1 



 



  −  



 = − b

y

a

x

x y (19a) 

( )

2

2

0

2

2

0 cos sin cos sin

g , 1 



 



 − φ + φ −  −  



 φ + φ − = − b

x y v

a

x y u

x y (19b) 

For this example with a1 = 0.5, b1 = 1.1, a2 = 0.4, b2 = 1.5, u0 = 0.3, v0 = 1.8, and φ = –30°, the 

solution of the simultaneous equations 0 f(x, y) = and 0 g(x, y) = are 





 



  =



 





0.672119

0.395809

1

1

y

x 



 



  =



 





 0.0965897

 0.498069

2

2

y

x (20a) 

This case is depicted in Figure 12 where a line segment connects the two solutions. The contact force 

direction angle is the perpendicular to the slope of the line 

 2

tan

2 1

1 2 1 π  +



 





−

− θ = −

x x

y y (20b) 

and is equal to 10.08° in the above case. 

In general, there will be two point solutions in two dimensions and a set of solutions comprising a 

closed loop in three dimensions. The two-dimensional case shown in Figure 12 defines a line connecting 

the two points of intersection. A normal vector to this line, located at the midpoint of the line, is one 

method of defining an effective contact force direction as an input prelude to the contact force model. The 

contact force model uses the overlap area or volume information to determine the magnitude of contact 

repulsion. 
NASA/TM—2010-216257 12 

Figure 12.—Contact detection of two particles 

described by analytical expressions. The contact 

normal direction is θ = 10.08°. 

In the three-dimensional case, the force direction of contacting particles is easily defined for spherical 

particles. Most DEM software, such as PFC3D (Ref. 5b) by Itasca, defines the force direction, in the case 

of simple spherical particles, as along the axis connecting contacting particle centers. In the case of more 

complex particles composed of clusters of overlapping spheres of different sizes, the force direction can 

still be modeled as along the axis connecting sphere centers of the two local spheres involved in the 

contact, assuming that the deformation effect due to contact can be isolated to only the two contacting 

spheres. In this case, the torque about the center of mass of the cluster must also be included in the 

equations of motion. 

The details of force direction are intimately coupled to the details of the particle shape model. For 

example, in the shape model used by Hogue (Ref. 16), particles are made up of triangular surfaces and 

nodes, where the points of the triangles are nodes. When two particles contact, a node of one will usually 

penetrate the triangular surface of another. The force direction is defined as the normal to the triangular 

surface. Other possibilities include node-to-node and triangle-to-triangle contact, which are treated as 

special cases. 

A common characteristic of most approaches to modeling contact force direction is the explicit or 

implicit need to know the normal to the surface of the particle at the point of contact. The gradient of the 

particle surface is a convenient means to quantify the surface normal of the continuous particle shape. For 

example, going back to the ellipsoid of Equation (3), the normal direction at any point on the surface is 

found from the gradient: 

f(x, y) 

g(x, y) 

y

Particle 1

(x2, y2) 

Contact 

Normal 

φ

1 00 0

00 0

cos sin

cos sin

uu v

vv u

−    φ − φ ⋅ =    φ + φ   

R

θ

Particle 2

(x1, y1) 

x
NASA/TM—2010-216257 13 

































=











  − 

 



  +  



  +  





































∂

∂

∂

∂

∂

∂

∇ =

2

2

2

2 2 2

2

2

2

f( , , ) 1

c

z

b

y

a

x

c

z

b

y

a

x

z

y

x

x y z (21) 

The difficulty of applying Equation (21) directly to the general case of arbitrary shaped particles is 

that it must be decided at what point on the contacting surface and associated normal direction is most 

relevant to the final contact force direction. A brute force method could involve integrating Equation (21) 

over the overlapping surface areas of both contacting particles to get an average normal direction, and 

thus a final force direction. 

On the other extreme of contact detection methods, the full discrete approach to contact detection can 

be performed using any particle shape since once the shape has been defined, the particle location is 

simply designated by a cluster of pixels. As shown in Figure 13, contact detection is the simultaneous 

enabling of pixel bits in the quantized space of the overlapping particles. The advantage of full discrete 

contact detection is that there is no penalty due to particle shape. However, the obvious disadvantage is 

resolution of the particle states (position and velocity) in that discrete space. 

The discrete version of the contact algorithm involves rectangular gridding of the particles so that any 

overlap can be detected and described by simultaneous occupancy of adjacent particles in one or more 

grids (or pixels for the two-dimensional case). Once a contact condition is detected by the overlap of one 

or more pixels associated with a particle pair, the direction of the contact normal can be computed by 

flagging the set of overlapping pixels. As is true for most mathematical models that substitute a simple 

solution for a more complex mathematical description, numerous approaches are possible. Below we 

derive (as an example) a method for calculating contact force direction that can be implemented as a 

computer algorithm. However we do not imply that this method is necessarily a method of choice by any 

specific DEM software vendor. 

The set of M-occupied pixels Xi can be transformed to a principal component direction by using the 

Hotelling transform, and then by taking the perpendicular direction, the contact normal direction can be 

defined. The procedure is described below. 

First, the center of mass mx of the set of overlapping pixels is calculated: 

1

1 M

x i M i=

m X =  (22) 

Next, the centered covariance matrix is computed: 

1

1 M

T T

x ii ii M i=

C XX mm = ⋅ −⋅  (23) 

The eigenvalues of Cx define the principal axis directions. Adding 90° to the angle that defines the 

first principal axis defines the normal direction, which is the contact force direction. 

Several commercial off-the-shelf (COTS) DEM software packages are available. In addition, several 

custom software packages have been developed, some of which may become COTS products in the 

future. In the final analysis, the choice of particle shape is predetermined by the details of specific DEM 

software and the particle shapes it supports. Some of these are listed in Table I. 
NASA/TM—2010-216257 14 

Figure 13(a).—Contact detection of two particles 

described by discrete particle representations 

with fine spatial grid. 

Figure 13(b).—Contact detection of two particles 

described by discrete particle representations 

with coarse spatial grid. 

TABLE I.—PARTICLE SHAPE CAPABILITIES OF SOME DEM SOFTWARE 

Software Supported particle shapes 

PFC2D/3D (Ref. 5) Spheres, rigid sphere clusters: overlapping and nonoverlapping 

EDEM (Ref. 6) Spheres, sphere clusters: overlapping and nonoverlapping 

BLOKS3D (Ref. 16) Polyhedra 

ROCKS3D Polyhedra 

MIMES (Ref. 17) Superquadric shapes 

YADE (Ref. 18) Spheres 

Polyhedra (work in progress) 

3.0 Analytical Representations of Particle Size Distribution 

The data-fitting procedure outlined in Bagnold (Ref. 20) provides valuable insight into the nature of 

the desert sand particle size distribution (PSD). The technique used by Bagnold, which was appropriate 

for the time period, assumes that the tools of choice are a pencil, graph paper, and a ruler. Fortunately, 

current computer software tools are beginning to provide a reliable means to fit arbitrary nonlinear 

functions to a data set. Taking this approach, it makes sense to define the PSD as a single formula, 

covering the entire particle size range.1

 

1

When dealing with particle size histograms and the numerous forms of size distributions and their moments, 

carefully keeping track of units reduces confusion and the risk of error. Also, it should be noted that histograms are 

summed, while distributions are integrated. 
NASA/TM—2010-216257 15 

3.1 Size Distribution by Mass 

The fractional mass of particles in a differential size increment dD, based on Bagnold’s work is 

( ) ( ) ( )

D

D D D D

N N D D s c ( )d d

/

0

/

0

0

γ γ γ +

= (24a) 

where the small slope s is always positive and coarse slope c is always negative. The PSD of 

Equation (24a) is normalized by setting the integral of N(D), for sizes from 0 to ∞, equal to 1: 

( ) 1

0

= 

∞

N D dD (24b) 

Note that in the case of Bagnold’s data, the integral is set equal to 100 percent. 

For the small-diameter case, D << D0, so that Equation (24a) reduces to 

( )s N D N0 D D0 ( ) = (25) 

If Equation (25) is plotted on a log-log plot, as described by Bagnold, a straight line graph of slope s 

results in 

0 0

0 0 log log log log

y sR sR

y N N s D s D

= + −

≡ = + − (26) 

where R ≡ log D . 

In the large diameter case, D >> D0, so that Equation (24) reduces to 

( )c N D N0 D D0 ( ) = (27) 

If Equation (27) is plotted on a log-log plot, a straight-line graph of slope c results, which is 

equivalent to 

0 0

0 0 log log log log

y cR cR

y N N c D c D

= + −

≡ = + − (28) 

Recall that c has been defined as negative value according to Bagnold. 

The exponent parameter γ in Equation (24) is introduced as a means to control the transition rate 

between the small-diameter positive slope regime and the large-diameter negative slope regime. 

 
NASA/TM—2010-216257 16 

Figure 14.—Particle size distribution fit to desert sand using 

Equation (24a). 

3.3 Particle Size Distribution 

Equation (24) describes a size distribution based on particle mass. This formula is an empirical fit to 

observed grain size distributions of desert sand (Fig. 14) based on measurements using a multilevel sieve. 

However, the fundamental size distribution of primary interest for an assemblage of particles is P(D), 

which may be expressed in units of [m–3 mm–1], where P(D) dD is the number of particles per unit volume 

of space within a differential size increment dD. The total number of particles of all sizes in a sample 

volume VS is then 



∞

=

0

NV VS P(D)dD (29) 

Assuming the packing density is homogeneous and repeatable, the total mass of the assemblage of 

particles in a volume VS is 



∞

= πρ

0

3

6

1 MV VS D P(D) dD (30) 

where ρ is the particle grain density, assuming all grains have equal density and are spherical. 

Equation (30) can be rewritten in terms of the volume density of the assemblage of particles by dividing 

both sides by VS: 



∞

ρ = πρ

0

3

6

1

V D P(D)dD (31) 

 
NASA/TM—2010-216257 17 

Now the mass size distribution described by Equation (24) can be expressed in terms of P(D): 

V

D P D

D P D D

D P D

N D

ρ

ρ = π

πρ

πρ =



∞

( )

( ) d

( ) ( )

3

6

1

0

3

6

1

3

6

1

(32) 

Equation (32) can now be inverted, thus providing a means to estimate the particle size distribution 

P(D) from the measured mass distribution N(D): 

( ) ( ) ( ) γ − γ + γ

⋅ π ρ

ρ =

π ρ

ρ =

s c

D D D D

N

D

N D

D

P D

V

V

0 0

0

3

6

1

3

6

1

/ /

( ) ( )

(33) 

where the parameters N0, s, c, γ, ρ, and ρV are determined by measurement. Note that the particle packing 

ratio is less than one ( ) ρV ρ < 1 , because of gaps between particles of an assemblage. 

3.4 FPA Data of JSC-1a 

Instruments that measure and count individual particles typically output a histogram count of particles 

in bins of equal intervals for the property of interest. This is the case with the fine particle analyzer (FPA) 

used to measure a sample of JSC-1a at the NASA Kennedy Space Center Granular Mechanics and 

Surface Systems Laboratory. An example of particle size counts from the first nine size bin intervals of an 

FPA measurement are shown in Table II. Each bin interval is defined by a ΔD =1/ 256 [mm], shown in 

the first column. The second column shows counts Ck of particles that fall into the size bin centered at the 

equivalent spherical diameter Dk where k is the bin index. 

To convert from a histogram of counts to a size distribution n(D) in units of [mm–1], the counts are 

divided by the bin interval, ΔD : 

D

C

n D n k

k k Δ ( ) ≈ = [mm–1] (34) 

The third column of Table II shows nk, the corresponding values of the number of particles in bin k, 

determined from the Ck and ΔD =1/ 256 [mm], using Equation (34). The nk are plotted in Figure 15, 

along with an empirical fit based on a sum of exponentials: 

1 2 1 2 ( ) D

D

D

D

n D D N e N e − − γ = + [mm–1] (35) 
NASA/TM—2010-216257 18 

TABLE II.—FPAa PARTICLE COUNT DATA FOR JSC-1a 

Particle size, 

Dk 

[mm] 

Particle count, 

Ck 

Count per size interval, 

nk = Ck/ΔD

[mm–1] 

0.003906 1232155 315431680 

.007813 629101 161049856 

.011719 676994 173310464 

0.015625 659910 168936960 

.019531 602111 154140416 

.023438 546315 139856640 

0.027344 496099 127001344 

.03125 448268 114756608 

.035156 403354 103258624 

∶ ∶ ∶ a

FPA refers to fine particle analyzer at NASA Kennedy 

Space Center. 

TABLE III.—FPAa MASS DATA FOR JSC-1a 

Particle size 

Dk, 

[mm] 

TL ratio, 

thickness/length 

Mass per interval, 

mk, 

[g⋅mm–1] 

0.003906 0.30649069 0.005223 

.007813 .51233492 .041281 

.011719 .54465904 .162795 

0.015625 0.56233665 0.392799 

.019531 .57190647 .716264 

.023438 .57940036 1.143148 

0.027344 0.58504106 1.670413 

.03125 .58960083 2.277133 

.035156 .59410737 2.948024 

⋮ ⋮ ⋮ a

FPA refers to fine particle analyzer at NASA Kennedy 

Space Center. 

Figure 15.—Fine particle analyzer (FPA) size distribution 

of JSC-1a, fitted to sum of exponentials. 

Particle mass data corresponding to the counts and diameters of Table II are shown in Table III. In the 

second column, L is defined as the longest chord across the particle, where a chord is defined as any line 

segment that crosses from one part of a particle's perimeter to another part of the same perimeter; and T is 

defined as the longest chord of all chords that are perpendicular to L. Typically in this context L denotes 

particle length and T, thickness. 

The size distribution is converted to mass distribution by multiplying by particle volume and particle 

density: 

( ) 2

(1 ) ( ) 3

6

1 n D R R

m D D TL + TL = π ρ [g⋅mm–1] (36) 

where the particle density ρ = 2.65 [g⋅cm–3]. The geometric factor RTL (1 + RTL ) / 2 in Equation (36) 

improves the particle volume estimate using a geometric correction factor to account for the nonspherical 

nature of JSC-1a. The TL ratio, RTL shown in Table III, is an output of the FPA data processing. Note 

that as RTL approaches 1, the geometric correction factor also approaches 1. 
NASA/TM—2010-216257 19 

Using the empirical n(D) function from Equation (35) and the equation for mass shown above, the 

FPA particle mass data is plotted in Figure 16. 

Going back to Bagnold’s method, Equation (24), the particle mass data from Table III can be fitted as 

shown in Figure 17. In this case the parameter N0 is expressed so that P(D) from Equation (29) is in units 

of [g⋅mm–1]: 

( ) ( ) γ





 



 +

=

γ γ s c D D D D

N

m D

/ /

( )

0 0

0 [g⋅mm–1] (37) 

A small-slope parameter s and coarse-slope parameter c is determine by this fit. The FPA particle 

count size distribution nk is plotted in Figure 18 along with the Bagnold’s fit, found by inverting 

Equation (36) with m(D) from Equation (37): 

( ) ( )( ) γ γ − γ

−

 + 

 



 +

π ρ

=





 



 + = π ρ

/

0

/

0

3

6

1

0

1

3

6

1

/ /

2

(1 )

( ) 2

(1 ) ( )

TL TL s c

TL TL

D D D D R R D

N

m D R Rn D D

[mm–1] (38) 

For comparison, a lognormal fit to the mass size distribution is given by 

( )

( ) D D 

e

D

m

m D

D

D

0

ln 2

0

0

2 ( )

2 2

0

π σ =

 σ 

 



 −

[g⋅mm–1] (39) 

Figure 16.—Fine particle analyzer (FPA) mass size distribution of JSC-1a, fitted to sum of 

exponentials, Equation (35). 
NASA/TM—2010-216257 20 

Figure 17.—Fine particle analyzer (FPA) mass size distribution of JSC-1a, fitted to Bagnold’s model. 

Figure 18.—Fine particle analyzer (FPA) size distribution of 

JSC-1a, fitted to Bagnold’s model using Equation (38). 
NASA/TM—2010-216257 21 

Figure 19.—Fine particle analyzer (FPA) mass size distribution of JSC-1a, fitted to lognormal distribution. 

The fit of Equation (39) to JSC-1a is shown in Figure 19. 

4.0 Concluding Remarks 

This report discusses topics that are important to the modeling of the interaction and dynamics of 

granular particles, and in particular, lunar stimulants. In Section 2.0, discrete element method (DEM) 

particle shapes were presented, where the choice of particle shape is motivated by many opposing 

tradeoffs, such as 

(1) The number of parameters describing the simulated particle surface (thus defining the interior 

versus exterior region). The goal is to attempt to closely match the shape of a real particle. 

(2) The total number of particles that can be incorporated in a simulation. 

(3) Computer processing resources needed: memory and total simulation time. 

(4) The difficulty versus ease of the contact detection algorithm and subsequent implementation of a 

force model. 

Most commercial and custom DEM software supports some type of complex particle shape in order 

to introduce asymmetry into the particle shape, so as to better simulate real granular particles interactions 

and behavior. 

In Section 3.0, a complete particle size distribution (PSD) formula is shown to fit the example desert 

sand data from Bagnold, using a coarse-grade exponent c and small-grade exponent s. An additional 

empirical parameter γ has been introduced to connect the two log-log linear fits in the small- and coarsegrade particle size regions. The key finding in this section is that Bagnold's model is a fair representation 
NASA/TM—2010-216257 22 

of the PSD except for particles below 20 μm, where it diverges somewhat. Also, the double exponential 

formulation is a fair representation of the size distribution, but the lognormal formulation is not very 

successful at matching the JSC-1a lunar stimulant used in this evaluation. Particle size data of JSC-1a 

obtained from a fine particle analyzer (FPA) at NASA Kennedy Space Center are also fitted to the PSD in 

order to extract equivalent small- and coarse-grade exponents. Note that the FPA data are uncorrelated 

quantitatively with other accepted PSD measuring tools, so there may be a bias in data that affects the 

quality of the PSD fit. 
NASA/TM—2010-216257 23 

Appendix A.—Microscope Image, Shadowgraphs, and 

SEM Spectrum of JSC-1a 

JSC-1a is a lunar regolith simulant that was developed by the NASA Johnson Space Center. It is a 

basaltic ash with a high glass content. Figure 20 is an optical microscope image of the larger particles in a 

sample of JSC-1a. Figure 21 presents shadowgraphs of JSC-1a stimulant showing the variety of shapes 

and sizes in the mixture. Figure 22 is an elemental analysis spectrum of a sample of JSC-1a using a Joel 

75000F scanning electron microscope (SEM) with energy-dispersive x-ray spectroscopy (EDS). Figure 23 

contains elemental maps of the sample in Figure 22, showing K-α emission lines for the detected 

elements (K lines result when an electron transitions to the innermost K shell from a 2p orbital of the L 

shell). 

Figure 20.—Microscope image of >600-μm particles of JCS-1a. 

 

Figure 21.—Shadowgraphs of JCS-1a. 
NASA/TM—2010-216257 24 

Figure 22.—Joel 75000F scanning electron microscopy (SEM) spectrum of JSC-1a. 

 
NASA/TM—2010-216257 25 

Figure 23.—Energy-dispersive x-ray spectroscopy maps, showing spatial mapping of elements corresponding to 

scanning electron microscopy (SEM) spectrum of Figure 22.

NASA/TM—2010-216257 27 

Appendix B.—Summary of Particle Shape Considerations 

Discrete element method (DEM) simulation software considers a granular material as being made up 

of a size distribution of idealized particle shapes. These shapes may be simple two-dimensional 

nonoverlapping or overlapping discs or ellipses, three-dimensional spheres or ellipsoids that overlap or do 

not overlap, or complex shapes such as polyhedrons or superquadrics. Table IV compares some of the 

pros and cons of these idealized particle types. 

TABLE IV.—COMPARISON OF DISCRETE ELEMENT METHOD (DEM) PARTICLE SHAPE TYPES 

Particle shape Pros Cons Comments and possible 

solutions to cons 

Incorporated in DEM 

software 

Sphere Simplest to 

implement; requires 

the minimum CPUa

resources, allowing a 

maximum number of 

simulated particles. 

Particle behavior is 

unrealistic because 

of simplified shape. 

Implements more complex 

force and friction models 

to match real particle 

behavior. 

PFC2D/PFC3D (Ref. 5) 

EDEM (Ref. 6) 

Yade (Ref. 19) 

Rigid 

nonoverlapping 

sphere cluster 

Relatively simple to 

implement. 

Approximates real 

particle shapes. 

In order to approximate 

a real particle shape, 

many spheres of 

various sizes are 

required, significantly 

increasing CPU 

resources. 

Decreases number of freestate parameters to only 

that needed to simulate a 

single rigid body. 

PFC2D/PFC3D (Ref. 5) 

EDEM (Ref. 6) 

Rigid 

overlapping 

sphere cluster 

Adds ability to better 

match real particle 

shapes. 

Again, in order to 

approximate a real 

particle shape, many 

spheres of various sizes 

are required. Some 

DEM codes (e.g., 

EDEM) can correctly 

account for total cluster 

mass density. 

Specifies a variable sphere 

mass density so that total 

cluster is the correct value. 

This increases complexity 

of sphere density 

specification. 

PFC2D (Ref. 5a) 

Ellipsoid Better matches real 

particle shapes. 

Contact force 

detection has some 

similarities to that of a 

sphere. CPU resources 

required are 

comparable to that of a 

sphere. 

Contact force detection 

is more complicated 

than for a sphere. The 

ability to match real 

particles is much better 

than a sphere, but is 

still limited. Contact 

detection is more 

difficult than for 

sphere. 

Ellipsoid is comparable to 

rigid sphere cluster as far 

as tradeoffs between CPU 

resources, number of 

possible particles, and 

similarity of simulated to 

real particle shape. 

Wang (Ref. 21) 

Superquadric More simulated 

particle parameters are 

provided in a closed 

form description, 

which result in more 

shapes. All of the 

particle parameters 

can be computed 

analytically. 

Similar difficulties 

encountered as with 

ellipsoids. When 

concave exponent 

values are used (n < 2), 

additional difficulties in 

contact detection and 

force model are 

present. 

Because purely analytic 

methods used for contact 

detection and force model 

can get messy, a gridded 

spatial model of particle 

and particle interactions 

may be a better choice. 

MIMES (Ref. 17) 

Asymmetric 

superquadric 

Similar to 

superquadric, but 

many more parameters 

can be added, 

providing a better 

match to real particle 

shapes. 

All problems described 

for superquadric 

particle shape are 

compounded by the 

inclusion of additional 

parameters. 

Even though an analytical 

approach to contact 

detection and force model 

is possible, it may not be 

practical in this case. A 

gridded numerical strategy 

may be the best approach. 

None found 
NASA/TM—2010-216257 28 

Particle shape Pros Cons Comments and possible 

solutions to cons 

Incorporated in DEM 

software 

Eight-quadrant 

superquadric 

A subset of 

asymmetric 

superquadric particle 

type, which provides a 

method to easily 

guarantee no surface 

discontinuities. 

Problems are similar to 

the asymmetric 

superquadric. 

Gridded representation of 

particle and particle 

interactions is probably the 

most practical approach. 

None found 

Poly-ellipsoid A subset of eightquadrant superquadric 

particle type, which 

guarantees no surface 

gradient 

discontinuities (i.e., no 

sharp ridges or edges). 

Problems with contact 

detection are more 

severe than with the 

ellipsoid, but less 

severe than with the 

general superquadric. 

A gridded approach may 

be best, but other pseudo 

analytical approaches have 

been developed which 

claim good success. 

Peters (Ref. 13) 

Polyhedra The polyhedron-based 

particle shape is the 

most general shape of 

those discussed in this 

summary. 

Analytical closed form 

contact detection 

methods are very 

difficult except in 

specific simple cases. 

Contact force 

determination on edges 

is also difficult and 

ambiguous. 

A gridded approach is 

almost a necessity. Sharp 

edges can be eliminated by 

performing a spatial 

convolution of the particle 

surface or volume with a 

sphere (Minkowski sum, 

Ref. 22). 

Nezami (Ref. 15) 

a

CPU refers to central processing unit.
NASA/TM—2010-216257 29 

Appendix C.—Eight-Quadrant Superquadric Particle Shape Examples 

Using Equation (8), the 12 particle shape parameters are randomly modified to demonstrate the 

resulting diversity of particle shapes. The figures in this appendix were created using Mathematica 

version 7.0. 

( )

( )

( )











θ ≤ θ < π

θ ≤ θ <

θ φ =









≤ θ < π

θ φ π ≤ φ < π

θ φ ≤ φ < π

θ φ =

≤ θ < π











θ φ ≤ φ <

θ φ − ≤ φ <

θ φ =

π

π

π π

π π

2

2 /

2

2

2 /

1

2 /

2

2 /

1

2

3

2

2 /

2

2 2

2 /

1

cos for

cos for 0

( , )

0

sin sin for 2

sin sin for 0

( , )

0

sin cos for

sin cos for

( , )

2

1

2

1

2

1

c 

c z

b 

b y

a 

a x

M

M

L

L

K

K

 (8) 
NASA/TM—2010-216257 30
NASA/TM—2010-216257 31
NASA/TM—2010-216257 32
NASA/TM—2010-216257 33

NASA/TM—2010-216257 35 

References 

 1. Cleary, Paul W.; and Sawley, Mark L.: DEM Modelling of Industrial Granular Flows: 3D Case 

Studies and the Effect of Particle Shape on Hopper Discharge. Appl. Math. Model., vol. 26, no. 2, 

2002, pp. 89–111. 

 2. Santamarina, J.C.; and Cho, G.C.: Soil Behaviour: The Role of Particle Shape. Adv. Geotech. Eng. 

Skempton Conf. Proc., 2004, pp. 604–617. 

 3. Zeghal, M.: Effect of Particle Shapes on the Resilient Behaviour of Aggregate Materials. Canadian 

Society for Civil Engineering Annual Conference, Victoria, 2001, pp. 1–5. 

 4. Peña, A.A.; García-Rojo, R.; and Herrmann, H.J.: Influence of Particle Shape on Sheared Dense 

Granular Media. Granular Matter, 2007, pp. 279–291. 

5a. Distinct Element Modeling for Micromechanical Analysis of Geomaterials and Particulate Systems 

in Two Dimensions. Itasca, Minneapolis, MN, 2010. http://www.itascacg.com/pfc2d/ Accessed 

October 12, 2010. 

5b. Distinct Element Modeling for Micromechanical Analysis of Geomaterials and Particulate Systems 

in Three Dimensions. Itasca, Minneapolis, MN, 2010. http://www.itascacg.com/pfc3d/ Accessed 

October 12, 2010. 

 6. EDEM Software. DEM Solutions, Edinburgh, Midlothian, 2010. http://www.demsolutions.com/software/edem-software.php Accessed October 12, 2010. 

 7. Ferellec, Jean-Francois; and McDowell, Glenn R.: A Method to Model Realistic Particle Shape and 

Inertia in DEM. Granular Matter, vol. 12, no. 5, 2010, pp. 459–467. 

 8. Ashmawy, Alaa K.; Sukumaran, Beena; and Vinh Hoang, V.: Evaluating the Influence of Particle 

Shape on Liquefaction Behavior Using Discrete Element Modeling. Proceedings of the 13th 

International Offshore and Polar Engineering Conference, Honolulu, HI, 2003, pp. 1089–1096. 

 9. Markauskas, D.: Discrete Element Modelling of Complex Axisymmetrical Particle Flow. 

Mechanika, vol. 62, no. 6, 2006, pp. 32–38. 

10. Vu-Quoc, L.; Zhang, X.; and Walton, O.R.: A 3–D Discrete Element Method for Dry Granular 

Flows of Ellipsoidal Particles. Computer Methods in Applied Mechanics and Engineering on 

Dynamics of Contact and Impact Problems, vol. 187, vols. 3–4, 2000, pp. 483–528. 

http://clesm.mae.ufl.edu/~vql/pdf/Vu-Quoc2000.pdf Accessed October 12, 2010. 

11. Preece, D.S.; and Perkins, Eric D.: Sand Production Modeling Using Superquadratic Discrete 

Elements and Coupling of Fluid Flow and Particle Motion. SAND99–0366C, 1999. 

12. Wolfram Mathematica 7. Wolfram Research, Champaign, IL, 2010. 

http://www.wolfram.com/products/mathematica/index.html Accessed October 12, 2010. 

13. Peters, John F., et al.: A Poly-Ellipsoid Particle for Non-Spherical Discrete Element Method. 

Engineering Computations, vol. 26, no. 6, 2009, pp. 645–657. 

14. Zhao, Dawei, et al.: Three-Dimensional Discrete Element Simulation for Granular Materials. Eng. 

Comput., vol. 23, no. 7, 2006, pp. 749–770. 

15. Nezami, Erfan G., et al.: Shortest Link Method for Contact Detection in Discrete Element Method. 

Int. J. Numer. Anal. Methods Geomech., vol. 30, no. 8, 2006, pp. 783–801. 

16. Hogue, Caroline: Shape Representation and Contact Detection for Discrete Element Simulations of 

Arbitrary Geometries. Eng. Comp., vol. 15, no. 3, 1998, pp. 374–390. 

17. Nezami, Erfan G., et al.: Simulation of Front End Loader Bucket—Soil Interaction Using Discrete 

Element Method. Int. J. Numer. Anal. Meth. Geomech., vol. 31, 2007, pp. 1147–1162. 

https://netfiles.uiuc.edu/hashash/www/web/Publications/Nezami_etal_2007.pdf Accessed 

October 12, 2010. 

18. Williams, John R., et al.: Dynamic Wave Propagation in Particulate Materials With Different Particle 

Shapes Using a Discrete Element Method. Proc. Eng. Mech., vol. 1, 1995, pp. 493–496. 
NASA/TM—2010-216257 36 

19. Yade. 2010. https://yade-dem.org/wiki/Yade Accessed October 12, 2010. 

20. Bagnold, Ralph Alger: The Physics of Blown Sand and Desert Dunes. Dover Publications, Methuen, 

London, 1954. 

21. Wang, Linbing; Park, Jin-Young; and Fu, Yanrong: Representation of Real Particles for DEM 

Simulation Using X-Ray Tomography. Constr. Build. Mater., vol. 21, no. 2, 2007, pp. 338–346. 

22. Alonso-Marroquín, Fernando; and Wang, Yucang: An Efficient Algorithm for Granular Dynamics 

Simulations With Complex-Shaped Objects. Granular Matter, vol. 11, no. 5, 2009, pp. 317–329. 

http://www.springerlink.com/content/w58x84mgj8k82083/fulltext.pdf Accessed October 12, 2010. 
REPORT DOCUMENTATION PAGE Form Approved 

OMB No. 0704-0188 

The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and maintaining the 

data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions for reducing this 

burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. 

Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid OMB 

control number. 

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS. 

1. REPORT DATE (DD-MM-YYYY)

01-12-2010 

2. REPORT TYPE

Technical Memorandum 

3. DATES COVERED (From - To)

4. TITLE AND SUBTITLE

A Review of Discrete Element Method (DEM) Particle Shapes and Size Distributions for 

Lunar Soil 

5a. CONTRACT NUMBER

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

6. AUTHOR(S)

Lane, John, E.; Metzger, Philip, T.; Wilkinson, R., Allen 

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

WBS 387498.04.01.01.02.03 

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration 

John H. Glenn Research Center at Lewis Field 

Cleveland, Ohio 44135-3191 

8. PERFORMING ORGANIZATION

 REPORT NUMBER

E-17243 

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration 

Washington, DC 20546-0001 

10. SPONSORING/MONITOR'S

 ACRONYM(S)

NASA 

11. SPONSORING/MONITORING

 REPORT NUMBER

NASA/TM-2010-216257 

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified-Unlimited 

Subject Categories: 31, 59, and 91 

Available electronically at http://gltrs.grc.nasa.gov 

This publication is available from the NASA Center for AeroSpace Information, 443-757-5802

13. SUPPLEMENTARY NOTES

14. ABSTRACT

As part of ongoing efforts to develop models of lunar soil mechanics, this report reviews two topics that are important to discrete element 

method (DEM) modeling the behavior of soils (such as lunar soils): (1) methods of modeling particle shapes and (2) analytical 

representations of particle size distribution. The choice of particle shape complexity is driven primarily by opposing tradeoffs with total 

number of particles, computer memory, and total simulation computer processing time. The choice is also dependent on available DEM 

software capabilities. For example, PFC2D/PFC3D and EDEM support clustering of spheres; MIMES incorporates superquadric particle 

shapes; and BLOKS3D provides polyhedra shapes. Most commercial and custom DEM software supports some type of complex particle 

shape beyond the standard sphere. Convex polyhedra, clusters of spheres and single parametric particle shapes such as the ellipsoid, polyellipsoid, and superquadric, are all motivated by the desire to introduce asymmetry into the particle shape, as well as edges and corners, in 

order to better simulate actual granular particle shapes and behavior. An empirical particle size distribution (PSD) formula is shown to fit 

desert sand data from Bagnold. Particle size data of JSC-1a obtained from a fine particle analyzer at the NASA Kennedy Space Center is 

also fitted to a similar empirical PSD function. 

15. SUBJECT TERMS

Geotechnical engineering; Particle size distribution; Lunar soil; Shapes; Discrete modeling 

16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF

 ABSTRACT

UU 

18. NUMBER

 OF 

 PAGES

35 

19a. NAME OF RESPONSIBLE PERSON

STI Help Desk (email:help@sti.nasa.gov) 

a. REPORT

U 

b. ABSTRACT

U 

c. THIS 

PAGE

U 

19b. TELEPHONE NUMBER (include area code)

443-757-5802 

Standard Form 298 (Rev. 8-98)

Prescribed by ANSI Std. Z39-18