Recent advances in dem modelling of tumbling mills

0892--6875(01)00145-5

Minerais Engineering, Vol. 14, No. 10, pp. 1295-1319,2001 © 2001 Published by Elsevier Science Ltd Ali rights reserved 0892-6875/01/$ - see front matter

RECENT ADVAN CES IN DEM MODELLING OF TUMBLING MILLS·

P.W.CLEARY CSIRO Division of Mathematical and Information Sciences, Private Bag 10, Clay ton, Vic, 3169. Australia E-mail: [email protected] (Received 8 May 2001; accepted 19 fuly 2001)

ABSTRACT

ln This paper we focus on recent developments in DEM modelling oftumbling mil/s. Examples of 3D models of SAG mills are presented along with detailed predictions of power draw, liner wear rates, liner stresses and energy spectra. Comparisons of simulation results with photographs of charge motion in a scale model SAG mill demonstrate a high degree of accuracy. The inclusion of particle breakage in the DEM is demonstrated as is the methodology for predicting breakage rates, mill throughput, and equilibrium particle size distribution of the charge. Final/y, the variation of power draw with liner shape throughout the lifter life cycle is illustrated. @ 2001 Published by Elsevier Science Ltd. Al/ rights reserved. Keywords Comminution; modelling; simulation; mineraI processing

INTRODUCTION Particle breakage is an essential component of mineraI processing. The aim is to reduce the particle sizes from those supplied to the much smaller ones required for further processing as efficiently as possible, producing maximum mill throughput and minimising operating costs. Such milling processes are believed to use less than 5% of the supplied energy for particle breakage. Opportunities exist to improve the performance and control of such mills using information obtained by computational modelling using the Discrete Element Method (DEM). Simulation of tumbling and other types of mills by DEM offers the opportunity of better understanding the internaI mill dynarnics and for developing improvements to mill design and operation that can lead to large increases in mill efficiency and throughput. DEM has been used for modelling a wide range of industrial applications, particularly in milling. Early work on baIl mills by Mishra and Rajamani (1992, 1994) has been followed by Cleary (1998a, 2oo1a) and Cleary and Sawley (1999). Similarly SAG mills were modelled by Rajamani and Mishra (1996) and subsequently by Cleary (2001b). Experimental validation has been performed for a centrifugaI mill (Cleary and Hoyer 2000) and charge motion predicted in a HICOM mill by Cleary and Sawley (1999). Other applications include vibrational segregation by Haff and Werner (1986), Ristow (1994) and Cleary (1998c), flows in rotatingdrums by Ristow (1994), Pôschel and Buchholtz (1995) and others, sampling by Robinson and Cleary (1999), filling of dragline excavators by Cleary (1998b, 2000), hopper flows by Hoist et .al. (1999), Potapov and Campbell (1996), Ristow (1994) and many more). • Presented at Comminution '01, Brisbane, Australia, March 2001

1295

P. W. Cleary

1296

THE DISCRETE ELEMENT METHOD DEM simulation involves following the motion of every particle in the flow and modelling each collision between the particles and between the particles and their environment (such as the mill liner). Industrial applications place heavy demands on the geometrical capabilities of DEM codes, particularly for threedimensional geometries. In our DEM code, any two-dimensional cross section can be constructed using an appropriate number of line segments, circular segments or discs. This gives essentially unrestricted geometric capability. In three-dimensions, boundary objects can be constructed by either extruding two dimensional cross sections or by importing triangular finite-element surface. meshes. Such meshes can be produced using any reasonable mesh generator from solid models generated in suitable CAD packages. This pro vides enormous flexibility in specifying three-dimensional environments with which the particles interact. The particles are currently modelled as spheres in three-dimensions and by discs or super-quadrics in twodimensions. Super-quadrics are defined by the equation

(1) where the power N determines the blockiness or angularity of the resulting particle and A determines the aspect ratio. Aspect ratios of up to 12: 1 and angularity factors of up to 20 can be easily used. Arbitrary size, shape and density distributions can be specified for the particles. They can either be added during the simulation as a stream of particles or can be placed in packed microstructures in user specified locations. The general DEM methodology and its variants are well established and are described in review articles by Campbell (1990), Barker (1994) and Walton (1994). Here we use a reasonably standard implementation that is described in more detail in Cleary (1998a, 1998c). Briefly, the particles are allowed to overlap and the amount of overlap dx, and normal Vn and tangential VI relative velocities determine the collisional forces via a contact force law. There are a number of possible contact force models available in the literature that approximate collision dynamics to various extents. We use a linear spring-dashpot model which is shown diagrammatically in Figure 1. For more complex models see Walton (1994) and Schtifer et. al. (1996). "Soft Particlc" contact model I1x is the partide overlap, kil and kt are the normal and tangential spring constants, vn and vt are the normal and tangential velocities, Cn and C, are the normal and tangential damping coefficients. iJ is the friction coefficient

Fig.1

Contact force law consisting of springs, dashpots and frictional sliders used for evaluating forces between interacting particles and boundaries.

The normal force (2)

consists of a linear 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 kn of the spring in the normal direction. Typically, average overlaps of 0.1-0.5% are desirable, requiring spring constants of the order of 107 N/m in two dimensions and 106 N/m in three dimensions. The spring constants

Recent advances in DEM modelling of tumbling mills

1297

required to achieve these overlap targets are lower in three dimensions than in two dimensions, since the force exerted on particles is proportional to their masses, and the mass of a spherical particle is much smaller than that of a cylinder of unit length. The normal damping coefficient en is chosen to give the required coefficient of restitution E (defined as the ratio of the post-collisional to pre-collisional normal component of the relative velocity), and is given in Cleary (1998c). The tangential force is given by

F;

=

min ~Fn ,kr fvrdt + Crvr } ,

(3)

where the vector force Fr and velocity Vr are defined in the plane tangent to the surface at the contact point. The integral term represents an incremental spring that stores energy from the relative tangential motion and models the elastic tangential deformation of the contacting surfaces, while the dashpot dissipates energy from the tangential motion and models the tangential plastic deformation of the contact. Depending on the history of the contact, it is possible for the spring to be loading in one direction and simultaneously unloading in the orthogonal direction. The total tangential force Fr is limited by the Coulomb frictional limit J.iFn , , at which point the surface contact shears and the particles begin to slide over each other. The discrete element algorithm has three main stages: 1.

2.

3.

