Topic7 ContactStressesAndDeformations

INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING
Int. J. Numer. Meth. Engng 2017; 110:186–200

Published online 22 September 2016 in Wiley Online Library (wileyonlinelibrary.com). DOI: 10.1002/nme.5361
Establishing stable time-steps for DEM simulations of

non-collinear planar collisions with linear contact laws
Shane J. Burns*,† and Kevin J. Hanley

Institute for Infrastructure and Environment, School of Engineering, The University of Edinburgh, Edinburgh EH9 3JL,
UK
SUMMARY
The discrete element method, developed by Cundall and Strack, typically uses some variations of the central
difference numerical integration scheme. However, like all explicit schemes, the scheme is only condition-
ally stable, with the stability determined by the size of the time-step. The current methods for estimating
appropriate discrete element method time-steps are based on many assumptions; therefore, large factors of
safety are usually applied to the time-step to ensure stability, which substantially increases the computa-
tional cost of a simulation. This work introduces a general framework for estimating critical time-steps for
any planar rigid body subject to linear damping and forcing. A numerical investigation of how system damp-
ing, coupled with non-collinear impact,
q affects the critical time-step is also presented. It is shown that the
critical time-step is proportional to m k
if a linear contact model is adopted, where m and k represent mass
and stiffness, respectively. The term which multiplies this factor is a function of known physical parame-
ters of the system. The stability of a system is independent of the initial conditions. © 2016 The Authors.
International Journal for Numerical Methods in Engineering Published by John Wiley & Sons Ltd.
Received 2 June 2016; Revised 17 August 2016; Accepted 23 August 2016
KEY WORDS: DEM; time-step; stability
1. INTRODUCTION
Contact mechanics primarily involves the analysis of the impulses and reaction forces that develop
as a result of two or more bodies colliding or being in contact, and the governing dynamics which
ensues [1–4]. The main objectives of contact mechanics are to develop and understand methods for
calculating the changes in motion, for the case of a collision between bodies [5], and to calculate
interaction forces for bodies already in contact [6]. The understanding of the mechanics of contact
has far-reaching applications in areas such as civil engineering, mechanical engineering, vehicle
collision analysis and sport science [7].
One such method for modelling contact problems is the discrete element method (DEM) [8].
DEM is typically used to describe the behaviour of bulk granular materials in which there is a large
number of particles and multiple simultaneous contacts [9]. The bodies are assumed to be rigid,
but a small amount of overlap at the contact region is permitted. The impact process is modelled
using a combination of springs, dampers and sliders. Using Newton’s second law together with
Euler’s second law, the contact process can be described mathematically as a system of second-order
differential equations (ODEs).
*Correspondence to: Shane J. Burns, Institute for Infrastructure and Environment, School of Engineering, The University
of Edinburgh, Edinburgh EH9 3JL, UK.
† E-mail: [email protected]
This is an open access article under the terms of the Creative Commons Attribution License, which permits use,
distribution and reproduction in any medium, provided the original work is properly cited.
© 2016 The Authors. International Journal for Numerical

Methods in Engineering Published by John Wiley & Sons Ltd.
ESTABLISHING STABLE TIME-STEPS FOR DEM SIMULATIONS 187
It is necessary to choose a numerical integration procedure to solve the equations of motion

