Procedia CIRP 31 ( 2015 ) 369 – 374
Available online at www.sciencedirect.com
2212-8271 © 2015 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license
(http://creativecommons.org/licenses/by-nc-nd/4.0/).
Peer-review under responsibility of the International Scientific Committee of the “15th Conference on Modelling of Machining Operations
doi: 10.1016/j.procir.2015.03.021
ScienceDirect
15th CIRP Conference on Modelling of Machining Operations
Discrete Element Modelling of Drag Finishing
Eckart Uhlmanna, Alexander Eulitza*, Arne Dethlefsa
aTechnische Universität Berlin, Institute for Machine Tools and Factory Management, Berlin, Germany
* Corresponding author. Tel.: +49-30-314-24963; fax: +49-30-314-25895. E-mail address: [email protected]
Abstract
Drag finishing is a machining process that is used to improve the surface topology of workpieces. Workpieces are moved through a bulk of
differently shaped abrasives, the so called media. Material removal is caused by the relative motion between workpiece and media. The
material removal rate is mainly depending on the contact intensity between workpiece and media. Up to now there is no viable way to
determine the intensity of single contacts empirically. However, a sound understanding of single contacts with respect to impact forces and
velocities could greatly improve process comprehension and reduce trial and error process design efforts. For that reason the movement of
media and workpiece is modelled using the Discrete Element Method (DEM). In this paper a comprehensive approach is presented covering
formulation, calibration, validation and utilization of the DEM. Media is considered as an aggregation of elastic particles that are subject to
contact, damping and gravitational forces causing particle movement. Geometric boundary conditions, i.e. workpiece and drag finishing bowl,
are implemented as elastic facets. Contact forces are calculated according to a non-linear, simplified Hertz-Mindlin contact force model. Energy
is dissipated by viscous damping and friction at contacts. Necessary parameters of the model are determined experimentally. The validation of
the model’s behaviour shows good agreement with experimental data. Finally the model is used to determine local contact intensities on the
workpiece surface and between particles. By analysing simulated contact forces, the formation of dominant contact chains between particles is
observed and investigated.
© 2015 The Authors. Published by Elsevier B.V.
Peer-review under responsibility of The International Scientific Committee of the “15th Conference on Modelling of Machining
Operations”.
Keywords: Finishing; Modelling; Discrete Element Method (DEM)
1. Introduction
Drag finishing is a versatile machining process that is used
to improve the surface topology of workpieces. These are
moved in a bulk of abrasives, the so called media. Caused by
the relative motion of workpiece and media, material of the
workpiece is removed. In contrast to other mass finishing
processes, like vibratory or centrifugal disc finishing,
workpieces are fixed by a clamping device and media is not
agitated in drag finishing. Industrial drag finishing processes
are usually carried out with a fixed track setup that guides the
workpieces through the media on a set or only slightly
adaptable trajectory. A major benefit of the fixed workpiece
setup is that workpieces cannot collide. Typical applications
of drag finishing include medical implants, gear components,
turbine blades or rounding of cutting tool edges [1,2,3]. A
newly developed and patented setup is robot guided drag
finishing [4], where an industrial 6-axes robot is used to guide
workpieces through media, resulting in versatile trajectories
and precise setting of workpiece velocities [5]. By increasing
workpiece velocities higher pressure than in vibratory
finishing processes can be achieved resulting in increasing
material removal rates and quicker processing. In addition the
defined trajectories of the workpieces allow for enhanced
machining of certain features at complex workpieces.
A lot of researchers have presented models for mass
finishing processes so far. Most of them are descriptive
models for vibratory finishing that correlate process
parameters with surface quality for certain workpieces.
MATSUNAGA [6] investigated contact mechanisms and
velocity fields in vibratory finishing, but did not establish a
correlation between processing results and measured or
observed process characteristics. Additionally HASHIMOTOS
findings [7] are of high relevance to the field of mass
© 2015 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license
(http://creativecommons.org/licenses/by-nc-nd/4.0/).
Peer-review under responsibility of the International Scientifi c Committee of the “15th Conference on Modelling of Machining Operations
370 Eckart Uhlmann et al. / Procedia CIRP 31 ( 2015 ) 369 – 374
finishing. Based on empirical studies he established basic
rules for mass finishing processes that were confirmed in
large parts by a recent investigation of UHLMANN ET AL. [8].
SPELT’S group [9,10,11,12,13] established a correlation
between process forces, contact mechanisms and reachable
surface quality in vibratory finishing. An important finding of
SPELT’S group is the identification of three major contact
modes that occur when spherical media is used in vibratory
finishing: single impact, free rolling of single spheres and
rolling of spheres in packages [13].
Nomenclature
Į empirical constant
cn damping coefficient in normal direction
ct damping coefficient in tangential direction
ds submersion depth
en coefficient of restitution in normal direction
et coefficient of restitution in tangential direction
E Young’s modulus
Fn normal force
Ft tangential force
G shear modulus
hf height of filling
kn contact stiffness in normal direction
kt contact stiffness in tangential direction
ij friction coefficient of material combination i, j
m mass
nb number of balls
Rb radius of balls
Ri radius of particle i
ȡi density of material i
įn displacement in tangential direction
įt displacement in normal direction
ǻIJ time step
tp process time
Ȟi Poisson’s ratio of material i
vt contact tangential velocity
vw workpiece velocity
The significance of contact forces and mechanisms in mass
finishing processes has been emphasized by many researchers
[9,14,15,16]. It has been shown that contact forces and impact
velocities have great influence on process outcome in terms of
material removal, roughness and hardness of the workpiece.
As a consequence there have been many investigations on
how to measure forces and impact velocities on the
workpiece. Most of them make use of in situ measurement of
contact forces. Costly measurement devices with high
sampling rates have to be used to record maximum impact
forces [9]. Furthermore a major handicap of this approach is
the dependence of force magnitudes on sensor surface
compliance [9,17]. CIAMPINI ET AL. [14] proposed to extract
normal impact velocity from measured forces instead. Based
on the determined impact velocities, energy imparted to
arbitrary workpieces can be inferred, disregarding
compliances. Nevertheless forces have to be measured for all
considered machine settings and workpiece geometries.
Directly obtaining impact velocities in mass finishing
processes is also challenging, as measuring devices have to be
placed appropriately in the bulk media. HASHEMNIA ET AL.
[18] have developed a laser displacement probe that is able to
measure impact and bulk flow velocities in vibrationally
fluidized beds. Other reported approaches concentrate on
measuring workpiece or bulk flow velocity, e. g. using spatial
and time-based analyses of photographs [19]. Instead of force
or velocity based approaches CIAMPINI ET AL. [11] used an
almen strip system, common in shot peening applications, to
characterise the effect of different process parameters on the
aggressiveness of a vibratory finishing process. Having
simulated the occurring impacts with an electromagnetic
apparatus they found that normal impacts are dominating in
vibratory finishing. Because in situ force and impact velocity
measurement is a complex task susceptible to disturbances in
the measuring chain, a different approach is presented, that
allows for modelling of contact intensities for arbitrary
workpiece geometries and materials using the Discrete
Element Method (DEM).
Up to now results from two different approaches towards a
deterministic numerical process model for mass finishing
have been reported:
x Modelling of the motion of spherical media with the DEM
[20,21,22,23,24]
x Modelling of the media’s velocity fields using
Computational Fluid Dynamics (CFD) for centrifugal disk
finishing [19] and vibratory finishing of immobilized
cylinders [25]
Notably, only CARIAPA [19] and UHLMANN ET AL. [20]
look at the three-dimensional movement of abrasive media in
their models. All other approaches consider mass finishing
processes with spherical steel or glass media or are restricted
to two dimensions. CARIAPA’S approach allows for
a qualitative prediction of material removal. The work of
UHLMANN ET AL. on a comprehensive process model [20]
combines an empirical geometry based model for surface
roughness prediction [8] with contact intensity simulation
using DEM. Using this comprehensive process model a
quantitative prediction of material removal and roughness
should be possible.
2. Discrete Element Method (DEM)
In the following a model to numerically simulate contacts
between abrasive media and workpieces using DEM will be
presented. The goal is to simulate the number, type and
intensity of contacts between workpieces and abrasive media
for finite areas of any workpiece. The DEM-model is used to
simulate a small-scale drag finishing process. Findings from
the model are utilized to broaden process comprehension.
The approach presented in this paper is based on the work
of UHLMANN ET AL. [8], in which bulk motion is modelled
using the open source framework Yade [26]. Media is
considered as a bulk of discrete, elastic, spherical particles
that overlap depending on contact forces, the so called soft-
particle approach. In contrast to the hard-particle approach,
that is based on instantaneous exchange of momentum when
contact occurs, the soft-particle approach allows for multiple
371
Eckart Uhlmann et al. / Procedia CIRP 31 ( 2015 ) 369 – 374
particle contacts, which are common for dense granular
media, at the cost of smaller time steps, resulting in higher
computation times [27]. Most soft-particle models make use
of an artificial increase in time step by allowing softer
interaction. This is usually done by setting particle stiffnesses
lower than physically observable [28].
Geometric boundary conditions, i. e. bowl and workpieces,
are implemented as facets. As a first step, spherical media is
considered. In that case only point-contacts occur, which
reduces the effort for collision detection and contact
calculation considerably. In order to simulate the number,
type and intensity of contacts for finite areas of any
workpiece, a suitable model for particle-particle and particle-
facet contacts has to be chosen. This is presented in the next
section.
3. Model formulation
In the presented approach, contact forces between two
bodies are calculated according to the non-linear viscoelastic
simplified Hertz-Mindlin contact force model, see Fig. 1. It is
generally accepted that the Hertz-Mindlin contact force model
is suitable to model contact between two elastic spheres or
one sphere and a half space [29,30].
Fig. 1. Spring dash-pot slider contact model for two spheres in contact.
The model is based on Hertzian contact equations [31] in
normal direction and the Mindlin no-slip approach [32,33] in
tangential direction. The non-linear correlation between
normal force Fn and normal displacement įn is given by
2/3
nnn kF
G
. (1)
Contact stiffness kn is expressed by
eqeqnREk 3
4
, (2)
where Eeq is the equivalent Young’s modulus expressed by
2
2
2
1
2
111
1
EEEeq
QQ
, (3)
with E1 and E2 being the Young’s moduli and Ȟ1 and Ȟ2
Poisson ratios of particles in contact. Req is the equivalent
radius of the contact which is given by
21
111
RRReq
, (4)
where R1 and R2 are radii of particles in contact. According to
MINDLIN [32] and MINDLIN AND DERESIEWICZ [33] the
relation between the tangential force Ft and tangential
displacement įt is depending on the normal displacement įn.
They proposed equations that constitute this relationship for
several contact configurations. If a contact surface is allowed
to slip, equations become complicated because the history of
the contact has to be considered. To reduce computing time a
non-slipping contact surface is assumed in this paper. I. e. the
equation for tangential contact stiffness kt is given by
n
m
m
eqt
G
Rk
G
Q
2
4, (5)
where Gm is the mean shear modulus and Ȟm is the mean
Poisson’s ratio of the particles in contact. In general shear
modulus G can be expressed by
)1(2
Q
E
G. (6)
Note that tangential contact stiffness kt is a function of the
normal displacement įn. In Yade, tangential contact force Ft is
calculated incrementally by
W
W
' tttt kFF v
1, , (7)
where Ft,IJ-1 is the tangential contact force at the previous time
step, vt is the contact tangential velocity and ǻIJ is the time
step.
3.1. Damping
Energy of the particulate media in the system is dissipated
through friction and viscous damping at the contacts. Force
transmitted in tangential direction between two particles is
limited by the Mohr-Coulomb friction law. It should be noted
that in the original DEM-formulation by CUNDALL AND
STRACK [34] a global damping was introduced. I. e. damping
forces were part of the equations of motions. This does not
represent a physical phenomenon but is an artificial numerical
damping. In this work non-linear viscous damping at the
contacts (also known as local damping), represented by the
dash-pots in Fig. 1, is used. It has been shown by ZHANG AND
WHITEN [35] that non-linear contact damping is more suitable
to resemble experimental results than linear contact damping.
Damping force in normal direction Fn,d is calculated by
nnndn vcF )(
,
G
, (8)
where the damping coefficient cn is a function of normal
displacement įn. TSUJI ET AL. [36] derived a formulation for
the damping coefficient cn
372 Eckart Uhlmann et al. / Procedia CIRP 31 ( 2015 ) 369 – 374
4/1
)( nneqnn kmc
GDG
, (9)
where Į is an empirical constant related to the coefficient of
restitution in normal direction en and meq is the equivalent
mass, calculated in analogous manner as (4). ANTYPOV AND
ELLIOTT [37] mapped the Hertzian contact model onto the
linear spring dash-pot model. They provided an analytical
relationship between the empirical constant Į from (9) and the
coefficient of restitution in normal direction en, i. e.
²)²(log
)log(5
)(
S
D
n
n
ne
e
e. (10)
Conforming to the work of TSUJI ET AL. [36] the damping
coefficient in tangential direction ct is assumed to be equal to
the damping coefficient in normal direction cn.
4. Model setup
A large number of abrasive particles are used in usual drag
finishing processes, see Fig. 2a. This results in immense
calculation times for setups modelled exactly according to
these processes. For that reason a smaller bowl and rod are
considered as a first step, see Fig. 2b. This reduces the
number of particles to 1422. Geometric boundary conditions
of the modelled setup can be found in Fig. 3.
bowl
workpiece
media
v
W
60 mm
v
w
workpiece
25 mm
a) b)
Fig. 2. Drag finishing setup (a) real process, (b) DEM-model.
Calibration approaches for modelling of granular material
can roughly be divided into two groups: physical and
phenomenological approaches. In physical calibration
approaches, direct correlations between model parameters and
experimental results are used to calibrate the model. Thus
model parameters are determined on a per-particle base
(micro-scale) according to well-known mechanical
relationships. Typical setups comprise sphere drop, tensile or
compression tests [38,39]. Because bulk materials often
consist of randomly shaped particles, physical calibration is
not always feasible. In such cases, phenomenological
calibration, which is focused on the macro-scale bulk
behaviour of particles, can be used [30]. Thereby, model
parameters are varied until certain macroscopic behaviour of
granular media is observed. Hereby particle stiffness is
usually set at considerably lower values than real stiffness,
resulting in larger time steps and hence faster
calculation [27,38]. By calibrating the bulk behaviour of soft
particles it is supposed that, disregarding unrealistic
stiffnesses, simulation results can be transferred to real
granular media. Examples for corresponding calibration
experiments are triaxial tests, shear tests and bulk
compression tests. If micro-scale behaviour, e. g. magnitude
of contact forces, is of interest, as it is the case for modelling
of drag finishing, physical calibration is necessary.
40 mm
10 mm
v
W
F
n
F
t
bowl
circlepath
workpiece
Fig. 3. Geometries of modelled small-scale process.
The input parameters for the contact model are Young’s
modulus Ei, Poisson’s ratio Ȟi, density ȡi, coefficient of
restitution in normal direction en and friction coefficient ij for
all materials i and material combinations i, j (i, j א Ȃ, i j
with number of materials Ȃ).
Young’s modulus Ec of the spherical ceramic-bond media
that was used for this paper was determined by considering
the media as a multiphase material. Using micro indentation
tests and photo analyses of microscope images of cross-
section polished media, Young’s moduli of the single phases
were calculated. The overall Young’s modulus of media was
estimated at Ec = 77 GPa, the overall Poisson ratio at
Ȟc = 0.19, respectively, both using the approaches presented
by HASHIN AND SHTRIKMAN [40] and ONDRACEK [41].
Friction coefficients ij for all material combinations were
obtained in a scratch test, i. e. the friction coefficient for
ceramic sphere scratching on another one cc = 0.34 and
ceramic sphere scratching on a bearing steel (1.3505) plate
cs = 0.28. The density of the ceramic media was determined
at ȡc = 2656 kg/m³, for bearing steel (1.3505) it was taken
from a data sheet with ȡs = 7610 kg/m³.
Damping behaviour of colliding particles is often
investigated using a sphere drop setup [29,39]. In this paper
the coefficient of restitution in normal direction en is
determined by iteratively comparing experimental with
simulated forces for varying coefficients of restitution in
normal direction en in simulation. This means, a sphere
(ceramic-bond media FSG 6MM BALLS, Walther Trowal)
which is released from a height of 10 mm hits a horizontal
plate (bearing steel 1.3505), that is mounted on a force sensor,
multiple times. For the calibration of the damping behaviour
the first three impacts are considered. The change in the
magnitude of the maximum impact forces and the time since
the first impact were used to fit the coefficient of restitution in
normal direction en. Doing so a coefficient of restitution in
normal direction for ceramic-steel contact en,cs = 0.68 was
estimated. Due to FSG media is not available as a plate the
coefficient of restitution in normal direction for
ceramic-ceramic contact en,cc could not be determined with the
described setup. Instead it is taken from similar investigations
done by KUMAR AND SATHYAN [39] who obtained en,cc = 0.86.
373
Eckart Uhlmann et al. / Procedia CIRP 31 ( 2015 ) 369 – 374
5. Validation
The model was verified using the described small-scale
bowl that was filled with ceramic media and an exemplary set
of process parameters. In the experimental setup, the total
force acting on the workpiece during a drag finishing process
was measured using a KISTLER model 9273 dynamometer.
Normal force Fn from the experiment is shown in Fig. 4a.
The mean experimental normal force is
Fn,exp,mean = -2.59 N ± 0.66 N (± indicates standard deviation).
The base level of the force is caused by static bulk pressure.
In addition there are some peculiar peaks of the normal force
Fn, e. g. at process times tp = 6.250; 10.055 s.
In the simulation the total force acting on the rod is
measured in a workpiece aligned coordinate system, i. e.
normal force Fn is always pointing in the direction of the
workpiece velocities vw and the tangential force Ft is
perpendicular to it, see Fig. 3. Fig. 4a shows the total
simulated forces in normal and tangential direction with
spherical ceramic media and a steel rod for one revolution of
the workpiece for a small scale process with geometries as
described in Fig. 3. For simplicity reasons bowl and rod are
made of the same material. This reduces the number of
occurring material combinations which results in less
parameters that have to be determined experimentally
beforehand. Mean simulated normal force is
Fn,sim,mean = -1.99 N ± 0.63 N, i. e. it deviates from the mean
experimental normal force by about 23.17 %. However,
simulated forces show comparable amplitudes and variation
as well as some distinct peaks as it was the case for
experimental forces. Using new means of contact intensity
analysis, the so called interaction network, it can be shown
that these peaks are caused by the formation of contact chains.
In Fig. 4b-c interaction networks after 4.5 s of process time
are shown. They are coloured according to the magnitude of
the normal contact force of the corresponding interaction,
ranging from green for no contact force to red for highest
overall contact force. The radius of each cylinder is also
scaled by normal contact force, i. e. the cylinder of the
interaction with the maximum contact force has a radius equal
to the mean radius of the balls Rb = 3.27 mm whereas
cylinders of interactions with lower contact forces have
smaller radii. It can be seen that contact forces in direction of
workpiece velocity vW are considerably higher than in other
directions. Furthermore in the region behind the workpiece
the smallest contact forces can be observed. Interestingly the
active zone inside the bulk reaches far ahead of the
workpiece. Additionally, Fig. 4b confirms that contact forces
tend to increase with depth, as expected because of the
increasing pressure caused by the weight of media.
The most important finding from Fig. 4b-c is the formation
of dominant contact chains between rod and bowl. Contact
chains with comparable orientations in xy-plane, see Fig. 4b,
were observed after various process times and typically go in
hand with distinct peaks in simulated normal forces. The
altitude in the bowl and orientation in xz-plane
varies, see Fig. 4c. This phenomenon will be subject to further
investigations, in which a correlation with known contact
mechanisms will be checked. First analyses of microscope
pictures reveal that on some workpiece surfaces deep
scratches can be found on regions pointing in the direction of
the bowls boundary. This is considered as a first indication for
presence of contact chains. It is noted that further efforts have
to be made using other process parameter sets for validation.
0 12.5003.125
process time t
P
6.250 s
0
-3.0
-6.0
simulation experiment
N
-4.5
normal force F
n
4.5
Model setup:
Submersion depth:
d
s
=50 mm
Workpiece velocity:
v
w
= 0.01 m/s
Radius of balls:
R
b
=3.27mm
į
b
=0.15mm
Number of balls:
n
b
=1422
Height of filling:
h
f
=70mm
Contact law and parameters:
Hertz-Mindlin no-slip,
viscous damping
Density:
ȡ
c
=2656 kg/m
3
ȡ
s
=7610 kg/m
3
Young‘s Modulus:
E
c
=77GPa
E
s
=210GPa
Friction coefficient:
cs
=0.28
cc
=0.34
Poisson‘s number:
Ȟ
c
=0.19
Ȟ
s
=0.30
Coefficient of restitution (COR):
e
n,cc
= e
t,cc
=0.86
e
n,cs
= e
t,cs
=0.68
a)
b)
v
w
X
Y
-1.2
0
normal force F
n
at
single interactions [N]
c) X
Z
Fig. 4. Simulation results for small-scale process (a) normal forces for
simulation and experiment (b, c) interaction networks.
374 Eckart Uhlmann et al. / Procedia CIRP 31 ( 2015 ) 369 – 374
6. Conclusion and outlook
After validation of the presented model, it was shown, how
such a model can be used to broaden the understanding of
process fundamentals. For the first time contact chains were
identified by the evaluation of interaction networks. Now it is
possible to determine the number, type and intensity of
contacts between workpieces and abrasive media for finite
areas of any workpiece. This will be subject to further
investigations. As a next step several tests consisting of single
and multi-contact, quasi-static and dynamic scenarios will be
carried out, in which forces, measured in experiments and
calculated in simulations will be compared. Furthermore the
laser displacement probe introduced by HASHEMNIA ET AL.
[18] could be used to validate DEM predictions of impact
velocities on the workpiece surface. Finally contact intensities
calculated with the DEM-model described in this paper
should be linked to the output of the material removal model,
mentioned above [8,20], thus creating a comprehensive
process model.
Acknowledgements
This work is supported by the German Research
Foundation (DFG), project UH 100/180-1.
References
[1] Gillespie L K. Mass Finishing Handbook. 1st ed. New York: Industrial
Press; 2007.
[2] Holzknecht E. Drag Finishing Displaces Manual And Robotic Grinding
& Polishinig. MFN Metal Finishing News 2014:42-4.
[3] Risse K. Einflüsse von Werkzeugdurchmesser und Schneidkanten-
verrundung beim Bohren mit Wendelbohrern in Stahl, Berichte aus der
Produktionstechnik, Aachen: Shaker; 2006;15.
[4] Uhlmann E, Hasper G, Mihotovic V, Fraunhofer Gesellschaft, TU Berlin,
Verfahren und Vorrichtung zum Gleitspanen eines Werkstücks; Patent
DE 10 2009 024 313; 2011.
[5] Uhlmann E, Dethlefs A. Polieren komplexer Bauteile. WB Werkstatt +
Betrieb 2011;6:2831.
[6] Matsunaga M, Hagiuda Y. Researches on barrel finishing. Report of the
Institute of Industrial Science, Tokyo: The University of Tokyo
1967;17:107-55.
[7] Hashimoto F. Modelling and Optimization of Vibratory Finishing
Process. Annals CIRP 1996;45:303-6.
[8] Uhlmann E, Dethlefs A, Eulitz A. Investigation into a geometry-based
model for surface roughness prediction in vibratory finishing processes.
Int J Adv Manuf Technol 2014;75:815-23.
[9] Wang S, Timsit RS; Spelt JK. Experimental investigation of vibratory
finishing of aluminum. Wear 2000;243:147-56.
[10] Ciampini D. Impact Velocity, Almen Strip Curvature and Residual Stress
Modelling in VibratoryFinishing. Toronto: University of Toronto; 2008.
[11] Ciampini D, Papini M, Spelt JK. Characterization of vibratory finishing
using the Almen system. Wear 2008;7-8:671-8.
[12] Ciampini D, Papini M, Spelt JK. Modeling the development of Almen
strip curvature in vibratory finishing. J Mat Proc Technol
2009;6:2923-39.
[13] Yabuki A, Baghbanan MR, Spelt JK. Contact forces and mechanisms in
a vibratory finisher. Wear 2002;252:635-43.
[14] Ciampini D, Papini M, Spelt JK. Impact velocity measurement of media
in a vibratory finisher. J Mat Proc Technol 2007;2-3:347-57.
[15] Kumar PP, Sathyan S. Simulation of 1D Abrasive Vibratory Finishing
Process. Adv Mat Res 2012;565:290-5.
[16] Baghbanan MR, Yabuki A, Timsit RS, Spelt JK. Tribological behavior of
aluminum alloys in a vibratory finishing process. Wear
2003;7-12:1369-79.
[17] Johnson KL, Contact Mechanics, Cambridge: Cambridge University
Press 1985:363-5.
[18] Hashemnia K, Mohajerani A, Spelt JK. Development of a laser
displacement probe to measure particle impact velocities in vibrationally
fluidized granular flows. Powder Technol 2013;235:940-52.
[19] Cariapa V, Park H, Kim J, Cheng C, Evaristo A. Development of a metal
removal model using spherical ceramic media in a centrifugal disk mass
finishing machine. Int J Adv Manuf Technol 2008;1-2:92-106.
[20] Uhlmann E, Dethlefs A, Eulitz A. Investigation of Material Removal and
Surface Topography Formation in Vibratory Finishing. Procedia CIRP
2014;14:25-30.
[21] Kawamoto Y, Kataoka M, Uekusa M, Terumichi Y. Behavior of Two
Kinds of Particles in Rotary Barrel with Planetary Rotation. Multibody
Sys Dyn 2004;3:187-207.
[22] Naeini SE, Spelt JK. Two-dimensional discrete element modeling of a
spherical steel media in a vibrating bed. Powder Technol
2009;195;2:83-90.
[23] Naeini SE, Spelt JK. Development of single-cell bulk circulation in
granular media in a vibrating bed. Powder Technol 2011;1:176-186.
[24] Takahashi Y, Kataoka M, Uekusa M, Terumichi Y. Behavior of Three
Kinds of Particles in Rotary Barrel with Planetary Rotation. Multibody
Sys Dyn 2005;2:195-209.
[25] Wan S, Sato T, Hartawan A. Vibratory Finishing of Immobilized
Cylinders. Adv Mat Res 2012;565:278-283.
[26] Šmilauer V, Catalano E, Chareyre B, Dorofeenko S, Duriez J, Gladky A
et al. Yade Reference Documentation. In: Šmilauer V, editor. Yade
Documentation, The Yade Project, 1st ed.; 2010;. http://yade-
dem.org/doc/.
[27] Hoomans BPB, Kuipers JAM, Briels WJ, van Swaaij WPM: Discrete
particle simulation of bubble and slug formation in a two-dimensional
gas-fluidised bed: A hard-sphere approach. Chem Eng Sci
1996;1:99-118.
[28] Yuu S, Abe T, Saitoh T, Umekage T. Three-dimensional numerical
simulation of the motion of particles discharging from a rectangular
hopper using distinct element method and comparison with experimental
data. Adv Powder Technol;1995;6;4:259-69.
[29] Thornton C, Cummins SJ, Cleary PW. An investigation of the
comparative behaviour of alternative contact force models during elastic
collisions. Powder Technol 2011;3:189-97.
[30] Johnstone M. Calibration of DEM models for granular materials using
bulk physical tests. Dissertation. Edinburgh: The University of
Edinburgh; 2010.
[31] Hertz H. Über die Berührung fester elastischer Körper. J für die Reine
und Angewandte Mathematik 1881:156-71.
[32] Mindlin RD. Compilance of Elastic Bodies in Contact. J Appl Mech
1949;20:259-68.
[33] Mindlin RD, Deresiewicz H. Elastic spheres in contact under varying
oblique forces. J Appl Mech 1953;20:327-44.
[34] Cundall PA, Strack ODL. A discrete numerical model for granular
assemblies. Géotechnique 1979;1:47-65.
[35] Zhang D, Whiten WJ. The calculation of contact forces between particles
using spring and damping models. Powder Technol 1996;1:59-64.
[36] Tsuji Y, Tanaka T, Ishida T. Lagrangian numerical simulation of plug
flow of cohesionless particles in a horizontal pipe. Powder Technol
1992;3:239-50.
[37] Antypov D, Elliott JA. On an analytical solution for the damped Hertzian
spring. EPL 2011;94:1-6.
[38] Grima AP, Wypych PW. Investigation into calibration of discrete
element model parameters for scale-up and validation of particlestructure
interactions under impact conditions. Powder Technol 2011;1:198-209.
[39] Kumar PP, Sathyan S. Simulation of 1D Abrasive Vibratory Finishing
Process. Adv Mat Res 2012;565:290-5.
[40] Hashin Z, Shtrikman S. A variational approach to the theory of the elastic
behaviour of multiphase materials. J Mech Phys Solids 1963;2:127-40.
[41] Ondracek G. Zum Zusammenhang zwischen Eigenschaften und
Gefügestruktur Mehrphasiger Werkstoffe. Z Werkstofftech
1978;9:96-100.
