Discrete Element Method

The discrete element method (DEM):

is intended for modeling events in which large numbers of discrete particles contact each other;
models each particle with a single-node element that has a rigid spherical shape, which may represent an individual grain, tablet, shot peen, or other simple body;
is a versatile tool for modeling particulate material behavior in pharmaceutical, chemical, food, ceramic, metallurgical, mining, and other industries; and
is not meant for modeling deformation of a continuum, but DEM can be used together with finite elements for modeling discrete particles interacting with deformable continua or other rigid bodies.

This page discusses:

Introduction
Applications
Strategies for Creating and Initializing a DEM Model
Strategy for Reducing Solution Noise
Time Incrementation Considerations for Non-Hertzian Contact
Automatic Time Incrementation for Hertzian Contact
Initial Conditions
Boundary Conditions
Loads
Elements
Constraints
Interactions
Output
Limitations
Input File Template
References

Introduction

The discrete element method (DEM) is an intuitive method in which discrete particles collide with each other and with other surfaces during an explicit dynamic simulation. Typically, each DEM particle represents a separate grain, tablet, shot peen, etc. DEM is not applicable to situations in which individual particles undergo complex deformation. Therefore, DEM is unlike, and conceptually simpler than, the smoothed particle hydrodynamic (SPH) method in which groups of particles collectively model a continuum body (see Smoothed Particle Hydrodynamics).

For example, DEM is well-suited for particle mixing applications, such as that shown in Figure 1. In this application DEM is used to model the initially separated blue and white particles, and rigid finite elements are used to model two mixing augers and the box-shaped container. The sequence of deformed plots in Figure 1 shows the particle response as the augers turn. DEM results for simulations such as this are often best viewed with animations. Another example of using DEM for a mixing application is described in Mixing of granular media in a drum mixer.

Each DEM particle is modeled with a single-node element of type PD3D. These elements are rigid spheres with specified radii. Nodes of PD3D elements have displacement and rotational degrees of freedom. Rotations of DEM particles can significantly influence contact interactions when friction is considered.

General contact definitions are easily extended to include interactions among DEM particles and interactions between DEM particles and finite-element-based (or analytical) surfaces. Large relative motion among particles is typical for DEM applications. Particle-to-particle interactions can involve like or unlike particles. Each particle can be involved in many contact interactions simultaneously. DEM particle interactions use finite contact stiffness, which introduces some compliance into the particle systems. For example, the contact stiffness can be specified such as to reflect the macroscopic stiffness of a packed granular material model with DEM.

For example, consider the interactions between the two spherical particles shown in Figure 2.

The three cases show two undeformed spheres just touching, two deformed spheres pushed toward one another with contact strictly enforced, and two rigid spheres pushed toward one another with some penetration. The distance between the centers of the spheres is the same for the cases shown in the middle and on the right in Figure 2. The physical behavior corresponds to the middle case. The case on the right corresponds to a DEM approximation. If the variable $δ$ is defined as

δ = r_{1} + r_{2} - d,

where $r_{1}$ and $r_{2}$ are the radii of the two spheres and d is the distance between the sphere centers, $δ = 0$ when the undeformed spheres are just touching and $δ > 0$ if the distance between the sphere centers is less than the combined radii. For the DEM approximation, $δ$ corresponds to the maximum penetration distance between the particles. You can improve the accuracy for some DEM applications by tuning the contact stiffness relationship (contact force F versus penetration $δ$ ) for DEM particles to reflect the Hertz contact solution (middle case in Figure 2). See Mixing of granular media in a drum mixer for further discussion of tuning the contact stiffness.

Applications

DEM is a versatile tool for modeling particulate material behavior in pharmaceutical, chemical, food, ceramic, metallurgical, mining, and other industries. DEM applications include the following categories:

Particle packing: involves processes such as pouring or deposition under gravity (such as sandpiling), vibration after deposition of particles, and compaction.
Particle flow: may occur under gravity only (as in the case of a hopper) or under gravity and other driving forces (such as for mixers and mills).
Particle-fluid interaction: occurs in transport of granular material within a fluid flow, during wavelike motion, and during fluidization (wherein fluid flows upwards through a bed of particles).

DEM can provide insight for many situations that are difficult to investigate with other computational methods or with physical experiments.

Strategies for Creating and Initializing a DEM Model