for the system. Most developments in the area of numerical integration have come from the finite
element method, which, like DEM, requires a numerical solution to a system of equations. The most
commonly used schemes are the central difference algorithm, Position–Verlet and Gear’s Predictor-
Corrector. A comprehensive numerical comparison of these time integration schemes, in terms of
accuracy, stability and performance, is given by Rougier et al. [10]. All of these time integration
schemes are of second-order accuracy and can be formulated using a variable or fixed time-step
t [11]. For the analysis in this work, only the constant time-step formulation will be considered.
Like all explicit numerical integration procedures, the scheme is only conditionally stable, with the
stability dependent on the size of t . It is well known that too small a time-step can slow down the
simulations unnecessarily and thereby increases the computational cost whereas too large a time-
step can lead to numerical instability. This type of instability is often manifested by a non-physical
generation of energy [12]. The methods being used at present to estimate appropriate DEM time-
steps are based on many assumptions [13], some of which lack physical justification. The different
approaches often lead to greatly differing estimates of the critical time-step. To account for the
numerous assumptions made in the estimation, large ‘factors of safety’ are applied to the calculated
time-steps, which then leads to the adoption of unnecessarily conservative time-steps.
There are two commonly used approaches for calculating the critical time-step tcri t . The first,
which calculates tcri t as a function of mass m and stiffness k, typically considers each degree of
freedom as a system of uncoupled point masses. The eigenvalues of the system matrix are deter-
mined and used in the criteria for calculating tcri t [14]. O’Sullivan and Bray extended this analysis
by drawing an analogy between the nodes of a finite element mesh and particles of a DEM simula-
tion [15]. They were able to calculate critical time-steps for regular packings of uniform disks (2D)
and spheres (3D) using this analogy. This work however has not been rigourously extended to 2D
or 3D collisions and has not been extended to allow for damped collisions.
The second method, the Rayleigh criterion, is based on the idea that energy cannot propagate from
a particle beyond its immediate neighbouring particles in a single time-step [16]. The assumption is
made that all energy transferred across a particulate system is due to Rayleigh waves [17], and thus,
the contributions of dilational and distortional waves can be neglected. The critical time-step is then
calculated using the theoretical expression for the Rayleigh wave velocity for the system in question.
Both approaches have been shown to calculate time-step values differing by orders of magnitude
for certain cases. This would lead to an unnecessarily conservative time-step if the wrong criterion
were used and would drastically slow down the simulations. An additional shortcoming of the cur-
rent approaches is that they assume perfect symmetry in the particles under consideration. Most real
particles are not perfect discs or spheres and have a complex, irregular geometry. Because particle
shape strongly influences the dynamics of granular systems, there has recently been growing interest
in adopting non-spherical particle shapes in DEM [18], using approaches such as potential particles
[19] or multi-sphere clumps [20]. At present, there are no guidelines available for users to choose
a time-step for DEM simulations of non-symmetric particles. It is clear therefore that the methods
being used currently to calculate suitable time-steps are inadequate, particularly if more complex,
non-symmetric geometries are incorporated. As DEM continues to grow in popularity, it is very
important to improve upon these approaches to ensure stability of simulations without unnecessary,
inefficient conservatism. The overall objective of this work is to improve upon existing guidelines
for choosing appropriate time-steps for DEM simulations using linear contact models. For this pur-
pose, we will introduce an alternative way of selecting a time-step, which builds upon the current
‘state of the art’ methodology, using a mathematically rigorous technique which is capable of incor-
porating any planar geometry together with system damping. This work acts as a guide and reference
for DEM users on how to select a stable time-step for more complex planar geometries, and it also
gives a mathematically exact analysis of planar collisions with linear damping and forces.
This article is organised as follows. In Section 2, we will discuss the current methodology for
estimating the critical time-step as a function of mass and stiffness, and we will further compare the
predicted critical time-step obtained using this technique with the natural frequency of oscillation
of the system. A general framework for calculating the natural frequency of oscillation of a non-
symmetric system is derived in Section 3 and is used to estimate critical time-steps for a number of
© 2016 The Authors. International Journal for Numerical Int. J. Numer. Meth. Engng 2017; 110:186–200
Methods in Engineering Published by John Wiley & Sons Ltd. DOI: 10.1002/nme
188 S.J. BURNS AND K.J. HANLEY
planar cases in Section 4. In Section 5, we carry out a numerical investigation of the changes in the
critical time-step under parameter variation. The paper concludes with a discussion in Section 6 that
provides an insight for engineers and other researchers working with DEM.
2. STABILITY AS A FUNCTION OF MASS AND STIFFNESS
In this section, we will introduce and discuss one of the most frequently used techniques for selecting
a stable time-step.
2.1. Stability using Rayleigh’s theorem

In this paper, when we refer to current methodologies, we are referring to the approach whereby the
critical time-step is calculated as a function of the mass and stiffness of the system particles. This
commonly used approach makes use of a corollary of Rayleigh’s theorem, the essential details of
which will be given here.
Consider the following second-order system of differential equations
M qR C C qP C Kq D F; (1)
where M is the mass matrix, C is the damping matrix, K is the stiffness matrix, F is the column
vector of external forces and q, qP and qR are the displacement, velocity and acceleration vectors,
respectively. Belytschko [14] derives the stability criterion for the discretised form of (1) by making
use of modal decomposition to reduce the system to a single degree of freedom system. The
maximum stable time-step tcri t is a function of the eigenvalues of the amplification matrix and is
given by
2
tcri t D ; (2)
!max
for a linear, undamped system [11, 15]. Belytschko [14] gives the following equation
p
!max D max ; (3)
where max is the maximum eigenvalue of
M 1 K: (4)
Equation (3) can be derived using an extension of Rayleigh’s bounding theorem [14], which relates
the eigenvalues of any two systems which are equivalent except for linear constraints. It is important
to note, however, that this condition was derived for a second-order system of the form
xR D x: (5)
Recasting this system as a system of first-order differential equations gives the following equivalent
relation
!max 2 D max : (6)
A major shortcoming of this analysis is that it requires that the modal equations of motion will
decouple. In order for this to be possible, it is necessary to impose certain restrictions on C , which
are often unphysical. The method employed in [14] assumes Rayleigh damping in that the damping
matrix C is directly defined as a linear combination of the mass and stiffness matrices. The main
consequence of this is that it may not always be correct to use (2) to calculate the critical time-step
for a system in which damping is present. This method has also not been developed to allow for
collisions of non-symmetric bodies.
2.2. Natural frequency of the system

The principle of natural frequency of a system can be used to theorise numerical stability for the
discretised form of the system [21]. The time at which the compression phase ends, tc , can be
calculated in terms of the natural frequency of oscillation of the system [1] and is given by

tc D ; (7)
2!
where ! is the natural frequency of oscillation with respect to the normal or tangential direction.
Equation (7) is valid only when the tangential and normal dynamics are decoupled, i.e., for collinear
collisions. If it is assumed that all forces and deformations occur around the contact point, it is
also intuitive to consider the natural frequency of oscillation at a frame of reference located at the
contact point [1]. Choosing a time-step greater than (7) could lead to particles passing through
each other without an impact being detected or could lead to an excessively large overlap between
particles, which, as a result, would give an unphysical reaction force leading to instability in the
system.
The time at which the compression period ends, given by (7), is almost the same as the critical
time-step predicted by (2). It makes sense physically that one should not choose a time-step greater
in value than the theoretical time for compression. As discussed in Section 2.1, a number of assump-
tions are necessary in order to derive (2). We propose using (7) as an alternative way of calculating
critical time-steps such that

