Discrete Element Method (DEM) is a numerical technique for predicting the behavior of bulk solids.

A bulk solid
is a large collection of solid particles; otherwise known as granular media. Some examples of granular media
flows include grains being moved through processing equipment, ore being passed through mining machinery,
and sand falling through an hourglass. Granular media flow can be quite complex as these flows are known to
exhibit solid-like, fluid–like, or a combination of both behaviors. For example, sand in an hourglass behaves like a
fluid while a stockpile of sand can have a solid-like stress-strain response.

DEM is a mesh-free method and

does not solve the continuum
equations of motion. Hence, no
stress-strain constitutive law for
the material is needed. Instead,
a stress-strain relationship can
be obtained as an output from
the DEM model. A general DEM
algorithm is shown in Figure 1.

The equations of motion for

every individual particle are
numerically integrated with
time. For this process, the total
force on a particle needs to be
known. The total force is the
resultant of contact forces
(between particles and with the
boundary) and body forces.
Typical body forces are gravity
(weight), fluid and other forces,
such as electrostatic, electro-
magnetic, and so on.

Presented in this document is a

description of the contact forces
used by Rocky, which is a
software package developed by
ESSS Rocky, S.L. Rocky is a very
actively evolving package with
many new models implemented
on a near monthly basis. This
document is a general guide for
the latest developments of
these Rocky models.
Figure 1: General DEM Algorithm Flowchart

Copyright  2017, ESSS Rocky, S.L. All Rights Reserved. Page 1

The main sections in this document include the following:
 Contact Forces in Rocky
 Heat Transfer in Rocky
 Simulation Timestep Calculations

The contact forces in any DEM code (including Rocky) consist of the following two parts:

 Normal-to-contact plane
 Tangential-to-contact plane

For spherical particles, the contact plane is perpendicular to the line that connects the centers of two spheres. In
the case of particle-to-boundary contact, the line connects the center of a sphere and the closest point of a
triangle. For non-round particles, the algorithm for determining a contact plane is more complex; it involves
calculating one of any of the following:

 the closest points of two particles, or

 the closest points of a particle and triangle, or
 the two points with the maximum overlap distance in the case of a physical contact

The contact plane is a plane perpendicular to the line connecting these closest points. The descriptions for the
normal and tangential force calculations used in Rocky, as well as the models used by these forces are presented
below.

The normal force model for a DEM simulation has two major requirements. First, the force has to be a repulsive
one. Second, the normal contact force model has to allow significant energy dissipation since a granular medium
is an extremely dissipative system. A number of models have been proposed for these purposes. The models
implemented in Rocky are discussed below.

The first elastic-plastic (repulsive and dissipation) normal contact model implemented in Rocky is a linear
hysteresis model. This model allows simulation of the plastic energy dissipation on a contact without introducing
the overhead of long simulation times. In addition, since no viscous damping term is used, the energy dissipation
is not dependent on the relative velocities of neighboring particles, making the energy dissipation insensitive to
other contacts. An additional advantage of this model is that compressible materials can be accurately modeled
due to the fact that the contact forces can be almost zero even at residual overlaps.

Copyright  2017, ESSS Rocky, S.L. All Rights Reserved. Page 2

The linear spring hysteresis model as implemented in Rocky can be described by the following set of equations:

Where:

 is the normal elastic-plastic contact force at the current time, t

 is the force at the previous time, , where is the timestep
 is the change in the contact normal overlap during the current time (it is assumed to be positive as
particles approach each other and negative when they separate)
 and are the normal overlap values at the previous and the current time
 and are the values of loading and unloading contact stiffnesses respectively
Note: The unloading force is limited by value of to ensure that force remains zero with
zero overlap

Let us discuss the definition of contact stiffnesses for loading and unloading parts of the normal contact force.
The loading and unloading stiffnesses are defined by the particle size, bulk Young’s Modulus, and by the
restitution coefficient of contacting materials, the last two of which the user inputs into Rocky. The coefficient of
restitution ε in Rocky is a measure of energy dissipation for the contacting pair of materials. The particle's
loading stiffness is calculated in the following manner:

Where:

 is the particle material bulk Young’s Modulus or Elastic Modulus (user input). A value of 1.0x108 Pa is
a reasonable number for practical particle sizes.
 is the boundary material Young’s Modulus (user input).
 is the particle size.

Copyright  2017, ESSS Rocky, S.L. All Rights Reserved. Page 3

For the contact of two particles, or of a particle with a boundary, the loading and unloading stiffnesses are
defined respectively as:

or par le par le onta t

or par le boundary onta t

Where subscripts _1 and _2 refer to particle 1 and particle 2 of two contacting particles.

The elastic-plastic normal contact model is appropriate for simulating non-adhesive and dry granular materials.
Unfortunately, quite often granular materials can have various moisture contents that can cause them to adhere
to themselves and other surfaces with which they come into contact. In order to capture this behavior, the
repulsive normal force has to be supplemented by the attractive normal force to accurately predict its flow
characteristics. For example, this is required for wet and sticky materials where the adhesion is from liquid
bridge forces or other mechanisms. The models implemented in Rocky for this purpose are discussed below.