A search grid is used to periodically build a near-neighbour interaction list that contains aIl the particle pairs that are likely to experience a collision in the short term. Using only particle pairs in this list reduces the force calculation to an O(M) operation, where M is the total number of particles. Industrial simulations with 500,000 particles are now possible in reasonable times on CUITent single processor workstations. The forces on each pair of colliding particles and/or boundary objects are evaluated in their local reference frame using the spring-dashpot interaction model (see Figure 1), and then transformed into the simulation frame of reference. Al! the forces and torques on the particles and objects are summed and the resulting equations of motion are integrated. Time integration is performed using a second-order predictor-coITector scheme and typically uses between 15 and 25 time steps to integrate accurately each collision. This leads to very small time steps (typically 10- 4 to 10- 6 s depending on the control!ing length and time scales for each application).

a=

Typical values of 0.3 and i = 0.75 are chosen for the coefficients of elasticity and friction respectively. It was shown in Cleary (1998a) that the power draw and charge behaviour was relatively insensitive to the choices of these two physical quantities, so the same values are used throughout al! the simulations in this paper. Quantitative predictions of boundary stresses, wear rates and distributions, collision forces, energy speCtra, power consumption, torques and flow rates, sampling statistics, mixing and segregation 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 Cleary (1998c). This DEM code has been used extensively for studying mill behaviour (see Cleary 1998a, 1999, 2000a, 2001a, 2oo1b and Cleary and Hoyer 2000).

DEM MODELS OF SAG MILLS The main limitation on DEM models is the number of particles that can be simulated in reasonable computational times. The greatest numbers of particles in the mill are in the smallest size ranges, which are not expected to affect the charge dynamics much (since their very smal! masses mean they have only marginal effects on large particles). The first approximation that can be made is therefore to truncate the size distribution at a specific cutoff size. Cleary (2001b) has showed that for a 5 m ball mil! truncation of the size distribution at between 1 and 1 mm has !ittle effect on mill power or charge motion. This suggests that a 20 mm cutoff for a 10 m SAG mill is reasonable. There are still too many particles above this size in a SAG mill to be modelled on CUITent computers, so a further approximation is required. The two alternatives used here are:

°

1298

P. W. Cleary

1.

Approximate the type of particle size distribution (PSD) by reducing the numbers of partic1es in the smaller fractions. This then allows the full mill to be modelled inc1uding the ends.

2.

Use the correct partic1e size distribution, but use only a periodic radial slice from the middle of the mill. Only part of the mill can be modelled in this way and axial transport cannot be examined.

Full SAG mill with approximate PSD Our first sag mill example is of a 10 m diameter and 4 m long SAG mill with flat end walls. It has 72 lifters with face angles of 5 degrees. The fill level is 30% (by volume) and the mill rotates at 75% of critical speed. The rocks have diameters uniformly distributed between a topsize of 250 mm and the 20 mm lower cutoff. The balls have diameters uniformly distributed between 50 and 200 mm. The full size mill simulation contains about 281,000 particles Figure 2 shows two snapshots of this SAG mill. In the upper frame the particles are coloured by diameter and in the lower frame by velocity. The particles near the mill shell rotate with the liner rotation and at moderate speed (green) from the toe position around and up to the shoulder position. The partic1es between the lifters are thrown in a high dilute cataracting stream which impacts on the liner around the 3 o'c1ock position. The impact region of the cataracting stream is well above the toe, even though the mill speed is only 75% of critical. This occurs because the lifters are close packed and have very steep face angles. The majority of the charge avalanches down the steep free surface to a point halfway across the mill and then flows more gradually across the remaining distance to the toe down a gently sloped surface. The upper plot shows a layer of fines (blue particles) at the top of the cataracting stream. These fines have segregated to the outside of the charge and predominantly lie between the lifters. The presence of the end walls on the mill give frictional support to the charge and help raise the shoulder position higher. The velocity colouring shows a large band of slowly moving and stationary particles in the middle of the charge. The particles near the lifters (green) move upwards at the peripheral mill velocity while the partic1es in the cascading zone flow down at a similar speed. This 1eads to high shear rates between the upward stream and the stationary zone and between the stationary zone and the cascading region, which is expected to lead to significant fines production via abrasion, chipping and attrition. The partic1es in the cataracting stream can be seen to change colour as they accelerate down, becoming red as they approach the toe of the charge. Slice model of a SAG mill with realistic PSD As mentioned above, another choice for the DEM model of the SAG mill is to use the correct PSD and to simulate only a narrow slice of the mill. Figure 4 shows the charge distribution in a 0.5 m section of a common design 36' SAG mill. The mill attributes are summarised in Table 1. There are again 72 rows of symmetric close packed steep face angle lifters (7 degree). The baIl and rock size distributions are given in Table 2. These are much more realistic than the ones on the previous example with the bulk of the rocks less than 50 mm and the bulk of the baIl bigger than 63 mm. This DEM model contains 184,000 particles with a minimum particle size of 25 mm. Using a standard 6 spring stiffness of 10 gives a timestep of 3.4xl0- 5 sand requires about 31 cpu hours 1 second of simulation on a 500 MHz EV6 alpha workstation. Charge motion

The overall flow pattern is quite similar to the previous SAG mill mode!, but the charge distribution is much clearer. This results from the larger proportion of fines in the mill giving more definition to the surfaces and the absence of the end walls leading to less disorder in the flow from end wall collisions. The charge surface has a very clear bi-linear shape, rising slowly from the toe to the middle of the charge and then steeply up to the shoulder. A substantial cataracting stream departs from the lifters above the shoulder

Recent advances in DEM modelling of tumbling mills

a)

b)

Fig. 2

Full size DEM model of a SAG mill a) particles coloured by radius (blue are fines and red are coarse) and b) particles coloured by velocity (blue are low speed are red are high speeds).

1299

P.

1300

w. C1eary

a)

b)

Fig.3

Slice model of a 36' SAG mill with a) particles coloured by size (blue are fines and red are coarse) and b) particles coloured by velocity (blue are stationary or slowly moving and red are faster particles).

Recent advances in DEM modelling of tumbling mills

1301

and falls in a well defined stream to impact on the liner above the toe at around the 3 o'clock position. The effects of radial segregation are very clear with a concentration of fine particles against the mill shell and in the upper part of the cataracting stream. The bulk of the lower part of the cataracting stream consists of larger rocks and balls. Segregation is use fuI in this context by ensuring that what balls are in the cataracting stream are at least on the shallower trajectories and lead to less damaging impacts with the liner. The colouring of the particles by velocity again reveals the slow moving nature of the middle of the charge and the high shear rates between this and the mill shell and the cascading flow above. The acceleration of the cataracting particles to more than 13 mis is also very clear.