t D : (8)
2!max
Equation (8) does not require any restrictions on the system damping other than linearity with
respect to the generalised coordinates. For non-collinear collisions, (8) will not be exact. However,
it serves as an adequate bound for estimating the critical time-step. It is also important to highlight
the physical meaning of the critical time-step, and the link between natural frequency of oscillation
and theoretical time for compression facilitates this. In the next section, we will present a general
framework for deriving the natural frequency of oscillation for any planar rigid body.
3. GENERAL FRAMEWORK
Consider the two bodies H and H 0 in free flight as shown in Figure 1(a) and in contact at point P as
shown in Figure 1(b). The position and rotation of the centre of mass G of body H can be described
in the n1 n2 plane by the coordinates q1 and q2 and the angle , and similarly, q10 and q20 are the
coordinates and 0 the rotation of the centre of mass G 0 of body H 0 (Figure 1(b)).
We let
q D .q1 ; q2 ; /T ; q 0 D .q10 ; q20 ; 0 /T ;
T T
qP D qP 1 ; qP 2 ; P ; qP 0 D qP 10 ; qP 20 ; P 0 ;
T T
qR D qR 1 ; qR 2 ; R and qR 0 D qR 10 ; qR 20 ; R 0
be the displacement, velocity and acceleration vectors of the centre of mass of H and H 0 , respec-
tively. In order to translate the contact forces to the centre-of-mass reference frame, we need to
consider

@qP 1 0 L sin./ @qP 1 0 L0 sin. 0 /
D and D ;
@q 0 1 L cos./ @q 0 0 1 L0 cos. 0 /
Figure 1. (a) A free-body force diagram for two planar rigid bodies in free flight. (b) Geometry of two planar
rigid bodies H and H 0 in contact at a point P .
where L is the distance from G to P and where L0 is the distance from G 0 to P . Further, we can
define e and e 0 as the eccentricity of collision such that
e D L cos./ and e 0 D L0 cos. 0 /:
Then, using Newton’s second law, the equations of motion for the two bodies at contact can be
written as

@qP T T
M qR D F T C (9)
@q
and
T
T @qP T
M 0 qR 0 D F 0 C 0 ; (10)
@q 0
where M and M 0 are, respectively, the mass matrices for H and H 0 given by
0 1 0 0 1
m 0 0 m 0 0
M D @ 0 m 0 A and M 0 D @ 0 m0 0 A ;
0 0 I 0 0 I0
and where the external forces, F and F 0 , which respectively act on H and H 0 , are given by
F D .F1 ; F2 ; R/ ;
and

F 0 D F10 ; F20 ; R0 :
The subscripts 1 and 2 respectively represent the components of the external forces acting in the
tangential and normal directions, and R and R0 are the external torques acting on the bodies as
shown in Figure 1(a). Similarly, we define and 0 as the forces generated at impact of each body
given by
D .T ; N / and 0 D .T ; N / :
It is important to note that this set-up is general and does not specify the mechanism that generates
the normal and tangential forces, N and T , respectively. In the context of the analysis presented
in this article, when we refer to external forces, we are referring to all forces which do not feature
in the contact force .
As discussed in Section 2.2, it is intuitive to use a reference frame located at the contact point
between the contacting rigid bodies. For this purpose, we consider (9) translated to the contact point
at coordinates qP so that

d qP P @qP @qP 1 @qP T T @qP 1 T
D qR D M C M F D ´1 T C f .F1 ; F2 ; R; q; q/ P (11)
dt @q @q @q @q
and
T
d qP P @qP @qP @qP @qP T T T
D 0 qR D 0 .M 0 /1 0 C
.M 0 /1 F 0 D .´0 /1 0 C f 0 F10 ;F20 ;R0; q 0; qP 0 :
dt @q @q @q 0 @q 0
(12)
This is the rate of change of the contact point velocities as a function of time, in which ´1 and
.´0 /1 are the symmetric matrices given by
T
1 @qP 1 @qP A B
´ D M D ;
@q @q B D
T
0 1 @qP @qP A0 B 0
.´ / D .M 0 /1 D ;
@q 0 @q 0 B 0 D0
where
1 L2 sin2 ./ L2 sin./ cos./ 1 L2 cos2 ./

AD C ; BD ; DD C
m I I m I
1 L0 2 sin2 . 0 / L0 2 sin. 0 / cos. 0 / 1 L0 2 cos2 . 0 /

A0 D C ; B0 D ; D0 D C :
m0 I0 I0 m0 I0
For the analysis of the system at contact, we do not need to consider terms which do not change
throughout the impact phase, and therefore, we can neglect the functions f and f 0 ([22]). This
framework is completely general and can be used with any planar geometry, provided the moment
of inertia is known, together with any contact law which is linear with respect to the generalised
coordinates.
3.1. Linear contact interaction

We model the contact process using linear springs and dampers acting in the normal and tangential
directions (see Figure 2). Let k1 and c1 be the spring and damping constants for the components
acting in the tangential direction and let k2 and c2 be the spring and damping constants for the
components acting in the normal direction. The contact forces and 0 are then given by