Particulate media often consist of randomly distributed grains of varying sizes. Generating an initial mesh for a DEM analysis can be challenging. A common strategy for DEM is to specify approximate initial positions of particles with some gaps between them and to allow the particles to settle into position under gravity loading during the first step. For example, this strategy is used for the mixing analysis shown in Figure 1: the augers are kept stationary during the first step in which the particles are allowed to settle, and the augers are turned during the second step to study the mixing behavior (which is the focus of the figure shown). The particle generator can be used to create DEM models (see Particle Generator).

Strategy for Reducing Solution Noise

The solution noise generated by numerous opening and closing contact conditions can be reduced by applying a small amount of mass proportional damping. For more information, see Alpha Damping.

Time Incrementation Considerations for Non-Hertzian Contact

DEM uses the explicit dynamics procedure type. In most cases Abaqus/Explicit automatically controls the time increment size, as discussed in Automatic Time Incrementation, based on stiffness and mass characteristics of the model. The relationship between the maximum stable time increment size, mass, and stiffness properties is complex. The stable time increment size tends to be proportional to the square root of mass and inversely proportional to the square root of stiffness. However, a stable time increment cannot be computed for each PD3D element because the particles are rigid, so you must specify a fixed time increment size for purely DEM analyses (see Fixed Time Incrementation).

Contact interactions among DEM particles can influence the appropriate time increment size. DEM analyses without tightly packed particles may simply call for a contact stiffness that is large enough to avoid significant penetrations, rather than a contact stiffness that is highly representative of physical stiffness characteristics of the particles (which are each modeled as rigid with DEM). If you do not specify the contact stiffness, Abaqus/Explicit assigns a default (penalty) contact stiffness based on the time increment size and mass/rotary inertia characteristics of the particles. This is undesirable since it is difficult to ensure that the time increment size is small enough to result in a sufficiently large penalty stiffness.

If you specify the non-Hertzian DEM contact stiffness, you must ensure that the time increment used for the analysis is small enough to avoid numerical instabilities. For dense three-dimensional packing of particles where each particle simultaneously contacts many particles, the numerical stability considerations are complex. A general guideline is that the time increment should not exceed $0.4 \sqrt{\frac{m}{k}}$ , where m and k represent the particle mass and contact stiffness, respectively. In some applications an even smaller time increment, such as $0.1 \sqrt{\frac{m}{k}}$ , may result in an improved solution.

If particle velocities become very large, the amount of incremental motion can influence the appropriate time increment size. Accurate resolution of particle motion sometimes requires specifying a smaller time increment than the maximum numerical stability time increment.

Automatic Time Incrementation for Hertzian Contact

The time increment size for a stable and accurate DEM analysis depends on several different factors, such as the underlying contact properties, the size of the particles, and the relative motion of the particles. Since these controlling factors are problem dependent and vary during the analysis, choosing an appropriate direct time increment for a DEM analysis can be difficult. When the Hertz or JKR-type pressure overclosure is specified, Abaqus/Explicit automatically controls the time increment size to achieve a stable and accurate solution.

The following criteria are used to control the time increment size:

Stability: based on particle mass, inertia, contact stiffness, and the number of contact interactions with other particles
Tracking accuracy: based on particle size and velocity
Duration of collision: restricted to a fraction of the duration of the collision between impacting particles
Rayleigh wave propagation: based on the time taken for a Rayleigh wave to travel from pole to pole of a particle
Separation distance for JKR model: limit the relative normal motion based on the separation distance for the JKR-type pressure overclosure
Tangential tracking accuracy: limit the relative tangential slip between particles. The slip depends on the particle rotation.

The chosen time increment is the smallest of the above six criterion. A scaling factor is associated with each of the above criterion.

Depending on the nature of the analysis, different criterion may control the time increment size at different times during the analysis. The stability criterion may dictate the time increment when particles are confined and subjected to compression, whereas for fast-moving impacting particles, the duration of collision criterion may control the time increment. The default scaling factor applied to each of the above criterion is based on numerical experimentation for low- to moderate-speed collisions. For high-speed impacts where the contact stiffness increases rapidly, the scaling factors related to the Rayleigh wave and collision duration may need to be reduced. You can specify values for the scaling factors using section controls.

The automatic time incrementation also works in conjunction with the particle generator provided the contact interactions specify the Hertz or JKR-type pressure overclosure type.

Input File Usage