One important point to note is that cohesion/adhesion of a bulk solid is a function of the stress. The linear force
adhesion contact model in Rocky captures this physical phenomenon by scaling the cohesion/adhesion with the
force of contact.

With proper calibration the constant force adhesive model can be used to simulate the behavior of adhesive
materials that do not exhibit stress consolidation effects (e.g., when adhesion is due to liquid bridge forces). This
is a simplest adhesion model in Rocky and includes the following two parameters:

 The minimum distance between particles or particle-boundary surfaces when the force is activated
 The value of this force is expressed as a fraction of the particle gravity force

If the force fraction is 1.0, that means the adhesive force will be equal to the gravity force applied to the
particle. In the case of contact between two particles of different mass, the smallest mass is considered for
gravity force calculations.

The constant adhesive force model as implemented in Rocky can be described by the following set of equations:

Copyright  2017, ESSS Rocky, S.L. All Rights Reserved. Page 4

Where:

 is the normal adhesive contact force at the current time, t

 is the contact normal overlap during the current time (it is assumed to be positive as particles
approach each other and negative when they separate)
 and are the mass of the particles in contact
 is the gravity acceleration
 is the force fraction

Like the constant force model, the linear force adhesive model requires calibration of the material behavior. This
model is appropriate for granular materials with stress consolidation. Like the constant force model, the linear
force model is also defined by two parameters:

 The minimum distance between particles or particle-boundary surfaces when the force is activated
 The ratio of the adhesive force stiffness to the contact loading stiffness

The linear adhesive model, in essence, behaves like an attractive linear spring; the value of this force is zero if
the contact distance is larger than the minimum adhesive contact force distance and increases proportional to

the difference between this distance and actual particles contact distance/overlap with the coefficient of
proportionality defined by and the above mentioned ratio.

The linear adhesive force model as implemented in Rocky can be described by the following set of equations:

Where:

 is the normal adhesive contact force at the current time, t

 is the contact normal overlap during the current time (it is assumed to be positive as particles
approach each other and negative when they separate)
 is the stiffness ratio
 is the values of loading contact stiffnesses

Figure 2 shows a force-overlap plot of a particle collision with a wall for both dry and adhesive impacts. It is
evident that the larger the contact overlap the larger the attractive force, thus allowing an accurate model to
represent granular material stress consolidation.

Copyright  2017, ESSS Rocky, S.L. All Rights Reserved. Page 5

 Figure 2: Force-Overlap Plot

The tangential force in Rocky is elastic-frictional. It is described by the following equations:

If no sliding is taking place for the contact, and

If sliding is taking place for the contact

Where:

 is the tangential contact force at the current time

 is the same force at the previous time
 is the tangential relative particles’ displacement during the current time
 and are the static and dynamic friction coefficients

The sliding is considered to be taking place on the contact the first time the tangential force exceeds the limit of
. Once the tangential force falls below the value of , the contact is considered non-sliding again.

Copyright  2017, ESSS Rocky, S.L. All Rights Reserved. Page 6

If the thermal model is activated in Rocky, an additional equation for the energy balance is solved along with the
equations governing the motion of the particle. In the current implementation, the particle temperature, , is
assumed to be uniform; i.e., no radial or circumferential temperature variation is admitted. This approximation
is reasonable for small or highly conductive particles.

The temperature variation of a particle can be obtained over time according to:

Where is the particle mass, is the particle specific heat and is the total particle heat transfer rate. This
heat transfer rate accounts for the heat transfer that occurs during the contact with other particles or walls, ,
and the heat transfer between particle and fluid phase, , if coupling with a CFD solver, according to:

The heat exchanged due to contact with other particles or walls is calculated according to:

Where is the number of collisions, is the heat exchange between particle and a particle or wall ,
calculated using a linear model as:

where is the contact conductance, calculated according to:

where is the effective contact thermal conductivity, calculated as:

The contact diameter, , is considered much smaller than the diameter of the contacting particles and and
is calculated according to:

Copyright  2017, ESSS Rocky, S.L. All Rights Reserved. Page 7

where is the contact normal force between particles and , is the geometric mean of the colliding
particles diameters and is the contact equivalent Young's modulus, calculated as:

where is the Poisson ratio of the particle or wall or, in the case of collision against wall, .

The timestep determines the amount of time required for the completion of the simulation. The simulation
timestep should be large enough to ensure reasonable completion time but small enough to ensure the
accuracy and stability of the simulation. The description of the timestep calculation procedure is described
below.

The timestep for the linear hysteresis normal force model is defined by material stiffness, effective mass of
contacting pair, and number of timesteps per loading cycle. The formula for calculating timestep is:

Where:

 is the effective mass between the contact

 is the minimum number of timesteps per loading cycle (user inputs this Loading N-Steps value into
Rocky)

The default value for in Rocky is 15. To ensure the accuracy of the calculations, it is recommended that this
Loading N-Steps value be at least 10.

Copyright  2017, ESSS Rocky, S.L. All Rights Reserved. Page 8