See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/350614607

A combination of density-based clustering method and DEM to numerically

investigate the breakage of bonded pharmaceutical granules in the ball
milling process

Article  in  Particuology · April 2021

DOI: 10.1016/j.partic.2021.03.008

CITATIONS READS

0 111

5 authors, including:

Alexander Krok Peter Peciar

Munster Technological University Slovak University of Technology in Bratislava
17 PUBLICATIONS   74 CITATIONS    48 PUBLICATIONS   59 CITATIONS   

SEE PROFILE SEE PROFILE

Keith Bryan Sandra Lenihan

Cork Institute of Technology Cork Institute of Technology
24 PUBLICATIONS   45 CITATIONS    20 PUBLICATIONS   99 CITATIONS   

SEE PROFILE SEE PROFILE

Some of the authors of this publication are also working on these related projects:

Hip Joint Design View project

Development of finite element model to predict the behaviour of bilayer tablets during uniaxial tensile test View project

All content following this page was uploaded by Alexander Krok on 13 April 2021.

The user has requested enhancement of the downloaded file.

 Particuology 58 (2021) 153–168

Contents lists available at ScienceDirect

Particuology
journal homepage: www.elsevier.com/locate/partic

A combination of density-based clustering method and DEM to

numerically investigate the breakage of bonded pharmaceutical
granules in the ball milling process
Alexander Krok a,∗ , Peter Peciar d , Kieran Coffey c , Keith Bryan b , Sandra Lenihan a
a
Department of Process, Energy and Transport Engineering, Munster Technological University, Bishopstown, Cork T12 P928, Ireland
b
Department of Mechanical, Biomedical and Manufacturing Engineering, Munster Technological University, Bishopstown, Cork T12 P928, Ireland
c
Pfizer Newbridge, Littleconnell, Newbridge, Co. Kildare W12 HX57, Ireland
d
Institute of Process Engineering, Faculty of Mechanical Engineering, Slovak University of Technology, Nam. Slobody 17, 81231 Bratislava, Slovakia

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

Article history: Ball milling is widely used in industry to mill particulate material. The primary purpose of this process
Received 20 October 2020 is to attain an appropriate product size with the least possible energy consumption. The process is also
Received in revised form 5 January 2021 extensively utilised in pharmaceuticals for the comminution of the excipients or drugs. Surprisingly, for
Accepted 18 March 2021
ball mill, little is known concerning the mechanism of size reduction. Traditional prediction approaches
Available online 3 April 2021
are not deemed useful to provide significant insights into the operation or facilitate radical step changes
in performance. Therefore, the discrete element method (DEM) as a computational modelling approach
Keywords:
has been used in this paper. In previous research, DEM has been applied to simulate breaking behaviour
Ball milling
Granular pharmaceutical lactose
through the impact energy of all ball collisions as the driving force for fracturing. However, the nature
Density-based clustering of pharmaceutical material fragmentation during ball milling is more complex. Suitable functional equa-
Discrete element method tions which link broken media and applied energy do not consider the collision of particulate media of
Breakage different shapes or collisions of particulate media (such as granules) with balls and rotating mill drum.
Mill rotation speed This could have a significant impact on fragmentation. Therefore, this paper aimed to investigate the frag-
ABAQUS mentation of bounded particles into DEM granules of different shape/size during the ball milling process.
A systematic study was undertaken to explore the effect of milling speed on breakage behaviour. Also,
in this study, a combination of a density-based clustering method and discrete element method was
employed to numerically investigate the number and size of the fragments generated during the ball
milling process over time. It was discovered that the collisions of the ball increased proportionally with
rotation speed until reaching the critical rotation speed. Consequently, results illustrate that with an
increase of rotation speed, the mill power increased correspondingly. The caratacting motion of mill
material together with balls was identified as the most effective regime regarding the fragmentation,
and fewer breakage events occurred for centrifugal motion. Higher quantities of the fines in each batch
were produced with increased milling speed with less quantities of grain fragments. Moreover, the rela-
tionship between the number of produced fragment and milling speed at the end of the process exhibited
a linear tendency.
© 2021 Chinese Society of Particuology and Institute of Process Engineering, Chinese Academy of
Sciences. Published by Elsevier B.V. All rights reserved.