Use the following options to modify the default scaling factors for one or more criterion for automatic time incrementation:

DISCRETE SECTION, CONTROLS=section_control_name
…
DISCRETE ELASTICITY
effective Young's Modulus, effective Poisson's ratio

SECTION CONTROLS, NAME=section_control_name
blank line
blank line
blank line
blank line
 $f_{s t a b i l i t y}$ ,  $f_{t r a c k i n g}$ ,  $f_{c o l l i s i o n}$ ,  $f_{R a y l e i g h}$ ,  $f_{s e p a r a t i o n_d i s t}$ ,  $f_{r o t a t i o n}$

STEP
DYNAMIC, EXPLICIT, ELEMENT BY ELEMENT
…

Initial Conditions

Initial conditions pertinent to mechanical analyses can be used in a discrete element method analysis. All of the initial conditions that are available for an explicit dynamic analysis are described in Initial Conditions.

Boundary Conditions

Boundary conditions are defined as described in Boundary Conditions. Boundary conditions are rarely applied on individual particles in DEM.

Loads

The loading types available for an explicit dynamic analysis are explained in About Loads. Gravity loads are very important for settling and particulate flow analysis in DEM. Concentrated loads are rarely applied on particles.

Elements

The discrete element method uses PD3D elements to model individual particles. These 1-node elements define individual grains of a particulate media, are spherical in shape, and are modeled as rigid (any compliance is built into the contact model). These particle elements use existing Abaqus functionality to reference element-related features such as initial conditions, distributed loads, and visualization. You can define these elements in a similar fashion as you would define point masses or rotary inertia. The coordinates of the node of a particle correspond to the center location of the physical grain of material. PD3D elements are assigned to a discrete section definition, where particle characteristics are specified. For more information, see Discrete Particle Elements.

A particle tracking box is established at the beginning of the analysis to define the rectangular region within which the particle search (finding all neighbors for all particles) is performed. You can modify the size of the box as discussed in Using Section Controls to Define the Particle Tracking Box for DEM and SPH Particles.

Input File Usage

Use the following options to define a discrete element medium:

ELEMENT, TYPE=PD3D, ELSET=particle_body
element number, node number

DISCRETE SECTION, ELSET=element_set_name

Constraints

Since PD3D elements are Lagrangian elements, their nodes can be involved in other features such as constraints. Although the PD3D element has a spherical shape, it is possible to model grains of complex shapes by clustering particles together, as illustrated in Figure 3. A cluster is a group of particles that are held together either rigidly or via compliant connections.

The particles in a cluster may overlap with each other. Contact forces that try to push apart overlapping particles of a cluster will exist unless contact exclusions are specified among particles of a cluster (see Specifying Contact Exclusions). These contact forces will have no effect on particles held together rigidly.

The particle-cluster approach may not replicate the precise geometry of actual grains. For example, the cluster shown in Figure 3 may approximate an ellipsoidal shape (indicated by the dashed line in the figure). More spherical particles of various sizes can be added to the cluster to obtain a closer approximation of the true shape.

Define BEAM-type multi-point constraints between nodes of a group of particles to create a rigid cluster. Clusters of overlapping particles that do not involve multi-point constraints may exhibit nonphysical behavior. For more information, see General Multi-Point Constraints.

Interactions

Contact is an essential ingredient for DEM analyses, as discussed above. General contact is used to define contact involving DEM particles. A DEM particle can be involved in multiple contact interactions simultaneously with

another particle with the same discrete section definition;
another particle with a different discrete section definition;
a surface based on finite elements; and
an analytical rigid surface.

Modeling contact between DEM particles requires that the particles be explicitly included in general contact as element-based surfaces using contact inclusions (see Element-Based Surface Definition). See About General Contact in Abaqus/Explicit for a discussion of general contact. By default, the particles are not part of the general contact domain, similar to other 1-node elements (such as point masses).

Contact stiffness for DEM is often intended to account for the physical stiffness characteristics of the particles because DEM models each particle as rigid; therefore, nondefault contact property assignments are common for DEM interactions.

Normal and Tangential Contact Forces

Figure 4 is a schematic representation of the contact stiffness and damping between two particles. The spring stiffness $K_{n}$ acts in the contact normal direction and may represent a simple linear or a nonlinear contact stiffness. The dashpot $C_{n}$ represents contact damping in the normal direction. The tangential spring stiffness $K_{t}$ along with the friction coefficient $μ$ represent friction between the particles. The dashpot $C_{t}$ represents contact damping in the tangential direction.

