Minerals Engineering 24 (2011) 181–187

Contents lists available at ScienceDirect

Minerals Engineering
journal homepage: www.elsevier.com/locate/mineng

CFD–DEM modelling of particle flow in IsaMills – Comparison between

simulations and PEPT measurements
C.T. Jayasundara a, R.Y. Yang a,*, B.Y. Guo a, A.B. Yu a, I. Govender b,c, A. Mainza b,
A. van der Westhuizen b, J. Rubenstein d
a
Lab for Simulation and Modelling of Particulate Systems, School of Materials Science and Engineering, University of New South Wales, Sydney, NSW 2052, Australia
b
Center for Minerals Research, Department of Chemical Engineering, University of Cape Town, P/Bag Rondebosch 7701, South Africa
c
Department of Physics, University of Cape Town, P/Bag Rondebosch 7701, South Africa
d
Xstrata Technology, Vancouver, BC, Canada V6B 2Z6

a r t i c l e i n f o a b s t r a c t

Article history: The IsaMill™ is a high speed stirred mill with a horizontal configuration that offers advantages such as
Available online 5 August 2010 energy efficiency and an inert grinding environment. A combined computational fluid dynamics (CFD)
and discrete element method (DEM) approach was developed to investigate the particle and fluid flows
Keywords: inside a simplified IsaMill™. The configuration of the mill was simpler than that of an actual IsaMill™ and
Computational fluid dynamics no feed flow or rotor was considered. The CFD–DEM model is a progression from earlier DEM only models
Discreet element modelling of ‘‘dry” systems which did not account for the fluid phase. The properties of flows at a macroscopically
Comminution
steady state, such as velocity field, distributions of particle velocity and acceleration in the radial direc-
Industrial minerals
Process optimisation
tion and power draw, were analysed. Detailed comparisons were carried out between the simulation
results and Positron Emission Particle Tracking (PEPT) measurements under similar conditions. The com-
parisons showed reasonable agreements, confirming that both techniques can capture the key features of
the flow. The discrepancies between simulated and measured results were discussed. The findings indi-
cated that the proposed model can be used to generate microdynamic information that is useful in lead-
ing to a better understanding of the underpinning physics of flow inside mills.
Ó 2010 Elsevier Ltd. All rights reserved.

1. Introduction In this light, a numerical model based on discrete element

method (DEM) was developed to simulate the dry particle flow
The IsaMill™ is increasingly being used in the minerals industry in a stirred mill (Jayasundara et al., 2006). Grinding performance
for grinding applications due to high energy efficiency and the in- in terms of impact energy was investigated for different operating
ert grinding environment. The mill is a type of horizontal stirred conditions. However, in a real system particle flow is affected by
mills with grinding discs fixed to a central shaft which rotates at the presence of slurry. Therefore, despite the fact that it adds great
a tip speed of around 20 m/s. During operation, the slurry feed is complexity to the simulation, it is necessary to include the slurry
introduced into one end and the ground product exits from the flow. A combined computational fluid dynamics (CFD) and DEM
other end. The high speed of the discs induces intense agitation approach was developed in the previous work to investigate slurry
of the media (usually ceramic but sometimes sand) which in turn flow in a stirred mill (Jayasundara et al., 2009). Slurry properties,
establishes the grinding action. Experimental and theoretical stud- such as density and viscosity, were varied to study their effects
ies on stirred mills have been reported in the past to understand on the flow properties in terms of flow velocity, power draw and
the grinding mechanisms (Gao and Forssberg, 1992; Berthiaux impact energy. The results showed that slurry density and viscosity
et al., 1996; Blecher et al., 1996; Kwade et al., 1996; Varinot have strong effects on flow patterns and impact energy inside the
et al., 1999; Eskin et al., 2005). They largely cover the macroscopic mill.
scale and the relationships established from these studies are However, the nature of the high speed stirred mill poses signif-
mostly empirical. A better understanding of the process can be icant practical difficulties in verifying the model. While experi-
achieved through detailed analysis of the microdynamic properties mental data such as overall flow pattern and power draw is
that are related to grinding media and slurry. comparable with simulation results, it offers little or no insight into
the microscopic details of flow. A major benefit of the model is that
it does describe the microscopic details but this advantage would
* Corresponding author. Tel.: +61 293856787; fax: +61 293855956. be undermined if it is difficult to confirm its predictions. A compar-
E-mail address: r.yang@unsw.edu.au (R.Y. Yang). ison of flow properties at particle scale with measured data is

