ON THE VALIDATION OF DEM AND FEM/DEM MODELS IN 2D AND 3D
Jiansheng Xiang, Imperial College, j.xiang@imperial.ac.uk
Antonio Munjiza, Queen Mary London University, a.munjiza@qmul.ac.uk
John-Paul Latham, Imperial College, j.p.latham@imperial.ac.uk
Romain Guises, Imperial College, romain.guises@imperial.ac.uk
ABSTRACT
As particulate systems evolve, sliding, rolling and collision contacts all produce forces that DEM
methods aim to predict. Verification of friction rarely takes high priority in validation studies even
though friction plays a very important role in applications and in mathematical models for numerical
simulation. We use analytical solutions for ‘block sliding’ to verify our tangential contact force
implementation of 2D FEM/DEM. Inspired by the kinetic art work “Liquid Reflections” by Liliane
Lijn, which consists of free balls responding within a rotating shallow dish, we use DEM to simulate
rolling, sliding and state-of-rest of spherical particles relative to horizontal and inclined, concave and
flat spinning platforms. Various material properties, initial and boundary conditions are set which
produce different trajectory regimes. Simulation output is found to be in excellent agreement when
compared with experimental results and analytical solutions. The more widespread use of
analytically solvable benchmark tests for DEM and FEM/DEM codes is recommended.
1. INTRODUCTION
Validation of our DEM community’s numerical simulation tools has taken many forms. Reasonably
accurate mathematical representations of the contact physics between particles and container walls
including friction effects have been implemented with varying sophistication in DEM codes. Li et al
(2005) produced a classic validation study. They performed experiments to complete the
determination of all material properties including friction coefficients applicable to glass and steel
sphere simulations so that ‘sandpile’ experiments with spherical particles could be numerically
simulated with DEM. In their DEM implementation, energy dissipation occurs through mechanisms
governed by accurately measurable material parameters, rather than invoking numerical damping
through un-measurable coefficients. Their results for angle of repose, when compared with the
equivalent experimentally produced angles are therefore all the more convincing. Theirs is one of
many cases of an emergent property of the granular system being used to compare and validate the
simulator’s accuracy. Other pseudo-static emergent properties include porosity, coordination number
statistics, fabric tensor statistics, force chain characteristics. The dynamics of granular systems in
flow, flow rates through apertures, velocity field statistics, and the special emergent properties
associated with tumbling mills, including power draw, all have been subject to experimental versus
DEM comparisons. (e.g. Cleary 2002). Unlike the work of Li et al (2005), Asmar et al (2002)
produced a numerical validation study. In their study, they set up eight ‘mathematical tests’ based on
artificial situations. Although they didn’t compare their results with that of experiments, their study
shows some simple simulations can help to verify code and reveal bugs.
Using FEM/DEM, Latham and Munjiza (2004) presented a comparison of cube-packing
experiments and their equivalent numerical simulation. They found sensitivity to initial conditions to
be even more prevalent for these non-spherical bodies, but emergent bulk behaviour in close
agreement.
One underused tool of validation for accuracy of contact physics might be to subject the motion
1
of simple systems involving few particles to deeper scrutiny. This is usually considered impractical
because there is so much sensitivity to initial conditions (chaos) in inter-particle mechanics,
especially for non-spherical shapes and increased particle numbers. The advantage of a rotating
container simulation is that there is a constant supply of energy to the system, the particle(s) can be
spatially constrained and there may be some interesting boundary conditions for which analytical
solutions exist or experiments can be devised with a reasonably high degree of parameter control.
Friction plays a very important role in most technical applications and from a numerical view
point, friction is an important part of the mathematical model. Verification of friction effects during
simulation should take high priority in the work of validation. This paper begins by addressing
sliding friction in FEM/DEM and rolling friction in DEM.
2. COMPUTER IMPLEMENTATION
2.1 FEM/DEM with sliding friction
In this paper, there are two groups of validation tests carried out based on two programs. One is
named Y code developed by Munjiza (2004) using a coupled discrete-finite Element Method
(FEM/DEM). Details of this method are presented in the references (Komodromos and Williams
2004, Munjiza & Owen 1995) and are further explained with reference to topics in computational
mechanics of discontinue in a Munjiza’s book (Munjiza 2004). Another one is called 3DSY which is
described in section 2.2.
Here, we develop further the FEM/DEM method by taking account of the sliding friction force.
The well-known classic Coulomb type friction is implemented and described as follows,
Ft = − Ft ,n − ηvt (1)
If Ft is bigger than the friction force obeying the Coulomb-type friction law, the particles slide
over each other and the tangential force is calculated by:
Ft = − µ Fn (2)
where is the coefficient of viscous dissipation, Fn is normal elastic force, and Ft,n is tangential
η
elastic force.
2.2 DEM with rolling friction
In the paper of Asmar et al (2002), eight numerical tests were set up to validate their code, named as
DMX. These tests constitute a standard benchmark. However, the effect of rolling friction was
neglected and not implemented in their model. Zhou et al (1999) investigated the mathematical
importance of rolling friction on the formation of a heap of spheres. They demonstrated that the
rolling friction plays a critical role in achieving physically or numerically stable results, and the
angle of repose increases with rolling friction coefficient and decreases with particle size. However,
it is perhaps surprised that to date validation of rolling friction rarely takes high priority in validation
studies.
A DEM code called 3DSY is developed by Xiang using a linear-spring and dashpot model. The
details of this model are presented in cited references (Tsuji et al 1993, Xiang 2004, Mustoe &
Miyata 2001). Beer and Johnson’s model (1976) which is described below is used to calculate the
rolling friction.
M = − µ r Fcnωˆ (3)
2
where M is the torque due to rolling friction, µ r is coefficient of rolling friction, Fn is the normal
contact force, and ω̂ is the direction of angular velocity.
3. VERIFICATION FOR FEM/DEM
The case of a rectangle which is paced on a horizontal plane with initial horizontal velocities is
carried out as benchmark test for 2D FEM/DEM with sliding friction. The stop distance L can be
theoretically derived and is given by the following equation,
2
vi (4)
L=
2µg
where is the coefficient of sliding friction, g is the gravity acceleration, and vi is the initial
horizontal velocity.
4
stop distance (theoretical analysis)
3.5 stop distance (numerical results with time step 1.0e-7s)
stop distance (numerical results with time step 1.0e-8s)
3
2
vi
L=
Stop distance L (m)
2.5
2 µg
2
1.5
0.5
0
0 1 2 3 4 5 6
Initial velocity v (m/s)
Figure 1 Stop distance plotted against initial velocity for a rectangle with frictional coefficient, µ=0.5
The FEM/DEM simulation of sliding rectangle is performed for a square rectangle with length,
l=0.05 m, density ρ=2650 kg/m3, frictional coefficient µ=0.5, and Young modulus E=1.0×109 Pa.
The rectangle is given with initial velocity vi then slows down due to the frictional effect with the
base. Finally it stops at a distance L. Two sets of simulation tests were repeated with different initial
velocity with two different time steps, ∆t=1.0×10-7 s and ∆t=1.0×10-8 s. Results are compared with
the analytical solution (see Fig 1). It shows that with a small time step, the numerical results are in
excellent agreement with the theoretical values. It is worth noting that with the larger of the two time
steps, the errors become significant. However, using the larger time step, the calculation of
FEM/DEM with zero friction is fairy stable. This suggests, in order to reduce the numerical error for
calculation of tangential forces, the smaller time step is required.
4. VALIDATION for DEM
Ehrlich and Tuszynski (1994) carried out several experimental tests to validate their theoretical
model. In their study they analyzed the motion of a ball on both a flat and a parabolically curved
rotating turntable. Compared with DEM, their model is quite simple. In their model it is assumed
that sliding friction force and rolling friction torque are constants. In our DEM simulation, the
contact forces are calculated using a spring, dashpot and slider (Xiang 2004).
3
Soodak and Tiersten (1996) further explained E&T’s experimental results and provided the
understanding of the effects of rolling friction and other ‘perturbations’, e.g. a tilt of the turntable. In
this paper we choose E&T’s experimental results of a flat turntable to validate our DEM code.
The works of art “Liquid Reflections” by Liliane Lijn, were made in 1966-68. They incorporate
a rotating circular platform and large acrylic spheres combined with spot lighting to produce an
intriguing example of 1960’s kinetic art that can be seen in many national modern art collections
(see Fig. 2 and the video “force fields” at http://www.lilianelijn.com/vid02.html).
Figure 2 “Liquid Reflections” by Liliane Lijn. In the words of the artist: “The movement of the balls on the surface
of the disc is governed by the laws of momentum, as well as centrifugal force and the pull of gravity induced by the
concavity of the disc. The balls also act as moving magnifying lenses, bringing to life now one area of the disc, now
another, with a strange lunar landscape of reflections and shadows”.
Drawing inspiration from Lijn’s experiment setup, we did several tests and found that if a ball
was released from rest, the ball tended to stay at a certain height and finally ran on a closed circular
orbit. With a certain angular velocity of the turntable, at a certain height, the contact force between
the ball and turntable can balance centrifugal force and gravity to maintain circular motion. The
stable height is derived below,
The critical speed for a banked curve with zero friction is v 2 = rg tan θ .
As for spherically curved turntable, v = ωr , r = R 2 − (R 2 − h 2 ) , and tan θ =
(
R2 − R2 − h2 )
, then
R−h
we can obtain
g
h = R − (5)
ω 2
According to equation (5), height increases with the increasing of angular velocity. When the
g
angular velocity is unlimited, the height equals R. If ω ≤ , h ≤ 0 . This is not physical correct, so
R
we make a criteria for the calculation and modify equation (5),
4
g g
R − , if ω ≥
h = ω 2 R
(6)
g
0, if ω <
R
In our DEM simulation, we examine spherical particle motions in flat and spherically curved
rotating platform geometries, beginning with behaviour for which analytical solutions can be derived
for the motion, progressing to more complex regimes.
In comparing with E&T’s experimental results and Lijn’s visual results, we set up the same set
of conditions. The physical properties of ball and geometry of the turntable are shown in Table 1.
Table 1. Turntable and ball properties used in the DEM models
Turntable
shape flat spherically curved
radius R (m) 0.2 20.0
Ball
shape spherical spherical
density ρ (kg m-3) 7800 7500
radius rs (mm) 6.35 3.00
frictional coefficient µ 0.15 0.5
coefficient of restitution e 0.8 0.8
coefficient of rolling friction µr 1.5×10-6 5.0×10-5
5
spring stiffness k 1.0×10 1.0×103
4.1 Flat turntable
Test 1
Several studies (Weltner 1979, Ehrlich and Tuszynski 1994) show that if the ball is given an initial
push, the ball runs on a stable circular orbit. The angular velocity is 2/7 of the turntable’s angular
velocity. This phenomenon is not related to initial position and radius of ball. In our simulation test,
before the ball is released on the rotating turntable it is given the initial spin velocity such that there
is no slipping between ball and rotating turntable. If the rolling friction is turned off, we get a stable
circular orbit (see Fig. 3). The angular velocity is 2.855 rad/s, 2/7 of turntable’s angular velocity.
Test 2
Figure 4 shows the case where the turntable is slightly tilted (0.18º) towards north with non-zero
rolling friction, the orbit radius grows in size, and its centre drift at fairly constant speed in the
direction of 225.00º which is the direction of kˆ × gˆ . (note: we define 0.00º represents east →, 90.00º
represents north ↑, 180.00º represents west ←, and 270.00º represents south ↓, see Fig. 4) k̂ is the
unit vector perpendicular to the surface of the turntable, directed from the surface to the ball centre.
ĝ is the unit vector in the direct of inclination. The numerical results agree with the experimental
results (Ehrlich and Tuszynski, 1994) and theoretical explanations (Soodak and Tiersten, 1996).
Test 3
When the turntable is given a slightly greater tilt degree, its orbit centre drifts at fairly constant
speed in the direction of 180.00º which is the direction of kˆ × gˆ . The simulation reproduces the
5
trajectory in Fig. 5 in E&T’s paper (Ehrlich and Tuszynski, 1994). When the rolling friction is
increased to 2.5×10-5 (see Fig. 6), the orbit centre drifts to the centre of turntable and orbit radius
grows in size. The output simulation is in good agreement with the numerical results in the paper
(Soodak and Tiersten, 1996).
Test 4
Fig. 7 shows a special trajectory when the orbit has a nearly fixed centre which grows slowly in size.
The drift is very small but in the direction of kˆ × gˆ . The simulation successfully reproduces the
trajectory of experimental results (Ehrlich and Tuszynski, 1994). Fig. 8 shows the trajectory when
the turntable tilt is changed from 0.17º to 0.00º. The numerical results in Fig 7 of E&T’s paper,
compare very well with our result in Fig. 8.
Test 5
Fig. 9 shows an interesting trajectory which is called, ‘Approach-to-the-centre trajectory’ (Ehrlich
and Tuszynski, 1994). The translational velocity of the ball was set to 0.00 m/s and the ball was
released under no slipping condition. Numerical results fit experimental results (Ehrlich and
Tuszynski, 1994) which shows that the ball approaches a certain point along a curved path (note that
this certain point is not exactly the centre of the turntable). The ball firstly reaches centreline of the
turntable, x axes in a certain angle. Then it approaches to the point very slowly after the inflection
point. After it reaches the point, the orbit radius grows in size and its centre drifts in the direction of
kˆ × gˆ .
N N
0.40 0.40
z (m)
z (m)
0.35 0.35
0.30 0.30
0.25 0.25
E W E
W 0.20 0.20
x (m)
x (m)
0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40
0.15 0.15
0.10 0.10
0.05 0.05
0.00 0.00
S S
Figure 3 A numerically computed trajectory of Figure 4 A numerically computed trajectory of
moving ball on flat turntable with tilt angle 0.00º, no moving ball on flat turntable with tilt angle 0.17º in
slipping initial condition, angular velocity of the direction of 135º, no slipping initial condition,
turntable 10.00 rad/s, initial velocity 0.0 m/s, initial angular velocity of turntable 10.00 rad/s, initial
radius 0.04 m, initial angle 180.00º, rolling friction velocity 0.107 m/s in 270.00º, initial radius 0.04 m,
0.00. initial angle 180.00º, and rolling friction 1.5×10 -6.
6
0.40
N 0.40 N
z (m)
z (m)
0.35
0.35
0.30
0.30
0.25
0.25
0.20 E W0.00 0.20 Ex (m)
W 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40
x (m)
0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40
0.15 0.15
0.10 0.10
0.05 0.05
0.00 0.00
S S
Figure 5 A numerically computed trajectory of Figure 6 A numerically computed trajectory of
moving ball on flat turntable with tilt angle 0.52º in moving ball on flat turntable with tilt angle 0.52º in
the direction of 90.00º, no slipping initial condition, the direction of 90.00º, no slipping initial condition,
angular velocity of turntable 17.60 rad/s, initial angular velocity of turntable 17.60 rad/s, initial
velocity 0.008 m/s in 90.00º, initial radius 0.01 m, and velocity 0.008 m/s in 90deg, initial radius 0.01 m, and
initial angle 0.00º, and rolling friction 1.5×10-6. initial angle 0.00º, and rolling friction 2.5×10-5.
N 0.40
N
0.40
z (m)
z (m)
0.35 0.35
0.30 0.30
0.25 0.25
W E E
0.20 W
x (m)
0.20
x (m)
0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40
0.15 0.15
0.10 0.10
0.05 0.05
0.00 0.00
S S
Figure 7 A numerically computed trajectory of moving Figure 8 A numerically computed trajectory of moving
ball on flat turntable with tilt angle 0.17º in the direction ball on flat turntable with tilt angle 0.00º in the
of 90.00º, no slipping initial condition, angular velocity direction of 90.00º, no slipping initial condition,
of turntable 17.60 rad/s, initial velocity 0.009 m/s in angular velocity of turntable 17.60 rad/s, initial
270.00º, initial radius 0.13 m, and initial angle 180.00º, velocity 0.009 m/s in 270.00º, initial radius 0.13 m,
and rolling friction 1.5×10-6. and initial angle 180.00º, and rolling friction 1.5×10-6.
7
0.40
N
0.25
z (m)
z (m)
0.24
0.35
0.23
0.22
0.30
0.21
0.20 x (m)
0.15 0.16 0.17 0.18 0.19 0.20 0.21 0.22 0.23 0.24 0.25
0.25
0.19
W E 0.18
0.20 x (m)
0.17
0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40
0.16
0.15
0.15
0.10
0.05
0.00
S
Figure 9 A numerically computed trajectory of moving ball on flat turntable with tilt angle 0.10º in the direction of
180.00º, no slipping initial condition, angular velocity of turntable 17.60 rad/s, initial velocity 0.00 m/s, initial radius
0.08 m, and initial angle 45.00º, and rolling friction 1.5×10-6.
4.2 Spherically curved turntable
21.0 36
20.8
32
20.6
28
20.4
24
20.2
z (m)
z (m)
20.0 20
19.8
16
19.6
12
19.4
8
19.2
19.0 4
19.0 19.2 19.4 19.6 19.8 20.0 20.2 20.4 20.6 20.8 21.0 4 8 12 16 20 24 28 32 36
x (m) x (m)
Figure 10 A numerically computed trajectory of Figure 11 A numerically computed trajectory of
moving ball on spherically curved turntable with tilt moving ball on spherically curved turntable with tilt
angle 0.00º, angular velocity of turntable 0.5 rad/s, angle 0.00º, angular velocity of turntable 0.9 rad/s,
initial radius 1 m, and initial angle 180.00º, and rolling initial radius 1 m, and initial angle 180.00º, and rolling
friction 5.0×10-5. The ball is released from static. friction 5.0×10-5. The ball is released from static.
A ball was released from static on a rotating of spherical i.e. constant curvature turntable. When
the turntable rotates slowly, the ball spirals to the centre of the turntable. Finally, it comes to rest at
the centre (see Fig. 10). When the turntable rotates fast, ball spiral outward and finally stays
8
at a certain height and on a stable circular orbit (see Fig. 11). Several tests were carried out with a
wide range of angular velocity of the turntable. Fig. 12 shows the comparison of numerical results of
stable heights with analytical solution according to equation (6). The figure shows when the angular
velocity is below 0.7 rad/s which is the critical value calculated by equation (6), the ball tends to
stay at the centre of turntable and remains stationary. When angular velocity is above 0.7 rad/s, the
stable height increases quickly with the increase of angular velocity. When angular velocity
increases further to more than 2 rad/s, the stable height increases slowly and is close to its maximum
value of 20 m.
20.0
15.0
Stable Heights (m)
10.0
5.0
Numerical results
Analytical solution according to Equation (6)
0.0
0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0
Angular velocity (rad/s)
Figure 12 Stable height to which ball will climb against angular velocity of turntable
5. CONCLUSION
The paper illustrates the extent to which our DEM and FEM/DEM codes agree with theoretical
solutions and suggests the possible need for further analytically solved benchmark tests for DEM
and FEM/DEM codes as an extension to those proposed by Asmar (2002). The potential exists to
extend this approach to tackle code validation using physical experiments conducted with system
collisions using more than one particle. Study of orbiting system geometries may also throw light on
validation for non-spherical shapes.
REFERENCES
1. B.N. Asmar, P.A. Langston, A. J. Matchett and J. K. Walters, 2002, Validation tests on a
distinct element model of vibrating cohesive particle systems, Computers & Chemical
Engineering, 26,785-802
2. F.P. Beer, E.R. Johnson, 1976, Mechanics for Engineers – Statics and Dynamics, MacGraw-
Hill, New York.
3. P.W. Cleary, R. Morrisson and S. Morrell, 2003. Comparison of DEM and experiment for a
scale model SAG mill International Journal of Mineral Processing, 68, 129-165
4. R. Ehrlich, J. Tuszynski, 1995, Ball on a rotating turntable: Comparison of theory and
experiment, American Journal of Physics 63, 351-359.
5. P.I. Komodromos, J.R. Williams, 2004, Dynamic simulation of multiple deformable bodies
9
using combined discrete and finite element methods, Engineering Computations, v 21, n 2, 431-
48
6. J.P. Latham, A. Munjiza. 2004, The modelling of particle systems with real shapes,
Philosophical Transactions of the Royal Society London, Series A (Mathematical, Physical and
Engineering Sciences), v 362, n 1822, 1953-72
7. G.G.W. Mustoe, M. Miyata, 2001, Material Flow Analyses of Noncircular-Shaped Granular
Media Using Discrete Element Methods, Journal of Engineering Mechanics, 2001.
8. Y. Li, Y. Xu, C. Thornton, 2005, A comparison of discrete element simulations and experiments
for ‘sandpiles’ composed of spherical particles. Powder Technology, 160, 219-228
9. A. Munjiza, D.R.J. Owen; N. Bicanic, 1995, Combined finite-discrete element method in
transient dynamics of fracturing solids, Engineering Computations (Swansea, Wales), v 12, n 2,
145-174
10. A. Munjiza, 2004, The Combined Finite-Discrete Element Method, John Wiley & Sons.
11. H. Soodak and M.S. Tiersten, 1996, Perturbation analysis of rolling friction on a turntable.
American Journal of Physics 64, p. 1130-1139
12. Y. Tsuji, T. Kawaguchi, T. Tanaka, 1993, Discrete particle simulation of two-dimensional
fluidized bed, Powder Technology 77, 79–87
13. K. Weltner, 1979, Stable circular orbits of freely moving balls on rotating discs, American
Journal of Physics 47, 984-986.
14. J.S. Xiang, 2004, Investigation of particle motion in dense phase pneumatic conveying,
Glasgow Caledonian University.
15. Y.C. Zhou, B.D. Wright, R.Y. Yang, B.H. Xu, A.B. Yu, 1999, Rolling friction in the dynamic
simulation of sandpile formation, Physica A, v269, n2, 536-553
10