Figure 4 shows that the tangential contact forces acting on particle surfaces cause moments at particle centers. Interactions involving DEM particles account for moment transfer across the interface.

Hertz Normal Contact between Particles

The Hertz elastic solution relating contact force, $F$ , to the approach distance, $δ,$ for two contacting spherical particles is:

F = \frac{4}{3} E^{*} \sqrt{R} \sqrt{δ^{3}},

where

R = \frac{R_{1} R_{2}}{R_{1} + R_{2}}

and

\frac{1}{E^{*}} = \frac{1 - ν_{1}^{2}}{E_{1}} + \frac{1 - ν_{2}^{2}}{E_{2}} .

$E_{1}$ , $ν_{1}$ and $E_{2}$ , $ν_{2}$ are the effective Young's modulus and Poisson's ratio of the two contacting particles, respectively. $R_{1}$ and $R_{2}$ are the radii of the two contacting particles, respectively. The normal contact stiffness is $K_{n} = 2 E^{*} \sqrt{R} \sqrt{δ}$ . You must specify the effective Young's modulus and Poisson's ratio for a contacting particle. In addition, you must specify the Hertz-type pressure overclosure.

$K_{m a x}$ is the limiting value of the normal contact stiffness, $K_{n}$ , beyond which the normal contact force increases linearly using the slope of the $K$ – $δ$ curve at $K_{m a x}$ .

Input File Usage

Use the following options to define Hertz-type normal contact behavior:

DISCRETE SECTION, ELSET=element_set_name
DISCRETE ELASTICITY
effective Young's Modulus, effective Poisson's ratio
...
SURFACE INTERACTION
SURFACE BEHAVIOR, PRESSURE-OVERCLOSURE=HERTZ
 $K_{m a x}$

Johnson-Kendall-Roberts Adhesive Normal Contact between Particles

The JKR model relates contact force, F, to the contact area, a, as follows:

F = \frac{4 E^{*} a^{3}}{3 R^{*}} - \sqrt{8 π Γ E^{*} a^{3}} .

The approach distance, $δ,$ for two contacting spherical particles is related to the contact area, a, as follows:

δ = \frac{a^{2}}{R^{*}} - \sqrt{\frac{2 π Γ a}{E^{*}}} .

In the above equations

R = \frac{R_{1} R_{2}}{R_{1} + R_{2}}

and

\frac{1}{E^{*}} = \frac{1 - ν_{1}^{2}}{E_{1}} + \frac{1 - ν_{2}^{2}}{E_{2}} .

In the above equations $Γ$ is the surface energy per unit area of the two contacting particles. Like nonadhesive Hertz contact $E_{1}$ , $ν_{1}$ and $E_{2}$ , $ν_{2}$ are the effective Young's modulus and Poisson's ratio of the two contacting particles, respectively. $R_{1}$ and $R_{2}$ are the radii of the two contacting particles, respectively. You must specify the effective Young's modulus and Poisson's ratio for a contacting particle. In addition, you must specify the JKR-type pressure overclosure.

$K_{m a x}$ is the limiting value of the normal contact stiffness, beyond which the normal contact force increases linearly.

Figure 5 shows the force-displacement curve for the JKR adhesion model. Adhesion between particles is triggered upon first contact. The pull-off force, $F_{c} = - \frac{3}{2} π Γ R^{*}$ , is the maximum tensile force. Particles continue to experience adhesion force even after physical separation occurs. The adhesion force between two particles goes to zero at a separation distance of

δ_{s e p a r a t i o n} = - \frac{3}{4} {(\frac{π^{2} Γ^{2} R^{*}}{E^{*}^{2}})}^{\frac{1}{3}} .

It can be seen from Figure 5 that adhesive forces are nonzero at $δ = 0$ and reduce to zero at $δ = δ_{F 0}$ . In some situations it may be desirable to have zero adhesive forces when $δ = 0$ . You can achieve this by applying a horizontal shift of $δ = δ_{F 0}$ to the force-displacement curve shown in Figure 5. Figure 6 shows the shifted curve. The pull-off force remains unchanged due to the shift, whereas the new separation distance $δ_{s e p a r a t i o n_s h i f t e d} = δ_{s e p a r a t i o n} + δ_{F 0}$ . Abaqus automatically computes the horizontal shift $δ_{F 0}$ for the “shifted” JKR type adhesive behavior.