PQ ; k2 qQ P 2 C c2 qPPQ 2
D k1 qQ P1 C c1 qP1 PQ ; k2 qQ P 2 c2 qPPQ 2 ;
and 0 D k1 qQ P1 c1 qP1
where the particle overlap is given by the relative displacement of the contact points and is given by
qQ P1 D qP0 qP :
Defining the relative acceleration of the contact points of H and H 0 as qRQ P WD qR P0 qR P , the
equations of motion now take the form

qP1RQ k1 qQ P1 c1 qQP P1
D Ó C Ó ; (13)
qPQR 2 k2 qQ P 2 c2 qQP P 2
Figure 2. Two planar rigid bodies H and H 0 with initial velocities V and V 0 in contact at a point P together
with a magnification of the contact region.
where

AO BO
Ó D ´1 C .´0 /1 D
BO DO
for
1 L2 sin2 ./ 1 L0 2 sin2 . 0 / L2 sin./ cos./ L0 2 sin. 0 / cos. 0 /
AO D C C 0C ; BO D ;
m I m I0 I I0
1 L2 cos2 ./ 1 L0 2 cos2 . 0 /

DO D C C 0C :
m I m I0
This general framework can be used for any planar rigid-body provided the moment of inertia,
together with the orientation terms given in A, B and D, are known. In the following section, we
will derive the general solution to (13) in order to determine the numerical stability for this system.
4. TIME-STEP CALCULATIONS
In this section, we will discuss the general solution to the ODE system derived in Section 3.1 and
discuss what features of the solution affect the stability of the discretised system. We will then
calculate analytical expressions for the critical time-step for a number of cases.
In order to solve (13), it is first necessary to recast the second-order system as a system of first-
order ODEs. We let x1 WD qQ P1 , x2 WD qPQ P1 , x3 WD qQ P 2 and x4 WD qPQ P 2 . Equation (13) can now be
recast as the following system
0 1 0 1
xP 1 x1
B 2C
xP B 2C
x
@ xP A D Q @ x A ; (14)
3 3
xP 4 x4
where the system matrix

0 1
0 1 0 0
B k1 AO c1 AO k2 BO c2 BO C
Q WD B
@ 0
C (15)
0 0 1 A
k1 BO c1 BO k2 DO c2 DO
subject to the initial conditions

T
x0 D qQ P1;0 ; qPQ P1;0 ; qQ P 2;0 ; qPQ P 2;0 ; : (16)
The general solution to (14) is determined by calculating the eigenvalues and eigenvectors of (15)
and using the initial condition vector (16) to give
0 1
x1 .t /
B x2 .t / C t t t t
@ x .t / A D ˛1 e 1 V1 C ˛2 e 2 V2 C ˛3 e 3 V3 C ˛4 e 4 V4 ; (17)
3
x4 .t /
where 1 , 2 , 3 , 4 , and V1 , V2 , V3 , V4 are the eigenvalues and corresponding eigenvectors of

(15), respectively. It is important to note that the eigenvalues and eigenvectors are not a function of
the initial conditions (16). The initial conditions form part of the constants ˛1 , ˛2 , ˛3 , ˛4 , which
scale the solution. As discussed in Section 2, the stability condition is only a function of the system
eigenvalues. This means that for linear contact interactions, the stability of the discretised system is
entirely independent of the initial conditions, and further that the natural frequency does not depend
on the initial velocity at contact, and furthermore, the compression time is independent of the initial
velocity.
4.1. Collinear case

For the special case of a collinear collision BO D 0; D 0 D 2 , the normal and tangential
dynamics decouple, and (8) becomes exact as a result. This represents the case in which the centre of
mass of each body can be joined by the line normal to the tangent plane passing through the contact
point P , for example, in the case of two spheres colliding. The eigenvalues of (15) are then given by
q
1
1 D c2 DO C c2 2 DO 2 C 4k2 DO (18)
2
q
1
2 D c2 DO c2 2 DO 2 C 4k2 DO (19)
2
q
1
3 D c1 AO C c1 2 AO2 C 4k1 AO (20)
2
q
1
4 D c1 AO c1 2 AO2 C 4k1 AO : (21)
2
Using (18), (20), (8) and (6) gives

0 1
B C
tcri t D min @ q ; q A: (22)
c2 DO C c2 2 DO 2 C 4k2 DO c1 AO C c1 2 AO2 C 4k1 AO
Note that, given the generality of this formulation, it is possible to have contact with initial relative
tangential velocity as well as relative normal velocity.
4.2. Zero damping

Consider the special case of zero damping (c1 D c2 D 0). The eigenvalues of (15) are now given by
r !
1 2
1 D p k2 DO C k1 AO C k2 DO k1 AO C 4k1 k2 BO 2 (23)
2
r !
1 2
2 D p k2 DO C k1 AO C k2 DO k1 AO C 4k1 k2 BO 2 (24)
2
r !
1 2
3 D p k2 DO C k1 AO k2 DO k1 AO C 4k1 k2 BO 2 (25)
2
r !
1 2
4 D p k2 DO C k1 AO O O O 2
k2 D k1 A C 4k1 k2 B : (26)
2
Using (23), (8) and (6) gives

tcri t D min p ; p : (27)
2 1 2 3
4.3. Two-disc collision