0892-6875/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved.
doi:10.1016/j.mineng.2010.07.011
182 C.T. Jayasundara et al. / Minerals Engineering 24 (2011) 181–187

therefore crucial in order to realize the full potential of the model. @ðqf euÞ
þ r  ðqf euuÞ  r  ðegruÞ
Positron Emission Particle Tracking (PEPT) has been increasingly @t
P
used in the particulate research (Jones and Bridgwater, 1997; Lau- Fd;i
rent et al., 2000). PEPT is a non-invasive technique to capture the ¼ erp   r  ðeqf u0 u0 Þ þ qf eg ð4Þ
DV
trajectory of a selected tracing particle based on which to deter-
mine the flow field of particles in a system. Previous studies have where u, u0 , g and p are, respectively, the fluid velocity, fluctuating
compared PEPT with DEM simulations for particle flow in bladed fluid velocity, fluid dynamic viscosity and pressure; e and DV are
cylindrical mixers (Stewart et al., 2001) and rotating drums (Yang the local porosity and volume of a computational cell. The shear
et al., 2003). While the comparisons were reasonably good, no stress transport model (SST) is used to solve the Reynolds’ stress
work has been done with the presence of fluid which requires term qf u0 u0 (ANSYS CFX-11 reference manual). Due to the disc
CFD–DEM modelling. This work therefore aims to perform a de- rotating motion, the simulation domain is specified with a rotating
tailed comparison of CFD–DEM simulations and PEPT measure- frame so that the holes on the disc and entire computational mesh
ments on the flow in a simplified IsaMill™, focusing on the flow are fixed in the rotating coordinates system without mesh deforma-
properties in terms of particle velocity field, distribution of velocity tion. Details of the force models in Eqs. (1)–(4) can be found from
and acceleration in the radial direction, and the power draw of the previous work (Jayasundara et al., 2009).
mill. The CFD–DEM coupling is achieved by combining an in-house
DEM code developed at UNSW with the commercial software
2. Simulation and experimental method CFXÒ. At each time step, based on fluid flow field, DEM calculates
particle-related information based on fluid flow field, such as the
Flow in a mill is composed of discrete particle flow simulated by positions and velocities of individual particles, to determine poros-
DEM and continuum fluid flow simulated by CFD. In the DEM mod- ity and volumetric particle–fluid interaction force in the individual
el, a particle surrounded by fluid has two types of motions, trans- computational cells. CFD then uses these data to determine the
lational and rotational motions, determined by Newton’s second fluid flow field which then yields the particle–fluid interaction
law of motion (Cundall and Strack, 1979), forces acting on individual particles. Incorporation of the resulting
forces into DEM produces information about the motion of individ-
dm i X n
mi ¼ Fd;i þ Fp;i þ mi g þ ðFij þ Fsij Þ ð1Þ ual particles for the next time step, and the process continues until
dt the system reaches its steady state, as determined by analyzing
velocity field and power consumption.
dxi X
Ii ¼ ðRi  Fsij  lr Ri jFnij jxi Þ ð2Þ A fully coupled CFD–DEM model requires the exchange of infor-
dt mation at each time step. However, the exchange can be less fre-
where mi, Ii, vi and xi are, respectively, the mass, moment of inertia, quent when the system reaches its steady state at the
and translational and rotational velocities of particle i. Ri is a vector macroscopic level. This is because, at a steady state, fluid flow is
running from the centre of the particle to the contact point with its less sensitive to the change of particle flow and vice versa. Such
magnitude equal to particle radius Ri. lr is the coefficient of rolling treatment, which only exchanges CFD and DEM information in a
friction. certain time, can significantly reduce the simulation time, making
In the CFD model, the fluid field solved from the continuity and the modelling of large, complicated systems such as IsaMill™ more
the Reynolds averaged Navier–Stokes equations based on the local feasible. As this work studied the flow at a macroscopically steady
mean variables over a computational cell (Tsuji et al., 1992; Xu and state, this method was therefore adopted in this work. Data was
Yu, 1997). The governing equations are given by: exchanged every 10° of the disc rotation and intermediate stages
were obtained by linear interpolation. More frequent data transfers
@e at smaller rotation angles were also tested but no significant differ-
þ r  ðeuÞ ¼ 0 ð3Þ
@t ence was observed.

Fig. 1. Geometry of the three disc model from different directions (dimensions in mm).
 C.T. Jayasundara et al. / Minerals Engineering 24 (2011) 181–187 183