TABLE 1 SAG mill parameters Mill diameter

10.86 m

Mill slice

0.5 m

Lifter face

7 degree

Lifter height

200 mm

Lifter width

150 mm

# Lifters

72

Fillievei

30% (by volume)

Ball charge

10% (by volume)

Mill speed

10 rpm (78% critical)

TABLE 2 Bail and rock size distributions for the slice SAG mill model Bali size distribution:

Rocks size distribution: (sg=2.8)

-125 + 88 mm: 30%

-200 + 140 mm: 8%

-88 + 63 mm: 40%

- 140 + 100 mm : 12%

-63 + 44 mm: 20%

-100 + 70 mm : 15%

-44+31mm:1O%

-70 + 50 mm: 15% - 50 + 35 mm : 20% - 35 + 25 mm : 30%

P. W. C1eary

1302

Powerdraw An important prediction of DEM models is the power draw resulting from the motion of the particulates in the mill. Note that this does not include power draw relating to the charge slurry or those due to mechanical los ses in the drive or bearings. Figure 4 shows the instantaneous power draw predicted by two different methods. The fluctuating black curve shows the power supplied by the mill to the particle charge. After sorne rapid fluctuations upon the startup of the mill motion this rapidly approaches it equilibrium value with a moderate peak around 3 s relating to lifting the bulk of the charge that has reached the toe after the first major slope collapse. After this the power draw is essentially steady. The smooth dark grey curve is a smoothed version of this power draw. The high frequency fluctuations produced by specific events in the mill, such as variations in the avalanching free surface, are less than 10% of the mean. There is also sorne evidence of a gentle longer term periodicity of around 4 s (approximately the mill rotation period) that may be an echo of the first major collapse and the peak at 3 s.

0

0

0 N

0

0

--.. Il) ...... ~

-

0

3:

......

.:::t. ~

CI)

0

a.

0

0

0 0

Il)

o

5

10

TIme (s) Fig. 4 Instantaneous and average power draw for a 36' SAG mill calculated in two ways. The power draw can also be calculated by integrating aIl the energy dissipation in both the normal and shear directions in aIl the collisions. This is shown by the fluctuating light grey curve. This shows very little energy dissipation between 1 and 2 s as the initially packed charge rotates rigidly up to the dynamic angle of repose. This leads to a major si ope failure and a massive avalanche producing a very sharp peak in the power dissipation at 3 s. As the charge seules in the toe, this power draw measure falls to just below the equilibrium value. The smooth light grey curve shows a smoothed version of this with the high frequency fluctuations removed. Between 5 sand 8 s there is a small difference between the two smoothed power draw measures. Such differences represent potential energy storage in the charge and indicates that the shoulder is moving slowly higher during this period. After 8 s the two power draw curves converge and are constant. This gives the prediction of the average mill power draw for these specific conditions. This figure demonstrates that this mill has not reached a suitable equilibrium until after 8 s of operation. This is around 2 mill revolutions and shows that the mill charge responds rapidly to changes in conditions. The equilibrium average power draw predicted for this 36' SAG mill is 778 kW / 0.5 m mill slice or 19.87 kW / tonne of charge. For a typical 6 m long SAG mill the total power draw produced by the charge is predicted to be 9,320 kW. This is very close to the power draw of around 10,000 kW that one might typically expect for such a SAG mill under these conditions.

1303

Recent advances in DEM modelling of tumbling mills

SAG mille eollisional energy speetra Prediction of collisional energy spectra allows the various contributions to the overall energy dissipation within the mill to be understood. Knowledge of the changes in the se spectra, with operational and geometric mill attributes, offers the opportunity to improve the energy efficiency of breakage and to reduce the energy wastefully dissipated such as in balI-liner collisions. Figure 5 shows the energy spectra for the 36' SAG mill being studied. In the DEM simulation the energy dissipation associated with each separate particle encounter can be calculated. The frequency distributions of these energy losses, both in the normal and shear directions, can be determined from the individual collision events. These spectra are then shown as collision frequency versus the energy level of those collisions. Figure 5a shows the spectra (both shear and normal) for aIl the collisions combined. The largest proportion of collisions (around 105 per second) are in the lowest energies of 0.01 J and below. For the normal energy dissipation, the frequency declines linearly (on this log-log sc ale) until there are only a few collisions above 100 J with a maximum energy of 300 J. The shear dissipation curve is noticeably curved with a larger proportion of collisions in the 0.1 to 5 J range. High shear energy events are less likely than high normal energy ones, but the highest energy dissipation of around 500 J occurs for the shear mechanism.

a)

b) Energy spectro: Ali collisions

Energy spectro: Rock-Rock

'b

'b Normal energy

~

0

--

.g §

...t"

0-

:i'5" u

'';re,..,r

toe.

e.(1f"~gy 10SS

>-

.

g

"

Normal energy 10••

"0

~:.,-

0 0 0

0-

8

~ 8 c

~

0

~

"0 u

-

~

ci

ci 0.01

0.1

1

10

100

0.01

0.1

Energy (J)

'b Normal energy 10..

(,

-

Sheor energy los9

Normal energy 108.

(,

."il-

~

·2" â u

~

~" ~

-

u

~

. c

"

-

â

-

ci 0.1

1

10

100

0.01

0.1

Energy (J)

li! 0

u

'b Normal energy 108.

è

-

Sheo( energy loss

§

0

.

0-

c

~

0

~

-

Normal energy 108.

~

~ c

...!"

ê

"0 u

ci

-

Sheor energy los9

§ 8 ~

ci

0.01

0.1

1 Energy (J)

Fig. 5

100

Energy spectro: Boll-Uner

'b

~"

10

f) Energy spectre: Boil-Boil

"

1 Energy (J)

e)

0-

IOS9

~ 8

0.01

~

Sheor enerqy

§

0-

8

ci

."

100

Energ)' spectro: Rock-Uner

'b

~

10

d) Energy .pectre: Rock-Boil

...!

1 Energy (J)

c)

~

Sheor er.e~JY loss

10

100

0.01

0.1

1

10

100

Energy (J)

Energy spectra for the slice model of the SAG mil! for, a) aIl collisions, b) rock-rock collisions, c) rock-bail collisions, d) rock-liner collisions, e) bali-bail collisions and f) bali-liner collisions.

1304

P. W. C1eary

