DEM Modelling of Comminution Devices
DEM Modelling of Comminution Devices
Paul Cleary*.
CSIRO Division of Mathematic01 and InfonnntionSciences. Priurrte Bag 10, Clayton South MDC, Clayton, Vic. 3169, Australia
SUMMARY
Particle size reduction, or milling, is an essential component of mineral processing and is important In other
industry sectors. This needs to be done as efficiently as possible, maximizing mill throughput while
minimizing operating costs. Such milling processes typically use only 1-5 per cent of the supplied energy for
particle breakage, which leaves room for improvement. In discrete element modelling (DEM)of granular
flows the trajectories, orientations and spins of all the particles and objects in the system are calculated and
their interactions with other particles and with their environment are predicted. It is necessary to simulate
particles of many different sizes and densities interacting with complex-shaped objects moving in different
ways. Particle flows in three types of mills; a 5 m ball mill, a 10 mm SAG mill and a 15 cm diameter
centrifugal mill are predicted Charge behaviour, torque and power draw are analysed for a range of rotation
rates from 50 to 130 per cent of the critical speed for the ball mill. Sensitivity of the results to material
properties and size distribution are examined. Radial size segregation is shown to occur and increases
strongly with mill speed. Charge motion and power consumption for the SAG mill are predicted. Compari-
son of simulated flow patterns for the centrifugal mill with high-speed experimental photographs reveals
close agreement. The limitations and restrictions of this type of DEM model are discussed in detail.
Copyright 02001 John Wiley & Sons, Ltd.
KEY WORDS: comminution; discrete element modelling; segregation; ball mill; sag mill; centrifugal mill
1. INTRODUCTION
D E M modelling has become accepted a s a powerful way t o model the flow of granular materials.
Review articles by Campbell [I] Barker [2] and Walton [3] describe the D E M methodology and
many of the findings made about granular flows. The extension to industrial applications has
occurred over the past several years with modelling of ball mills [4-91 dragline excavators [lo]
and mineral sampling from conveyor belts [I 11.
In this paper, previous results of ball mill modelling are extended and the methodology is also
applied to a 10 m diameter SAG mill. Finally, comparison of simulated charge motion in
a centrifugal mill with high-speed photographs of an experimental mill demonstrate that D E M
can produce accurate predictions for such industrial particle flows.
Correspondence to: Paul Cleary, CSIRO Division of Mathematical and InformationSciences, Private Bag 10, Clayton
South MDC, Clayton, Vic, 3169, Australia
' E-mail: [Link]@[Link]
An essential feature of discrete element simulation is that collisional interactions of particles with
each other and with their environment are detected and modelled using a contact force law.
Equations of motion are then solved for the particle motions and for the motion of any boundary
objects with which the particles interact.
For mills of varying types, it is important to be able to model the complex boundary geometry
of the liner, including any lifters. In our two-dimensional code, boundary objects are constructed
from segments which can be line segments, circular segments or discs. Nearly arbitrary shaped
two-dimensional mill cross-sections can be constructed. Arbitrary rigid-body and surface
motions for these objects can also be specified.
In this paper, the particles are modelled as disks with varying size distributions. Non-circular
particles can also be modelled with our DEM code as super-quadrics xn + (y/A)" = s", where the
power n determines the blockiness of the particle and A determining the aspect ratio. Aspect
ratios of up to 12: 1 and blockiness factors of up to 20 can be used.
In the simulation, particles are allowed to overlap and the amount of overlap Ax, and normal
v, and tangential v, relative velocities determine the collisional forces via a contact force law.
There are a range of possible contact force models available in the literature that approximate the
collision dynamics to various extents. We use a linear spring-dashpot model. For more complex
models see References [3, 121. Figure 1 shows diagramatically the collisional force model used
here. The normal force
consists of a spring to provide the repulsive force and a dashpot to dissipate a proportion of the
relative kinetic energy. The maximum overlap between particles is determined by the stiffness k of
the spring in the normal direction. Typically average overlaps of 0.1-1.0 per cent are desirable,
requiring spring constants of the order of lo6-10' N/m. The normal damping coefficient C, is
chosen to give a required coefficient of restitution E (defined as the ratio of the post-collisional to
pre-collisional normal component of the relative velocity)
Normal Force -
Fn= knAx + C,vn
Tangential Force F, =min( pF, , k, v,dt + C,v, )
Figure 1. The contact force model involves a spring and a dashpot in the normal direction and an
incrementing spring and dashpot limited by the sliding friction in the tangential direction.
Copyright 02001 John Wiley & Sons. Ltd. Int. J . Numer. Anal. Meth. Geomech., 2001; 25:83-105
MODELLING COMMINUTION DEVICES USING DEM
where
In ( 4 -
and mij = mimj
+
?=- J-) mi mi
is the reduced mass of particles i and j with masses mi and mi, respectively. This arises from the
analytic solution of the normal collision equation (1) for two such particles. Each particle can
have a different coefficient of restitution (and corresponding damping constant), but we generally
use the same values for all the particles because little experimental data on particle to particle
variations in E is available.
The tangential force is given by
where the integral of the tangential velocity v, over the collision behaves as an incremental spring
that stores energy from the relative tangential motions and represents the elastic tangential
deformation of the contacting surfaces. The dashpot dissipates energy from the tangential motion
and models the tangential plastic deformation of the contact. The total tangential force (given by
the sum of the elastic and plastic components) is limited by the Coulomb frictional limit at which
point the surface contact shears and the particles begin to slide over each other.
The contact force model requires both the coefficient of elasticity and friction for the particles
to be supplied. Such quantities are in practice very difficult to obtain experimentally and vary
from particle to particle and during the flow. Typical values are therefore chosen and the
sensitivity of the numerical results of these choices is examined. - -
where xi, ui and Fii are the position, velocity and collisional forces on particle i, and Oi and oi
are the particle orientation and spin produced by the moments M i j . Here g is the gravity
force.
The integration scheme is a second-order predictor-corrector. Between 20 and 50 timesteps are
required to accurately integrate each collision. This necessitates very small timesteps (typically
10-3.-10-6s depending on the controlling length and time scales of each problem). If the
integration is accurate then the coefficient of restitution from which the damping coefficient was
chosen should be recovered.
Copyright 0 2001 John Wiley & Sons, Ltd. Int. J . Numer. Anal. Meth. Geomech., 2001; 2583-105
86 [Link]
Quantitative predictions about wear rates and distributions, collision forces, dynamic loads on
boundaries, power consumption, torques and flow rates, sampling statistics, mixing and segrega-
tion rates and many other quantities can be made from the information available in
DEM simulations. For more details on the simulation method and on the data analysis see
Reference [13].
Ball mills are used in mineral processing to help liberate minerals. They consist of a rotating
cylinder up to 5 m in diameter. The rock generally comes from a crusher or perhaps a semi-
autogenous (SAG)mill and is fed into one end of the mill. Grinding media consisting of steel balls,
typically with diameters from 50 to 200 mm, are already present in the mill along with rock that is
still too large to pass out through the extraction grate at the opposite end. Particles migrate
slowly along the length of the mill while circulating rapidly in the plane orthogonal to the mill
axis. The motion in this plane is assisted by regularly spaced lifter bars attached to the mill shell.
Replaceable liner plates are bolted to the shell between the lifters.
A typical 5 m diameter ball mill consumes around 3-4 MW of power and has a very low energy
efficiency. Significant economic and environmental benefits could be obtained by improving the
efficiency even modestly. There are also large costs involved in replacing the liners (commonly
made from expensive wear-resistant cast molybdenum stainless steels) arising both from the liner
replacement cost and from lost production. Further benefits can be obtained through higher
downstream recoveries if the exit particle size distributions can be made closer to the optimum for
the subsequent flotation processes.
Here we simulate a ball mill with 23 lifter bars spaced circumferentially around the mill shell.
Wear plates protect the mill shell between lifters. The ball mill rotates clockwise at various
constant fractions N of the critical speed of 19.5 rpm, at which the particles begin to centrifuge.
The bulk of the charge in such a mill is rock of various sizes with a smaller proportion of steel
balls. Previously, Cleary [9] has shown that the inclusion of rocks in the simulationsis important,
since predictions based on a charge consisting entirely of balls predict around 10 per cent higher
power consumption. Here we use a realistic charge consisting of steel balls with diameters
uniformly distributed between 50 and IOOmm and rocks with diameters uniformly distributed
between 5 and 50 mm. The rock frequency is 10 times higher than the ball frequency. The charge
therefore consists of 57 per cent (by volume) rock. This charge composition is used for most of the
ball mill predictions made in this paper. The torque required to maintain the specified rotation
rate and the resulting energy consumption are predicted. Segregation in the charge is also
evaluated.
We use default values of 0.3 for the coefficient of restitution E and 0.75 for the friction coefficient p.
Copyright ((9 2001 John Wiley & Sons, Ltd. Int. J . Numer. Anal. M r t h . Gromech., 2001; 2583-105
MODELLING COMMINUTION DEVICES USING DEM
Figure 2. Charge shape predicted for a ball mill filled with rocks and steel balls for: (a) N = 60
per cent, (b) N = 70 per cent, (c) N = 80 per cent and (d) N = 90 per cent. The particles are
shaded according to their speed.
and the charge fluidizes. Most of the material then flows down the steep free surface into the
horizontal toe region. This is generally described as the cascading part of the flow. A small
amount of material in the shoulder region is trapped between the lifters and is ballistically ejected
into the mill cavity. This material drops onto the charge in the toe region and sometimes onto the
unprotected liner of the mill. This part of the flow is generally termed cataracting.
The particles are shaded according to their speed with light grey being low speed and dark grey
being high. The fastest particles are the cataracting ones just as they impact the toe of the charge.
The next fastest group of particles are the cascading or avalanching particles composing the free
surface. They reach a peak speed just below the centre of the mill and slow before arrival at the
toe. The next fastest group of particles are the ones rotating with the mill near the shell. There is
a semi-circular band (indicated by the light grey) of particles that move quite slowly. At the very
centre of this band is the centre of the circulation pattern for which the particle speed is zero.
Changes in the shading indicate shear flow that leads to particle size reduction by both abrasion
Copyright (02001 John Wiley & Sons, Ltd. Inr. J . Numer. Anul. Merh. Geomech., 2001; 2583-105
88 [Link]
and attrition. Such shear regions are clearly visible between the faster moving particles around
the outside of the mill and the slowly moving interior particles and also between the cascading
surface particles and the slow-moving interior. The largest amounts of shear (indicated by the
sharpest changes in grey scale) occur in toe region where material abruptly changes its motion
from cascading down the surface to being pulled under the bulk of the charge by the lifters and
swept to the left.
At N = 70 per cent, (Figure 2(b))the shoulder position moves slightly higher and the toe moves
slightly lower. The inflexion point on the free surface moves much closer to the centre of the mill
as the free surface profile becomes almost bi-linear. The amount of cataracting material has
increased significantly and it possesses a clear banded structure. Each band is slightly curved, is
almost vertical buttipped to the right a bit at the top and leads to the front face of one of the lifters
above. Each lifter leaving the shoulder region acts somewhat like a filled bucket. As the lifter rises
and its angle changes, the particles trapped inside gradually pour out with the last particles
departing as the lifter reaches vertical. Each band corresponds to the material that poured from
the space in front of the lifter to which it leads. The impacts of cataracting material with the
charge and mill liner in the toe region are therefore not continuous in time, but are bunched
together and produce a cyclic variation in the impact forces in this region. For some rotation
rates, spectral analysis of the instantaneous power draw detects a signal with a frequency given by
the product of the mill frequency and the number of lifters which corresponds to the one that such
cyclic impacts should produce. The distribution of particle speeds, and therefore shear, within the
charge is little changed from the lower speed. The relative speed of the particles moving around
the shell with the lifters compared to the speed of the cascading free surface particles has increased
sIightIy.
Similar changes continue to occur when the mill speed is increased to N = 80 per cent (see
Figure 2(c)). The amount of cataracting material has again increased strongly and their trajecto-
ries take particles close to the lifters above. The band structure of these particles is more clearly
visible. The shoulder again moves slightly higher and the toe marginally lower. The speed of the
particles moving with the lifters and shell at the bottom again increases relative to the speed of the
cascading particles. The semi-circular band of slow-movingparticles has become better defined as
the shear to either side becomes more concentrated. Near the shoulder region the centrifugal force
causes the free surface to become less sharply defined as this balances the gravitational force.
By N = 90 per cent (Figure 2(d)) the outermost cataracting particles barely move away from
the shell and collide with the backs of the lifters slightly above the centerline of the mill. This
would increase wear on the lifters and is very undesirable. The fluidization of the upper cascading
free surface by the centrifugal force is now strong and the concentration of the shear into two
semi-circular bands from the toe to the shoulder continues. For higher rotation rates particles
become permanently trapped between the lifters and are then centrifuged for N > 100 per cent.
The depth of the centrifuged layer increases with N.
Copyright G 2001 John Wiley & Sons, Lid. Int. J . Numcr. Anal. Merh. Geomech.. 2001; 2583-105
MODELLING COMMINUTION DEVICES USING DEM
Figure 3. (a) Average torque and power as a function of mill speed N and (b) power draw for particles with
different restitution and friction coefficients. The fill level was 50 per cent.
seconds gives an invariant measure of the power draw that can be meaningfully compared to
laboratory and field measurements.
Figure 3(a) shows the variation of the average torque and power with mill speed for the ball mill
described above with a 6ll level of 50 per cent. The torque increases slowly until the peak is
attained around N = 80 per cent. As the amount of cataracting and then centrifuging increases,
the balance of the charge improves and the torque required to maintain the asymmetric charge
position declines steadily. The power is given by the product of this torque curve with the mill
speed. It has a much rounder peak centred on N = 105 per cent before dropping sharply for
higher rotation rates. For a mill operating at a typical 80 per cent of critical, the power predicted
by the model is 3 x lo5 W/per m of mill length. For a typical 7 m long mill, the power required is
therefore predicted to be 2.1 MW. This compares well with the peak power of 3.3 MW that is
typical for a mill of this size.
3.2.1. Sensitivity to materialparameters. The friction and restitution coefficientsfor the charge are
not well known and vary between mines and within mines and are also likely to be modified by
the presence of water in the mill. We therefore examine the sensitivity of the power prediction to
variations in the material parameters. Figure 3b shows the power draw as a function of mill speed
for a range of material properties. For this analysis we consider only the steel balls. The base case
uses E = 0.3 and p = 0.75 and is shown by the circles. .Increasing the coefficient of restitution to
0.5 has only a modest effect on the power draw curve. There is essentially no change at low speeds
N < 75 per cent. For intermediate speeds (at which the mills are most commonly operated) and
high speeds there is a reduction in power consumption between 2 and 6 per cent. Increasing E to
0.8 produces even smaller changes in the power draw.
It may seem surprising that although the change in restitution co.:fficient from 0.3 to 0.5
reduces the energy dissipation per collision from 91 to 75 per cent it produces a much smaller
reduction in the overall power consumption. This is easily explained by understanding that in all
the non-cataracting part of the charge, the particles are densely packed and these particles
undergo very rapid series of collisions (on millisecond time scales) with their neighbours. It does
Copyright 0 2001 John Wiley & Sons, Ltd. Int. J . Numer. Anal. Meth. Geomech., 2001; 2583-105
90 P. CLEARY
not matter whether 2, 3 or 4 collisions in rapid sequence are required to dissipate all relative
motion between the particles. All the relative motion is dissipated much faster than any time scale
for macroscopic movement of the charge mass. For the cataracting material this is not true. These
particles are more sparsely spaced and can move reasonable distances between collisions. For
N < 75 per cent there is little cataracting material so the power draw is effectively independent of
c. For higher rotation rates the amount of cataracting material (for which the power consumption
is sensitive to E ) is a small proportion of the total charge. Therefore it is reasonable to expect
a reduction in power draw for the less dissipative materials, but proportionally much less than
one may have estimated based on the energy losses for one collision.
Reducing the friction coefficient from 0.75 to 0.25 (a reasonable value considering the lubricat-
ing effects of the interstitial fluid), we again find only a small change in the power draw below 75
per cent of the critical speed. One may have expected large decreases in the friction to significantly
decrease the power consumption, but this change does not much affect the dynamics of the charge
in the lower left part of the mill. Figure 2 demonstrates the distribution of shear or sliding.
Changing the friction coefficient cannot alter either the speed of the particles moving with the mill
liner or the speed of the stationary eye of the charge circulation. The friction coefficient does
determine how strongly concentrated this shear is, but it cannot change the total amount of shear
or sliding. The overall energy dissipation is determined by the total amount of sliding, not by its
distribution. In the avalanching layer, lower friction can reduce power consumption, but these
particles are again a minority of those in the charge and the particle contacts are generally not
long duration sliding ones. Lower friction therefore leads to only a small reduction in the energy
consumption.
At higher sub-critical speeds and at super-criticalspeeds power consumption actually rises with
the decreases in friction. Again this may be unexpected if only isolated collisions were considered.
However, the lower friction increases the mobility of loose particles near the surface of the charge
and particularly in the centrifuged layer that appears at high speeds. These more mobile particles
then dissipate more energy via their normal contacts. A decrease in the friction therefore leads to
a small decrease in the dissipation in the tangential direction, but leads to a much larger increase
in the dissipation in the normal direction due to increased particle mobility. At N = 130 per cent
the power consumption is actually doubled.
3.2.2. Eflect o f j l l level. Next, we examine the effectof changes in the fill level of charge in the mill
on the power draw. Figure 4 shows the power draw as a function of mill speed for fill levels
between 10 and 50 per cent. The same 5 m ball mill and the same type of charge (consisting of
both balls and-rocks with a bottom size of 5 m).
Examining the power consumption predictions for the five fill levels shows that over most of
the rotation speed range (Figure qa)), a decrease in the fill level leads to a decrease in the power
consumption. This should be expected, since all the material in the mill contributes to power
consumption to one degree or another. So any reduction in fill level for all but the highest
rotation rates should reduce the power consumption.
More importantly, the relative size of the decrease in power consumption increases as the fill
level becomes smaller. The largest decrease is by a factor of one-half and occurs between fill levels
20 and 10 per cent. These changes can be more easily understood by considering the distribution
of the charge within the mill (shown in Figure 5 for these five fill levels). When decreasing the fill
level from 50 to 40 per cent, the dominant change in the charge structure is the lowering of the
location of the avalanching free surface. There is also a marginal increase in the amount
Copyright 0 2001 John Wiley & Sons, Ltd. lnt. J . Numer. Anal. Meth. Geomech., 2001; 25:83-105
MODELLING COMMINUTION DEVICES USING DEM
I ' I - I . 1 '
-
40% bad
0-30Xbod
-
+ 2M bod
.-1Mbad
- -
I I I- I - I I I I
60 80 -100 120 60 80 100 120
(a) Mill Rotation (X critical) (b) Mill Rotation (X critical)
Figure 4. For ball mill fill levels between 10 and 50 per cent, (a) power draw and
(b) power consumption/tonne.
cataracting material because some particles which would have collided with the 50 per cent
height free surface remain ballistic in the enlarged space above the charge. The centre of mass
(COM) of the charge is below and to the left of the centre of the mill. Crudely, the torque
required to move the charge is the product of the weight of the charge and the horizontal
distance between the centre of mass of the charge and the centre of the [Link] is the torque
required to maintain the charge mass off centre. The power consumed is the product of this
torque and the rotation rate. When the fill level is reduced from 50 to 40 per cent the COM of the
charge moves only marginally further to the left. This is because the material removed from the
charge in the 50 per cent case is very close to the centre of the mill and has made very little
contribution to the charge being off centre and therefore very little contribution to the overall
power consumption.
Reducing the fill level from 40 to 30 per cent again lowers the free surface level. The shoulder
and toe locations retreat only slightly indicating that the smaller amount of material is spread out
more thinly over only slightly smaller a fraction of the mill circumference. As before the cascading
material that has been removed is the material closest to the centre of the mill and makes the
lowest contribution to power consumption. Comparing the material removed when changing the
fill level from 40 to 30 per cent with that removed when going from 50 to 40 per cent, we find that
horizontal distance of the centre of mass of this removal material from the centre of the mill is
more than double the size for the second removal than it was for the first. The amount of material
removed is the same in both cases. So the more than doubling of the moment arm of the removed
material more than doubles the reduction in power consumption.
This behaviour continues with each subsequent reduction of the fill level, with the COM of the
remaining charge approaching the shell of the mill and each removal making a proportionally
larger reduction in the power consumption.
Examining the charge distributions shown in Figure 5 more carefully, we observe several
trends:
. The free surface profile changes from being bi-linear to being increasing curved.
,... . ,
Copyright 02001 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech., 2001; 2583-105
.:
,
[Link]
Figure 5. Charge shape predicted for ball mill rotating at 80 per cent of critical and filled with rocks and
balls with fill levels: (a) V = 50 per cent, (b) V = 40 per cent, (c) V = 30 per cent, (d) V = 20 pel: cent and
(e) V = 10 per cent. The particles are shaded according to their speed.
. The shoulder position moves lower but only slowly and the toe moves to the left but again only
slowly; Therefore the proportion of the mill circumference covered by the charge decreases
much more slowly than the charge volume, leading to a charge mass that is progressively
thinner.
Copyright 02001 John Wiley & Sons, Ltd. Int. J . Numer. Anal. Meth. Geomech., 2001; 2583-105
MODELLING COMMINUTION DEVICES USING DEM 93
The velocity of particles near the shell remains that of the mill shell and the velocity of the
cascading particles decreases only gradually with the slow lowering in the shoulder location.
The combination of decreasing charge depth with relatively constant velocity difference across
the charge leads to steadily increasing shear rates within the charge and should therefore lead to
more efficient grinding. The increase in the shear rates can be seen in Figure 5 by the increasing
gradients in the shading of the particles. In frame d, the fast-moving particles at the bottom and
top of the charge are very close to the curved central region of slow-moving or stationary
particles.
For a fill level of around 10 per cent there are almost no cascading particles. The charge depth is
barely deeper than the depth of the lifter and almost all the particles move in the cataracting
regime.
The power consumption per tonne (known either as power density or specific power consump-
tion) is more interesting than the total power draw and is shown in Figure 4b. This shows that
(except for the highest rotation rate) there is a considerable increase in specific power consump-
tion, by between 15 and 20 per cent, when the fill level declines from 50 to 40 per cent. So although
the total power consumed declines the amount of power consumed per unit mass actually
increases. The rate of particle grinding is actually dependent on the local rate of energy
dissipation or the specific power rather than on the total power consumed. This demonstrates
that more efficient grinding will occur for a lower total energy consumpiion with a 40 per cent full
mill, rather than one jammed to near or above capacity at 50 per cent. This clearly illustrates the
well-known observation that an overfilled mill suffers from reduced grinding efficiency.
Decreasing the fill level to 30 per cent again increases the specific power consumption leading
again to an increase in grinding efficiency. For all three fill' levels the actual peak power
consumption occurs at 105 per cent, although the shapes of the power curves near these peaks
suggest a very mild tendency for higher power consumption to increase slightly for higher
rotation rates. Reducing the fill level to 20 per cent produces a significant change in the structure
of the specific power curve. The curve has developed two peaks. The major one has moved to
a significantly lower rotation rate at N = 85-90 per cent of critical and a minor peak at 120 per
cent. This is the first time that we have observed the location of the peak power consumption
being sensitive to any simulation variable. Previously, it has been shown that the peak is
insensitive to material properties of the charge and in Reference [9] it was shown not to depend
on the particle size distribution of the charge. The maximum specific power consumption at 20
per cent fill level is slightly higher than the peak for the 30 per cent fill level. For rotation rates of
N 2 98 per cent it appears that mill efficiency will be higher for a 30 per cent fill and for lower
rotation rates it will be more efficient at lower fill levels.
Decreasing the fill level to 10 per cent produces a specific power curve whose peak has
increased in magnitude and whose location has moved down to N = 70 per cent. So for N < 82
per cent the grinding rate is largerfor a fill level of 10 per cent and for higher rotation rates is
higher for a fill level of 20 per cent.
Note that operating a mill at a 10 per cent fill level and a speed of 70 per cent would have two
serious drawbacks:
1. Most of the material for this fill level is cataracting and so most of the breakage results from
high-energy impacts with the liner leading to serious wear problems.
2. Although the grinding rate is maximized, the total amount of material in the mill is small, so
the mill output is relatively low.
Copyright (:I 2001 John Wiley & Sons, Ltd. Int. J . Numer. Anul. Meth. Geomech., 2001; 2563-105
94 [Link]
Real mills are operated at slightly higher rotation rates and target fill levels between 25 and 35
per cent. This gives somewhat lower specific power consumptions (between 13 and 14 kW/tonne
depending on the precise fill level) but much higher mill throughput. Optimizing mill performance
is about choosing fill levels and rotation rates that maximize the specific power consumption and
the mill throughput. In an environment of continuously changing charge, this is not a simple task.
These DEM predictions suggest that for rotation rates of N 75 per cent the maximum
specific power consumption (and therefore grinding rate) should occur for mill levels between 20
and 30 per cent. This seems to be consistent with industrial practice.
3.2.3. Eflect of ball or rock densiv. In the previous simulations the steel balls and the rocks
were given the same-density of 4000 kg/m3. Now we examine the changes in the flow and the
power consumption produced by increasing the ball density to its proper value of 7800 kg/m3.
Figure 6 shows the power draw for same configuration as used earlier (with fill level 50 per cent)
for the two cases of ball density. The increase in ball density increases the weight of the charge and
therefore the power draw, as one would expect. The specific power consumption is very close to
identical for the two cases. A very small reduction for higher rotation rates (N 2 80 per cent is
probably not significant). No differences are detectable in the charge shape or flow either. This
suggests strongly that the density of specific parts of the charge (either ball or rock) has negligible
effect on either the flow in the ball mill or the specific power consumption. The use of steel balls
simply increases the overall power consumption by increasing the weight of the charge. There is
no obvious reason why this should produce benefits in terms of improved grinding. It would seem
beneficial in terms of power consumption to minimize ball use, subject to the other requirements
for ball inclusion in mills.
Steel balls are generally used because they have good wear properties and maintain the number
of particles in the mill that are in the range required for attrition. If this requirement can be
satisfied for a particular ore by the larger incoming ore particles, then there are potential benefits
to lowering the loading of balls and running the mill more like an A G mill than a SAG mill. This
will depend entirely upon the specific details in each mill. These predictions only indicate that
A -
N , , , I , , , -
C
C
0
.;:
-3 2 -
ln
f
(U m -
P
1
U
\
L LO-
o - Ball density 4000 kg/m3
(U n - Ball density 7800 kg/m3 -
3
n
Figure 6. (a) power draw and (b) power consumption/tonne when the balls in the ball mill have
densities of 4000 and 7800 kg/m3.
Copyright .C 2001 John Wiley & Sons,Ltd. Int. J . Numer. Anal. Meth. Gcomech., 2001; 2583-105
MODELLING COMMINUTION DEVICES USING DEM 95
minimizing the use of steel will have benefits in reduced power consumption provided other
constraints on the particle size distribution can be met in other ways. They also indicate that the
effect of heavier balls in the cataracting stream breaking particles in the toe region better by
higher-energy impacts is negligible.
Figure 7. Segregation state of the charge after 60 s of motion for: (a) N = 50 per cent, (b) N = 70 per cent,
(c) N = 90 per cent and (d) N = 110 per cent.
Copyright (02001 John Wiley & Sons, Ltd. Inr. J . Numer. Anal. Merh. Geomech., 2001;25:83-105
96 P. CLEARY
The segregation state can be measured by calculating the coefficient of variation of the
distribution of local average diameters (see Reference [13]). For N = 50 per cent, this quantity q is
mostly slightly larger than the randomly mixed level, confirming the small amount of segregation
observed above. For N = 70 per cent, q remains below the randomly mixed level indicating that
there is no discernible segregation. The time series of q is shown in Figure 8 for N = 90 and 110
per cent. The upper dashed line is the fully segregated limit whilst the lower curve is the randomly
mixed limit. The central curve is q and shows the amount of segregation. In both cases the charge
segregates rapidly reaching its asymptotic level after around 20 s. The simulations are continued
to 240 s in order to demonstrate that the amount of segregation remains unchanged. The average
segregation level is 13 per cent for N = 90 and 40 per cent for N = 110 per cent and shows that
the amount of segregation increases rapidly with mill speeds above 80 per cent.
There seem to be two basic mechanisms producing size segregation in mills. At low speeds the
segregation is produced by percolation of fines through the surface avalanching layers. This
causes large particles to migrate to the outside of the charge (as seen in Figure 7(a)). Within
a deforming or flowing granular media the second mechanism causes large particles to move
against the acceleration field and smaller particles to move in its direction (this is the same
mechanism that causes segregation of particles in a vibrating container - the muesli effect [6]).
In a mill at higher speeds the centrifugal force is comparable to gravity and is directed radially
outwards from the centre of the mill. This force causes the large particles to move into the centre
of the charge and the fines to the outside (as seen in Figure 7(c) and (d)). The two mechanisms
cause particles to move in opposite directions, so for intermediate rotation rates these two
mechanisms cancel and leave a randomly mixed charge (Figure 7(b)). -
Segregation will have a detrimental effect on the mill performance. Since the steel balls are
confined to the slow-moving low shear central region and there are few rocks remaining there, the
principal ball milling attrition mechanism of small rocks being nipped between balls is much
reduced. The secondary mechanism of balls in the cataracting stream falling long distances and
crushing rocks with high-energy impacts in the toe region is largely inhibited.
Copyright C;)
2001 John Wiley & Sons, Ltd. Inr. J . Numer. Anul. Merh. [Link]., 2001; 25:83-105
MODELLING COMMINUTION DEVICES USING DEM
Semi-autogenous grinding (SAG) mills are physically similar to ball mills but much larger and are
fed with rock with a wide size range, generally straight from a crusher, and have a small load of
steel balls inside. The grinding mechanisms in SAG mills are quite different to those in ball mills
and they are notoriously difficult to control. The highly variable material entering the SAG mill
means that the grinding efficiencies for different size fractions are continually changing. These
variations can lead to the mill over filling and choking or to over grinding and producing low
charge levels that lead to damage as the steel balls impact the liner rather than rock. The only
control parameters generally available to manage the mill are the feed and water flow rates and
sometimes the mill speed. Substantial benefits are potentially available to mill operators if the
behaviour of SAG mills can be understood.
We model a 10 m SAG mill with 72 lifter spaced equally around its circumference. These lifters
are about 200 mm high. The critical speed for such a mill is 13.38 rpm. The charge is composed of
rocks with size range 10-100 mm and density 2500 kg/m3 and steel balls with sizes 50-100 mm
and density 7800 kg/m3. The relative frequency of rock to ball is 10: 1. The proportion of steel
balls is therefore about 16 per cent by volume. The mill is loaded to 40 per cent.
Figure 9 shows the charge distribution for four different rotation rates. At N = 60 per cent of
critical speed, the amount of cataracting material is small and it falls onto the upper cascading
slope of the charge far short of the toe. The cascading material moqes much faster than the
material adjacent to the liner. Comparing this profile to that ofthe ball mill for the same fraction
of critical speed (Figure 2(a)), one fmds that the free surface is much more sharply bi-linear and the
cataracting material much less important. The semi-circular low velocity/stationary zone is
similar in shape and location for both the ball and SAG mills.
Increasing the rotation rate produces similar changes for the SAG mill as for the ball mill. The
free surface becomes progressively more sharply bi-linear. The toe and shoulder positions move
only marginally. The amount of cataracting material increases and moves through a progress-
ively larger fraction of the mill. The rotation rate at which the cataracting material begins to
impact the liner directly is around 70 per cent for the ball mill and nearly 80 per cent for the SAG
mill. Significant buildup of material at the back of the lifter is observed above the level of the toe
for a speed of 90 per cent. The patterns of change in the velocity distributions within the charge
are also similar with the shear steadily increasing with rotation rate. Broadly, the doubling of the
mill diameter and the increase in the number of lifters has had only marginal effects on the flow
pattern within the mill.
Figure 10 shows the power consumption and specific power consumption for both the ball and
SAG mills. Total power consumption in the SAG mill is predicted to be significantly higher for
the SAG mill. A typical 7 m long mill of this diameter operating dry and at a speed of 80 per cent
critical would consume about 9.1 MW of power. This is a consistent with real mill power
consumption.
The largest contribution to the greater power consumption of the SAG mill relative to the ball
mill is due to the four-fold increase in the volume of the charge due to the doubling of the mill
diameter while maintaining the same fill level. The specific power consumption (Figure 10(b))
removes this effect and shows that the increase in mill diameter has lead to a real increase of 45
per cent in the power consumed per tonne of charge. The grinding rate is proportional to this
energy which is dissipated by inter-particle collisions. This demonstrates why the industry trend
of increasing the mill diameter is expected to lead to improved grinding rates. Essentially the balls
Copyright 0 2001 John Wiley & Sons, Ltd. I n ( . J . Numer. Anul. Meth. Geomech., 2001; 2583-105
P. CLEARY
Figure 9. Charge shape predicted for a 10 m SAG mill with 72 lifters for: (a) N = 60 per cent, (b) N = 70 per
cent, (c) N = 80 per cent and (d) N = 90 per cent. Particles are shaded according to their speed.
A - 10 m YG mill with 72 l i e m
I
Figure 10. Comparison of a 10 m SAG mill with a 5 m ball mill for mill speeds between 50 per cent and 130
per cent of critical speed, (a) power consumption (b) specific power consumption.
Copyright 0 2001 John Wiley &: Sons: Ltd. I n t . 3. Numcr. A n d . Meth. Gcomeclt., 2001; 2583-105
MODELLING COMMINLJTION DEVICES USING DEM 99
are raised higher in the larger mill and therefore have more potential energy to dissipate via
collision and produce more particle breakage.
Centrifugal mills are used for fine grinding or milling. These are very important processes in
industries ranging from mineral processing, where final particle sizes sometimes need to be
reduced to 1 pm or lower in order. to liberate key mineral components, to pharmaceuticals where
key constituents need to be reduced to similar sizes for blending with other powders to produce
pharmaceutical products.
A centrifugal mill is a cylindrical vessel that is attached to a rotating arm that revolves at high
speed. The mill vessel is also able to rotate at a controlled speed around its centre of mass at the
end of the arm. This allows a range of motions that generate extremely high centrifugal forces
(with rotation rates between 400 to 1000 rpm) thereby producing very fine grinding. Previously
Hoyer [14] has reported the results of high-speed photography which was used to record the
behaviour of the charge in such a mill under a range of conditions.
This is a particularly good validation test problem because for the case where the mill arm and
chamber rotation are equal and opposite (so the chamber maintainsthe same orientation with the
ground) there are two distinct flow regimes and a well-deiined bifurcation between them:
..For higher fill levels the charge has a constant profile that simply rotates with the mill rotation.
For lower fill levels the charge adopts a pointed three-sided structure that flops around the
inside of the mill spraying particles through the middle of the mill vessel.
The flow transition was observed experimentally [14] to occur at 34 per cent for a ratio of
diameter of mill arm to chamber diameter of 0.4.
The constancy of the charge profile for higher fill levels, the bifurcation leading to the change in
flow behaviour at lower fill levels with its characteristic motion and actual steady charge shapes
provides three key elements for the simulations to reproduce.
The test mill used in the experiments had a chamber that was 150 mm in diameter and 150 mm
deep (in the axial direction) with four equi-spaced lifter bars. The particles used were plastic beads
that were essentially somewhat rounded cylinders mean particle diameter of about 3 mm.
Figure ll(a) shows the experiment and simulation for a 75 per cent fill level, when the toe
charge is near the top lifter. The simulation and experimental profiles are almost indistinguish-
able, with the charge just separating from the bottom lifter and the contact point of the toe with
the mill liner being the same distance from the upper lifter. No variation of the simulation profile
is observable with time in the same way as there is negligible change in the experimental profile.
For the 75 per cent fill level the flow is very steady and closely matches the experiment.
For the 50 per cent fill level, the experimental and simulation results are compared in Figure
1l(b). The underlying profile of simulated charge is very similar (ignoring the few outlying
particles) to the experimental one. In each case, the charge mass has just separated from the lifter
on the right and similar separations of charge from lifter are observable. The amount of curvature
on the charge mass is also quite similar. This experimental photograph shows that the flow in the
mill is not completely two dimensional. The particles adjacent to the perspex are quite clearly
defined and the simulation profile near the toe is very similar to this. But to the left of the clearly
defined particles in the photograph is blurred mid-grey region, corresponding to particles further
Copyright 6 2001 John Wiley & Sons, Ltd. Int. J . Nunter. Anal. Meth. Geomech., 2001; 25:83-105
[Link]
Figure 11. Charge motion in a centrifugal mill rotating at 695 rpm for (a) and (b) 75 per cent
fill level, (c) and (d) 50 per cent fill level and (e) and (f) 25 per cent fill level. Experimental
photographs are on the left and simulations are on the right.
into the mill chamber. This front face of the charge rotates moderately in advance of the material
at the perspex mill end. The lagging of the profile at the end is due to the additional frictional
effects of the perspex end wall. The DEM prediction again demonstrates invariancy of the charge
profile, which simply rotates with the mill rotation. The amount of variation from the experi-
mental profile, although small, seems to increase as the fill level in the mill declines.
Hoyer [14] showed that there is a change in the dynamic behaviour of the charge for fill levels
around 34 per cent for GID = 0.4. Figure ll(c) compares the profiles of the experiment and the
Copyright 02001 John Wiley & Sons, Ltd. Int. J . Numer. Anal. Meth. Geomech., 2001; 25:83-105
MODELLING COMMINUTION DEVICES USING DEM 101
simulation at the same point in the mill rotation. The experiment shows the three-sided structure
that is characteristic of these low fill levels. There are also quite a number of outlying particles
being sprayed around. The simulation shows the same characteristic shape with the particles at
the top of the mill occupying the same parts of the mill vessel. In the simulation, this three-pointed
charge mass is observed to 'flop' around the inside of the mill chamber in the same way as
occurred in the experiments, confirming the DEM simulations' ability to predict the bifurcation
between flow regimes.
The DEM implementation used here contains several approximations and restrictions. It is
important to understand their impact so that the usefulness of the results can be interpreted.
The modelling described here has been in two dimensions. This raises a question about how
applicable two-dimensional results are to particle flows in real three-dimensional systems. The
answer to this is determined largely by the geometry of the objects with which the particles
interact and the information that is to be extracted from the model. In general a particle flow in
a geometry that is inherently two dimensional and which is far away from walls in the third
dimension can be modelled with reasonable accuracy p i n g a 2D DEM model. This limits the
application to flows such as those described in-this paper. The ball, SAG and centrifugal mill
models all correspond to flows in cross-sections through the middle of each mill. This is
a reasonable approximation because of the axial symmetry of the mill chambers and because the
flow in the axial direction is much slower than in the radial or circumferential directions. In the
ball and SAG mills, water is flushed through the charge to remove the smaller particles through
an end wall grate. Generally this flow velocity is small compared to the circumferential speed
produced by the high mill rotation rate, so the out of plane motions are relatively unimportant
over the time scales that are involved in the simulations presented here. If the information to be
calculated from the model involved simulation of times that were comparable to the residence
times of particles in the mill or the fluid flow velocity was increased substantially, then the axial
particle motion would be important and a full 3D model would be required. For the centrifugal
mill there is no axial flow, so the 2D model is valid over a much wider range of conditions.
Three-dimensional simulation of such mills is possible, but the benefit of this depends on what
information one seeks to predict. One can model either the entire mill (which is a huge
computational problem) and truncate the particle size distribution at unreasonably large dia-
meters or one can simulate a slice of the mill and use periodic boundaries in the axial direction.
This second approach produces feasible size computations, but the additional reality in the model
is not actually as great as one might expect. The major out of plane (axial) motions in such mills
are produced by fluid flow from the feed end to the grate end. This cannot be captured by the 3D
slice model since the flow rates through the periodic boundaries for each size fraction are not
known nor are the composition gradients within the mill. The only real improvement obtained by
the slice model is in the disorder of the three-dimensional microstructure. The 2D circular
particles used here actually represent cylinders, all having the same length which is independent of
the particle diameter. The 2D model can be interpreted as a 3D slice model where each 2D disc
corresponds to a column of identical particles in the axial direction. All the identical particles in
each particle column are constrained to move together in the radial/circumferential plane.
Essentially the only difference between the 2D and 3D slice models for these mill applications is
Copyright :o 2001 John Wiley & Sons, Ltd. Int. J . Numer. Anal. M e t l ~ Geomech.,
. 2001; 25:83-105
102 P. CLEARY
that in 3D the particles in these columns are allowed to move independently and become
disordered. This has various effects and it depends on what information one seeks as to whether
this is important.
From a dynamical perspective, the 2D simulations are able to .reproduce the quantitative
features of 3D particles in a 2D geometry surprisingly well. Popken and Cleary [15], using a 2D
DEM and a 3D kinetic scheme, showed that the density and energy profiles across the channel in
a granular Couette flow were quite close. This demonstrated that the 2D simulations, in an
inherently 2D geometry, satisfactorily predicted the 3D energy and density distributions. Of more
importance is the effect of the disorder of the particles in 3D on the strength of the microstructure.
This is a real but modest effect. However, it is expected that neglecting the real particle shapes by
approximating them as spheres is much more important. This will be described in a future paper.
Ultimately the adequacy or otherwise of 2D simulations for these applications depends on the
information being calculated. Predictions, such as power and charge motion (which are presented
here), depend entirely on the ability of the simulation to satisfactorily predict the overall locations
of the particles. This is done well by the 2D model. Predictions that rely on the specific details of
the local microstructure and the detailed motion of individual particles, such as breakage rates,
will be much more sensitive to the particle shapes, the dimensionality of the modelling, the
inclusion of an interstitial fluid and the form of the contact force law. Similarly, the segregation
rates predicted in 3D will be quantitatively different to those predicted in 2D because of the very
strong effect that local microstructure has on percolation rates.
Similar considerations apply to the contact force laws. There are a number of alternatives
proposed in the literature capturing various aspects of the physics of particle contact to varying
degrees. These are well summarized by Schifer et al. [12]. In general, these laws model Hertzian
deformation of two spheres under various normal and oblique conditions with a variety of
frictional laws. These are then applied in very varied circumstances from dilute rapid granular
flows to static and quasi-static granular arrays. Two specific problems relating to the use of the
soft particle models for multi-body encounters are described in Reference [12]. There is little data
available on the form that the contact force should actually take in many of these regimes and
particularly in the transition between regimes. This is strongly complicated by neglecting the
particle shape effects. The linear spring-dashpot model represents essentially the lowest level of
physics that can be included to obtain reasonably realistic particle motions. This is demonstrated
by the closeness of the simulation and experimental results for the centrifugal mill presented
above. At present there are insufficient good data for large-scale flows to discriminate between the
predictions of the bulk flow produced by different contact models. It would not be surprising if,
when such data become available, that many or all of the existing contact force laws will need to
be improved in some way.
Finally, it is worth considering the effect of the choice of the spring stiffnesses. These are the
only numerical parameters in this particular DEM model. In general the stiffnesses are chosen to
be the same for all the particles. For most real materials the Young's moduli range from 10' to
lo", whereas numerically we typically use values around lo7. This means that the simulation
particles are between 2 and 4 orders of magnitude softer than the real ones. The only effects of this
softening of the particles is to increase the amount of deformation at the contbtt (given by the
overlap) and to increase the collisional time scale by 1-2 orders of magnitude. The bulk effects of
increasing the particle deformations are to potentially increase the number of contacts per
particle as the microstructure deforms and for the particle solid fraction to increase. Popken and
Cleary [15] explicitly calculated the 2d correlation function for particle contacts modelled using
Copyright 0 2001 John Wiley & Sons, Ltd. Int. J . Numer. Anal. Meth. Geomech., 2001; 2583-105
MODELLING COMMINUTION DEVICES USING DEM 103
DEM for a range of spring stiffnesses. For stiffnesses that produced average overlaps of < 1 per
cent the correlation function was found to be negligibly different to that obtained for rigid
non-deforming spheres for the full range of solid fractions. For overlaps around 2.5 per cent the
correlation function was close to that of rigid particles only for solid fractions below 0.2 and
diverged moderately for higher solid fractions as the particle deformation changed the actual
collision dynamics. For overlaps of 8 per cent and greater, the correlation function was always
smaller than that for the rigid case, with extreme divergence at higher solid fractions. This showed
that an average overlap criteria of 1 per cent was a threshold below which the particle dynamics
behaved indistinguishably from those of rigid spheres (with instantaneous collisions), but for
more deformable particles the dynamics rapidly diverged.
The 1-2 order of magnitude increase in the collision time scale is also not important in these
applications because the collision time scales still remain much smaller than the other important
time scales in the problem.-For the ball mill, typical collision time scales using our standard spring
stiffness, are of the order of 0.1 ms. The mill rotation period is of the order of 4 s whilst the time
scales for structural changes to occur to the microstructure lie between these. If the particles
become softer though, then the collision time scale increases further and begins to overlap the
other dynamical time scales, at which point the simulation dynamics begin to diverge from the
real ones.
It is important to understand these considerations since they limit the range of application of
the model. Models such as this (2D, circular particles) caa be used to calculate some quantitative
information for a good range of real two-dimensional geometries, but they should not be used
beyond these limits. The boundary of applicability depends not only on the geometry and the
nature of the particles, but on the type of information being predicted.
7. CONCLUSION
Discrete element modelling of three mill applications has been performed and a range
of quantitative predictions have been made. For the ball and SAG mills the power consump-
tions predicted are consistent with the observed power consumption of real mills with these
dimensions.
Power consumption versus rotation rate curves were predicted for several variations of the ball
mill. The sensitivity of the power curve to material properties of the charge was found, in general,
to be low. Variations were typically around 5 per cent or less for normal operating conditions.
For super-critical speeds, low friction (such as that produced by water in the mill) increases the
power draw considerably.
The effect of the mill fill level on the power curve was analysed. Equal reductions in the fill level
produced increasingly large reductions in power consumption, since the material removed was
always that which was closest to the mill centre. Conversely, the specific power consumption (a
first-order predictor of the grinding rate) was found to rise steadily with decreasing fill level. The
rotation rate at which peak power consumption occurred was similar for fill levels between 30
and 50 per cent but dropped rapidly for lower fill levels: ,Based on these DEM simulations, the
highest grinding efficiency is expected for fill levels between 20 and 30 per cent. Higher grinding
rates can be produced for lower fill levels, but the resulting mill throughput is too low to be useful.
The effect of the density of the steel balls was analysed. The increase in the weight of the charge
increased the power consumption proportionally, but the specific power consumption was
Copyright 02001 John Wiley & Sons, Ltd. Inr. J . Numer. Anal. Merh. Geomech., 2001; 25:83-105
104 P. CLEARY
unchanged. This means that the use of denser material does not effect the local dissipation of
energy and is unlikely to have large effects on the grinding rate of the rock.
Two mechanisms were identified for size segregation in mills. One operates at low speed and
concentrates large particles around the outside of the charge. The other operates at high speeds
and causes large particles to accumulate in the centre of the charge. At intermediate speeds these
mechanisms cancel. The segregation rates are rapid and almost complete radial segregation can
be achieved for a charge of rocks and balls in 95 s for N = 110 per cent.
The primary differences between a SAG mill and a ball mill are in diameter and the number of
lifters around the circumference. The greater mill diameter of the SAG mill allows it to hold
a larger amount of material and therefore increases the power consumption. In addition, the
higher position of the shoulder leads to more energetic particle collisions and leads to a 45 per
cent increase in the specific power consumption, suggesting a comparably greater grinding rate
for this larger mill. Otherwise. the charge distribution and its flow pattern are quite similar to
those of the ball mill and their variation with mill speed are also broadly similar.
The centrifugal mill predictions were compared with high-speed photographs in order to
validate the simulation predictions. At high fill levels (> 34 per cent) the simulated charge
demonstrated a steady profile that rotated with the mill chamber, in the same way as that
demonstrated by the experiments. The actual charge profile is indistinguishablefrom the experi-
ment for a fill level of 75 per cent and is very close for 50 per cent. For fill levels below 34 per cent,
the simulations displayed an equivalently observed flow bifurcatiori leading to a chaotic flow with
a three-pointed structure that flopped around the inside of the mill. The predicted motion was
qualitatively very similar to that shown in the experiment.
The closeness of the predictions of the charge profile and behaviour to that of experiments for
the centrifugal mill and the closeness of the power predictions for the ball and SAG mills to those
observed for full size mills suggest that DEM models are able to capture much of the essential
dynamics of these important industrial particle flows.
REFERENCES
1. Campbell CS. Rapid granular flows. Annual Review of Fluid Mechanics. 1990; 22:57-92.
2. Barker GC. Computer simulations of granular materials. In Granular Matter: An Interdisciplinary Approuch, Mehta
A (ed.). Springer: New York, 1994.
3. Walton OR. Numerical simulation of inelastic frictional particle-particle interaction. In Particulate Two-Phase Flow,
chapter 25, Roco MC (ed.). 1994; 884-911.
4. Mishra BK, Rajamani RJ. The discrete element method for the simulation of ball mills. Applied Mathematical
Modelling 1992; 16598-604.
5. Mishra BK, Rajamani RK. Simulation of charge motion in ball mills. Part 1: experimental verifications. International
Journal of Mineral Processing 1994; 40:171-1 86.
6. Cleary PW. Modelling industrial granular flows. In Proceedings of 1st Australian Engineering Mathematics Confer-
ence, Melbourne, 11-13th July, Eastson AK, Steiner J M (eds). Studentlitteratur, 1994; 169-177.
7. Rajamani RK, Mishra BK. Dynamics of ball and rock charge in SAG mills. Proceedings of SAG 1996. Department of
Mining and Mineral Process Engineering, University of British Columbia, 1996.
8. Cleary PW. Discrete element modelling industrial granular flows. Proceedings oF 2nd Australian Engineering
Mathematics Conference, Yuen WY, Broadbridge P, Steiner J M (eds). Institute of Engineers, Australia. 1996;
301-308.
9. Cleary PW. Predictingcharge motion, power draw, segregation, wear and particle breakage in ball mills using discrete
element methods. Mineruls Engineering, 1998; 11:1061-1080.
10. Cleary PW. The filling of dragline buckets. Murhemarical Engineering in Industry 1998; 7:l-24.
11. Robinson GK, Cleary PW. The conditions for sampling of particulate materials to be unbiased-investigation using
granular flow modelling. Minerals Engineering. 1999; 12: 1 101-1 118.
Copyright 6 2001 John Wiley & Sons. Ltd. Int. J . Numer. Anul. Meth. Geomech., 2001; 25:83-105
MODELLING COMMlNUnON DEVICES USING DEM 105
12. Schier J, Dippel S, Wolf DE. Force schemes in simulation of granular material. Journal de Physique Institute of
France 1996; 65.
13. Cleary PW. Discrete element modelling of industrial granular flow applications. TASK. Quarter1.v-Scientific Bulletin
1998; 2: 385-416.
14. Hoyer DI. Particle trajectories and charge shapes in centrifugal mills. Proceedings ofthe International Conference on
Recent Advances in Mineral Sciences and Technology, MINTEK, South Africa, 1984; 401-409.
15. Popken L, Cleary PW. Comparison of kinetic theory and-discrete element schemes for modelling granular Couette
flows. Journal of Computationnl Physics 1999; 1551-25.
Copyright 02001 John Wiley & Sons, Ltd. Int. J . Numer. Anal. Meth. Geomech., 2001; 25:83-105