A 20 l mill with three discs was simulated as this setup allows used, i.e. glass beads (GB) and ceramic beads (CB), and the mill
the study of both flow near wall and flow in the middle. Fig. 1 was operated at different loadings and speeds. The values of the
shows the schematic representation of the mill used for this study. parameters considered in the simulations are listed in Table 1.
The values of the parameters considered in the simulations are PEPT experiments were carried out by UCT under similar condi-
listed in Table 1. The disc spacing, shaft diameter, and internal shell tions. In the experiments, the glass and ceramic tracers were made
diameter closely match those of the 20 l pilot scale mill. Note this radioactive. Gamma rays were then picked up by the two gamma
simplified system (and thus the model) does not represent a real ray detectors on each side of the mill to define a line on which
IsaMill™ since it does not include feed flow, a rotor, a full set of the tracer was located. Thousands of these lines were recorded
grinding discs and a separator. Nevertheless, the current system every second and the tracer position was determined using a trian-
provides a base system to validate the model before simulating gulation routine. The behavior of the whole mill can then be de-
more complicated systems. Two types of grinding media were duced from the behavior of the single particle if the tracking was
done over a sufficiently long period of time. More details of the
PEPT method used in this work is presented in a separate paper
Table 1 (van der Westhuizen et al., 2010).
Physical parameters used in the simulations.

Parameter Values 3. Results and discussion

3 3 3
Media density, q (kg/m ) 2.5  10 (GB), 3.7  10 (CB)
Particle–particle sliding friction coefficient 0.1 Fig. 2 shows the velocity fields of the fluid at different mill load-
Particle–wall sliding friction coefficient 0.3 ings. The particles are glass beads and the rotation speed
Young’s modulus, Y (N/m2) 2.0  107 X = 1200 rpm. In the radial direction (Fig. 2a), both loadings have
Restitution coefficient, e 0.7 (GB), 0.5 (CB)
maximum fluid velocities near the hole positions. 40% loading,
Poisson’s ratio, r
~ 0.29
Rotating speed, X (rpm) 800, 1200
however, has higher velocity than 80% loading due to smaller num-
Volumetric media loading (%), J 40%, 60%, 80% ber of particles in the mill. In the axial direction (Fig. 2b), the fluid
Total number of particles, N 112,000, 168,600, 224,800 velocity has a symmetric circulating field and its magnitude is
Fluid density, qf (kg/m3) 1000 higher near the discs. Such circulating flow was also observed in
Fluid viscosity, g (Pa s) 0.001
the previous CFD study (Blecher et al., 1996), which is due to the

(a) (b)
Fig. 2. Fluid velocity for 40% (top) and 80% (bottom) mill loadings at different cross sections: (a) radial plane (section Y–Y) and (b) axial plane (section X–X) when
X = 1200 rpm.
184 C.T. Jayasundara et al. / Minerals Engineering 24 (2011) 181–187

fluid being driven outward by the rotating discs and then diverted loading increasing from 40% to 80% at X = 1200 rpm. As one might
at the outer wall towards the shaft. expect, decreasing rotation speed from 1200 rpm (Fig. 4c) to
Fig. 3 shows the corresponding particle flows. The particles 800 rpm results in decreased particle velocity.
move towards the mill shell in the radial direction at J = 40% Fig. 5 shows the measured and simulated particle velocity in
(Fig. 3a) but are more uniformly distributed at J = 80%. In the axial the radial direction. Particle velocity gradually increases from
direction, 40% loading shows a symmetric particle distribution the shaft to the disc hole, then decreases again towards the mill
with fewer particles in the mid region and more particles near shell. Fig. 5a shows that the highest velocities for the glass beads
the end walls. As mill loading increases to 80%, the particles are are generally observed within the region of the disc hole, although
more uniformly distributed and more particles can be observed sometimes they shift towards the outer edge of the holes. There is
in between discs. a small change in velocity from the tip of the disc to the shell. The
Fig. 4 shows the comparison of simulations and PEPT measure- particle velocity decreases with increasing loading, consistent
ments on the particle velocity. To calculate the velocity distribu- with the observation from Fig. 4. With increase rotation speed,
tion, the mill was divided into a set of square cells of cross the particle velocity increases as expected due to the high energy
section area 4  4 mm2 and a length same as the mill length, then transfer from the discs. Fig. 5b shows the particle velocity in the
the velocity of the particles within each cell was averaged. As the radial plane for ceramic beads. Compared to the glass beads, the
single particle in the PEPT experiments cannot generate sufficient ceramic beads have lower velocity due to the higher density.
data to isolate the analysis to smaller sections of the system, The simulation results are overall comparable with PEPT data
Fig. 4 actually shows the velocity profiles from the entire length although the curves from the experiments are flatter. The discrep-
of the mill projected onto the end view of the mill. In other words, ancy may be attributed to slight difference in operation condi-
it is the entire data set collapsed onto one plane. The overall com- tions, limited data resolution in experiments and experimental
parison on velocity distributions is reasonable. Both experiments and numerical errors.
and simulations show a ring of high velocity near the disc holes Fig. 6 shows the acceleration of the particles in the radial plane.
which gradually decreases towards the shell. The simulations show In the experiments, the acceleration of the tracer was obtained by
more clearly defined rings compared to experimental results, calculating the change in velocity with respect to time at each loca-
which is expected considering measurement error and greater var- tion. The average acceleration in each of the 4  4 mm2 bins was
iation in the process. A gradual decrease in velocity can be seen for then calculated. The plots show that the acceleration generally