Figure 5b shows the spectra for just the rock-rock collisions. The shear spectra is essentially the same as for the overall one, but the normal spectra is intriguingly different. There is little difference up to about 3 J, but after this there are noticeably lower collision rates with the higher energies. In particular, the normal spectra is now lower than the shear one. The spectra for the rock-ball collisions are both similar to the overall spectra but with reduced frequencies around 2x104 at low energies instead on 105 • This reflects the lower proportion of balls in the charge compared to rock. The highest energy dissipation levels are also lower, with few collisions dissipating more that 200 J. The rock-liner spectra are quite different though. The high energy part of the spectra is similar to previous ones, but there is a pronounced and non-smooth bend in the spectra at around 2 J with significantly lower collision rates for the lower energies. This reflects the fact that most of the rock-liner collisions occur for particles in the cataracting stream and these all tend to collide with the liner with high energy. The small numbers of low energy collisions result from particles in the cascading stream striking the liner in the toe region and from secondary impacts from slowed particles from the cataracting stream. The baIl-baIl spectra are shown in Figure 5e and are again quite different. These curves have pronounced curvature with the more shear dissipation events in the mid-range energies and most of the high energy events being impact (normal) ones. The collision rates for low energies at 2xl03 is still quite high and demonstrates that significant energy is dissipated in baIl-baIl collisions that at best waste energy and at worst lead to accelerated wear rates for the balls, and at the higher impact energies, sorne baIl breakage. Figure 5f shows the balI-liner spectra. This is most instructive. The spectra are quite flat with re1atively small numbers oflow energy collisions (around 20 collisions/second at the 0.01 J level) but with significant rates still at the highest energy levels. In particular, high energy impacts (normal dissipation events) with energies between 100 and 350 J occur for balI-liner with about the same frequency as experienced for the overall charge. This explains the absence of most of the high energy impacts from the rock-rock spectra when compared to the overall spectra that was observed above. Most of the high energy tail of overall spectra is caused by bali-liner collisions. The occurrence of such events can be confirmed by observing the trajectories of the cataracting stream in the DEM animations and observing that many of these particles are in fact balls. This demonstrates that for this mill much of the high impact energies introduced by the use of a reasonably high charge is wasted on producing liner damage rather than doing use fui breakage. Liner stress and wear distributions Stress and wear data is collected on a high resolution triangular mesh that covers the millliner. The average element length is 20 mm, which is smaller than the size of the smallest particles in the simulation, giving very good spatial information about the stress and wear distributions. The contribution of every particle collision with the mill liner is accumulated in the triangular elements in which the collisions occur. This raw data is then aggregated across the depth of the mill slice and spatially smoothed to remove noise on the scale of the particles. Figure 6 shows a section of the liner coloured by the stress components and sorne of the wear measures. Figure 7 shows the stress distributions along the liner as a line plot. The distance is measured from the bottom of the front face of the lifter. Each section of the lifter/liner stress is separated by a vertical dashed line. The first section shows the stress along the front face, then along the lifter top, the lifter back and the last section shows the stresses on the liner plate. The normal and shear stresses have very similar distributions, but the normal stress is around three times higher than the shear stress. The stress is zero at the base of the lifter (since no particles can make contact with this part of the lifter) and rises steadily with distance up the front of the lifter to a peak occurring near the top corner. This reflects the force transmitted to the front of the lifter as they push into and then lift the charge. The highest stresses for this new liner actually occur exactly at the corners, but these are not shown in Figure 7 since the magnitudes are much higher and any amount of rounding of the corner by wear will substantially reduce these peaks. The stress along the top of the lifter is relatively constant and about half the level of at the top of the face. The stress on the back of the lifter is low at about half the level of the top of the lifter. This reflects the fact that only a small proportion of the charge is supported by the backs of the lifters. The stress on the liner plate is higher and increases as the front of the next lifter is approached. Note that there are modest size but weIl defined peaks near the front and back of each lifter. These occur 18 mm from the corners and correspond closely to the average particle radius. This coherent peak occurs because a

Recent advances in DEM modelling of tumbling mills

1305

particle al ways sits in each of the corners while the charge is in contact and the average location of this contact is deterrnined by the average particle size in the mill. Our ability to capture such features demonstrates the spatial accuracy that is now possible to obtain using DEM simulation. If one were to use a charge consisting of a number of discrete sizes (as opposed to the continuo us sizes used within each sieve size) then one could also resolve up to three more coherent peaks along the liner and up the front face with magnitudes decreasing with distance from the corners.

Fig.6

Stress and wear distributions on the SAG liner, a) normal stress, b) shear stress, c) impact damage and d) abrasion damage. In each case red represents high magnitudes, green is moderate and blue is low.

Normal stress (newton/metreI\2) ; ; ;

!

;

cij

i

!

! !

'~!

11!

, !

i i

!

i

!

i -

o

200

ln

-0

400 600 Distance a10ng surface (mm)

800

newton/metrel\2 i1 1

1 1

!

1

! 1

1

~i '

: ! ------1-

o Fig. 7

200

!

i 1

Irl

-----0

400 600 Distance along surface (mm)

800

Normal and shear stresses along the lifter and liner plate. The distance is measured from the base of the lifter front. Each region (lifter front, top, rear and lifter plate) is separated by a verticalline.

P. W. C1eary

1306

To predict the impact damage on the liner we use two measures. The first is the energy dissipated in the normal direction in collisions between the particles and the liner. The second is a measure of excess kinetic energy of impact. Low speed collisions ( < 0.1 mis) which are large in number but of limited importance for damage make no contribution to the damage estimate, but high speed collisions do much more damage because of the quadratic dependence on speed. These impact damage measures are shown in Figure 8. The distributions are quite different to the normal stress distributions. In particular, the wear is substantially higher across the entire top of the lifter with peaks near the corners. There is little impact damage on either front or back faces since they are protected from the cataracting stream by the steep face angles and close lifter spacing. There is sorne moderately higher wear on the upper part of the leading face that is produced when the lifters smash into the charge in the toe region. The wear on the liner plate though is meaningfully higher with a maximum in the middle. This damage is produced by the penetration of the cataracting stream (and particularly balls) between the lifters where they pound the liner. The middle of the liner plate is clearly the most exposed to impact and so has the highest wear rate predicted.

Normal work 1500

!l

1 i V " :

1000

1! 1 1

! ! !

500

1

!

i!

i1

i

o~~

o

i

"'\ __~~~-~'1__-~~~~__~~~~________~~b

200

400

600

800

Distance along surface (mm)

e 400

1

1 \ ) \ ;

300

~

i

i

1

1

1

1

I

i

I i1 i ! i !

o-

---- ----1----

o

Fig. 8

200

-l---~.........--"-,/

400 600 Distance along surface (mm)