We will consider an example system of a two-disc rigid collision to allow for a direct comparison
with the current state of the art. For the case of a collision between discs of radii r and r 0 and masses
m and m0 , we have that BO D 0; D 0 D 2 . Considering the specific case analysed in [15] of an
undamped pair of discs (c1 D c2 D 0), (22) reduces to
0 1 0 1
B C B C
tcri t D min @ q ; q A D min @ q ; q A: (28)
mCm 0 mCm 0
2 k2 DO 2 k1 AO 2 k2 mm0 2 k1 3 mm0
For a two-disc collision of discs of equal mass m and equal spring constants k, (28) further reduces to
r r
m m
tcri t D p Ñ 0:64 : (29)
2 6 k k
q
m
This value is slightly more conservative than the estimate predicted in [15]. For the case of
k
q
close-packed configurations, as also considered in [15] and [16], the 0:64 m k
estimate reported here
is less conservative. This is due to the fact that a close-packed configuration will tend to increase
the natural frequency of the system. The estimate predicted by (28) shows that when the spring
constants are equal, the tangential frequency dictates the system stability. In practice, however, the
normal spring stiffness is typically three times greater than the tangential spring stiffness.
4.4. Full system

The eigenvalues of (15) can be calculated analytically; however, the expression for each eigenvalue
contains a large number of terms, making it non-intuitive to give here. We instead consider a numer-
ical investigation to examine the calculated stable time-step for a variety of system parameters as
shown in Section 5.
5. NUMERICAL RESULTS
Here, we will focus on the numerical analysis of how the critical time-steps calculated in Section 4
are affected by the system parameters. The purpose of this is threefold. Firstly, we want to show
how system damping affects the critical time-step. Secondly, we want to investigate the effect the
eccentricity of collision has on the critical time-step, and thirdly, we will exploit the fact that our
analysis allows for time-step calculations for any planar rigid body to investigate the effect of shape
on the critical time-step.
5.1. Undamped
In this section, we will consider the critical time-steps for a number of undamped cases. The purpose
of this is to enable a comparison with the damped case in order to correctly interpret the effect of
damping on the system. To determine the effect of eccentricity on the critical time-step, we will use
the eccentricity parameters e and e 0 defined in Section 3. For this purpose, we define a value for the
moment of inertia of each body and vary the eccentricity as shown in Figure 3. This acts as a rota-
tion about the contact point P of each body away from the collinear configuration. In Figure 3(a),
we show different contours of the spring constant k1 for zero damping in the system. For each
value of k1 , the critical time-step decreases with increasing eccentricity, showing that the system
Figure 3. The critical time-step t crit as a function of eccentricity with (a) varying contours of the spring
constant in the tangential direction k1 and (b) varying contours of the spring constant in the normal direction
k2 . In both figures, we use the initial conditions m1 D m2 D 0:003 kg, L D 0:006 m, I D Ip D 5:4
108 kg m2 and c1 D c2 D 0 kg s1 . In (a), k2 D 5 106 kg s1 . In (b), k1 D 1 107 kg s1 .
m
Figure 4. The critical time-step t crit as a function of mass ratio m 0 for a rigid disc–rigid disc collision
(a) varying contours of the total mass in the system and (b) varying contours of the spring constant in
the normal direction k2 for a collinear collision. In (a), we use the initial conditions D 1500 kg m3 ,
k1 D 4 107 kg s1 and k2 D 5 107 kg s1 . In (b), we use the initial conditions D 1500 kg m3 ,
m C m0 D 0:0005 kg and k1 D 5 106 kg s1 .
becomes less stable the further we move away from a collinear configuration. In Figure 3(b), we
show different contours of the spring constant k2 for zero damping in the system. As for Figure 3(a),
we see the same effect of decreasing stability with increasing eccentricity. It is interesting to note
6 1
how the critical time-step changes with increasing
6 1
p k1 . When we increase k1 from 1 10 kg s to
2 10 kg s in Figure 3(a), we do not obtain a 2 decrease in tcri t as expected. However, p when
we increase k1 from k1 D 2 106 kg s1 to k1 D 4 106 kg s1 , we observe the 2 decrease in
tcri t . The reason for this ‘nonlinearity’ is due to the fact that the contour k1 D 1 106 kg s1
corresponds to the eigenvalue 1 , but the contours k1 D 2 106 kg s1 , k1 D 4 106 kg s1
correspond to the eigenvalue 3 given by (23) and (25), respectively.
In Figure 5, we investigate the effect of particle size on the critical time-step. We consider a rigid
disc–rigid disc collision for systems with constant mass and investigate how the stability is affected
m 0
by decreasing the ratio m 0 , where m 6 m . In Figure 4(a), we consider different contours of the total
system mass, and in Figure 4(b), we consider varying contours of the spring constant k2 . In both
m
cases, it is clear that as the ratio m 0 7! 0, tcri t 7! 0. In Figure 4(a), we see that as the system mass
Figure 5. Schematic illustration of the rigid-body orientations used for the analysis shown in Figure 6 and
10. In (a), a rigid rectangular plate–rigid rectangular plate collision is shown, in which an isolated point
contact as opposed to a surface contact is assumed. In (b), a rigid rectangular plate–rigid disc collision is
considered, and in (c), a rigid disc–rigid disc collision is shown.
Figure 6. The variation of critical time-step t crit with spring constant ratio k 1
k2
for a rectangular plate–
rectangular plate collision in (a), and a rectangular plate–rectangular plate collision (e D e 0 D 4 104 m),
a rectangular plate–disc collision (e D e 0 D 0 m) and a disc–disc collision in (b). The rectangular plate
has dimensions a b, with a D 0:001 m and b D 0:00075 m with mPlate D mDisc D 0:003 kg, rDisc D
5 104 m. In both (a) and (b), we use the initial conditions k2 D 8 107 kg s1 , c1 D c2 D 0 kg s1 .
is increased, the stability also increases, consistent with previous findings. Similarly, in Figure 4(b),
we see that as k2 is increased, the stability decreases, as expected.
To examine the effect that particle shape and inertia has on the critical time-step, we consider two
types of laminar, namely, a rigid rectangular plate and a rigid disc (see Figure 5). In Figure 6(a),
we consider a rigid rectangular plate–rigid rectangular plate collision as a function of tangential
to normal spring constant ratio kk12 for varying contours of eccentricity of collision. We see that
the collinear case is the most stable, as expected, and that the stability decreases as the ratio kk12
increases. In Figure 6(b), we consider three cases, namely, rigid rectangular plate–rigid rectangular
plate collision, rigid rectangular plate–rigid disc collision and a rigid disc–rigid disc collision. We
see that the non-collinear rigid rectangular plate–rigid rectangular plate collision is the least stable
and that the collinear rigid rectangular plate–rigid disc collision is the most stable. This is perhaps a
surprising result in that it is a more stable configuration than the disc–disc collision. Similar to 6(a),
we see that the stability decreases as the ratio kk12 increases. In both cases, we see that for certain
conditions, tcri t remains constant with respect to kk12 . This corresponds to the transition between
the eigenvalues 1 and 3 given by (23) and (25). A similar transition was observed in the different
contours in Figure 3(a).
Figure 7. The critical time-step t crit as a function of eccentricity with (a) varying contours of the spring
constant in the tangential direction k1 and (b) varying contours of the spring constant in the normal direction
k2 . For both figures, we use the initial conditions m1 D m2 D 0:003 kg, L D 0:006 m, I D Ip D
5:4 108 kg m2 . For (a), k2 D 5 106 kg s1 , c1 D c2 D 0 kg s1 . For (b), k1 D 1 107 kg s1 ,
c1 D c2 D 250 kg s1 .
Figure 8. The critical time-step t crit as a function of damping in the normal direction with (a) varying
contours of the damping in the tangential direction for a non-collinear collision .e D e 0 D 0:0005 m/ and (b)
varying contours of the spring constant in the normal direction k2 for a collinear collision. For both figures,
we use the initial conditions m1 D m2 D 0:003 kg, I D Ip D 3:79 108 kg m2 , k1 D 4 107 kg s1 .
In (a), k2 D 5 107 kg s1 . In (b), c1 D 0 kg s1 .
5.2. Damped
In Figure 7, we present the same analysis as in Figure 3 but with system damping included. For each
value, the critical time-step decreases with increasing eccentricity and increasing damping, showing
that the system becomes less stable with increased damping and eccentricity. In Figure 7(b), we see
the same trend. In both cases, we see that the critical time-step is almost four times smaller than in
the undamped case, showing the significant effect system damping has on the stability.
In Figure 8, we investigate the effect that damping in the normal direction has on the critical
time-step. Figure 8(a) shows different contours for the tangential damping constant, and Figure 8(b)
shows different contours for the normal spring constant. The critical time-step decreases with
increased damping and increased normal force. In Figure 8(b), the critical time-step is constant until
c2 Ñ 390 kg s1 for the contour k2 D 1 106 kg s1 . This corresponds to a transition between the
eigenvalues in Section 4.4. Once a large enough damping force has been reached, one of the sys-
tem eigenvalues dominates and dictates the system stability. The result that the system becomes less
stable with increased damping is consistent with the findings in [23].
To examine the effect that particle shape and inertia have on the critical time-step for the damped
case, we consider two different numerical investigations as shown in Figures 9 and 10. In Figure 9,
m
Figure 9. The critical time-step t crit as a function of mass ratio m 0 for a rigid disc–rigid disc collision
(a) varying contours of the total mass in the system (b) varying contours of the spring constant in the normal
direction k2 for a collinear collision. In (a), we use the initial conditions D 1500 kg m3 , k1 D 4
107 kg s1 and k2 D 5 107 kg s1 . In (b), we use the initial conditions D 1500 kg m3 , k1 D 5
106 kg s1 and c1 D c2 D 50 kg s1 .
Figure 10. The critical time-step t crit as a function of (a) damping constant ratio cc12 and (b) forcing
constant ratio k1
k2
for a rectangular plate–rectangular plate collision (e D 0:0004 m), rectangular plate–disc
collision and a disc–disc collision. The rectangular plate has dimensions a b, with a D 0:001 m and
b D 0:00075 m with mPlate D mDisc D 0:003 kg, rDisc D 5 104 m. In (a), we use the initial conditions
k1 D k2 D 5 107 kg s1 and c2 D 500 kg s1 . In (b), we use the initial conditions k2 D 8 107 kg s1
and c1 D c2 D 250 kg s1 .
we consider a similar set-up to that considered in Figure 5 in Section 5.1. Figure 9(a) considers the
critical time-step as function of the mass ratio for varying contours of system damping. Figure 9(b)
considers the critical time-step as function of the mass ratio for varying contours of k2 for constant
non-zero damping. In both cases, the same trends as in Figure 5 are observed with the heightened
effect because of the presence of damping. This shows that when there is a collision between two
discs with one significantly larger than the other, then it is necessary to adopt a very small time-step
to ensure numerical stability, particularly when damping is included.
Figure 10 considers the effect that the shape of the rigid body has on the critical time-step for
the damped case. In this analysis, three different cases are considered, namely, a rigid rectangular
plate–rigid rectangular plate collision, rigid rectangular plate–rigid disc collision and a rigid disc–
rigid disc collision. In Figure 10(a), the critical time-step is calculated as a function of the damping
constant ratio cc12 while in Figure 10(b), the critical time-step is calculated as a function of the forcing
constant ratio kk12 . The critical time-step decreases with increasing ratios of damping or forcing. Both
Figure 10(a) and (b) shows the change in the critical time-step with increasing shape complexity,
with the rigid rectangular plate–rigid disc collision being the most stable, as in the undamped case in
Figure 6(b). In Figure 10(a), we see that for certain conditions tcri t remain constant with respect
to cc12 . This corresponds to the transition between the eigenvalues 1 and 3 given by (23) and
(25), similar to the trend observed in Figure 6(a) and (b). However, we do not see this behaviour
in Figure 10(b). The presence of damping, for this configuration, has the effect that 3 is always
greater than 1 in this region of parameter space.
6. DISCUSSION AND CONCLUSIONS
In the framework of using system eigenvalues to determine numerical simulation stability, an exten-
sion to current methodologies has been proposed for taking into account system damping together
with complexity in particle geometry. This method is based on a general framework, which anal-
yses the contact phase in terms of contact point equations of motion. This general framework is
compatible with any planar rigid body subject to a linear contact law.
A comparison between the classical critical time-step calculation, based on a single degree of
freedom, undamped system, and the theoretical expression for the compression time of the con-
tact phase is made and suggested as an alternative way for calculating the critical time-step. The
framework presented in this work however can be used with the current methodology at the cost of
possible inaccuracy because of the damping assumptions.
For all q
of the cases considered here, the critical time-step calculations were found to be propor-
tional to m k
, consistent with current methodologies. In each instance, this factor was scaled by
some combination of the symmetric matrix terms A, O BO and D. O These terms represent the mass dis-
q
tribution of the system. It is well known that in practice, it is necessary to multiply mk
by a ‘factor
of safety’ in order to achieve a stable numerical simulation. This work has given a physical mean-
ing to this safety factor in terms of physical parameters of the system, albeit just for the general,
linear planar cases presented here. We have also shown that, for the linear cases analysed here, the
stability of the system is independent of the initial conditions.
We performed a numerical investigation of how the critical time-step varies as a function of the
system damping together with the eccentricity of collision. This analysis is completely novel and
provides useful insights for engineers and other researchers who want to perform stable numerical
simulations for non-symmetric systems with damping. The analytical results presented here have
also confirmed previous observations in the literature showing how simulations become less stable
with increased system damping. We have also shown that there are regions in parameter space in
which the critical time-step remains constant with respect to increasing spring constant and damping
constant. This is a novel result and is something which is of great use to researchers considering
DEM simulations.
This work, in its current form, is only applicable for two contacting particles. The framework
presented here, however, is general and can be extended to incorporate a close packing of DEM
particles: a configuration which is often encountered in many DEM applications. As mentioned
previously, the effect of increasing the number of particles will be to increase the natural frequency
of the system and thus reduce the critical time-step. We expect that there may also be some coupling
between the eccentricity of each body, which may yield some interesting results.
This work is the first step in an attempt to calculate critical time-steps for a number of general
configurations. We anticipate that when this analysis is extended to nonlinear interaction laws, the
factor of safety will be a function of the initial conditions, particularly the initial contact relative
velocity.
7. ACKNOWLEDGEMENTS
This work was supported by the EPSRC grant EP/N004477/1. We also thank Mark Cook and
Andreas Piskopakis from DEM Solutions for their useful discussions.
REFERENCES
1. Stronge WJ. Impact Mechanics. Cambridge University Press: Cambridge, 2000.