Fig. 3. Particle flow for 40% (top) and 80% (bottom) mill loadings at different cross sections: (a) radial plane (section Y–Y) and (b) axial plane (section X–X) when
X = 1200 rpm.
 C.T. Jayasundara et al. / Minerals Engineering 24 (2011) 181–187 185

Fig. 4. Comparison of the spatial velocity distribution for glass beads obtained from numerical simulations (top) and PEPT experiments (bottom): (a) 40%, 1200 rpm; (b) 60%,
1200 rpm; (c) 80%, 1200 rpm and (d) 80%, 800 rpm.

5
Disc
4
3
Hole
2
1 40%, 1200rpm
4

2
Velocity (m/s)

1
60%, 1200rpm
4 2
Disc
3

2 1
Hole
Velocity / (m/s)

1
80%, 1200rpm
80%, 800rpm
3
2
2

1
1
80%, 800rpm
80%, 1200rpm
0
0
20 40 60 80 100 120 20 40 60 80 100 120
r (mm) r (mm)

PEPT SIMULATION PEPT SIMULATION

(a) (b)
Fig. 5. Experimental and simulation velocity profiles in the radial plane under different mill speeds and mill loadings: (a) glass beads and (b) ceramic beads.

increases with increasing speed, decreasing loading or density. The beyond on the hole position (r > 75 mm). In the region closer to the
simulation results are often comparable to PEPT ones in the region shaft, however, there are discrepancies between two results. For
186 C.T. Jayasundara et al. / Minerals Engineering 24 (2011) 181–187

300
Disc
200

100 Hole
40%, 1200rpm

200
Acceleration (m/s)

100

60%, 1200rpm
200 40
Disc
150 30

100 20

Acceleration (m/s2)
Hole
50 80%, 1200rpm 10
80%, 800rpm
80

80 60

40
40
20 80%, 1200rpm
80%, 800rpm
0 0
20 40 60 80 100 120 20 40 60 80 100 120
r (mm) r (mm)
PEPT SIMULATION PEPT SIMULATION

(a) (b)
Fig. 6. Experimental and simulation acceleration profiles in the radial direction under different mill speeds and mill loadings: (a) glass beads and (b) ceramic beads.

6000 PEPT A - 40%, 1200 rpm (GB)

Simulation
Power draw (W)

B - 80%, 800 rpm (GB)

4000 C - 80%, 1200 rpm (GB)

D - 80%, 800 rpm (CB)

2000 E - 80%, 1200 rpm (CB)

GB - glass beads

CB - ceramic beads
0
A B C D E

Fig. 7. Power draw comparison for glass and ceramic beads at different operating conditions.

example, Fig. 6a shows that the measured acceleration for the glass power, 10% energy loss has been accounted when calculating the
beads has a monotonic decrease with radial direction at J = 40% experimental net power. In the simulations, torque on the discs
and X = 1200 rpm while the simulated one has a peak in the hole is caused by fluid and particles. At a particular time, each contact
region. Upward tails near the shaft can be observed in the simu- between a particle with the rotating shaft or discs produces a tor-
lated acceleration for 80% loading whereas the measured curve que on the mill which is the product of the contact force and dis-
has a flat curve or a drop in that region (Fig. 6a). The reason again tance between the contact point and centre line of the drum. The
could be due the lack of data resolution in experiments as individual torque is summed up to give the total torque which,
explained in the previous section. Since the velocity values are multiplied by the angular mill velocity, gives the power draw of
not exactly the same for both systems, the difference is further the mill at that particular time. Due to intense nature of the inter-
exemplified in the calculation of acceleration. Furthermore, the actions between particles and discs, fluctuation occurs in power
uncertainty in the interpretation of acceleration from velocity as draw. Averaging the power draw over a certain time gives a rela-
used in PEPT may also contribute to the difference. tively invariant value. The agreements are quite reasonable and
Fig. 7 shows the comparison of power draw for glass beads and the differences are well within the numerical and experimental
ceramic beads. Due to practical difficulty of measuring no load errors.
 C.T. Jayasundara et al. / Minerals Engineering 24 (2011) 181–187 187