Impact wear rate along the lifter and liner plate as given by the normal work and an impact damage measure.

The abrasive wear is also estimated using two measures. The first is the shear work which is the energy dissipated by sliding (tangential) interactions between particles and the liner. The second uses the kinetic energy of each collision with the inclusion of a strong angular dependence. This takes account of the fact that collisions at around 22° produce significantly more scouring/abrasion damage than do particles sliding directly along the boundary and normal impacts. Figure 9 shows the abrasive wear distributions. These are again both quite different to the stress and impact damage distributions shown above. The highest abrasive wear occurs on the front face of the lifter with the wear rate increasing with height up the lifter. This wear will lead to steadily increasing lifter face angle (as one would expect). There is also significant abrasive erosion from the top surface of the lifter which one would expect to lead to steadily decreasing lifter height. The abrasion on the back of the lifter and the liner plate predicted by the shear work is relatively even and around one third of the magnitude on the lifter top. The second abrasion measure gives similar predictions to the shear work for the front and top faces but predicts much lower wear on the liner and the rear face. This lower prediction arises from the much lower weighting given to the many low speed contacts sliding directly along these surfaces. They dissipate a reasonable amount of energy, but this is likely to be an overesti mate , as shown by the second measure. The second measure suggests that the wear will decrease with height down the back face which is reasonable since most oblique impacts will occur during the filling of the space between lifters as they pass through the toe region. There is a peak in the peak abrasion in the middle of the liner for sirnilar reasons.

1307

Recent advances in DEM modelIing of tumbling mills

The rate of normal work is around double the rate of shear work. This might cause one to conclude that the dominant erosion mechanism is from impact rather than abrasion. The actual erosion rates, however, must also be dependent on the material properties of the liner and its resistance to impact and abrasion damage. A high quality steel might be expected to be very resistant to impact and erode predominantly by abrasion. Our DEM code currently has the capability to evolve the shape of the liner in accordance with the wear rates predicted. However, it is not clear which of these wear rates or which combinations one should use in order to obtain the best quantitative predictions of the wear behaviour. Ultimately, it will be a combination of an impact and an abrasion measure weighted by the resistance of the liner material to each damage mechanism.

Shear work ·oules)

i

~ !

i ~!

! ! !

i

i! ! ,

,

,

!

!

i

ii

li

!~o

i

i

i

o

1!

j

200

400

600

800

Distance along surface (mm)

Abrasion

140r----------,------,---------,-------------,

~:1

120 100

1

1

i

i

80

j i ! !

60

!! i

i

! , !

j

-1

o

200

! ! ,

-

i

1·

400

600

800

Distance along surface (mm)

Fig. 9

Shear wear rate along the lifter and liner plate as given by the shear work and an angular dependent abrasion measure.

V ALIDA TION OF SAG CHARGE MOTION: COMPARISON WITH EXPERIMENT Comparison of predicted charge distributions with those from experiments performed by JKMRC (Ward, 1992) allows evaluation of the accuracy of the various DEM models for the charge motion in a sc ale model SAG mill. The mill used in the JKMRC experiments has an internai diameter of 592 mm and was fitted with Noranda lifters. Figure 10 shows a snapshot from a two-dimensional DEM simulation of this mill showing the charge and the lifters. The experimental charge consisted of basait with sg = 2.84. The full particle size distribution from the experiments from 2 mm to 22.4 mm was replicated in the simulations to give as precise a match as was possible. During the experiments photographs were taken using a 1/25 s time exposure. Each particle then appears as a streak. To facilitate comparison a similar process was used for the DEM with particle trajectories drawn for 1/25 sand with thicknesses matching the particle size. Figure Il shows the experiment and matching full 3D DEM prediction for a fill level of 35% and a speed of 75% critical. The charge distribution and flow patterns are very close, with good agreement on the shoulder, vortex center and free surface shapes and reasonable agreement around the toe region. The DEM shoulder angles lie in the range 45-50° (measured anti-clockwise from the 3 pm location) and are close to the experimentally observed values of 50-56°. Similarly the toe angles from the full 3D DEM simulation of 226-231 ° covers the experimental range of 232-238°. It should also be appreciated that the shoulder and toe angles do not provide a precise description of the charge position. Repeated measurements of the same photographs by the original experimenter (Ward, 1992) gave standard deviations of 1.5 - 2°.

1308

P. W. Cleary

Fig.l0 Conventional picture of DEM simulation of validation case. Tables 3 and 4 show the variation of shoulder and toe angles with the type of DEM model used. The model types considered here are two-dimensional with circular particles, two-dimensional with non-circular particles, a three dimensional slice model and a full three dimensional model including the end wall effects. AlI four produce qualitatively similar flow patterns, but with important differences that can be summarised: • The 2D circular model consistently under-predicts the shoulder and toe locations by around 10°. This occurs because the circular particle microstructure is weaker than it should be and fails more easily producing a more slumped charge. • The inclusion of particle shape increases the shear strength of the particle microstructure leading to higher shoulder locations (by around 3 -5 degrees). • The 3D slice model includes the locking effects of particles in 3D which also increases the strength of the microstructure. This has a much larger effect on the predicted location of the toe (increasing it by 6-8 degrees for lower mill speeds) • The full 3D model also includes the very strong frictional and supporting end wall effects. The motion of the walls helps lift the charge higher. The predictions of the full 3D model are quantitatively very close to those of the experiments. It is expected that the combination of particle shape effects with the full 3D model will allow extremely close matching of the data. A more detailed study and comparison is in progress with other fill levels and a range of lifter profiles being considered.

TABLE 3 Shoulder angles for the experiments and different DEM models for a 35% tillievei

Speed

Experiment

2D-C

2D-NC

3D slice

Full 3D

65%

46.4-50.4

37.1-42.6

41.3-45.4

35.1-41.7

41.0-44.0

75%

50.3-55.5

40.8-44.3

45.4-49.6

39.7-45.5

44.5-49.3

95%

68.6-72.4

58.6-70.0

59.2-70.5

65.9-72.6

68.0-77.6

Recent advances in DEM modelling of tumbling mills

1309

TABLE 4 Toe angles for the experiments and different DEM models for a 35% fillievei

Speed

Exp

2D-C

2D-NC

3D slice

Full 3D

65%

228.8-232.2

217.8-219.9

219.4-222.0

225.4-228.4

228.7-231.0

75%

232.2-237.6

218.8-220.9

219.2-223.1

224.7-228.9

226.7-231.4

95%

206.1-223.5

218.4-223.5

219.0-227.4

207.5-227.4

198.7-232.3