2. Stronge WJ. Rigid body collisions with friction. Mathematical and Physical Sciences 1990; 431(1881):169–181.
3. Brogliato B. Nonsmooth Mechanics. Springer-Verlag: London, 2007.
4. Burns SJ, Piiroinen PT. The complexity of a basic impact mapping for rigid bodies with impact and friction. Journal
of Regular and Chaotic Dynamics 2014; 19(1):20–36.
5. Brach RM. Mechanical Impact Dynamics, Rigid Body Collisions. John Wiley & Sons: New York, 2007.
6. Burns SJ, Piiroinen PT. Simulation and long-term behaviour of unconstrained planar rigid bodies with impact and
friction. International Journal of Nonlinear Mechanics 2015; 77:312–324.
7. Osorio G, Di Bernado M, Santini S. Chattering and complex behaviour of a cam-follower system. Proceedings of
ENOC 2005 European Conference on Nonlinear Dynamics, August 2005.
8. Cundall PA, Strack ODL. A discrete numerical model for granular assemblies. Géotechnique 1979; 29(1):47–65.
9. Huang X, O’Sullivan C, Hanley KJ, Kwok CY. Discrete-element method analysis of the state parameter. Géotech-
nique 2014; 64(12):954–965.
10. Rougier E, Munjiza A, John NWM. Numerical comparison of some explicit time integration schemes used in
DEM, FEM/DEM and molecular dynamics. International Journal for Numerical Methods in Engineering 2004; 61:
856–879.
11. Bathe KJ, Wilson EJ. Numerical Methods in Finite Element Analysis. Prentice Hall, INC: New Jersey, 1976.
12. Belytschko T, Liu WK, Moran B, Elkhodary KI. Nonlinear Finite Elements for Continua and Structures (2nd ed.)
Wiley: New York, 2013.
13. O’Sullivan C. Particulate Discrete Element Modelling: A Geomechanics Perspective (1st ed.) Taylor & Francis:
UK, 2011.
14. Belytschko T. An overview of semidiscretization and time integration procedures. In Computational Methods for
Transient Analysis, Belytschko T, Hughes TJR (eds.), Computational Methods in Mechanics Series. North Holland:
New York, 1983.
15. O’Sullivan C. Selecting a suitable timestep for discrete element simulations that use the central difference algorithm
time integration scheme. Engineering Computations 2004; 21(2/3/4):278–303.
16. Tavarez FA, Plesha ME. Discrete element method for modelling solid and particulate materials. International Journal
for Numerical Methods in Engineering 2007; 70(4):379–404.
17. Li Y, Xu Y, Thornton C. A comparison of discrete element simulations and experiments for ‘sandpiles’ composed
of spherical particles. Powder Technology 2005; 160(3):219–228.
18. Lu G, Third JR, Müller CR. Discrete element models for non-spherical particle systems: from theoretical develop-
ments to applications. Chemical Engineering Science 2015; 127:425–465.
19. Houlsby GT. Potential particles: a method for modeling non-circular particles in DEM. Computers and Geotechnics
2009; 36(6):953–959.
20. Favier JF, Abbaspour-Fard MH, Kremmer M, Raji A. Shape representation of axi-symmetrical, non-spherical
particles in discrete element simulation using multi-element model particles. Engineering Computations 1999; 16(4):
467–480.
21. Burns SJ, Hanley KJ. Improving estimates of critical time-steps for discrete element simulations. Proceedings of the
7th International Conference on Discrete Element Methods: Springer, Singapore, 2016.
22. Nordmark A, Dankowicz H, Champneys A. Discontinuity-induced bifurcation in systems with impacts and friction:
discontinuities in the impact law. International Journal of Non-Linear Mechanics 2009; 44:1011–1023.
23. Fraige FY, Langston PA. Integration schemes and damping algorithms in distinct element models. Advanced Powder
Technology 2004; 15(2):227–245.