4. Conclusions Cundall, P.A., Strack, O.D.L., 1979. A discrete numerical model for granular
assemblies. Geotechnique 29, 47.
Eskin, D., Zhupanska, O., Hamey, R., Moudgil, B., Scarlett, B., 2005.
A CFD–DEM model was developed to simulate particle and fluid Microhydrodynamics of stirred media milling. Powder Technology 156, 95–102.
flow in a simplified stirred mill. The flow properties of different Gao, M.-W., Forssberg, E., 1992. A study of the effect of parameters in stirred ball
milling. International Journal of Mineral Processing 37, 45–49.
particles, mill loadings and rotation speeds and grinding media
Jayasundara, C.T., Yang, R.Y., Yu, A.B., Curry, D., 2006. Discrete particle simulation of
were investigated and compared with PEPT measurements under particle flow in IsaMill. Industrial and Engineering Chemistry Research 45,
similar conditions. The overall comparison showed reasonable 6349–6359.
Jayasundara, C.T., Yang, R.Y., Guo, B.Y., Yu, A.B., Rubenstein, J., 2009. Effect of slurry
agreements, indicating the consistency of the two techniques.
properties on particle motion in IsaMills. Minerals Engineering 22, 886–892.
The discrepancies are attributed to many factors, including differ- Jones, J.R., Bridgwater, J., 1997. Using positron emission particle tracking to model a
ent operation conditions, relatively low sampling rate in PEPT than mixer as a series of interconnected tanks. The Institute of Chemical Engineers,
in the simulations and result interpretation. The comparison con- 409–412.
Kwade, A., Blecher, L., Schwedes, J., 1996. Motion and stress intensity of grinding
firmed that validity of the CFD–DEM model which can generate beads in a stirred media mill. Part 2: stress intensity and its effect on
microdynamic information for a better understanding of the flow comminution. Powder Technology 86, 69–76.
inside mills. Laurent, B., Bridgwater, J., Parker, D.J., 2000. Motion in a particle bed agitated by a
single blade. AIChE Journal 42, 3351–3363.
Stewart, R.L., Bridgwater, J., Zhou, Y.C., Yu, A.B., 2001. Simulated and measured flow
Acknowledgements of granules in a bladed mixer – a detailed comparison. Chemical Engineering
Science 56, 5457–5471.
Tsuji, Y., Tanaka, T., Ishida, T., 1992. Lagrangian numerical simulation of plug flow of
Authors (C T J, R Y Y and A B Y) are grateful to the Australia Re- cohesionless particles in a horizontal pipe. Powder Technology 71, 239–250.
search Council (ARC) and Xstrata Technology for the financial sup- van der Westhuizen, A., Govender, I., Mainza, A., Rubenstein, J., Tracking the motion
port for this work. of media particles inside an IsaMill using PEPT. Comminution 10. Cape Town,
South Africa, 2010.
Varinot, C., Berthiaux, H., Dodds, J., 1999. Prediction of the product size distribution
References in associations of stirred bead mills. Powder Technology 105, 228–235.
Xu, B.H., Yu, A.B., 1997. Numerical simulation of the gas–solid flow in a fluidized
Berthiaux, H., Heitzmann, D., Dodds, J.A., 1996. Validation of a model of a stirred bed by combining discrete particle method with computational fluid dynamics.
beads mill by comparing results obtained in batch and continuous mode Chemical Engineering Science 52, 2785–2809.
grinding. International Journal of Mineral Processing 44, 653–661. Yang, R.Y., Zou, R.P., Yu, A.B., 2003. Microdynamic analysis of particle flow in a
Blecher, L., Kwade, A., Schwedes, J., 1996. Motion and stress intensity of grinding horizontal rotating drum. Powder Technology 130, 138–146.
beads in a stirred media mill. Part 1: energy density distribution and motion of
single grinding beads. Powder Technology 86, 59–68.