Fig. Il Comparison of a full 3D DEM model with experiment for a scale model SAG mill for a fill level of 35% (by volume) and a speed of 75% critical a) experiment and b) matching DEM prediction.

DEPENDENCES OF POWER DRA WIN BALL MILLS As described earlier, predictions of power draw in mills can be made using DEM models. This allows the variation of power draw with various mill, charge and operational parameters to be studied. Figure 12a shows the predicted variation power draw for fillieveis between 10% and 50% for a 5 m diameter ball mil!. Significant changes in power draw are observed. In particular the location of the peak power moves sharply to lower speeds as the fillievei decreases. Figure 12b shows the change in power draw throughout the lifter wear cycle for four mill speeds. At low mill speeds the power draw is highest for new lifters and declines gradually as the lifters wear. At 80% critical speed, the power draw actually increases until the peak power is obtained when the face angle has worn to 60 degrees. This trend continues, with the power draw rising continuously as the lifters wear when the mill speed is 100% of critical. A comprehensive analysis of the dependences ofball mill power are given in Cleary (1998a, 2001b).

BREAKAGE IN DEM MILLS Breakage is obviously the most important prediction that we would want to be able to make using DEM. Here we de scribe our first generation of breakage models that have been included in our DEM mode!. Continuously feed mills can now be mode lIed with the fill rate set to match the rate of fines production and maintain a constant fill leve!. The equilibrium particle size distribution of the charge can be predicted as can the mill throughput. The breakage models take account of both high energy impact breakage and compressive breakage under the toe of the charge.

P. W. Cleary

1310

"c

~

c

,

Il)

0

.....

~

'-' III

III III

e E ~

III

e E

0

.,

0>

e .s=

,

U

U

L-

~

a..

'0....

60

80

~ 100

120

80

t.4i11 Rotation (% criticol) a)

-+

~

• - lQ1: fil

Il)

N=70"

~

x-2œfll 11- 1~ fil

e>

0

N

I? o .s=

L-

60

20

Ufter face ongle (degrees) h)

Fig.I2

Predictions of power draw for a 5 m ball mill, a) variation with fill level and b) throughout the lifter wear cycle.

Stress chains To begin, it is worthwhile demonstrating the importance of stress chains in granular materials. Figure 13 shows a stationary box that is partially filled with particles having a modest size distribution. The particles are shaded according to their total instantaneous force, with dark grey being high and light grey showing low force. To make the stress chain even clearer we connect the centers with lines whose widths depend on the magnitude of the inter-particle force. This is a traditional DEM way of showing stress chains. This example highlights several important facts: 1.

2. 3. 4. 5.

Unlike in fluids, the forces are far from isotropic and do not vary purely with depth The force is concentrated in stress chains Most of the particles experience only low stress with most of the weight carried by the small number of particles in the main stress chains Large particles tend to be in the stress chains because their larger areas provide more opportunity for contacts to transmit force to them. Stress chains foc us compressive force on larger particles deep in the granular material.

This last point is particularly important for tumbling mills, since compressive failure of large particles under the toe region is potentially an important mechanism for particle breakage. Figure 14 shows a 5 m baIl mill with particles similarly shaded by the total compressive force applied to them. A collection of major stress chains can again be seen to support most of the dynamic weight of the charge. These stress chains commence at the liner and extend up into the charge. They tend to have a curved shape and change on a rapid timescale (millisecond) over which there are only slight deformations of the charge microstructure. Many of the particles in the stress chains are large ones that one might expect to be crushed by the focusing of the compressive stresses on them. A tumbling mill seems to essentially act as a vast array of randomly configured crushing rolls.

Simple breakage in DEM Previously particle breakage has been modelled in DEM by constructing the macro-particles from very large numbers of small sub-particles. Thornton (1996) used sphericaVdisc sub-particles in this way to model breakage of agglomerates and Potapov and Campbell (1994, 1996) used triangular/tetrahedral subparticles to model both impact and compression breakage. This allows very detailed analysis of particle breakage under a range of conditions. Unfortunately the computational co st of this approach is linearly

Recent advances in DEM modelling of tumbling mills

1311

dependent on the number of sub-elements. In an industrial scale simulation this would increase the cost by factors of 103 to 104 , depending on the accuracy that one wished to obtain for the breakage. Given that such 3D industrial scale mill simulations using 105 to 106 particles have only become in the last few years, the additional cost is prohibitive. An alternative approach is to use a rule base to identify the conditions under which particles would be expected to fracture and then to replace the fracturing parent particle with an appropriate assembly of daughter particles packed into the space occupied by the parent. This is a trivial operation in a DEM since complete information is known about aIl particles and complete control of their behaviour is possible.

Fig. 13 Stress chains in a stationary box.

Fig. 14 Stress chains in a tumbling mil!.

P. W. Cleary

1312

We use this approach to including particle breakage into DEM so that wc can apply this machinery to complete mill models of the type discussed earlier in this paper. Our mIe base includes two particle breakage mechanisms: 1. 2.

High energy impacts: Particles are broken if the cumulative energy loss during a collision exceeds a size dependent breakage energy Compressive fracture: Particles are broken if the total instantaneous force applied to a particle exceeds a specified size dependent compressive stress threshold.

The actual mies used here are relatively cmde and are presently the subject of extensive work. At this stage we wish only to demonstrate the methodology and the additional critical predictions that it allows DEM to make for mills. Currently the parent and daughter fragments are circular (2D) / spherical (3D) and tbeir sizes are chosen using geometrical criteria to optimally pack the space occupied by the parent particles. A minimum fragment size is nominated that defines the size of fines which are not resolved by the DEM simulation. Figure 15 shows a simple high speed two particle collision in three dimensions. The parent particles in a) collide and dissipate enough energy in the normal direction to fracture and are replaced by the daughter fragments. The cleavage plane between the primary fragments is chosen to pass through the contact point of the particles in the collision. The fragments are then free to move independently according to the residual velocities of the parent particles and are free to collide and break furtber if there is sufficient energy to do so. In future the fragment distributions will be taken from drop weight or similar calibration tests. a)

b)

c)

d)

e)

f)

• Fig. 15 Simple breakage of two particles ina high speed collision. Impact breakage in a rotating box A simple demonstration of the breakage model involves taking a partially filled box (in two dimensions) and rotating it. Figure 16 shows such an example. The box is initially haif full of particles with sizes between 80 and 100 mm. The particles are shaded according to their size. After 0.69 s the particles are about to avalanche for the first time. The solid initial packing of the particles means that the slope needs to be close to vertical before failure occurs. Therefore, the leading particles in the avalanche have an