Input File Usage

Use the following options to define JKR-type normal contact behavior:

DISCRETE SECTION, ELSET=element_set_name
DISCRETE ELASTICITY
effective Young's modulus, effective Poisson's ratio
...
SURFACE INTERACTION
SURFACE BEHAVIOR, PRESSURE-OVERCLOSURE=JKR
Surface energy per unit area,  $K_{m a x}$

Use the following options to define the shifted JKR-type normal contact behavior:

DISCRETE SECTION, ELSET=element_set_name
DISCRETE ELASTICITY
effective Young's modulus, effective Poisson's ratio
...
SURFACE INTERACTION
SURFACE BEHAVIOR, PRESSURE-OVERCLOSURE=JKR
Surface energy per unit area,  $K_{m a x}$ , SHIFTED

Output

No element output is available for PD3D elements. The nodal output includes all output variables generally available in Abaqus/Explicit analyses (see Abaqus/Explicit Output Variable Identifiers).

Limitations

Discrete element method analyses are subject to the following limitations:

Volume average output for stress, strain, and other similar continuum element output is not available for DEM analysis.
Only a spherical shape is supported for PD3D elements.
PD3D elements cannot be part of a rigid body definition.
Models with rigid clusters using PD3D elements and BEAM-type multi-point constraints cannot be run in domain parallel mode.
It is not possible to specify cohesive or thermal contact between PD3D elements or between PD3D elements and other elements.
Rolling friction is ignored for contact between PD3D elements or between PD3D elements and other elements.
User-defined surface interaction is not supported for contact between PD3D elements.

DEM computations are distributed across parallel domains except when multiple discrete sections are defined with different alpha damping parameters (which will degrade parallel scalability). DEM analyses are subject to the following limitations if multiple CPUs are used:

Contact output is not supported for DEM secondary nodes.
Energy history output other than for the whole model is not supported.
Dynamic load balancing cannot be activated.
If any DEM particles participate in general contact, all DEM particles must be included in the general contact definition.
At least 10,000 DEM particles per domain is suggested to achieve good scalability.
A significant increase in memory usage may be needed if a large number of CPUs is used.

Input File Template

The following example illustrates a discrete element method analysis:

HEADING
…
ELEMENT, TYPE=PD3D, ELSET=name
Element number, node number
…
DISCRETE SECTION, ELSET=name
Particle radius
**
INITIAL CONDITIONS, TYPE=VELOCITY
Data lines to define velocity initial conditions
NSET, NSET=name, ELSET=name
SURFACE, NAME=name
,
**
SURFACE INTERACTION
**
STEP
DYNAMIC, EXPLICIT
DLOAD
Data lines to define gravity load
**
CONTACT
CONTACT INCLUSIONS
CONTACT PROPERTY ASSIGNMENT
**
CONTACT CONTROLS ASSIGNMENT
OUTPUT, FIELD
END STEP

References

Cundall, P. A., and O. D. Strack, “A Distinct Element Method for Granular Assemblies,” Geotechnique, vol. 29, pp. 47–65, 1979.
Johnson, K. L., K. Kendall, and A. D. Roberts, “Surface Energy and the Contact of Elastic Solids,” Proceedings of the Royal Society of London, vol. 324, pp. 301–313, 1971.
Munjiza, A., and K. R. F. Andrews, “NBS Contact Detection Algorithm for Bodies of Similar Size,” International Journal for Numerical Methods in Engineering, vol. 43, pp. 131–149, 1998.
O'Sullivan, C., and J. D. Bray, “Selecting a Suitable Time Step for Discrete Element Simulations that Use the Central Difference Time Integration Scheme,” Engineering Computations, vol. 21(2/3/4), pp. 278–303, 2004.
Zhu, H. P., Z. Y. Zhou, R. Y. Yang, and A. B. Yu, “Discrete Particle Simulation of Particulate Systems: A Review of Major Applications and Findings,” Chemical Engineering Science, vol. 63, pp. 5728–5770, 2008.
Zhu, H. P., Z. Y. Zhou, R. Y. Yang, and A. B. Yu, “Discrete Particle Simulation of Particulate Systems: Theoretical Developments,” Chemical Engineering Science, vol. 62, pp. 3378–3396, 2007.