Introduction et al., 2015; Toker et al., 2016; Weeber & Bakker, 1988). This is due to
process simplicity and effectiveness during application. Equipment
Ball milling belongs to an essential unit operation for many for this process generally represents a hollow cylindrical drum with
industries (such as mineral, chemical, food and pharmaceutical) lifters rotating around its axis. At the same time, this drum is par-
widely used to mill particulate material (Carmen et al., 2019; Loh tially filled with milling media (metallic or non-metallic balls) and
particulate material (powders, granules or agglomerates) and the
aim is to reduce original size of the particles, granules or agglomer-
ates to smaller quantities. Depending on the loading conditions,
∗ Corresponding author. the type of material and associated properties, the integrity of
E-mail address: alexander.krok@cit.ie (A. Krok).

https://doi.org/10.1016/j.partic.2021.03.008
1674-2001/© 2021 Chinese Society of Particuology and Institute of Process Engineering, Chinese Academy of Sciences. Published by Elsevier B.V. All rights reserved.
A. Krok et al. Particuology 58 (2021) 153–168

material may be affected due to brittle deformation or fatigue frac- element method (DEM) as a computational modelling approach has
turing (Cutt et al., 1987; Gong et al., 2015). It often encounters the been growing in capability over the last decade and this method
requirement that mill product must have the prescribed fraction is now able to provide significant insights into the operation of
distribution or specific surface area. Less often, it is requested that existing and future mills. For example, axial transport and grate
the fragments in mill product have the same shape in the whole discharge of particular material during milling process have been
batch. Subsequently, fraction distribution and specific surface area investigated by using DEM (Cleary, 2006). Moreover, centrifugal,
of particulate material affect critical physical properties such as, e.g. tumbling, planetary and stirred milling process has been simu-
solubility, the reaction rate in chemical reactions, surface activity, lated and experimentally validated (Burmelster & Kwade, 2003;
mechanisms of interparticle bonds (Castellanos, 2005; Cho & Sohn, Hacifazlioglu & Korkmasz, 2020; Mishra & Rajamani, 1992). Bbosa
2016; Mosharraf & Nystrom, 1995; Mangwandi et al., 2010). et al. (2016) employed the DEM and the Positron Emission Particle
The primary purpose of size reduction is to attain an appro- Tracking (PEPT) method to investigate the charge motion of glass
priate product size with the least possible energy consumption beads of different sizes. The authors have shown that the energy
(Fuerstenau & Abouzeid, 2002). Therefore, the ball milling pro- of the cascading charge represents significant contribution, and
cess must be adequately designed, and equipment is required to therefore must be incorporated in net mill power determination.
operate at optimum conditions as even a slight improvement in Bian et al. (2017) used the combination of DEM and experimen-
the efficiency will provide an opportunity for economic benefit to tal approach to investigate the effect of lifter geometry on the ball
the industry. There are several variables to control this process. milling process. Qualitatively and quantitatively behaviour of parti-
Nevertheless, variables such as rotation speed and properties of cles was validated. It was observed that with a low height of lifters,
balls and milling material have the most significant influence on the efficiency of the balls in milling decreased.
the final quality of the product. It is also necessary to consider the Moreover, streams of mill particles and cataracting motion of
optimum filling level of both milling material and balls (Deniz and particles are mainly affected by a number of lifters and number
Onur, 2002). However, the number of lifters and their associated of balls. It has been investigated that with increasing the number
geometry also significantly affects the motion of milling material of lifters and the balls, the particle streams become increasingly
and draw of balls. dense and the number of cataracting particles increases (Orozco
Moreover, many mathematical equations based on experimen- et al., 2019; Xu et al., 2019). The influence of contact parameters in
tal work were developed to understand the correlation between DEM on charge motion and power draw was studied by Boemer and
mill power consumption and product quality. Mainly, these equa- Ponthot (2017). The authors explored that damper and tangential
tions were categorised into two groups. The first group was spring in contact model has no significant influence on power draw
focused on the relationship between energy requirements and size and charge motion. It was observed that with increasing normal
reduction of material (Bond, 1961). These equations usually are stiffness, the charge motion and the power draw are also rela-
evaluating the grindability and power requirement to drive the tively non-sensitive. Cleary and Morrison (2011) investigated the
mill. The second group was related to power and equipment mill behaviour of powder with three different sizes of the fraction (1.18
size (Herbst & Fuerstenau, 1980). Both methods are not influenced mm–2.8 mm) during milling process with fill level in ranging from
by the operating parameters and show a proportional relationship 0% to 150%. The powder mixture was trapped around the shoul-
between the rate of breakage and net mill power. In contrast, there ders of the balls, for low fill level. By increasing the fill level up to
are few theories to improve the correlation between the energy for 100%, the material started to flow along the cascading region, and
the milling and size reduction of particles, based on a specific prop- the mobility significantly increases. For higher fill level, the powder
erty modulus of material (Rittinger, 1867). By considering the initial surrounded the balls, and direct ball-ball collisions in the mill were
filling level of materials into equipment and the final product size reduced. It was also observed that for this level, the powder was
distribution, many different breakage frameworks have also been insufficiently mixed.
developed (Charles, 1957; Stamboliadis, 2002). Besides, coupled DEM and PBM method significantly improved
Additionally, research based on the optimisation of energy product sizes prediction from the milling process, by considering
consumption in milling has been undertaken by using phenomeno- dissipated energy obtained from DEM and selected breakage func-
logical grinding kinetics models and population balance modelling tion of particle determinates from experiments (Datta & Rajamani,
(PBM) approach (Austin et al., 1984; Nomura et al., 1991). A con- 2002). This previous approach assists to understand the particle
stant value of breakage rate and the breakage function, which is flow of material inside mills, but there are difficulties in interpret-
independent on the milling rate provided the fundamental size- ing the microscopic behaviour during processing. Also, it is not clear
mass balance equation for a fully mixed batch in milling operations which form of the equation to calculate energy is more suitable.
(Deniz, 2013). Nevertheless, it was proven by (Bozkurt & Özgür, Nevertheless, previous studies do not focus on mechanistic
2007; Celik, 1988) that the parameters of breakage function must modelling of material fragmentation during the ball milling pro-
be experimentally estimated in PBM, which is a time-consuming cess. Even DEM has been applied to simulate breaking through
process. Likewise, a back-calculation method which minimises the the impact energy of all ball collisions as the driven force for frac-
error between predicted and experimental size distribution was turing (Cleary & Morrison, 2011; Datta & Rajamani 2002; Tuzcu
also widely employed (Austin et al., 2007). In PBM, it was postulated & Rajamani, 2011), the nature of material fragmentation during
that for long milling times, the ball mills exhibited accelerating in ball milling is more complex. Suitable functional equations which
the mill drum, which affects the breakage rate (Celik, 1988) and this link broken particular media and applied energy were proven to be
was determined by using G–H method (Rajamani & Guo, 1992). useful (Cleary & Morrison, 2011; Morrison & Cleary, 2004; Vogel &
Moreover, a phenomenological theory was proposed to explain Peukert, 2004). However, these equations do not consider the colli-
multi-particle interactions in ball milling (Bilgili & Scarlett, 2005). sion of particulate media of different shapes and size or collisions of
These approaches mentioned above provide a reliable solution particulate media with balls and rotated mill drum. This could have
to predict particle size distribution for ball mills. However, some a significant impact on a better understanding of fragmentation
degree of error has been shown in predictions with a variation of during milling. Therefore, the aim of this paper was to investigate
operation conditions and feed characteristics. the fragmentation of bounded particles into DEM granules of dif-
Besides, traditional mathematical equations or PBM are deemed ferent shape and size during the ball milling process. A systematic
not very useful in the fundamental understanding of the process or study was hence performed to explore the effect of milling speed on
facilitate radical step changes in performance. Therefore, discrete breakage behaviour of DEM granules during the processing. More-

154
A. Krok et al. Particuology 58 (2021) 153–168

Fig. 1. Scanning electron micrographs of granule lactose with magnitudes (a) ×500; (b) ×1000.

Fig. 2. Granule size distribution of lactose chracterised by particle size analyser.

over, a combination of a density-based clustering method and DEM

were used to quantify the quality of the final product from milling
over time.
Fig. 3. Graphical illustration of sphericity and roundness.

Experimental
culture, powder and food or pharmaceutical engineering. It was
Material recognised the difficulty of measuring 3D surface areas of particles
in the past (Krumbein, 1941; Wadell, 1933) and therefore some
Unlubricated and un-sieved granule lactose (LFA Tablet press, uncertainty exists which is the most effective shape in describing
USA) was used in this study. This material is a common excipient the particle form. During the last decade, alternative definitions
used in the pharmaceutical industry. Scanning electron micro- of sphericity (Rodriguez et al., 2012) have been created such as
graphs of this granular lactose is shown in Fig. 1. The true density area sphericity, diameter sphericity, circle sphericity, perimeter
was measured using a Helium Pycnometer (AccuPyc II 1340, sphericity or width to length sphericity (Altuhafi et al., 2013;
Micromeritics, UK) and bulk density characterised by graduated Mitchell & Soga, 2005). Also, advances in optical image process-
cylinder method (World Health Organization, 2012) yeilding val- ing methods yielded digitisation of particle projections and created
ues of 1525.3 kg/m3 and 510 kg/m3 , respectively. The initial relative potential to automate the shape characterisation procedure.
density of granular lactose is 0.33 (−), which corresponds to the Additionally, methods such as Fourier analysis method
packing density. A particle size analyser (Malvern Morphology, (Bowman et al., 2001; Wang et al., 2005; Wettimuny & Penumadu,
Malvern Panalytical, UK) was employed to measure the shape of 2004), angularity index method (Sukumaran & Ashmawy, 2001;
the granules and also granule size distribution (Fig. 2), yielding par- Tutumluer & Pan, 2008), and a fractal technique (Arasan et al.,
ticle sizes d10 , d50 , d90 of granule lactose 202.7 ␮m; 800.5 ␮m and 2011) were used to calculate sphericity and roundness as more eas-
1530 ␮m, respectively. ily computational methods. However, a combination of statistical
locally weighted regression smoothing approach, together with a
Experimental granule shape characterisation K-fold cross-validation method seems to be most effective (Zheng
& Hryciw, 2015). Therefore, this approach is also used in this paper
In order to describe the granule shape in detail, there are many to calculate the most characteristic shape of granular lactose. Zheng
terms, quantities and definitions used in the literature (Mitchell and Hryciw (2015) found that the definition of sphericity, according
& Soga, 2005; Wadell, 1933). Sphericity and roundness are two to Krumbein and Sloss (1951), could be the most practical and easy
most common shapes denoted in the characterisation of particulate to determinate. Hence, in this study, it was assumed (Krumbein &
material, and this shape terminology is used in many disciplines, Sloss, 1951) that sphericity (ϕ) is defined as the ratio of the diameter
such as geotechnical engineering, soil science and mineral, agri- (Dint ) to diameter (Dcir ) as it is shown in Eq. (1) and Fig. 3.

155
A. Krok et al. Particuology 58 (2021) 153–168

Additionally, roundness (ω) is indicating the sharpness of par-

ticle corners. Using 2D projections of a certain number of granules
(N), the roundness was defined as the ratio of the average radius of
curvature of the corners of a particle (ri ) to the radius of the max-
imum inscribed circle (rinit ) (Bareither et al., 2008; Cabalar et al.,
2013; Chapuis, 2012; Cho and Sohn, 2016; Mitchell & Soga, 2005;
Santamarina & Cho, 2004; Shin & Santamarina, 2013; Wadell, 1935)
as it is shown in Eq. (2) and Fig. 3.

Dint
ϕ= (1)
Dcir
 N 

ri /N
i=1
ω= (2)
rint

An image-based method was used to quantify the shape of

the granules. One hundred and fifty systematically capture image Fig. 4. Experimental sphericity and roundness of granular lactose.
frames of selected scan areas were collected in the form of 2D pro-
jections using commercial Morphologi 3D software. Each image
frame detects one irregular particle within a size range of 10–1550
␮m. Subsequently, MATLAB code developed by Zheng and Hryciw
(2015) was used in this study to calculate sphericity and roundness
of the particles.
For the calculation of sphericity and roundness of each consid-
ered 2D projection, the outline of each particle picture was first
discretised. During discretisation, multiple numbers of outer points
were created in each image followed by determining the mini-
mum circumscribing circle with a diameter of Dcir . The minimum Fig. 5. Representative granules with (a) SHAPE A; (b) SHAPE B; (c) SHAPE C used in
five outer points were connected. If other points were within this this study.
circle, then this was considered as the minimum circumscribing
circle. Otherwise, the point which lies furthest outside of the cir-
cle, and therefore a new minimum circle using another set of the the third type of granule with “rounded” shape (denoted by yel-
minimum number of outer points need to be found. The proce- low colour), it was between 0.6–0.85 for sphericity and 0.25–0.6
dure was repeated until no points lie outside the circle. A similar for roundness. Based on revised experimental data (Fig. 4), three
procedure was used for determining the minimum circumscrib- examples with elongated (Fig. 5a), blocky (Fig. 5b) and rounded
ing circle with a diameter of Dinit . Deviations of circles with Dcir (Fig. 5c) shape of granules were selected for further investigation
and Dinit and surface of particles from its mean plane were sub- of the ball milling process.
sequently used to characterise roughness at a given scale. In this
study, surface roughness through computational geometry was
addressed by non-parametric fitting technique LOESS. Fitting strat- Mill configuration
egy without any prior specification of a functional relationship
between the points is the main advantage of this technique. Dur- The lab-scale ball mill was defined as a steel cylindrical shell
ing this procedure, each data point of the circle with Dcir and Dinit with an inner diameter of 50 mm and a depth of 10 mm, respec-
was replaced with the smoothed value determined by a locally tively. The shell drum was fitted with eight lifters, each with
weighted regression. Hence, each connected outer point and its corresponding heights of 1.2 mm and width of 0.8 mm. It was
nearest neighbouring points over a span distance were used in shown that a smaller number of balls is more efficient since there
the process. The selection of an appropriate span was the key is less over-grinding as well as in a reduction in power consump-
aspect of the fitting procedure to ensure that the LOESS curve is tion (Orozco et al., 2019). Ten steel balls with a density of 7800
as close as possible to the mean surface. Later, K-fold trial-error kg/m3 and a diameter of 1 mm were chosen. The mill configuration
cross-validation techniques determined the optimal span value. used in this study is illustrated in Fig. 6. For steel drum and balls, it
The best-fit circle was found by minimising the sum of the squares was defined that Young’s modulus is 200 GPa, and Poisson’s ratio is
of the distances between the points and the circle. To complete 0.3 (−). Individual zones such as empty zone, dead zone, abrasion
Wadell’s roundness calculation, it was also important to fit corners zone and impact zone often observed in drum space during mate-
of the particle surface with appropriate circles. In this step, Fourier rial milling are illustrated in Fig. 6(b). Additionally, cascading and
analysis, angularity index and the fractal technique were applied cataracting motion of material in this figure is also demonstrated.
to measure curvatures over the entire particle outline. The sharper Lab ball milling equipment involves, at least, billions of particles
corners with larger curvature required longer segments, whereas and these macro-scale problems, cannot be effectively modelled
flatter corners required line segments. yet due to the limitation of computer resources currently available.
Sphericity-roundness characterisation of lactose can be seen in Therefore, the discrete element model assumed to have been sim-
Fig. 4. From this figure, it is shown, that for the first type of gran- plified. The DEM model has the same configuration as in lab scale
ule shape (denoted by red colour) the sphericity variate between device. Moreover, this model is also assumed to have the same load-
0.75–0.93 and roundness between 0.2–1.0. This type of granules can ing/boundary conditions (such as rotation speed and filling volume)
be considered as “elongated”. For the second type of granules with as in the real scale milling process. In this study, only the size of the
“blocky” shape (denoted by blue colour) the sphericity and round- drum was reduced. The intention of this study was to present how
ness dispersion are narrower in ranges 0.65–0.92 and 0.25–0.6. For the combination of clustering method and DEM can be useful in

156
A. Krok et al. Particuology 58 (2021) 153–168

Fig. 6. (a) Schematic ball mill configuration used in simulation; (b) Schematic view of material motion in the drum.

analysing the mill’s final quality at the end of the process and not
validating experimental data.

Numerical

Theoretical principles and aspects of the discrete element

method (DEM) as a modern numerical approach were described in
many studies (Cundall & Strack, 1979; Ghaboussi & Barbosa, 1990;
Mishra & Rajamani, 1992; Wu et al., 2003). At the same time, this
method was proven to be considerably useful in different areas. In
this paper, DEM was employed to study the ball milling of pharma-
ceutical granules of different sizes and with non-spherical shapes.
Fig. 7. (a) Contact model diagram of Hertz-Mindlin; (b) Model breakage force
In ABAQUS commercial software, each particle was modelled as (Fbound ) defined by the multi-point connector.
a PD3D single element surrounded by rigid spheres of specified
radii during an explicit dynamic simulation and these particles col-
provides a spring repulsive force and dashpot to dissipate a portion
lided with each other but also with rigid finite-element-based drum
of the relative kinetic energy.
surface. PD3D element has been defined with displacement and
rotational degrees of freedom. Also, general contact procedure with 4 ∗ ∗
Kn = E R ın (5)
contact inclusions was used for particles, and it was considered that 3

particles in the current study could be involved in multiple con- Kt = 8G∗ R∗ ın (6)
tact interactions simultaneously. ABAQUS 2018 software has some
limitation regarding DEM modelling. Particularly, python scripts The relative normal and tangential velocity and ϑtrel ) (ϑnrel
which were required to be incorporated into ABAQUS, do not allow behaves as incremental spring and correspond to elastic normal
to define coefficients of restitution. Additionally, rolling friction is and tangential deformation. On the other hand, the dashpot repre-
ignored for the contact between the elements. Conversely, the coef- sents normal and tangential plastic deformation of the contacting
ficient of rolling friction has been widely discussed, as evidenced surfaces.


by several publications dealing with this topic (Chen et al., 2017;

lnε 5 
Hlosta et al., 2020a,b; Markauskas et al., 2015). When spherical Cn = 2  2E ∗ R∗ ın m∗ (7)
6
particles are used, rolling friction should be included, but this is ln2 ε + 2
unnecessary when working with non-spherical particles. 
lnε 5
Moreover, the non-linear Hertz-Mindlin contact model was Ct = 2  Kt m∗ (8)
6
used for simulating interactions between particles. This model ln2 ε + 2
belongs to the standard model in DEM. The normal and tangential
R1 R2
forces (Eq. (3) and Eq. (4)) are calculated based on the Hertz con- R∗ = (9)
R1 + R2
tact theory (Hertz, 1942) and Mindlin & Deresiewicz work (Mindlin,
1949; Mindlin and Deresiewicz, 1953), while the contact between 1 − 12 1 − 22
the particles is represented by the spring and damper (Fig. 7a). E∗ = + (10)
E1 E2
Force-displacement (F-␦) response using the Hertz-Mindlin
Fn = −Kn ın + Cn ϑnrel (3)
model in the tangential and normal direction depend on loading
  history and can be calculated by employing Eqs. (3–10). Where
Ft = min Fn , K t ıt + Ct ϑtrel (4) E1 , ␯1 and E2 , ␯2 are Young’s modulus and Poisson’s ratio, and R1
and R2 represent the radii of the two contacting particles and G*
The spring and damper are used to represent the elastic, and is equivalent shear modulus, respectively. ın and ıt represent the
dissipative contacts (Eqs. (5–8)) between the particles and friction displacement of a particle in a normal and tangential direction.
coefficient () with sliding block represent the friction between From these equations,
√ √ normal contact stiffness can be derived as
the particles limited by Coulomb friction criteria (Ft ≤ Fn ). This follow: Sn = 2E ∗ R∗ ı, while this parameter represents the slope

157
A. Krok et al. Particuology 58 (2021) 153–168

of the force-displacement curve. Sn beyond which the normal con- were sufficient to ensure that there is minimum porosity in each
tact force (Fn ) increases linearly can be marked as Smax , and this granule. Additionally, it was found that with increasing the number
parameter represents the maximum value of the normal contact of particles, the computation time increased proportionally, while
stiffness, which must be defined in the ABAQUS model. granules’ response behaviour during milling changed only slightly.
Only one type of pharmaceutical excipient was considered in The drum should not be underfilled or overfilled with balls and
this study with Young’s modulus of 3 GPa and Poisson’s ratio of 0.2 milling material. Overfilling tends to generate fines at the top of the
(−) for all particles contain in granules (Perkins et al., 2007). Nor- drum, which reduced breakage impact. If there is a low filling level,
mal contact stiffness, Sn, is varies with displacement and highly then the ball-to-ball contact against the breakage rate. Based on
depends on the properties of granules but also on the correspond- previous experiences (Austin et al. 1984; Bond, 1958; Napier-Munn
ing size of the granules. Three different granules’ sizes have been et al., 1999), the milling material should not exceed 45% of the drum
considered for each simulation, and therefore three different values volume. For grate discharged mills material should occupy less than
should be implemented. However, commercial ABAQUS software 50% of drum volume. In this study, it was considered that the gran-
for the Hertz contact model allowed the incorporation of only the ules and balls fill 20% of drum volume. There is a preferred ratio of
maximum value for the average size of the granule. In this case,
≈ 0.4 (Eq. 12) and the volume fraction of voids between 0.6 and
the Smax employed had a value of 3.313 kN/m. More importantly, it 1.1 for balls and occupied granules (Shoji et al. 1982). Nevertheless,
was essential to ensure that the time increment used for the analy- pharmaceutical excipients are more breakable with much smaller
sis was small enough to avoid numerical instabilities. For example, granules, and therefore a value of
≈ 0.1 as it was postulated in a
contact interactions among particles can impact the appropriate previous study (Zhang et al. 2015; 2018).
time increment size, and incremental motion can be influenced
drastically if particle velocities become very large. The stable time Vg

= (12)
increment size tends to be proportional to the square root of mass Vb ε
and inversely proportional to the square root of stiffness:
 Where Vg (m3 ) is the volume of mill occupied by the granules; Vb
m∗ (m3 ) is the volume of the mill occupied by balls and ε (−) rep-
t = Q (11)
k resent the voids between these two medium respectively. Three
different sizes (0.6 mm; 0.8 mm and 1.2 mm) of each type of DEM
In general, Q should not exceed between 0.1–0.4 while m∗ and granules have been incorporated in this study, which represented
krepresents the reduce particle mass and contact stiffness, respec- 1220 granules with 21960 particles for SHAPE A, 915 granules with
tively. Minimum step time was usually less than 1 ␮s and each 16470 particles for SHAPE B and 610 granules with 10980 parti-
simulation take at least 48 h. The experimental determination of cles for SHAPE C for every batch filled into the milling drum before
the internal friction defining as interactions between the particles initiation of the ball milling process. Size of the granules has been
themselves is often challenging, and measurement is even more chosen based on experimental granule sized distribution (as pre-
difficult for grains of a different shapes. From previous observation sented in Fig. 2). The sizes of the particles in DEM granules represent
(Krok et al. 2016; Krok and Wu, 2019), granular lactose used in this the typical size of the monohydrate lactose particles in pharma-
study behave as very brittle material with poor flowability. There- ceutical applications (Bonakdar & Ghadiri, 2018; Krok et al., 2016;
fore, it was decided to use a higher value rather than lower value. Verheezen et al., 2004). Considering previous studies (Guenette
Additionally, stainless steel millers are permanently controlled and et al. 2009; Huang et al. 2013; Mangwandi et al., 2010), it was
clean in the pharmaceutical industry regarding cleaning validation. shown in experimental particle size distribution data that typical
It was assumed that the friction coefficient between the wall of size of monohydrate lactose particles is usually in the range 10–120
the drum and granules is low due the regular cleaning practices ␮m. Therefore, the particles’ size with 0.03; 0.037 and 0.05 mm in
adopted by the pharmaceutical industry. Besides, the balls made granule, belongs to d10 –d60 of this material’s typical particle size.
by chrome steel were considered with a small degree of roughness, In most cases, roll compaction and the shear milling process are
and it was assumed that wall friction of the balls was slightly higher. used as a part of dry granulation for milling of excipients such as
In this model, the inter-particle friction coefficient was defined as lactose. It can be noted that there are only a few studies in rela-
0.4 (−). The friction between the drum and balls as 0.2 (−) but also tion to ball milling of pharmaceutical excipient. Lab scale milling
friction between the drum and milled material was also considered is relevant as often it is necessary to work with a small amount of
(with wall friction coefficient 0.1 (−). material or when roll compactor is not available. However, advance
modelling approaches such as DEM can improve understanding
DEM granule shape characterisation how this process could be used effectively in the pharmaceutical
industry.
Software ABAQUS allows several types of constraints to be The ball milling of granules is a mechanical process, during
defined depending on the type of application. For example, multi- which the granular material is a fracture with consequent disrup-
point constraints connecting nodes with available components of tion of integrity. By disconnecting a certain number of granules
relative motion and linear or non-linear force versus displacement during this process, a larger number of smaller fineness particles are
behaviour can be defined by beam-type multi-point connectors formed. Therefore, it was considered in this study that connected
with specifying breakage force value. Thus, the granules have been particles into the granules by beam-type multi-point constraints
defined as a grain of complex shapes by bonding spherical parti- can be damaged. Moreover, it was stated in previous studies (Cho
cles together using compliant connectors. In particular, overlapping et al., 2006; Goh et al., 2018; Krok et al., 2014; Kumar et al., 2014)
particles (Wang et al., 2006) were connected by beam-type multi- that there is a direct correlation between the size and shape of
point constraints instead of adhesion integration. It is noted that granules and their corresponding strength. Additionally, granules’
this approach may not replicate the precise geometry of actual bonded properties significantly affect breakage behaviour, which
grains. Eighteen overlapping particles were used to distinguish directly impacts the milling process (Metzger & Glasser, 2012).
each type of granules to obtain a closer approximation of the exact Bonded particles in granule were connected via linear breakage
shape, as shown in Table 1. Three representative DEM granules force-displacement type of connectors. The maximum value of
(SHAPE A; SHAPE B and SHAPE C) were defined, while DEM gran- breakage force and displacement was defined individually between
ules’ shape is shown in Table1. It was determined that 18 particles the particles for each size of the granules. The approach for estima-

158
A. Krok et al. Particuology 58 (2021) 153–168

Table 1
Characterisation of DEM granules.

DEM granules Name Sphericity (−) Roughness (−) Number of Size of granules Size of particles
particles (−) (mm) (mm)

0.6 0.03
SHAPE A 0.851 0.505 18 0.8 0.037
1.2 0.05

0.6 0.03
SHAPE B 0.794 0.437 18 0.8 0.037
1.2 0.05

0.6 0.03
SHAPE C 0.752 0.373 18 0.8 0.037
1.2 0.05

tion of breakage force for different size of the granules is described it would be possible to measure irregular dry granules’ breakage
in Section “Estimation of breakage force for DEM granule”. force.
Subsequently, the whole agglomerate as one united, while the
fragmentation is investigated via impact breakage test (Liu et al.,
2010; Potapov & Campbell, 2001; Thornton et al., 1997; Wu et al.,
Estimation of breakage force for DEM granule 2008). It was observed that rearrangement of the particles in
agglomerate and the level of adhesion force plays a crucial role
The size reduction process is one of the most critical processes in the results. In the above papers, auto-adhesive particles were
in raw material transformation. Reasonable knowledge of breakage specified within surface energy, and also relatively weak agglom-
behaviour at the single-particle scale can improve the understating erates with low adhesion interaction were tested rather than strong
of breakage at the process scale. Therefore, extensive studies have agglomerates.
been undertaken to estimate breakage force for wet granules (Fu The impact breakage test was employed to estimate the crit-
et al. 2004; Lian et al. 1998; Reynolds et al. 2005; Wu et al., 2003). ical force for damaging of granules in this study, similarly as in
It is expected that dry granules, in contrast to wet granules, will (Metzger & Glasser, 2012). In this aforementioned study, a regu-
exhibit different breakage behaviour due to the different nature lar cubic lattice compromised of 5 particles in each direction were
of the particle bonding forces. In general, detailed studies in rela- used to represent agglomerates. Simultaneously, the effect of the
tion to dry granules’ measurement are lacking. In most cases, it particle and bond properties on breakage behaviour of agglom-
is possible to prepare relatively big and almost regular shape wet erates was investigated. Moreover, they investigated the effect
granules; however, dry granules are usually small, much stronger of impact velocity on the breakage of the agglomerates. It was
and irregular shaped. Hence, the irregular shape of dry granules presented that the stiffness and strength ratios had little impact
would complicate this study because it is impossible to ensure on the results, but with increasing the bond and contact radius
that granule will drop on the same fracture plate. Even a range of of particles, the strength of the bonds increased. It was shown
experimental systems are accessible to characterise the breakage that with increasing impact velocity, the breakage of agglomer-
of particles or granules; for dry granules, it is difficult to determi- ates increased for both low and high bond strengths. However,
nate the specifics due to the small length of the particle and times the number of fragments varied. Particularly, for high strength,
scale of fracture behaviour. Additionally, it is incredibly challenging a smaller number of fragments were generated. Moreover, they
to experimentally measure energy dissipation processes (Mishra & investigated that an agglomerate is stronger with 1000 small parti-
Thornton, 2001; Nisoli et al. 2004; Refahi et al., 2010; Salman & cles than an agglomerate with 125 bigger particles. This is justified
Gorham, 1997; Tavares & King, 1998). as agglomerates need to overcome a higher number of bond inter-
In general, disruption of the equilibrium of internal binding action. Tavares and King (1998) investigated the deformation and
forces between the particles in granules is the physical essence fracture of single particles by using the Ultrafast Load Cell UFLC.
of comminution. Particles emphasised by the external forces are The measurements were performed under impact loading simi-
subjected to interatomic attractive van der Waals and repulsive lar to the conditions in typical industrial ball mills. Particularly,
Born forces. Hence, both forces act simultaneously and create a cer- the effect of material type, particle size and particle shape on
tain balanced equilibrium. However, with increasing the external the three fundamental fracture characteristics of brittle materials
forces, the distance between the atoms is changing. When a certain were investigated. Twenty different mineral materials were exam-
distance is exceeded, the interaction of the force field between the ined with the variation of shape factors between 0.363–0.524 and
atoms in granules is disturbed, and new surfaces are formed in the densities between 1040–5240 kg/m3 . Linear correlation between
form of fracture surfaces. There are two most popular theoretical fracture energy and particle strength was also explored. They illus-
paradigms to understand the failure of granular materials. Rumpf trated for gelena, spharelite, quartz corundum the linear correlation
(1962) consider that granule fails by the simultaneous rapture of between particle stiffness and cumulative distribution of particles.
all the bonds along a fracture plane. Alternatively, Kenderll (1988) In contrast, geological materials have an inherent resistance to
postulated that a granule failed through crack nucleation and prop- mechanical size reduction impact. Additionally, rocks and ores had
agation. Both approaches are mainly applicable to wet granules. higher fracture energies, than pure minerals. Tavares (2009) also
The well-known theoretical prediction method has been developed demonstrated that particle size significantly affects the damage
based on experimental data for wet granules. Studies for dry gran- accumulation coefficient, while particle shape has only neglected
ules have been made only by the DEM method (Metzger & Glasser, influence. He suggested that this coefficient being higher for com-
2012). Currently there is no direct experimental device by which plex microstructures material, but lower for more brittle material.

159
A. Krok et al. Particuology 58 (2021) 153–168

Fig. 8. Fragmentation of granule during impact breakage test.

In this paper, the strategy for estimation of breakage force of

granules was deviated slightly from above. It was assumed that the
ribbon’s strength prepared by roll compaction process is the same Fig. 9. Estimated breakage force for three different shapes of granules with three
across the whole volume. Even from an experimental perspective, different sizes.
the ribbon strength is only slightly smaller close to the edge as it was
explored in previous studies (Krok and Wu, 2019; Miquélez-Morán
et al. 2009). Consequently, granules which were produced from
the ribbons would have similar strength independent of size. In
other words, it is assumed that all considered granules should have
similar strength independently on size and the breakage force per
granule is varies due to different the fracturing surface. Although
there is no experimental evidence of this behaviour, the results can
be justified and may be observed with the right set of experiments.
Through an alteration of desired impact velocity, the breakage force
which pushed overlapping particles apart during the dropping of
granules has existed unless the contact exclusions reached spec-
ified cut-off critical value. Graphical representation of breakage
force defined by the multi-point connector is shown in Fig. 7(b).
Notably, each type of granule (SHAPE A; SHAPE B and SHAPE C) with
three different sizes (0.6 mm, 0.8 mm and 1.2 mm) was subjected
by this impact measurement. Each granule of different strength was
Fig. 10. Number of fragments varying with impact speed for three different shapes
placed at the top of the cylinder in hight of 50 mm, and then granule of granules of three different sizes.
dropped towards the flat plate with eight different impact speeds
(0.46 m/s, 0.5 m/s, 0.54 m/s, 0.58 m/s, 0.62 m/s, 0.66 m/s, 0.7 m/s,
0.74 m/s). A typical DEM simulation of impact test is illustrated in Density-based clustering
Fig. 8.
The systematic trial-error approach was applied to estimate crit- Unsupervised machine learning clustering method (CM) belong
ical breakage force for all individual DEM granules. It was assumed to a popular technique in data science (Xu & Tian, 2015). The task
that the granules are impacting at the high velocity, whereas there of this method is to group a set of objects in such a way that
is minimal breakage at low impact speed (Ning et al., 1997). Hence, objects in the same group are similar to each other than to those
it was suggested, that for an impact speed of 0.46 m/s, the granule in other groups. Cluster analysis itself is not representing one spe-
will not break, and it will remain as one unit. For other applications, cific algorithm. The suitable clustering algorithm and parameter
this initial value can be easily changed. On the other hand, with neg- settings depend on the individual data set and intention to use
ligible change in breakage force value, the same granule will break the results (Saxena et al., 2017). The density clustering method
with an impact speed of 0.5 m/s. For example: with breakage force was used in this study to identify individual fragments during
6.3089 mN, the SHAPE A granule (0.6 mm) was not broken with an the ball milling process (Nanda and Panda, 2015; Sander et al.,
impact speed of 0.46 m/s, but it breaks with an impact speed of 0.5 1998).
m/s. All estimated values of the breakage force for each granule are DBSCAN (Density-Based Spatial Clustering of Application with
shown in Fig. 9. It was explored, that size of the granules signifi- Noise) algorithm goes through the positions of each particle within
cantly effects the value of estimated breakage force compare to the the fragments and its counting the number of particles nearby. For
granule shape or rearrangement of the particles when only slight this purpose, the MATLAB code developed by Sander et al. (1998)
changes have been observed. was used in this study. Typical identification of fragments and
Three main types of breakage have occurred during impact test: their sizes is presented in Figs. 11 and 12. Rearrangement of frag-
(1) no-breakage; (2) intermediate breakage and (3) complete dis- ments was taken from the impact breakage test, previously shown
integration when all original bonds between the particles were in Fig. 8. As an example, for granule with SHAPE A (size of 0.6
broken. No-breakage and intermediate breakage results yielded a mm) and impact speed of 0.5 m/s, eight fragments were created
range of shapes and size distributions. A strong correlation between (Fig. 11) and also three different sizes of fragments were identi-
the impact speed and a number of fragments has been found as fied (Fig. 12). MinPts and Eps are two critical parameters which
it is shown in Fig. 10, and results are in agreement with previ- need to be defined for DBSCAN method. It was assumed that the
ous work (Kafui & Thornton, 2000; Moreno et al., 2003; Samimi minimum number of particles which can represent individual frag-
et al., 2004). Small granules brakes into primary particles for most ments is equal to one (MinPts = 1) and a distance used to locate
impact speeds. Conversely, the number of particles for big granules the points in the neighbourhood between the bounded particle
increased proportionally. was defined as 0.5 mm (Eps = 0.5 mm). Each breakage event was

160
A. Krok et al. Particuology 58 (2021) 153–168

Fig. 11. Identification fragment numbers at the end of the impact breakage test by Fig. 12. Identification of fragment sizes at the end of the impact test by using
using DBSCAN method. DBSCAN method.

recorded into *txt file for further post-processing and data was later granule shape characterisation”. In Fig. 13, the rotation speed was
analysed. set up to 95 rpm and the total simulation time was 5 s. Moreover,
ten steel balls with a size of 1 mm have been filled into each batch.
Results and discussion Over the last decade, many aspects have been included in ball
milling simulation, which improved mechanistic understanding of
Ball milling simulation the process. However, mechanical fracturing was rarely consid-
ered. Limited research on the effect of rotation speed on breakage
Typical behaviour of a ball milling simulation with time is pre- behaviour in ball milling exists (Bian et al., 2017). From experi-
sented in Fig. 13. Each DEM model incorporates a batch of granules ments in previous studies, it was noted that during the milling, a
with three different shapes (SHAPE A, SHAPE B and SHAPE C) and small amount of the granules aggregated together during the whole
three different sizes (0.6 mm, 0.8 mm, and 1.2 mm) and each batch process for the specific operating condition. In contrast, most of
includes 2745 granules which represent 49410 particles. All details the granule fragments during milling are reducing from their orig-
pertaining to granule configurations are presented in Section “DEM inal size, and many of them also break into individual particles.

Fig. 13. Typical distribution of mill material in the drum for a DEM simulation of ball milling at a rotational speed of 95 rpm with time (a) 0.2 s; (b) 0.75 s; (c) 2.25 s and (d)
5 s.

161
A. Krok et al. Particuology 58 (2021) 153–168

process. Therefore, this paper is mainly focused on the effect of

rotating speed on the final product quality by incorporating the
mechanistic fragmentation. Explicitly, five different rotating speed
has been used in this study (47 rpm, 95 rpm, 143 rpm, 190 rpm
and 238 rpm). It can be noted that the rotation speed used in this
study reflects the industrial processing condition. Often the indi-
vidual rotation speed is expressed in percentage of critical speed,
while the ball mill critical speed represents the speed where the
balls will not fall from its position in drum into the bottom of the
shell, as a result of centrifugal forces. In other words, the critical
speed for a ball milling was defined as the rotational speed where
centrifugal forces equal gravitational forces at the mill shell’s inside
surface.
Many previous studies (Bbosa et al., 2016; Brandao et al., 2020;
Hlosta et al., 2020a,b) demonstrate that with a variation of rotation
speed, three primary mechanisms of material motion (slipping;
cascading and cataracting) in the drum are usually detected which
can be classified based on Froude number:
Fig. 14. Typical torque behaviour with time during ball milling.
ω2 R
Fr = (13)
g
Therefore, fracturing and mixing mechanisms were induced in DEM
simulations presented in this paper. Where ω (rad/s) is drum rotation speed, R (m) is the radius of the
The initial stage of the DEM simulation is illustrated in Fig. 13(a). drum and g (kg/m3 ) is the gravitational acceleration. For a better
The separation of granules based on their size in simulation time understanding of the relationship between rotation speed and final
equal to zero (t = 0.2 s) can be easily recognised. Fig. 13(b) represents quality of the product, the velocity distribution of milling material
the propagation stage of the milling. It was observed that the level (Fig. 15), the individual trajectory of the balls in the drum, and evo-
of mixing increases with simulation time = 0.75 s, even though most lution of the power over time are presented (Figs. 16 and 17). These
of the granules were remaining as one unit, and where the number plots are giving detail about transverse bed motion of filled mate-
of breaking bonds between the particles only slightly increased. rials but also information about energy consumption and potential
With time = 2.25 s, the mixing level reaches steady-state condition, collision of the balls. The results regarding the number and size of
and from this time, the same pattern of granule mixing in the drum the fragments in the final product generated due to the variation of
was observed (Fig. 13c). Steady-state condition for the generation milling speed also allowed to identified effectiveness of the process
of new fragments due to the breaking of granules was not achieved (Figs. 18 and 19).
until time = 3.3 s. Fig. 13(d) presents the final product quality at the
end of the simulation (time = 5 s). At this stage, both mixing and Velocity distribution
fragmentation behaviour of original granules reached steady-state Low rotation speed (>47 rpm) was not considered in this study
conditions. as it was expected that the material would exhibit slipping motion
Adapted DBSCAN method was employed to identify the number which is not relevant for mixing or ball milling processing (Bbosa
and size of the fragments in milling material. Details concerning et al., 2016; Hlosta et al., 2020a,b). In contrary, milling with a rota-
identification were previously outlined in Section “Density-based tion speed of 47 rpm (19.7% of critical speed) as in Fig. 15(a), was
clustering”. Original batch with 2745 granules was milled into 2494 classified with a rolling motion, which is a subtype of cascading
primary particles and 7946 new fragments with size between 0.03 mode. The balls and milled material in this figure are coloured by
mm and 1.2 mm. Moreover, the torque of the rotating drum was their size. In general, the rolling motion belongs to the effective
calculated because this parameter can yield unbiased information regime for the purpose of the mixing. However, for a rolling regime
about the energy aspect of the process (Fig. 14). The torque imme- with a Froude number of 0.064 (−), mixing occurred insufficient
diately increased to 4 N.mm during filling of the granules into the (poor mixing), and uniform flow of the granules was registered
drum and oscillation of this parameter was negligible until all gran- mainly only on the free-flowing zone (Fig. 6b) located at the top
ules were settled and rearranged into the drum. After simulation of the bed (red colour in Fig. 15a). Mostly, the granules were sliding
time of 1 s, the drum started to rotate, and granules together with along the drum wall (blue colour in Fig. 15a), and the abrasion zone
steel balls started to move. From this time, the torque started to was very dense.
oscillate due to fragmentation and balls collisions. After a simula- The rolling-cascading motion was saturated, when milling
tion time of 4 s, the drum ceased to rotate and subsequently torque speed increased slightly from 47 rpm to 95 rpm (39.9% of criti-
oscillation drastically reduced until the end of the process. Practi- cal speed; Froude number is 0.2547) as it is shown in Fig. 15(b).
cally, the lower part of the drum was affected caused by collisions It was noted that this milling mode was also not adequate for ball
between fractured granules and balls, and the typical distribution milling, and it leads to little in the way of fracturing and mixing.
of contact pressure is shown in Fig. 15(b) at the simulation time of The free-flowing zone slightly increased, and the density of par-
2.25 s. ticle in abrasion zone decreased, which is reflected in the volume
expansion of free flow granules (Fig. 15b). In this case, the dead and
The effect of milling rotation speed on the fragmentation of impact zones are neglected.
granules Cascading motion in Fig. 15(c) with a Froude number of 0.5734
(−) was achieved for ball milling with a rotation speed of 143 rpm
In a production environment, optimum process condition (such (60% of critical speed). The mixing of granules started to be more
as filling level; properties of milling material and balls; the geome- evident. Lifters worked more effective as outlined in Fig. 6(a) and
try of the drum and balls or milling speed) needs to be considered enable redistribution of granules inside the drum. The higher free-
with regards to energy consumption and particle size distribution. flowing volume of materials in a drum space was indicated, as
Rotation speed should be taken into consideration to optimise the illustrated in Fig. 15(c).

162
A. Krok et al. Particuology 58 (2021) 153–168

Fig. 15. The velocity distribution of flowing granules in the drum during ball milling with varying rotation speed (a) 47 rpm; (b) 95 rpm; (c) 143 rpm; (d) 190 rpm and (e)
238 rpm.

By increasing the rotation speed from 143 rpm to 190 rpm, the lution of ball trajectory with time viewed along the mill’s axis with
transition from non-effective to effective mode was accomplished. varying rotational speeds. It should be noted it is a cross-sectional
The transition between these two aforementioned regimes is diffi- view across the X-Y direction. Varying colours in this figure indicate
cult to ascertain due to highly dynamic behaviour. With a Froude the location of the balls in a drum over time. The dark blue colour
number of 1.019 (−) and rotation speed of 190 rpm (79.8% of critical demonstrates the position of the balls at the beginning of milling,
speed), cataracting motion of material was achieved. This rotation and orange colour expresses the position of the balls at the end of
speed for the milling process can be indicated as the most effective the process. Trajectories indicate where the balls go, but it does not
speed in this study. This mode was observed to be optimal. Fig. 15(d) yield information regarding the nature of the collisions, which may
presents the intensive movement of material in the whole volume cause a break of granules. Previous studies demonstrate that there
of the drum, where lifters have a significant influence on the redis- is a direct correlation between the trajectories and collisions of the
tribution of material. The granules and fragments near the lifters balls (Cleary & Morrison, 2011; Jonsén et al., 2014).
move upwards at the peripheral mill velocity and then flow down The weakly developed trajectory of moving balls in the drum
at similar speed into the impact zone which leads to high shear was observed for rolling and rolling-cascade motion regimes with
rates and significant fines production due to abrasion, chipping, the rotations speed of 47 rpm and 95 rpm. The balls remain near
and attrition. The empty zone observed is minimal (occupied vol- the centre of charge circulation, and slow grinding occurs due to
ume is reduced significantly) relative to the empty zones observed low ball inertia. The ball trajectories were similar with results pre-
in Fig. 15(a–c). In addition, the dead zone is practically non-existent sented by Cleary and Morrison (2011), where the collisions of the
which was evident in Fig. 15(a–c). ball raised proportionally with rotation speed (Fig. 16a and b). How-
Conversely, when a mill is running at critical speed, the drum ever, there are no significant discrepancies between regime with Fr
is rotating too fast, and the speed is reaching centrifuging value. = 0.064 (−) and Fr = 0.2547 (−). In both cases, the lifters have an
This centrifuging motion is presented in Fig. 15(e) when the Froude only negligible impact on final mill product due to the insufficient
number was 1.5928 (−), and the rotation speed of 238 rpm reaches a ball collisions.
critical value. In this particular case, the radial segregation of mov- Although the angle of the ball’s trajectory does not change sig-
ing material was quite significant. The material close to the mill nificantly for the rotation speed of 143 rpm, the balls contributed to
shell consists mostly of fines, while large fragments and balls were collisions more effectively (Fig. 16c). The balls usually collide with
concentrated in the centre of the charge. The results supporting each other, but also with lifters, the drum and milling granules
the evidence of radial segregation were also observed by Powell in both vertical and horizontal directions. Especially for cascad-
and Nurick (1996). ing motion with Fr = 0.5734 (−), the intensity of travelling balls
in the horizontal direction was denser as seen in Fig. 16(c), which
Ball collisions demonstrates that there was significant mobility in this direction.
Several DEM approaches have been developed to analyse mill Therefore, a higher quantity of ball’s kinetic energy was transmitted
power in terms of ball collisions (Bbosa et al., 2016; Cleary, 2006; to the granules, which have a positive impact on fragmentation.
Mishra and Rajamani, 1992). In general, the accumulated sum of The trajectories in Fig. 16(d and e) indicated that the impact of
forces from the weight of the colliding balls on the mill geometry ball collisions was substantial for both cataracting and centrifugal
was determined to calculate power consumption. It was high- motion regimes with Fr = 1.019 (−) and Fr = 1.5928 (−) and the balls
lighted that proper identification of ball motion is imperative to travel a longer distance in the drum. This is followed by acceler-
provide precise power predictions. Fig. 16 demonstrates the evo- ated collisions. With each circulation of the drum, the balls and mill

163
A. Krok et al. Particuology 58 (2021) 153–168

Fig. 16. Representative trajectory of the balls in the drum during milling with varying rotation speed (a) 47 rpm; (b) 95 rpm; (c) 143 rpm; (d) 190 rpm and (e) 238 rpm.

granules were transported upwards. As the shoulder of the drum

was approached, the mill material together with balls fall under
gravity and collide with lifters, granules and neighbouring balls.
Transmitted energy from the balls to the granules also depends on
operational parameters. Although the collision behaviour of balls
for centrifugal motion (Fig. 16e) was very similar to cataracting
motion, the granules and balls were driven to the outside of the
mill due to centrifugal forces, and therefore less breakage events
occurred.

Power consumption
The power required to operate the milling process is closely
related to the torque and angular velocity of the rotating drum.
Particularly, reaction force (Fr ) from the centre of rotation drum
multiply by the radius of the drum gives torque (T) which together
with drum rotation speed define the mill power (P) as follows:
Fig. 17. The mill power evolution overtime during the ball milling process for the
following speeds (a) 47 rpm; (b) 95 rpm; (c) 143 rpm; (d) 190 rpm and (e) 238 rpm.
P = Fr Rω (14)

In general, commercial ABAQUS software generate robust

databased of output data for further post-processing at the end of of the drum was calculated. As a result, oscillation and the integral
each simulation. Hence, a parameter such as reaction force from value of mill power over time were possible to investigate, as shown
a structural reference point in the centre of the drum from this in Fig. 17. Note that the calculated mill power overtime is not relat-
databased was extracted. Subsequently, the torque in the centre ing to the mechanical losses in the drive. For all rotation speeds, the

164
A. Krok et al. Particuology 58 (2021) 153–168

Fig. 19. Number and size fragments variation in the final product over different
milling speed.

tively, and for 190 rpm approximately 21760 new fragments were
generated every step time until the end of the process (t = 5 s). For
Fig. 18. Generation of fragments during the milling process with time for five dif- 238 rpm due to the undesirable consequence of centrifugal forces,
ferent rotation speeds. the production of new fragments decreased to 20010, even rotation
speed increased.
power significantly increased between 0–0.1 s, which is referring to Developed DEM models facilitate the prediction of the size dis-
the feeding of granules and balls into the drum. The material started tribution of new fragments in the final product with a variation of
to rotate from 1 s of simulation time. From this point, mill power rotation speed. The relationship between the size and number of
oscillated until the drum stopped to move (4 s). After some rapid fragments is shown in Fig. 19. In general, it can be concluded that
fluctuations upon the start-up of the mill motion, the mill pow- sufficient mechanical energy through the collisions during milling
der reached its equilibrium value, and the power was essentially was generated and most of the granules reduced their size due to
steady. Results indicate that with an increase of rotation speed, the fracturing. The simulations show that granules experience differ-
mill power also increased, which agrees with earlier experimental ent impact numbers due to the variation of rotation milling speed.
observations (McIvor, 1983; Sastri & Rao, 1997). It was determined It was found that ball milling affected the final state of the product
that for 47 rpm, the power was around 3.17 W, for 95 rpm it was 6.3 in terms of their finesses and complex structure of consisting frag-
W and for 143 rpm, 190 rpm and 238 rpm the mill power oscillated ments. The initial small number of large granules, being subjected
around 9.82 W, 12.63 W and 16.91 W. Besides, amplitude and fre- to increasing mechanical energy, broke to become more numer-
quency of oscillation with rotation speed rise proportionally, which ous and smaller in size. The number of primary particles number
is in line with others (Cleary et al., 2018; Govender et al., 2015; increased appreciably. The size distribution of fragments in final
Hlungwani et al., 2003). Moderated peaks over time were related products was significantly modified by rotation speed.
to lifting and collapsing of the material inside of the drum during In contrast, the energy level required to fracture certain size
milling. of granules may reach a level that is not provided by the miller.
Therefore, the granules remain unchanged during further milling.
Fragmentation of granules As milling proceeds material breaks down into a few large particles
The bonded granules of different shape and size were assem- and several fine particles with relatively few of intermediate size.
bled by primary particles as outlined in Section “DEM granule The size of the larger particles is related to the size reduction pro-
shape characterisation”. While each bond was modelled as a set of cess, and the size of the finer particles is dependent on the structure
multi-point connector’s elastic-springs distributed over a centre of of the material. About 0.62% and 0.64%, of the original granules, exit
a particle that transfers both translation and rotation motion. Each the mill without impacting, for the rotation speed of 47 rpm and
bonded particle remains independent and resisted tension pulling 95 rpm. For 143 rpm, 190 rpm and 238 rpm it was 0.53%, 0.27% and
until it reaches a cut-off value of breakage force due to impact. 0.31%. The simulations also have shown that the granules also expe-
Moreover, complex constitutive equations were not required to rience a large number of impacts. The quantity of fines in product
describe post-failure flow as the broken granules followed the stan- with rotation speed increased between 43–190 rpm. For 43 rpm,
dard DEM interaction rules. Number and size of the fragments 95 rpm, 143 rpm and 190 rpm it was 26.10%, 34.58%, 39.87% and
produced from original overlapping particles in granules during 49.34% respectively. Due to the centrifugal motion of mill mate-
the milling process were effectively identified by using DBSCAN rial with the rotation speed of 238 rpm, only 44.52% of fines were
method (more detail is stated in Section “Density-based cluster- generated.
ing”).
Fig. 18 illustrates the generation of new fragments in the drum Conclusion
due to milling of granules over simulation time of 5 s, for five dif-
ferent rotation speeds. The original batch of 2745 granules started A systematic study on DEM modelling of ball milling process
to be subject of fragmentation due to different nature of collisions was covered in this study. The aim was to reduce the original size
between the balls, mill material, the drum, and lifters. Gradual frac- of pharmaceutical granules to much smaller quantities. Granular
turing was found with the pharmaceutical excipient. Generation of lactose was experimentally characterised. A combination of sta-
new fragments for all five different rotation speeds was progres- tistical locally weighted regression smoothing approach, together
sively increased until 3.5 s. After that, a steady-state was reached, with a K-fold cross-validation method, was employed to iden-
and a stable number of fragments were generated. From 3.2 s of tify the sphericity and roundness of the granules. Based on the
simulation time, the milling process with the rotation speed of 47 experimental data, three types of the granules with elongated,
rpm, produced 6395 new fragments; for 95 rpm and 143 rpm of blocky and rounded shape were selected for further investigation
rotation speed, it was 10440 and 14420 new fragments respec- of the ball milling process. Commercial software ABAQUS was used

165
A. Krok et al. Particuology 58 (2021) 153–168

to numerically simulate the ball milling process, while standard Austin, L. G., Julianelli, K., Sampaio de Souza, A., & Schneider, C. L. (2007).
Hertz-Mindlin contact model rules the interaction between parti- Simulation wet ball milling of iron ore at Carajas, Brazil. International Journal of
Mineral Processing, 84, 157–171.
cles. The granules with three different non-spherical shapes and Austin, L. G., Klimpel, R. R., & Luckie, P. T. (1984). Process engineering of size
three different sizes were defined as a grain of complex shapes by reduction: Ball milling. New York: SME–AIME.
bonding spherical particles together using compliant connectors. Bareither, C. A., Edil, T. B., Benson, C. H., & Mickelson, D. M. (2008). Geological and
physical factors affecting the friction angle of compacted sands. Journal of
Overlapping particles were connected by beam-type multi-point Geotechnical and Geoenvironmental Engineering, 134, 1476–1489.
constraints instead of adhesion integration. It was considered in Bbosa, L. S., Govender, I., & Mainza, A. (2016). Development of a novel
this study that connected particles into the granules can be dam- methodology to determine mill power draw. International Journal of Mineral
Processing, 149, 94–103.
aged by impact force or collisions into smaller fragments or primary
Bian, X., Wang, G., Wang, H., Wang, S., & Lv, W. (2017). Effect of lifters and mill
particles. Hence, the systematic trial-error approach was applied speed on particle behaviour, torque and power consumption of a tumbling ball
to estimate critical breakage force for all DEM granules. It was mill: Experimental study and DEM simulation. Minerals Engineering, 105,
22–35.
explored, that size of the granules significantly affects the value
Bilgili, E., & Scarlett, B. (2005). Population balance modeling of non-linear effects in
of breakage force compare to the granule shape or rearrangement milling processes. Powder Technology, 153, 59–71.
of the particles when only slight changes have been observed. A Boemer, D., & Ponthot, J. P. (2017). DEM modeling of ball mills with experimental
strong correlation between the impact speed and a number of frag- validation: Influence of contact parameters on charge motion and power draw.
Computational Particle Mechanics, 4, 53–67.
ments has been found, and results seem to agree with previous Bonakdar, T., & Ghadiri, M. (2018). Analysis of pin milling of pharmaceutical
papers. DBSCAN algorithm was adapted to perform the identifica- materials. International Journal of Pharmaceutics, 552, 394–400.
tion of the clustering fragments effectively. This paper was mainly Bond, F. C. (1958). Grinding ball size selection. Mining Engineering, 10, 592–595.
Bond, F. (1961). Crushing and grinding calculation part I. British Chemical
focused on the effect of rotating speed on the final quality of the Engineering, 378–548.
product incorporating the mechanistic fragmentation. It was deter- Bowman, E. T., Soga, K., & Drummond, W. (2001). Particle shape characterisation
mined based on Froud number that material reaches rolling motion using Fourier descriptor analysis. Géotechnique, 6, 545–554.
Bozkurt, V., & Özgür, I. (2007). Dry grinding kinetics of colemanite. Powder
for 47 rpm, rolling-cascading motion for 95 rpm, cascading motion Technology, 176, 88–92.
for 143 rpm, cataracting motion for 190 rpm and critical centrifu- Brandao, R. J., Lima, R. M., Santos, R. L., Duarte, C. R., & Barrozo, M. A. S. (2020).
gal motion for 238 rpm in the drum during milling. For the rolling Experimental study and DEM analysis of granular segregation in a rotating
drum. Powder Technology, 364, 1–12.
regime, insufficient mixing occurred, and uniform flow of the gran-
Burmelster, C. F., & Kwade, A. (2003). Process engineering with planetary ball mills.
ules was registered mainly only on free-flowing layers located at Chemical Society Reviews, 42, 7660–7667.
the top of the bed. Neither rolling-cascading mode was not ade- Cabalar, A. F., Dulundu, K., & Tuncay, K. (2013). Strength of various sands in triaxial
and cyclic direct shear tests. Engineering Geology, 156, 92–102.
quate for ball milling, which leads to not effective fracturing and
Carmen, C. P., Fernandez-Prieto, S., & De Borggraeve, W. M. (2019). Ball milling: A
mixing. The zone with active layers slightly increased with cascad- green technology for the preparation and functionalisation of nanocellulose
ing motion, which reflected in the volume expansion of free flow derivative. Nanoscale Advances, Royal Society of Chemistry, 1, 937–947.
granules and the mixing of granules started to be more profound. Castellanos, A. (2005). The relationship between attractive interparticle forces and
bulk behaviour in dry and uncharged fine powders. Advance in Physics, Taylor &
Moreover, the lifters work more effectively and allow redistribution Francis, 54, 263–376.
of granules inside of the drum space. Rotation speed with 190 rpm Celik, M. S. (1988). Acceleration of breakage rates of anthracite during grinding in a
indicated to be as the most effective speed in this study. This mode ball mill. Powder Technology, 54, 227–233.
Chapuis, R. P. (2012). Estimating the in situ porosity of sandy soils sampled in
was observed to be optimal, where the most size reduction was boreholes. Engineering Geology, 19, 57–64.
taking place. In contrary, with increasing speed above this value, Charles, R. J. (1957). Energy–size reduction relationships in comminution. Trans.
a centrifugal motion was achieved, which have undesirable conse- AIME, 208, 80–88.
Chen, H., Xiao, Y., Liu, Y., & Shi, Y. (2017). Effect of Young’s modulus on DEM results
quences on the production of new fragments. Effectiveness of the regarding transverse mixing of particles within a rotating drum. Powder
rotation speed on the process directly correlated with the colli- Technology, 318, 507–517.
sions of the ball. It was discovered that with increasing the rotation Cho, J., & Sohn, H. Y. (2016). Effect of particle shape and size distribution on the
overall fluid-solid reaction rates of particle assemblages. The Canadian Journal
speed up to 190 rpm, the collisions of the ball raised proportion-
of Chemical Engineering, 94, 1516–1523.
ally, which demonstrate the development of the trajectory of the Cho, G. C., Dodds, J., & Santamarina, C. (2006). Particle shape effects on packing
balls. However, the calculated power input increased with rotation density, stiffness and strength: Natural and crushed sands. Journal of
Geotechnical and Geoenvironmental Engineering, 132, 591–602.
speed.
Cleary, P. W. (2006). Axial transport in dry ball mills. Applied Mathematical
Modelling, 30, 1343–1355.
Cleary, P. W., & Morrison, R. D. (2011). Understanding fine ore breakage in a
Declaration of Competing Interest laboratory scale ball mill using DEM. Mineral Engineering, 24, 352–366.
Cleary, P. W., Morrison, R. D., & Delaney, G. W. (2018). Incremental damage and
The authors report no declarations of interest. particle size reduction in a pilot SAG mill: DEM breakage method extension
and validation. Minerals Engineering, 128, 56–68.
Cundall, P. A., & Strack, O. D. L. (1979). A discrete numerical model for granular
assemblies. Géotechnique, 1, 47–65.
Acknowledgements Cutt, T., Fell, J. T., Rue, P. J., & Spring, M. S. (1987). Granulation and compaction of a
model system. II. Stress relaxation. International Journal of Pharmaceutics, 39,
This work was supported by the Career-FIT Fellowships, funded 157–161.
Datta, A., & Rajamani, R. K. (2002). A direct approach of modelling batch grinding in
through European Union’s Horizon 2020 research and innovation ball mills using population balance principles and impact energy distribution.
programme under the Marie Skłodowska-Curie grant agreement International Journal of Mineral Processing, 64, 181–200.
No. 713654. Part of the work was also supported by ACCORD (ITMS Deniz, V. (2013). Effects of mill speed on kinetic breakage parameters of four
different particulate pumices. Particulate Science and Technology: An
project code:313021X329), funded through the European Regional
InternationalJournal, 31(2), 101–108. https://doi.org/10.1080/02726351.2012.
Development Fund. 658903
Deniz, V., & Onur, T. (2002). Investigation of the breakage kinetics of pumice
samples as dependent on powder filling in a ball mill. International Journal of
References Mineral Processing, 67, 71–78.
Fu, J., Adams, M. J., Reynolds, G. K., Salman, A. D., & Hounslow, M. J. (2004). Impact
Altuhafi, F., O’Sullivan, C., & Cavarretta, I. (2013). Analysis of animage-based deformation and rebound of wet granules. Powder Technology, 140, 248–257.
method to quantify the size and shape of sand particles. Journal of Geotechnical Fuerstenau, D. W., & Abouzeid, A. Z. M. (2002). The energy efficiency of ball milling
and Geoenvironmental Engineering, 8, 1290–1307. ASCE139. in comminution. International Journal of Mineral Processing, 67, 161–185.
Arasan, S., Akbulut, S., & Hasiloglu, A. S. (2011). The relationship between the Ghaboussi, J., & Barbosa, R. (1990). Three-dimensional discrete element method for
fractal dimension and shape properties of particles. Journal of Civil Engineering, granular materials. Numerical and Analytical Methods in Geomechanis, 14,
15, 1219–1225. 451–472.

166
A. Krok et al. Particuology 58 (2021) 153–168

Goh, H. P., Heng, P. W. S., & Liew, C. V. (2018). Comparative evaluation of powder Mitchell, J. K., & Soga, K. (2005). Fundamentals of soil behavior (3rd edition). New
flow parameters with reference to particle size and shape. International Journal York, USA: Wiley.
of Pharmaceutics, 547, 133–141. Moreno, R., Ghadiri, M., & Antony, S. J. (2003). Effect of the impact angle on the
Gong, X., Chang, S. Y., Yeboah, F. O., Paul, S., Perumalla, S. R., Shi, L., Sun, W. J., Zhou, breakage of agglomerates: A numerical study using DEM. Powder Technology,
Q., & Sun, C. C. (2015). Dependence of tablet brittleness on tensile strength and 130, 132–137.
porosity. International Journal of Pharmaceutics, 493, 208–213. Morrison, R., & Cleary, P. W. (2004). Using DEM to model ore breakage within a
Govender, N., Rajamani, R. K., Kok, S., & Wilke, D. N. (2015). Discrete element pilot scale sag mill. Mineral Engineering, 17, 1117–1124.
simulation of mill charge in 3D using the BLAZE-DEM GPU framework. Mosharraf, M., & Nystrom, C. (1995). The effect of particle size and shape on the
Minerals Engineering, 79, 152–168. surface specific dissolution rate of microsized practically insoluble drugs.
Guenette, E., Barrett, A., Kraus, D., Brody, R., Harding, L., & Magee, G. (2009). International Journal of Pharmaceutics, 22, 35–47.
Understanding the effect of lactose particle size on the properties of DPI Nanda, S. J., & Panda, G. (2015). Design of computationally efficient density-based
formulations using experimental design. International Journal of Pharmaceutics, clustering algorithms. Data & Knowledge Engineering, 95, 23–38.
380(1–2), 80–88. Napier-Munn, T. J., Morrell, S., Morrison, R. D., & Kojovic, T. (1999). Mineral
Hacifazlioglu, H., & Korkmasz, A. V. (2020). Performance comparison of stirred comminution circuits their operation and op-timisation. JKMRC monograph series
media mill and ball (BOND) mill in bauxite grinding. Particulate Science and in mining and mineral processing. Brisbane: University of Queensland (2).
Technology, 38, 404–408. Ning, Z., Boerefiin, R., Ghadiri, M., & Thornton, C. (1997). Effect of particle size and
Herbst, J. A., & Fuerstenau, D. W. (1980). Scale-up procedure for continuous bond strength on impact breakage of weak agglomerates. Powders & Grains,
grinding mill design using population balance models. International Journal of Balkema, Rotterdam, 97.
Mineral Processing, 7, 1–31. Nisoli, A., Doherty, M. F., & Malone, M. F. (2004). Feasible regions for step-growth
Hertz, H. (1942). On the contact of elastic solids. Journal Reine Und Angewandte melt poly-condensation systems. Industrial & Engineering Chemistry Research,
Mathematik, 92, 156–171. 43(2), 428–440.
Hlosta, J., Jezerska, L., Rozbroj, J., Zurovec, D., Necas, J., & Zegzulka, J. (2020a). DEM Nomura, S., Hosoda, K., & Tanaka, T. (1991). An analysis of the selection function
investigation of the influence of particulate properties and operating for mills using balls as grinding media. Powder Technology, 68, 1–12.
conditions on the mixing process in rotary drums: Part 2 - process validation Orozco, L. F., Nquyen, D. H., Delenne, J. Y., Sornay, P., & Radjai, F. (2019).
and experimental study. Processes, 8, 184. Discrete-element simulation of comminution in rotating drums: Effects of
Hlosta, J., Jezerska, L., Rozbroj, J., Zurovec, D., Necas, J., & Zegzulka, J. (2020b). DEM grinding media. Powder Technology, 362, 157–167.
investigation of the influence of particulate properties and operating Perkins, M., Ebbens, S. J., Hayes, S., Roberts, C. J., Madden, C. E., Luk, S. Y., & Patel, N.
conditions on the mixing process in rotary drums: Part 1 – determination of (2007). Elastic modulus measurements from individual lactose particles using
the dem parameters and calibration process. Processes, 8(2), 222. atomic force microscopy. International Journal of Pharamceutics, 332, 168–
Hlungwani, O., Rikhotso, J., Dong, H., & Moys, M. H. (2003). Further validation of 175.
DEM modeling of milling: Effects of liner profile and mill speed. Powder Potapov, A. V., & Campbell, C. S. (2001). Parametric dependence of particle
Technology, 16, 993–998. breakage mechanisms. Powder Technology, 120(3), 164–174.
Huang, W., Shi, Y., Kongtong, C. W., Sun, F. Y., & Li, Y. (2013). Using spray-dried Powell, M., & Nurick, G. (1996). A study of charge motion in Rotary Mills. Part 3 –
lactose monohydrate in wet granulation method for a low-dose oral analysis of results. Mineral Engineering, 9, 399–418.
formulation of a paliperidone derivative. Powder Technology, 246, 379–394. Rajamani, R. K., & Guo, D. (1992). Acceleration and deceleration of breakage rates
Jonsén, P., Palsson, B. I., Stener, J. F., & Haggblad, H. A. (2014). A novel method for in wet ball mills. International Journal of Mineral Processing, 34, 103–118.
modelling of interactions between pulp, charge and mill structure in tubling Refahi, A., Aghazadeh Mohandesi, J., & Rezai, B. (2010). Discrete element modeling
mills. Minerals Engineering, 63, 65–72. for predicting breakage behavior and fracture energy of a single particle in a
Kafui, K. D., & Thornton, C. (2000). Numerical simulations of impact breakage of a jaw crusher. International Journal of Mineral Processing, 94(1–2), 83–91.
spherical crystalline agglomerate. Powder Technology, 109, 113–132. Reynolds, G. K., Fu, J. S., Cheong, Y. S., Hounslow, M. J., & Salman, A. D. (2005).
Kenderll, K. (1988). Agglomerate strength. Powder Metallurgy, 31, 28–31. Breakage in granulation: A review. Powder Technology, 60, 3969–
Krok, A., & Wu, C. Y. (2019). Evolutions of temperature and density during roll 3992.
compaction of a pharmaceutical excipient. International Journal of Rittinger, P. R. (1867). Lehrbuch der Aufbereitungskunde. Berlin: Ernst and Korn.
Pharmaceutics, 572, Article 118822. Rodriguez, J. M., Johansson, J. M. A., & Edeskar, T. (2012). Particle shape
Krok, A., Mirtic, A., Reynolds, G. K., Schiano, S., Roberts, R., & Wu, C. Y. (2016). An determination by two-dimensional image analysis in geotechnical
experimental investigation of temperature rise during compaction of engineering. Proceedings of Nordic Conference on Soil Mechanics and Geotechnics,
pharmaceutical powders. International Journal of Pharmaceutics, 513, 97–108. Copenhagen, Denmark: Danish Geotechnical Society.
Krok, A., Peciar, M., & Fekete, R. (2014). Numerical investigation into the influence Rumpf, H. (1962). The strength of granules and agglomerates. In W. A. Knepper
of the punch shape on the mechanical behaviour of pharmaceutical powders (Ed.), Agglomeration (pp. 379–418). New York: John Wiley and sons.
during compaction. Particuology, 16, 116–131. Salman, A. D., & Gorham, D. A. (1997). The fracture of glass spheres under impact
Krumbein, W. C. (1941). Measurement and geological significance of shape and loading. Powders & Grains, 97.
roundness of sedimentary particles. Journal of Sedimentary Research, 11, 64–72. Samimi, A., Moreno, R., & Ghadiri, M. (2004). Analysis of impact damage of
Krumbein, W. C., & Sloss, L. L. (1951). Stratigraphy and Sedimentation. Science, agglomerates: Effect of impact angle. Powder Technology, 143–144, 97–109.
114(2961), 335–336. https://doi.org/10.1126/science.114.2961.335 Sander, J., Ester, M., Kriegel, H. E., & Xu, X. (1998). Density-Based Clustering in
Kumar, A., Vercruysse, J., Bellandi, G., Germaey, K. V., Vervaet, C., Remon, J. P., Beer, Spatial Databases: The Algorithm GDBSCAN and Its Applications. Data Minig
T. D., & Nopends, I. (2014). Experimental investigation of granule size and and Knowledge Discovery, 2, 169–194.
shape dynamics in twin-screw granulation. International Journal of Santamarina, J. C., & Cho, G. C. (2004). Soil behaviour: the role of particle shape. In
Pharmaceutics, 475, 485–495. In advances in geotechnical engineering, The Skempton conference. pp. 604–617.
Lian, G., Thornton, C., & Adams, M. J. (1998). Discrete particle simulation of London, UK: Thomas Telford.
agglomerate impact coalescence. Chemical Engineering Science, 53, 3381–3391. Sastri, S. R. S., Rao, K. K. Simple empirical equations for prediction of mill power
Liu, L., Kafui, K. D., & Thornton, C. (2010). Impact breakage of spherical, cuboidal draw. Trans. Instn Min. Metall. (Sect. C: MineralProcess. Extr. Metall.), 106
and cylindrical agglomerates. Powder Technology, 199(2), 189–196. (1997) C91– 4.
Loh, Z. H., Sarmanta, A. K., & Heng, P. W. S. (2015). Overview of milling techniques Saxena, A., Prasad, M., Gupta, A., Bharill, N., Patel, O. P., Er, A. T., Weiping, D., &
for improving the solubility of poorly water-soluble drugs. Asian Journal of Chin-Teng, L. (2017). A review of clustering techniques and developments.
Pharmaceutical Sciences, 10, 255–274. Neurocomputing, 267, 664–681.
Mangwandi, C., Adams, M. J., Hounslow, M. J., & Salman, A. D. (2010). Effect of Shin, H., & Santamarina, J. C. (2013). The role of particle angularity on the
impeller speed on mechanical and dissolution properties of high-shear mechanical behaviour of granular mixtures. Journal of Geotechnical and
granules. Chemical Engineering Journal, 164(2–3), 305–315. Geoenviromental Engineering, 139, 353–355.
Markauskas, D., Ramírez-Gómez, Á., Kacianauskas, R., & Zdancevicius, E. (2015). Shoji, K., Austin, L. G., Smaila, F., Brame, K., & Luckie, P. T. (1982). Powder
Maize grain shape approaches for DEM modelling. Computers and Electronics in Technology, 3, 121.
Agriculture, 118, 247–258. Stamboliadis, E. Th. (2002). A contribution to the relationship of energy and
McIvor, R. E. (1983). Effects of speed and liner configuration on ball mill particle size in the comminution of brittle particulate materials. Minerals
performance. Mineral Engineering, 617–620. Engineering, 15(10), 707–713. https://doi.org/10.1016/S0892-6875(02)00185-1
Metzger, M. J., & Glasser, B. J. (2012). Numerical investigation of the greakage of Sukumaran, B., & Ashmawy, A. K. (2001). Quantitative characterisation of
bonded agglomerates during impact. Powder Technology, 217, 304–314. geometry of discrete particles. Géotechnique, 51, 171–179.
Mindlin, R. D. (1949). Compliance of elastic bodies in contact. Journal of Applied Tavares, L. M. (2009). Analysis of particle fracture by repeated stressing as damage
Mechanics, 16, 259–268. accumulation. Powder Technology, 190, 327–339.
Mindlin, RD, Deresiewicz, H. Elastic spheres in contact under varying oblique Tavares, L. M., & King, R. P. (1998). Single-particle fracture under impact loading.
forces ASME, September 1953: 327-344. International Journal of Mineral Processing, 54(1), 1–28.
Miquélez-Morán, A. M., Wu, C. Y., Dong, H., & Seville, P. K. (2009). Characterisation Thornton, C., Ciomocos, M. T., Yin, K. K., & Adams, M. J. (1997). Fracture of
of density distributions in roller-compacted ribbons using micro-indentation particulate solids. Powders & Grains, 97. Balkema, Rotterdam.
and X-ray micro-computed tomography. Eur. J. Pharm. Biopharm., 72, 173–182. Toker, O. S., Zorlucan, F. T., Konar, N., Daglioglu, O., Sagdic, O., & Sener, D. (2016).
Mishra, B. K., & Rajamani, R. K. (1992). The discrete element method for the Investigating the effect of production process of ball mill refiner on some
simulation of ball mills. Applied Mathematical Modelling, 16, 598–604. physical quality parameters of compound chocolate: Response surface
Mishra, B. K., & Thornton, C. (2001). Impact breakage of particle agglomerates. methodology approach. International Journal of Food Science & Technology, 52,
International Journal of Mineral Processing, 61(4), 225–239. 788–799.

167
A. Krok et al. Particuology 58 (2021) 153–168

Tutumluer, E., & Pan, T. (2008). Aggregate morphology effecting strength and Wettimuny, R., & Penumadu, D. (2004). Application of Fourier analysis to digital
permanent deformation behaviour of unbound aggregate materials. Journal of imaging for particle shape analysis. Journal of Computational Civil Engineering,
Material Civil Engineering, 20, 617–627. 18, 2–9.
Tuzcu, T., & Rajamani, R. K. (2011). Modeling breakage rates in mills with impact Wu, C. Y., Cock, A. C. F., Gilia, O. T., & Thompson, D. A. (2003). Experimental and
energy spectra and ultra fast load cell data. Mineral Engineering, 24, 252–260. numerical investigations of powder transfer. Powder Technology, 138, 216–228.
Verheezen, J. J. A. M., van der Voort Maarschalk, K., Faassen, F., & Vromans, H. Wu, C. Y., Thornton, C., & Li, L. Y. (2008). Rebound behaviour of spheres during
(2004). Milling of agglomerates in an impact mill. International Journal of elastic-plastic oblique impact. International Journal of Modern Physics B, 22,
Pharmaceutics, 278, 165–172. 1095–1102.
Vogel, L., & Peukert, W. (2004). Determination of material properties relevant to Xu, D., & Tian, Y. (2015). A comprehensive survey of clustering algorithms. Annals
grinding by particable lascale milling tests. International Journal of Minerals of Data Science, 2, 165–193.
Process, 74S, S329–S338. Xu, L., Luo, K., Zhao, Y., Fan, J., & Cen, K. (2019). Influence of particle shape on liner
Wadell, H. (1933). Sphericity and roundness of rock particles. The Journal of wear in tubling mills: A DEM study. Powder Technology, 350, 26–35.
Geology, 41, 310–331. Zhang, L., Jia, Y., He, H., Yin, J., Chen, R., Zhang, C., Shen, W., & Wang, X. (2018).
Wadell, H. (1935). Volume shape, and roundness of quartz particles. The Journal of Multiple factor analysis on preparation of cellulose nanofiber by ball milling
Geology, 43, 250–280. from softwood pulp. BioResources, 13(2), 2397–2410.
Wang, L., Wang, X., Mohammad, L., & Abadie, C. (2005). A unified method to Zhang, L., Tsuzuki, T., & Wang, X. (2015). Preparation of cellulose nanofiber from
quantify affregate shape angularity and texture using Fourier analysis. Journal softwood pulp by ball milling. Cellulose, 22(3), 1729–1741.
of Material Civil Engineering, 17, 498–504. Zheng, J., & Hryciw, R. D. (2015). Traditional soil particle sphericity, roundness and
Wang, X., Yu, J., Lv, F., Wang, Y., & Fu, H. (2006). A multi-sphere based modelling surface roughness by computational geometry. Géotechnique, 6, 494–506.
method for maise grain assemblies. Advanced Powder Technology.
Weeber, A. W., & Bakker, H. (1988). Amorphization by ball milling. A review.
Physica B: Condensed Matter, 153, 93–135.

168

View publication stats