l313

Recent advances in DEM modelling of tumbling mi Ils

unimpeded faH of most of the height of the box. In this example we use an impact breakage rule and we have reduced the breakage energies substantiaHy so that significant breakage occurs quickly so that we can demonstrate the behaviour. We use a 2 mm minimum fragment size (i.e. we are attempting to resolve ail fragments larger than 2 mm).

=

0.69 s

TIme = 0.87 s

TIme..

1.13 s

TIme

=

1.17 s

TIme

=

1.39 s

TIme

=

3.00 s

TIme

TIme

=

2.00 s

=

1.00 s

TIme =

1.258

=

4.00 s

TIme

TIme

Fig. 16 Impact breakage demonstration in a rotating two-dimensional box.

1314

P. W. Cleary

At 0.87 s the leading particle has hit the bottom of the box with sufticient energy to break. The fragments are shaded a light grey. At 1 s there have been a number of fracture events and the particles flow towards the left bottom corner of the box. The smaller daughter fragments rapidly fall through the particle microstructure and settle in the bottom corner. By 1.25 s there is already a significant bed of fines accumulated from the fractures produced by the tirst and largest avalanche. After 2 s, the bottom of the box is filled predominantly with smaller particles with some of the original particles embedded. By this time, most of the breakage events are concentrated in the toe region where the avalanching material strikes the wall. Occasionally a large particle is broken there, but mostly it is midsized particles that are seen to be trapped between a large incoming particle and either the wall or another large particle. This crushing of mid-sized particles is consistent with experimental observations. By 4 s, a significant mass of finer particles has been produced and more than half the original particles have been fractured. In such simulations the breakage rates of different size fractions can be measured and studied. In this case we see that the breakage rates of the original particles declines steadily with time with no breakage occurring after about 10 s. By this stage the large 'sea' offines that surrounds aIl the particles insulate the remaining large particles and pre vent them falling unimpeded for long enough to acquire sufficient energy to break. This demonstrates the importance of efficient removal of fines from tumbling mills and the penalty in terms of reduced breakage of large particles if they are not. Figure 17 shows the three-dimensional analogue of the previous example. The parent particles are again between 80-100 mm and fragments down to a minimum size of 10 mm are resolved. The dynamics are quite similar with many fines generated quickly and the se fines concentrating in the lower section of the cube. The large amounts of fines produces very fluid like motion. Breakage in a continuously filled mill Finally, we apply the breakage model to a tumbling mill. The mill details are not particularly important since this type of simulation can be performed for any size or style of tumbling mill in two or three dimensions. For demonstration purposes we use the same 5 m baIl mill as used in previous studies (Cleary 1998a, 2001b) rotating now at 70% cri tic a!. Figure 18 shows the operation of this mill model and the particles are shaded by their size. The charge initially consists of only 125 mm steel balls (dark grey). A stream of coarse feed particles (80-100 mm) is added continuously from a location just above the toe. A target charge mass of 35 tonnes / m of axial length is chosen. This produces a volumetric fill level of around 45%. The situation after 0.9 revs is shown in Figure 18 with the feed stream on the right. The feed material occupies a band just inside the liner with the steel balls occupying the region above. Breakage by both impact and compression is enabled, so these feed particles break when the flow circumstances dictate. The minimum fragment size is set to 10 mm to correspond to an exit grate size with this dimension. This means that fragments that should be produced that are smaller than this are automatically removed from the mi Il. This effectively means that sub-grate size particles have a very rapid axial transport rate and depart the mill immediately. So when the DEM simulation is setup in this way, this means that the rate of mass loss during breakage is equal to the rate of product generation. This then allows us to directly predict the solids discharge rate from the mill using DEM simulation. In addition, the feed rate is set equal to this discharge rate to maintain a constant fiIllevel in the mill mode!. By 1.9 revs there is a thick layer of feed particles surrounding a steel ball core. The diffusion rate of rocks into this region is quite slow. By 2.8 revs the miIl has reached the target J'ill level and the feed rate is changed to equal the solids discharge rate. Diffusion of rock particles into the core of the charge continues and smaIler particles (produced by breakage events) can now be observed in the cataracting stream. At 6.5 revs there is now plenty of rock evident in the core of the charge and there has been significant evolution in the charge size distribution. By 14 revs the charge flow pattern and distribution is close to equilibrium, but the size distribution continues to evolve becoming slowly finer. This can be seen by the graduaI lightening of the grey that shows the rock component of the charge. By 140 revs the charge PSD has stabilised and the miIl is in equilibrium. Since the charge PSD evolves dynamically from the feed PSD we are able to predict equilibrium rock PSD as weIl as the miIl throughput. In the next generation of this model we expect to be able to also make a

Recent advances in DEM modelling of tumbling mills

1315

plausible prediction of the product PSD. The discharge rate from the miU predicted by this model is critically dependent on the size dependent breakage energies and forces specified for the particles (as the y should be). This model is not limited to predicting steady state operating performance. The transient effect of changes in miII controls, rock hardness or feed rate can also be evaluated.

Fig. 17 Impact breakage in a rotating cube. FinaIly, we can record the location and type of aIl breakage events in the continuously filled mill. Figure 19 shows the distribution for the miII demonstration in Figure 18. High energy impacts produce breakage events (shown by the dark grey dots) that are predominantly near the surface of the charge where both the cataracting and avalanching streams of particle impact on the toe. High force (but low energy) compression produces breakage events (shown by Iight grey +s) lower in the charge under the toe and strongly clustered

1316

P. W. Cleary

around the lifters where high stresses are generated as the lifters accelerate the charge. The substantial difference in these distributions demonstrates that aIl breakage and size reduction mechanisms will need to be included in the model. Certainly one must progress beyond a focus on the effect of high speed balls in the cataracting stream. This is only a smaIl part of a much larger and more complex breakage story that needs to be understood if quantitatively accurate predictions are to be made.

Fig. 18 A 2D mill showing filling, breakage and radial fines transport.

Recent advances in DEM modelling of tumbling mills

1317

Fig. 19 Location and types of breakage events in the continuously fiIled mill.

CONCLUSIONS Recent developments in the modelling of the motions within tumbling mills have been described in this paper. They faIl into four major categories: 1.

2.

3.

4.