Topic7 ContactStressesAndDeformations

Uploaded by

Copyright:

Available Formats

Topic7 ContactStressesAndDeformations

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Topic7 ContactStressesAndDeformations

Uploaded by

Copyright:

Available Formats

INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING

Int. J. Numer. Meth. Engng 2017; 110:186–200

Establishing stable time-steps for DEM simulations of

Shane J. Burns*,† and Kevin J. Hanley

Received 2 June 2016; Revised 17 August 2016; Accepted 23 August 2016

KEY WORDS: DEM; time-step; stability

© 2016 The Authors. International Journal for Numerical

It is necessary to choose a numerical integration procedure to solve the equations of motion

2. STABILITY AS A FUNCTION OF MASS AND STIFFNESS

2.1. Stability using Rayleigh’s theorem

where max is the maximum eigenvalue of

!max 2 D max : (6)

2.2. Natural frequency of the system

q D .q1 ; q2 ; /T ; q 0 D .q10 ; q20 ;  0 /T ;

e D L cos./ and e 0 D L0 cos. 0 /:

 D .T ; N / and  0 D .T ; N / :

1 L2 sin2 ./ L2 sin./ cos./ 1 L2 cos2 ./

1 L0 2 sin2 . 0 / L0 2 sin. 0 / cos. 0 / 1 L0 2 cos2 . 0 /

3.1. Linear contact interaction

1 L2 cos2 ./ 1 L0 2 cos2 . 0 /

where the system matrix

subject to the initial conditions

where 1 , 2 , 3 , 4 , and V1 , V2 , V3 , V4 are the eigenvalues and corresponding eigenvectors of

4.1. Collinear case

Using (18), (20), (8) and (6) gives

4.2. Zero damping

4.3. Two-disc collision

4.4. Full system

6. DISCUSSION AND CONCLUSIONS

1. Stronge WJ. Impact Mechanics. Cambridge University Press: Cambridge, 2000.

You might also like

q D .q1 ; q2 ; /T ; q 0 D .q10 ; q20 ; 0 /T ;

e D L cos./ and e 0 D L0 cos. 0 /:

D .T ; N / and 0 D .T ; N / :

1 L2 sin2 ./ L2 sin./ cos./ 1 L2 cos2 ./

1 L0 2 sin2 . 0 / L0 2 sin. 0 / cos. 0 / 1 L0 2 cos2 . 0 /

1 L2 cos2 ./ 1 L0 2 cos2 . 0 /