Increasingly realistic industrial scale models of baIl and SAG mills are now possible. In particular full SAG mill modelling is now feasible if the exact charge PSD is not used or 3D slice models can be used when an accurate PSD (with lower cutoffs of 20--25 mm) is required Increasingly detailed and complete predictions of mill behaviour are possible. In particular, high resolution predictions of liner stresses and wear rates and detailed collisional energy spectra (classified by collision type or particle size) allow a much clearer understanding of where the energy consumed in the mill is going and what effect it has both on the charge and on the liner. These are in addition to traditional predictions of power consumption and charge motion. Comparison with quality experimental data allows us to understand the accuracy of various types of DEM models and the penalties associated with the various assumptions made by these models. Such validation exercises are cri tic al in establishing the accuracy of existing models and improving their accuracy in the future. Here we have demonstrated that full 3D models representing the entire charge PSD are able to predict charge profiles and shoulder and toe positions with good accuracy. The inclusion of the new particle breakage model in conjunction with the continuously feed miIl model described here in principle aIlows the mill throughput and the charge PSD to be predicted directly from the DEM model for the first time. In future we also expect that the product PSD will also be able to be predicted. It is important to understand that the quantitative accuracy of these predictions will be criticaIly dependent on the accuracy of the experimental data that is used in the breakage model (size dependent breakage energies and compressive stresses) and the measured size and shape distributions of daughter fragments for each of the parent size classes. Accuracy will also depend strongly on the correctness of the size reduction mechanisms identified and included in the DEM model. Size reduction from attrition and chipping is also important and will need to be implemented as weIl. There are also important issues relating to size dependent axial transport rates and the effect of subsequent breakage of particles that are too smaIl to resolve during their transport from the mill. This initial model is only the first step in the long development of DEM as a complete predictor of miIl breakage performance.

With increasing completeness of the information that can be predicted by DEM along with the increasing accuracy of prediction and the increasing size of the models that are feasible to compute, DEM is now

1318

P. W. Cleary

approaching the point where it can genuinely be used as a tool to both better understand mills and to be assist in their design and optimization.

ACKNOWLEDGEMENTS My thanks go to: • • • •

Rob Morrison, Walter Valery and David Royston for their advice on specifying the generic 36' SAG mil!. JKMRC and particularly to Rob Morrison for making their scale model data available for the validation study described here and for many helpful discussions. Geoff Robinson for his assistance with the spatial smoothing of the stress and wear data. François Debroux from our CFD group whose assistance with pre- and post-processing was invaluable.

REFERENCES Barker, G. C., Computer simulations of granular materials, in: Granular Matter: An interdisciplinary approach, Ed. Anita Mehta, 1994, Springer-Verlag, NY. Campbell, C. S., Rapid Granular Flows, Annual Rev. Fluid Mech., 1990,22,57-92. Cleary, P. W., Predicting charge motion, power draw, segregation, wear and particle breakage in bail mills using discrete element methods", Minerais Engineering, 1998a, 11, 1061-1080. Cleary, P. W., The Filling of Dragline Buckets, Mathematical Engineering in Industry, 1998b, 7, 1-24. Cleary, P. W., Discrete Element Modelling of Industrial Granular Flow Applications, TASK. QuarterlyScientific Bulletin, 1998c, 2, 385-416. Cleary, P. W., and Sawley, M. L., Three dimensional modelling of industrial granular flows", in Proc. 2nd Int. Conf on CFD in the Minerais and Pro cess Industries, Eds. M. P. Schwarz, M. R. Davidson, A. K. Easton, P. J. Witt, and M. L. Sawley, 1999, pp. 95-100. Cleary, P. W., DEM simulation of industrial particle flows: Case studies of dragline excavators, mixing in tumblers and centrifugai mills, Powder Technology, 2000,109,83-104. Cleary, P. W., Modelling Comminution Devices using DEM, Int. J. for Numer. Anal. Meth. Geomechan., 2001a,25,83-105. Cleary, P. W., Charge behaviour and power consumption in bail mills: Sensitivity to mill operating conditions, liner geometry and charge composition", to appear: Int. J. Min. Processing, 2oo1b. Cleary, P. W., and Hoyer, D., Centrifugai mill charge motion and power draw: comparison of DEM predictions with experiment, Int. J. Min. Proc., 2000,59,131-148. Haff, P. K., and Werner, B. T., Powder Tech., 1986,48,239. Holst, J. M., Rotter, J. M., Ooi, J. Y., and Rong, G. H., Numerical modelling of silo filling. II: Discrete element analysis, J. Eng. Mech., 1999,125,94-110. Mishra, B. K., and Rajamani, R. J., The discrete element method for the simulation of ball mills, App. Math. Modelling, 1992,16,598-604. Mishra, B. K., and Rajamani, R. K., Simulation of charge motion in bail mills. Part 1: experimental verifications, Int. J. Mineral Processing, 1994,40,171-186. Poschel, T., and Buchholtz V., Comp1ex flow of granular material in a rotating cylinder, Chaos, Solitons and Fractals, 1995,5,1901-1912 Potapov, A. V., and Campbell, C. S., Computer Simulation of Hopper Flows, Physics of Fluids, 1996,8, 2884-2894. Potapov, A. V., and Campbell, C. S., Computer simulation of impact-induced particle breakage, Powder Tech., 1994,81,207-216. Potapov, A. V, and Campbell, C. S., A three dimensional simulation of brittle solid fracture, Int. J. Mod. Phys., 1996, C7, 717-729. Rajamani, R. K., and Mishra, B. K., Dynamics of ball and rock charge in sag mills, Proc. SAG 1996, 1996, Department of Mining and Mineral Process Engineering, University of British Columbia. Ristow, G. H., 1994, Granular Dynamics: A review about recent Molecu1ar Dynamics Simulations, Annual Rev. of Comp. Phys., 1, 275-308.

1319

Recent advances in DEM modelling of tumbling mills

Robinson, G. K., and Cleary, P. W., The conditions for sampling of particulate materials to be unbiased Investigation using granular flow modelling", Minerais Engineering, 1999,12,1101-1118. Schiifer, J., Dippel, S., and Wolf, D. E., Force schemes in simulation of granular material, J. Physique 1 France, 1996,6,5. Thornton, c., Numerical simulation of the impact fracture and fragmentation of agglomerates, J. Phys. D., 1996,29,424-435. Ward, B.K,. Optimisation of the design of lifter profiles in grinding mills, B.Sc. Metallurgy (Hons) Thesis, University of Queensland (unpublished), 1992. Walton, O. R., Numerical simulation of inelastic frictional partic\e-partic\e interaction, Chapter 25 in: Particulate two-phase flow, ed. M. C. Roco, 1994, pp. 884-911.

Correspondence on papers published [email protected]

In

MineraIs Engineering

lS

invited bye-mail to