Verification of Polyhedral DEM with Laboratory Grinding Mill Experiments

Abstract

The simulation of grinding mills with the discrete element method (DEM) has been advancing. First, it emerged as a method for studying charge motion with spherical balls and predicting the power draw of the mill. Subsequently, studies on liner wear, charge motion with ellipsoidal and polyhedral shaped particles simulated with three dimensional DEM followed. Further, the impact energy spectra computed in the DEM algorithm is leading to the development of models for the breakage of brittle particles in mills. The core elements in such simulations are the shape of particles in the mill charge and the power draw of the mill due to operating variables. To advance the field, we present a set of experimental data and the corresponding DEM validation results for a 90 × 13 cm mill. The DEM algorithm uses the volume-overlap method which is more realistic for multifaceted irregular particle collisions. Further, we use the scanned shape of the rock media and multifaceted spherical shape for the grinding media to represent as close as possible the actual charge in the mill. First, we present DEM validation for spherical grinding media-only experiments, rock-only experiments, and a mixture of spherical grinding media and rocks, as well as aluminum cubes only to represent the theme of particle shape. Finally, a discussion of the contact mechanics parameters in the four modes of experiments is given. Since the feed ore to plant scale mills can vary in shape, mill simulations with scanned shape of typical particles are the future for more accurate results.

1. Introduction

In hard rock mining and metal production, the grinding unit operation follows stages of crushing mined ore. Naturally, occurring ores contain metal-bearing minerals. The purpose of grinding is to liberate metal-bearing mineral grains from the host particles. Next, the valuable mineral particles are recovered in subsequent processes. Since the grade of metal values in ores is in the range of one percent or less, thousands of tons of ore are milled every hour. As a result, the grinding mills are massive in size requiring a few megawatts of power to operate. Being such an energy consuming operation, grinding mills have been studied by researchers around the world to improve energy efficiency. Since the mills are massive in size, the tumbling environment inside the mill is very harsh, defying any attempt at quantifying the charge dynamics with sensors. However, laboratory-scale experiments have been done tactfully to provide information on power draw and charge profiles of material inside the mill. One of the critical needs in milling is to understand charge flow inside the cylindrical mill shell. Since the internal forces produced inside the charge mass are the cause for particle size reduction, this study has fallen under the realm of the numerical methods known as the Discrete Element Method (DEM) (Tavares et al., 2020). It provides a direct numerical simulation of the motion of thousands of particles tumbling inside the cylindrical mill shell, and besides, the forces acting on the particles can be used to model particle breakage. Djodjeric et al. (2006) studied the competency of large pieces of rocks in an AG mill with DEM. They presented charge pressure due to normal forces acting on the bed of rock to account for the breakage or non-breakage of rocks.

1.1 Grinding mills

Almost every mining plant operation employs either a Semi-autogenous (SAG) or Autogenous (AG) mill, followed by ball mills. Since SAG and AG mills are capable of performing the size reduction of one or two-stage crushing size reduction circuits, they are found extensively in many mineral ore processing operations. The product of SAG or AG mills is represented by size, called “transfer size”, which is in the range of few millimeters. Material contained in these sizes is suitable for further size reduction in ball mills to a size of 100–200 μm. Typical SAG mills are 9–12 m in diameter. The top size of feed ore to SAG/AG mills is around 10–35 cm. These mills usually operate with a screen deck at the discharge end to recycle the coarse portion of the material back into the mill.

AG mills take advantage of the weaker competence of the ore to break them. The ore particles impact directly against the mill shell and break readily. This mode of breakage is called self-breakage, and only a few ore deposits around the world are amenable for breakage in AG mills. In SAG mills, in addition to self-breakage, the grinding balls due to their greater density add to ball-to-rock impact breakage. Besides these two major modes of impacts, the grinding action includes attrition grinding or chipping of sharp edges. The entire system of grinding modes is influenced by design factors such as liners, lifters, grate opening, operating conditions (mill filling and speed), and ore hardness.

1.2 DEM simulation of grinding mills

It is apparent from the foregoing that real-world particles are irregular in shape. However, early numerical DEM studies began with spherical shapes, partly because they are the easiest to model collisions numerically (Mishra and Rajamani, 1990; 1992; Cleary, 2001; Cleary et al., 2003; Morrison and Cleary, 2004). Although, researchers in the DEM field recognized the importance of particle shape, yet the cost of computational time-restricted them to basic spherical shape and its close variants. The earliest of these were shapes generated with superquadrics algebra. During milling, particles with high angularity lose their sharp edge readily mainly due to chipping and abrasion, not to mention these sharp edges are the weakest in the modes of breakage. Hence, superquadrics in a preliminary sense were well adapted to simulate such rounded particles. Rothenburg and Bathurst (1991) presented a DEM algorithm for contact detection and collision of elliptical particles in the simulation of micromechanics of granular assemblies. Hopkins et al. (1991) used polygonal shape particles to create a two-dimensional of rocks and sea ice blocks to simulate the formation of multilayers floes. Potapov and Campbell (1996) created a three-dimensional technique to simulate the brittle solid fracture by “gluing” together convex polyhedral particles. Realizing the need for non-spherical shapes in DEM, Latham et al. (2008) published a library of three-dimensional shapes. Cleary (1998) was the first author to introduce superquadrics for the DEM simulation of grinding mills. This author (Cleary, 1998) described the simulation of a 5-m diameter ball mill in detail including particle segregation, collisional force distributions, and liner wear. Delaney et al. (2010) went a step further using superquadrics in the simulation of particle breakage in comminution processes: a set of smaller particles were packed inside a superellipsoid, and the superellipsoid was made to break upon experiencing a certain energy limit releasing the smaller particles representing the breakage of the larger ellipsoid. Besides the spherical grinding media in ball mills, a shape termed “cylpebs” was developed to produce improved media packing in the tumbling mill charge. Kiangi et al. (2013) used DEM simulation spherical and multi-faceted polyhedron shapes to predict the load and power draw in a pilot-scale mill filled with spheres and cylpebs. Xu et al. (2019) using the superellipsoid particle model studied the effect of spherical and cuboidal particles on the wear of a mill. Realistic DEM simulations of plant size grinding mills require particle filling in the millions. Govender et al. (2015) showed how graphic processing unit (GPU) computation can speed up 3D DEM simulation of mill charge. The framework for such high-performance computing is described in Govender et al. (2016). This was followed by the introduction of irregular shaped particles such as scaled hexagonal prisms, truncated tetrahedra, cubes, and other shapes in GPU-based DEM computing. Govender et al. (2018a) presented the grinding mill simulation of regular polyhedral particles such as cube, tetrahedra, hexagonal prism, Biluna, and mixed polyhedral particles. This study investigated the role of particle angularity and aspect ratio on power draw. In particular, the normal and shear energy dissipation components on the power draw were shown to differ depending on the shape of the particle. Particle shape and angularity contribute to interlocking and shear slippage between particles (Govender et al., 2018b) for the study of particle mixing. The next logical choice for irregular particles was a cluster of spheres or clumped spheres. Like superquadrics, this object too was amenable to easier collision detection (Latham et al., 2008), since at the time of contact only a pair of spheres in the outermost periphery of colliding objects are involved. Hence, the problem reverts to spherical particle collision. This cluster of spheres particles model is adopted in the commercial software EDEM (2011).

There have been numerous studies of DEM applied to the simulation of grinding mills. Some have even applied smoothed particle hydrodynamics (Cleary et al., 2006; 2008) to describe the motion of ore particle slurry along with the motion of larger DEM particles. The single most advance in DEM modeling of grinding mills is the use of collision energy spectra (number of collisions of an energy class versus collision energy) to predict the breakage rate of brittle ore particles (Tavares and Carvalho, 2009). Tavares and Carvalho (2009) described the breakage rate of coarse particles taking into account direct impact breakage of particles by grinding media as well as the weakening of particles by repeated collisions. In a review paper on this subject, Tavares (2017) describes the evolution of microscopic particle breakage rate models from the 1980s to recent advances such as the detailed models of breakage probabilities as a function of impact energy, body and surface breakage, mass of powder caught in each collision, weakening or damage of particles due to repeated collisions, and the incorporation of these models in DEM. In this manuscript, we present a DEM simulation of non-spherical or polyhedral shaped particles tumbling in a laboratory-scale mill. In particular, the scanned polyhedral shapes of rock particles are used in DEM simulation. This study demonstrates that mill power draw can be modeled directly with a mixture of spherical and polyhedral particles. Laboratory-scale measurements of power and captured images of charge motion are compared with DEM simulation done with Blaze-DEM GPU framework software. The method of modeling irregular particles is presented followed by experimental and simulated results. The novelty of this study lies in the fact that the introduction of real rock shapes mixed in with spheres allows improved charge motion representation in DEM simulation of grinding mills. We show the actual scanned polyhedra shapes can lead to accurate prediction of the power draw and charge motion in the tumbling mill. This work paves the way forward for the simulation of plant scale SAG/AG mills using actual large rock shapes passing through these mills.

2. DEM Computing

The simulation of few thousands of particles in DEM is relatively less cumbersome and can be done readily in single-core central processing units (CPU). For hundreds of thousands of particles, parallel computing on multi-core CPUs can be done in a reasonable time frame. However, when millions of particles are simulated, it takes several days of computing for a mere 1 second of real-time on CPUs. Hence, the next step in the development of the DEM framework is the use of graphic processor unit (GPU) computing that has become readily available in the last ten years and which is going to increase in popularity due to many software advances. Today’s GPU offers 24–48 GB of memory and hence handling millions of particles is not a problem. However, DEM requires more severe memory transactions compared to mere numerical computations. These transactions become an order of magnitude higher when computing with polyhedral particles with the number of faces as large as thirty-two. Hence, DEM algorithms must be extremely efficient even with GPU computing.

2.1 Blaze-DEM framework

Blaze-DEM is a GPU-based computational framework for the discrete element methodology (Govender et al., 2016). The evolution of this code has been described in previous publications and we present a few details about the simulation results here. The first feature is that the simulation object geometry can be represented using typical CAD geometry (STL format). This feature is especially useful in importing, for simulation, mill geometry with liners, lifters, grate plates, and pulp lifters in STL format. The code implements “world” objects, to indicate the process vessel under simulation, in this case, a primitive cylinder. Besides, the world object accommodates STL files depicted in triangulated geometry. Next, volume objects in the code can accommodate liners, lifters, grate plates, and pulp lifters. Finally, the material object includes steel spheres and rock particles of different shapes and solid densities. These objects are represented as convex polyhedral shapes joined by planes.

Contact detection is a key feature of DEM. The usual approach is to compute the distance between the center of mass of two colliding objects comparing it to the distance between the two foremost colliding surfaces. Blaze-DEM uses a unique approach named “ray-tracing” in which independent rays (computational threads) are projected from each vertex of polyhedron or center of mass of spheres to check if there is contact with the other surfaces. The ray-tracing algorithm is common in computer graphics. This approach is very convenient to detect contact between two polyhedra since polyhedron is just a collection of surfaces. More details of contact detection involving broad phase, narrow phase, and detailed phase of contact detection are found in Govender et al. (2018).

2.2 Contact mechanics

The mechanics used in the collision between two polyhedra is much different from that commonly used in a pair of colliding spheres. With spheres, the virtual penetration distance is the primary quantity from which forces are estimated in various ways. Herewith polyhedra, the forces are calculated from overlapping volumes (Govender et al., 2018). The advantage of contact volume is that both the direction and magnitude of forces can be resolved with an energy-conserving contact scheme (Govender et al., 2019). The contact scheme is illustrated with a cube and a polyhedral particle colliding with a larger cube, shown in Fig. 1. The intersecting volume is a convex hull formed by the vertices at the intersections between the cube’s edges and polyhedron surfaces.

Fig. 1

Intersecting vertices of overlap and volume overlap between cubes and polyhedron.

Computing the volume of a convex hull, with a computationally efficient algorithm is not a trivial task. The details of volume computing are given in Govender et al. (2018a; 2018b; 2018c). First, the bounding surfaces and their respective normal are determined. Using Green’s theorem, the triple integral of the volume is computed with double integration as a summation over the computed areas.

The simulation shown here involves a few hundreds of polyhedral rock particles and multifaceted steel spheres. As detailed above, even for a few hundred particles, contact detection, and computing intersecting volume becomes an intensive computational task. Besides, the simulations were carried out for eight complete revolutions of the mill cylinder. These computations took a few hours of computing time with an Intel Core i7-6850K 3.6 GHz with an NVIDIA GeForce RTX 2060 GPU.

Normal Contact: The well-known linear viscoelastic spring dashpot model is found to reproduce grinding mill results (Rajamani et al., 2000). The spring force that stores energy and the dissipative force that dissipates energy are given by Eqn. (1) and Eqn. (2), respectively.

  

$${{\mathit{F}}_{\text{n}}}^{\text{spring}}=(K_{\text{n}}\mathrm{\Delta }{v}^{\frac{1}{3}})\widehat{\mathit{n}}$$(1)

  

$${{\mathit{F}}_{\text{n}}}^{\text{dash}}=C_{\text{n}}(V_{\text{R}}\widehat{\mathit{n}})\hspace{0.17em}\widehat{\mathit{n}}$$(2)

where K_n is the spring stiffness ( $\frac{\text{m}}{\text{s}}$), $\hat{\mathit{n}}$ is the normal direction along which the force acts, C_n is the damping coefficient ( $\frac{\text{Ns}}{\text{m}}$), ν is the intersecting volume and V_R is the relative velocity between contacting polyhedra. In the simulation results presented here following refinement of spring constant K_n and dashpot constant C_n were used.

  

$$K_{\text{n}}=\frac{m_{\text{eff}}}{t_{\text{contact}}^2}\text{ln}{(\varepsilon )}^2+{\pi }^2$$(3)

  

$$C_{\text{n}}=\frac{2\text{ln}(\varepsilon )\sqrt{K_{\text{n}}{m}_{\text{eff}}}}{\sqrt{\text{ln}{(\varepsilon )}^2+{\pi }^2}}$$(4)

where ε is the coefficient of restitution, and m_eff is the effective mass of colliding bodies, which is given by Eqn. (5).

  

$$m_{\text{eff}}={\left(\frac{1}{m_1}+\frac{1}{m_2}\right)}^{-1}$$(5)

The common practice in contact force calculations is to ‘break contact’ or terminate contact calculations when the normal force changes direction and is turned inward toward the contact. However, by limiting the contact time, conservation of energy is maintained. Contact time is usually chosen as multiples of computational time steps depending on the properties of the pair of materials in contact.

Tangential contact: Although there are many variations of the spring-dashpot model for tangential contact, the best one that produces consistent mill simulation results in the history-dependent tangential contact proposed by Cundall and Strack (1979). The current time step tangential force is calculated by Eqn. (6) and Eqn. (7).

  

$${{\mathit{F}}_{\text{t}}}^{\text{spring}}=-K_{\text{t}}(V_{\text{t}}\mathrm{\Delta }t)+{{\mathit{F}}_{\text{t}}}^{\prime }$$(6)

  

$${\mathit{F}}_{\text{t}}^{\text{dashpot}}=\text{min\hspace{0.17em}}\left(C_{\text{t}}{V}_{\text{t}},\mu {\mathit{F}}_{\text{n}}^{\text{spring}+\text{dashpot}}\right)+{{\mathit{F}}_{\text{t}}}^{\prime }$$(7)

where F_t′ is the tangential force existing in the contact at the previous time step, but projected onto the current tangential plane and V_t is the relative tangential velocity. It should be mentioned implementing history dependence in GPU architecture is a formidable task since contact detection among millions of particles is done with a hashing algorithm. It is nearly impossible to synchronize a particular previous instant contact with the current one in hashing algorithms.

Power Calculation: The verification of DEM prediction for grinding mills relies on mill power prediction. In addition, a comparison of the experimental charge profiles with simulated profiles offers additional proof. Unfortunately, measurement of stresses on mill shells or lifters is only possible under very particular conditions in the laboratory and pilot-scale mills. The total energy consumed in a collision is simply the energy dissipated in the numerous collision events. Therefore, the power is calculated by summing all the dissipated energy over one full revolution of the mill, given by Eqn. (8).

  

$$E={\sum }_i^k({\mathit{F}}^{\text{dashpot}}\mathrm{\Delta }x)$$(8)

The F^dashpot force defined in Eqn. (2) and Eqn. (7), respectively, and Δx refers to the distance over which the dissipative force acts in each collision. In this intersecting volume approaching Δx in the value $\mathrm{\Delta }{v}^{\frac{1}{3}}$. Finally, the estimated power is given by Eqn. (9).

  

$$P=\frac{E}{t_{\text{rev}}}$$(9)

where t_rev is the time for one revolution. It is customary to run the simulation to about eight revolutions and take the average power over the last few revolutions.

3. Experiment work and Blaze-DEM validation

In this paper, we validate the volume-based polyhedral contact model with a history-dependent tangential model via experiments in a tumbling mill. First, we show the validation for the mill running with spherical steel balls, next with rocks, and finally with a mixture of steel balls and rocks. In addition, an experiment with aluminum cubes is also included for model validations.

3.1 Experimental setup

We used a laboratory-scale slice of a mill fitted with a clear plexiglass front plate for capturing the motion of the charge with a camera. The mill is 90 cm in diameter and 13 cm in internal length (15 cm including front and back plates). The mill is fitted with eight square lifters with a cross-section of 4 × 4 cm. The grinding media were 5.08 cm steel balls of density $7859\hspace{0.17em}\frac{\text{kg}}{{\text{m}}^3}$. The irregular ore particles had a length dimension in the range of 3 and 5 cm and a height dimension in the range of 3 and 4 cm. The mill is fitted with a torque sensor on the drive shaft attached to the cylinder to measure the power drawn by the charge. Thus, the energy losses in the electric motor, are avoided.

Ore particles: Rocks from a local quarry were used in these experiments. These rocks were chosen for their high hardness because they did not break during experiments retaining the integrity of their shape. To characterize the shape of the rocks, four rocks were chosen randomly, shown in Fig. 2. Although it is possible to map the three-dimensional shape with laser scanners, we needed to have fewer faces in the polyhedral representation. Otherwise, the simulation would take unusually long compute times. The rocks were placed in a three-dimensional graduated cardboard axis frame. First, the rocks were triangulated with a pen as well as possible. Next, the x, y, and z coordinates were measured with a ruler and the axis frame. Finally, the polyhedral particle was reconstructed in CAD software. We only expect to get an approximate polyhedron representation of the rock’s shape. Fig. 2 shows the reconstructed polyhedral shapes. It is our view that a much higher shape resolution is not necessary for tumbling mill predictions. These shapes are referred to as PolyA, PolyB, PolyC, PolyD, PolyE, and PolyF in the simulations. The measured density of the rock was $2636\hspace{0.17em}\frac{\text{kg}}{{\text{m}}^3}$. Grinding balls: The diameter of the balls was 2 in and only one single size was used in all the experiments.

Fig. 2

Four random rock particles a, b, c, and d on the left, and their polyhedral representation on the right.

4. Experimental results and DEM validation

As a prelude to experimental results in milling, it is useful to understand the general picture of the charge shape inside the mill. Due to the presence of lifters, the charge in contact with the mill shell moves at the same velocity as the mill shell except for slippage of charge in between the lifters. The particles leaving from the shell at the shoulder of the mill executed a parabolic trajectory. These trajectories meet the mill shell on the other side of the mill. To describe different locations on the mill shell, we use the ubiquitous mathematical notation depicted (on Fig. 3) of angle increasing counterclockwise from the x-axis, reaching 90° at the positive y-axis, and gradually increasing to 180° at the negative side of the x-axis followed by 270° at the negative y-axis, where the mill’s center is the origin for the x-axis and y-axis frames. At a higher percentage of critical speeds, the particles leaving the shoulder may strike in the vicinity of 180° or even higher at 150°.

Fig. 3

Description of the definition of shoulder and toe in the charge profile.

The fully cascading charge shape has a bi-linear shape for the free surface above the charge. The free surface starts in between 210° and 240° marks, termed “toe”, rising at a shallow angle, to the inflection point near the center of the mill, and then the charge rises steeply to the shoulder position, which is around 30° to 60° mark. This profile is entirely directed by the mill filling volume. Between the mill shell and charge surface, there is a layer traveling upward followed by layers cascading downward. Between these two layers, there are a few layers of particles that are nearly stationary.

Next, we show a set of experimental results and DEM simulation comparisons for mill charge flow and power draw. The experiments were run with balls only, rocks only, and then a mixture of balls and rocks. Finally, an experimental run with aluminum cubes is shown. The DEM simulation requires 10 sets of parameters for the respective contact mechanics between balls, rocks, and the mill shell. In each of the experiments shown here a table of the contact parameter is also included. A discussion of this parameter set is given at the end of this section.

Case 8: This case contains 190 steel spherical balls which contribute to 25 % volume filling of the mill for a total mass of 97.4 kg. The density of the steel sphere was set to $7800\hspace{0.17em}\frac{\text{kg}}{{\text{m}}^3}$. The mill speed was set at $31.6\hspace{0.17em}\frac{\text{rev}}{\text{min}}$. Table 1 gives the DEM parameters used in the simulation. Fig. 4(a) shows the image of the charge velocity profile for the DEM simulation, and Fig. 4(b) shows the image of the charge in the mill. With the mill filling being all steel spheres, the two charge profiles match very well at the toe, shoulder, and belly positions. There were a few balls cataracting in the actual experiment. The balls are hitting above the toe region and between 180° to 270°. This characteristic of the profile is also seen in the DEM simulation, which is a clear indication of the simulation accuracy. The experimental power is 1762 ± 30 W, whereas the DEM simulation predicts 1609 ± 40 W.

Table 1

Simulation parameter for Case 8, P2P indicates particle-particle and P2W indicates particle-wall settings.

Parameter	Symbol	Value	Unit
Normal Particle Stiffness	Kn	4 × 106	N/m
Tangent Particle Stiffness	Kt	3 × 106	N/m
Damping Coefficient	C	0.01
Coefficient of Restitution (P2P)	ɛ	0.85
Static Friction (P2P)	μs	0.7
Kinetic Friction (P2P)	μk	0.7
Coefficient of Restitution (P2W)		0.7
Static Friction (P2W)		0.6
Kinetic Friction (P2W)		0.6
Time Step	Δt	1 × 10−4	s

Fig. 4

(a) Simulation velocity profile and (b) experimental results for Case 8 (25 % balls, 31.8 rev/min).

We now present the simulation results for the polyhedral (rock) particles alone or in a mixture with spherical steel balls, comparing the power draw and charge profiles. The experimental mill reaches a steady-state in two full revolutions of the mill and hence all the results presented here are at a steady-state.

Case 11: This case is a mixture with 354 rock particles (12 %) and 99 5 cm steel balls (13 %) for a total mass of 66.3 kg. The mill speed was set at $25.1\hspace{0.17em}\frac{\text{rev}}{\text{min}}$. Table 2 gives the DEM parameters used in the simulation. Fig. 5(a) shows the image of the charge and Fig. 5(b) depicts the DEM profile colored by particle type. In this case, the spherical particle too was modeled as a 32-equisided polyhedron. Due to repeated impacts, the steel spheres become non-round over time. Hence, we capture the shape of the sphere that looks like a soccer ball with the 32-sided polyhedron. Such representations of the shapes, adds to the accuracy of the particle trajectories. The two profiles agree well regarding the toe position, shoulder position, and the belly of the charge. Also, few particles cataracting is well matched between experiment and DEM simulation. Due to a very high degree of angularity in the rock particles, the charge seems to lock up with each other and hence we see a steeper shoulder position than normal. The DEM simulation captures this behavior. The DEM power draw of 1,010 W ± 5 compares well with the experimental value of 1105 W.

Table 2

Simulation parameter for Case 11, P2P indicates particle-particle and P2W indicates particle-wall settings. * Rock/Spheres values.

Parameter	Symbol	Value	Unit
Normal Particle Stiffness	Kn	*2/4(× 106)	N/m
Tangent Particle Stiffness	Kt	*1/3(× 106)	N/m
Damping Coefficient	C	0.001/0.1
Coefficient of Restitution (P2P)	ɛ	0.4/0.85
Static Friction (P2P)	μs	0.2
Kinetic Friction (P2P)	μk	0.2
Coefficient of Restitution (P2W)		0.7
Static Friction (P2W)		0.6
Kinetic Friction (P2W)		0.6
Time Step	Δt	1 × 10−5	s

Fig. 5

(a) Experimental profile and (b) DEM results for Case 11 (12 % rocks, 13 % balls, 25.1 rev/min).

Case 1: In these experiments, a low-level filling of rocks (15 %) made up of 434 pieces was used. The mill speed was set at $27.4\hspace{0.17em}\frac{\text{rev}}{\text{min}}$. Fig. 6(a) shows the image of the experimental mill charge and Fig. 6(b) is the corresponding DEM simulation snapshot. Once again, the two charge profiles agreed quite well. In both the experimental and simulation, the charge is confined between 240° mark (toe) and 30° marks (shoulder). Further, Fig. 6(a) shows five to eight rock particles that are carried up by the square lifter, above the shoulder of the cascading charge. These rock particles then follow the parabolic trajectory downward. We see exactly similar transport by the lifter in the simulation too. Although not perfect, the nearly close charge profile structure between the experiment and DEM simulation, in this case, shows that the friction part of the contact model is bearing out the trends. This case stands in contrast to mill operating with only spherical steel spheres, wherein DEM simulation has produced picture-perfect results coinciding with experimental work (Venugopal (2001), Govender et al. (2015)). Finally, the DEM power 482 W ± 30 compares well with the experimental value of 424 W.

Fig. 6

(a) Experimental and (b) DEM result for Case 1 (15 % rocks, $27.4\hspace{0.17em}\frac{\text{rev}}{\text{min}}$ ).

Case 20: In this experiment, a very high filling made up of 1031 rock particles (35 % filling) was run to test the ability of DEM simulation to predict higher power draw, which in turn involves numerous polyhedra collisions. The primary purpose of these experiments work is to verify the volume of contact mechanics used in the Blaze-DEM code. Table 3 gives the DEM parameters used in the simulation. Fig. 7(a) shows the experimental mill snapshot and Fig. 7(b) shows the snapshot from DEM simulation. Due to the high load, the shoulder of the charge is almost at the 80° mark. The speed of the mill too was very high set at $33.6\hspace{0.17em}\frac{\text{rev}}{\text{min}}$ or 75 % of mill critical speed. The simulation tracks the overall structure of the central body of the charge including the toe and shoulder. However, the simulated charge is slightly expanded compared with the experimental frame. This is attributed to the representation of the rock particles with only four polyhedral shapes depicted in Fig. 2. A more comprehensive sample from the 1031 particles, could have provided more accurate modeling of the charge. Yet, in milling, we are interested mainly in the cascading and cataracting charge motion and mill power. Hence, the results presented here verify the ability of the polyhedral contact algorithms at best. The computed power draw of 796 ± 10 W compares well with the experimental power draw of 856 W.

Table 3

Simulation parameters for Case 1, P2P indicates particle-particle and P2W indicates particle-wall settings.

Parameter	Symbol	Value	Unit
Normal Particle Stiffness	Kn	2 × 106	N/m
Tangent Particle Stiffness	Kt	1 × 106	N/m
Damping Coefficient	C	0.01
Coefficient of Restitution (P2P)	ɛ	0.4
Static Friction (P2P)	μs	0.5
Kinetic Friction (P2P)	μk	0.2
Coefficient of Restitution (P2W)		0.7
Static Friction (P2W)		0.6
Kinetic Friction (P2W)		0.6
Time Step	Δt	5 × 10−5	s

Fig. 7

(a) Experimental and (b) DEM charge profile for Case 20: 35 % rocks, 33.6 rev/min.

Case 4C: A mill charge of cubes is an obvious choice for DEM verification since the entire mass is made up of cubes. In previous experiments, although we used irregular rock particles, we had to model only a few rock shapes as polyhedrons for the sake of expediency. Hence, there is a difference between experimental results and simulation results due to the variations in the rock shape. In this experiment, the entire charge was made up of 25.4 mm aluminum cubes, with a density of $1800\hspace{0.17em}\frac{\text{kg}}{{\text{m}}^3}$. These cubes were identical in shape since they were cut from a rod of square cross-section (25.4 × 25.4 mm). Table 4 gives the DEM parameters used in the simulation. Fig. 8(a) shows the mill operating with 884 aluminum cubes at a mill speed of $29.2\hspace{0.17em}\frac{\text{rev}}{\text{min}}$. There is a much closer resemblance between the experimental charge profile and the simulated profile (Fig. 8(b)), simply because all the particles are identical in shape. The shoulder and toe position agree very well. The cataracting stream from the shoulder also shows good agreement. Similarly, the cubes striking the toe region too is well predicted. The experimental power is 711 ± 11 W, whereas the DEM simulation predicts 752 ± 72 W.

Table 4

Simulation parameters for Cases 4c and 6c, P2P indicates the particle-particle and P2W indicates particle-wall settings.

Parameter	Symbol	Value	Unit
Normal Particle Stiffness	Kn	6 × 106	N/m
Tangent Particle Stiffness	Kt	3 × 106	N/m
Damping Coefficient	C	0.01
Coefficient of Restitution (P2P)	ɛ	0.8
Static Friction (P2P)	μs	0.6
Kinetic Friction (P2P)	μk	0.6
Coefficient of Restitution (P2W)		0.7
Static Friction (P2W)		0.6
Kinetic Friction (P2W)		0.6
Time Step	Δt	5 × 10−5	s

Fig. 8

(a) Experimental and (b) DEM charge profiles for Case 4c: 25 % filling with 25 mm aluminum cubes at $29.2\hspace{0.17em}\frac{\text{rev}}{\text{min}}$ .

Case 6C: in this experiment, the number of cubes was increased to 1061, making up a mill volume filling of 30 %. Here too, we see excellent agreement between DEM polyhedral model and experimental profiles, as shown in Fig. 9. The computed power draw of 752 ± 37 W agrees with experimental values of 740 ± 11 W.

Fig. 9

(a) Experimental and (b) DEM charge profiles for Case 6C: 30 % filling with 25-mm aluminum cubes at $29.2\hspace{0.17em}\frac{\text{rev}}{\text{min}}$.

4.1 Contact mechanics parameters

The parameters used in the simulation deserve an explanation here in the milling study. The four main parameters are the two spring constants and the two dashpot constants, normal and tangential, respectively. The spring constants mainly control the time of contact, i.e number of time steps the particle pair interpenetrate before disengaging. The dashpot constant controls the energy dissipated (computed as contact force reduction) during the entire duration of contact between a pair of particles. A reasonable choice of these four parameters produces acceptable simulation results with DEM. A point to note here is that, in a majority of cases published in the literature, DEM simulation target has been bulk charge motion simulation as inflow discharge from a silo, motion over a conveyor chute, mixing vessels, rock-scooping buckets, and many others. In all these simulations principally the metric in bulk particle velocity, for the simulation to resemble measured experiments. In fact, in most such cases, approximate agreement of bulk flow is shown between simulation and DEM results.

Milling, on the other hand, requires correct prediction of charge motion as well as mill power draw, so the four parameters mentioned earlier, have to reach a compromise between accurately predicting charge motion and power. In this compromise, the dashpot constants play a critical role. In addition, the spring constants determine the energy dissipated in each collision since they determine the duration of the collision event. Suffice it to say, there is no unique set of contact parameters that will universally describe a system. In any particular system, a set of parameters is arrived at through trial and error carried over numerous experiments and simulations.

Table 5 shows the contact parameters used in the simulation shown in this manuscript. They range from balls only, rocks only to a combination of rock and balls followed by aluminum cubes. In this particular set of experiments, the spring constant (normal/tangential) for balls only (4 × 10⁶, 3 × 10⁶) and rocks only (2 × 10⁶, 1 × 10⁶) seems to predict well for the mixture of balls and rocks (4 × 10⁶/2 × 10⁶, 3 × 10⁶/1 × 10⁶). The same trend is noticed with the coefficient of restitution values. The static and kinetic friction values trigger minor changes in the energy of collisions. This level of perfect agreement between rock, ball, and rock+ball parameters is fortuitous and may not hold good for all kinds of variations in mill operations. Nevertheless, the magnitude of the constants shown in Table 5 describes laboratory milling operation both in two-dimensional and three-dimensional simulations. This order of magnitude in parameters when carried on to large mill simulations (9 m diameter × 3.5 m length) also produces reasonable power results. In summary, the parameter set in DEM is still an ongoing open question.

Table 5

Summary of simulation parameters used during the simulation, where (b) ball only, (r) rock only, (b&r) mixed ball and rock, and (a) aluminum cubes.

Case #	6, (b)	1, (r)	20, (r)	11, (b&r)	4C, (a)	6C, (a)
ɛ	0.85	0.4	0.4	0.85/0.4	0.8	0.8
μs	0.7	0.7	0.7	0.2	0.8	0.6
μk	0.7	0.7	0.7	0.2	0.8	0.6
Kn(106)	4	2	2	4/2	6	6
Kt(106)	3	1	1	3/1	3	3
% mill filling	15	15	35	25	25	30
$\text{Speed\hspace{0.17em}}\frac{\text{rev}}{\text{min}}$	31.6	27.4	33.7	25.1	29.2	29

5. Conclusion

It is shown that volume-of-contact contact mechanics is able to successfully predict the motion of charge in a laboratory-scale mill. The polyhedral particle model helps in modeling rock and ball assembly that is commonly encountered in mineral ore grinding plant operations. The actual rock shapes, especially in the size range of 0.2–0.4 m can be scanned and modeled as irregular polyhedra. In addition, GPU-based codes such as Blaze-DEM enables one to study mill discharge through grate slots and pulp lifters, simply because such codes can handle millions of particles.

References

•  Cleary  P.W., Recent advances in DEM modelling of tumbling mills, Minerals Engineering, 14 (2001) 1295–1319. DOI: 10.1016/s0892-6875(01)00145-5
•  Cleary  P.W.,  Morrisson  R.,  Morrell  S., Comparison of DEM and experiment for a scale model SAG mill, International Journal of Mineral Processing, 68 (2003) 129–165. DOI: 10.1016/s0301-7516(02)00065-0
•  Cleary  P.W., Predicting charge motion, power draw, segregation and wear in ball mills using discrete element methods, Minerals Engineering, 11 (1998) 1061–1080. DOI: 10.1016/s0892-6875(98)00093-4
•  Cleary  P.W.,  Sinnott  M.,  Morrison  R., Prediction of slurry transport in SAG mills using SPH fluid flow in a dynamic DEM based porous media, Minerals Engineering, 19 (2006) 1517–1527. DOI: 10.1016/j.mineng.2006.08.018
•  Cleary  P.W.,  Sinnott  M.D.,  Morrison  R.D., DEM prediction of particle flows in grinding processes, International Journal for Numerical Methods in Fluids, 58 (2008) 319–353. DOI: 10.1002/fld.1728
•  Cundall  P.A.,  Strack  O.D.L., A discrete numerical model for granular assemblies, Geotechnique, (1979) 47–65. DOI: 10.1680/geot.1979.29.1.47
•  Delaney  G.W.,  Cleary  P.W.,  Sinnott  M.D.,  Morrison  R.D., Novel application of DEM to modelling comminution processes, IOP Conference Series: Materials Science and Engineering, 10 (2010) 012099. DOI: 10.1088/1757-899x/10/1/012099
•  Djordjevic  N.,  Morrison  R.,  Loveday  B.,  Cleary  P., Modelling comminution patterns within a pilot scale AG/SAG mill, Minerals Engineering, 19 (2006) 1505–1516. DOI: 10.1016/j.mineng.2006.07.003
• EDEM, User Guide, DEM solution Ltd., Edinburgh, Scotland, UK Copyright © (2011). <http://tm.spbstu.ru/images/2/28/EDEM2.4_user_guide.pdf> accessed 30.07.2021.
•  Govender  N.,  Rajamani  R.K.,  Kok  S.,  Wilke  D.N., Discrete element simulation of mill charge in 3D using the BLAZE-DEM GPU framework, Minerals Engineering, 79 (2015) 152–168. DOI: 10.1016/j.mineng.2015.05.010
•  Govender  N.,  Wilke  D.N.,  Kok  S., Blaze-DEMGPU: Modular high performance DEM framework for the GPU architecture, SoftwareX, 5 (2016) 62–66. DOI: 10.1016/j.softx.2016.04.004
•  Govender  N.,  Rajamani  R.,  Wilke  D.N.,  Wu  C.-Y.,  Khinast  J.,  Glasser  B.J., Effect of particle shape in grinding mills using a GPU based DEM code, Minerals Engineering, 129 (2018a) 71–84. DOI: 10.1016/j.mineng.2018.09.019
•  Govender  N.,  Wilke  D.N.,  Pizette  P.,  Abriak  N.-E., A study of shape non-uniformity and poly-dispersity in hopper discharge of spherical and polyhedral particle systems using the Blaze-DEM GPU code, Applied Mathematics and Computation, 319 (2018b) 318–336. DOI: 10.1016/j.amc.2017.03.037
•  Govender  N.,  Wilke  D.N.,  Wu  C.-Y.,  Tuzun  U.,  Kureck  H., A numerical investigation into the effect of angular particle shape on blast furnace burden topography and percolation using a GPU solved discrete element model, Chemical Engineering Science, 204 (2019) 9–26. DOI: 10.1016/j.ces.2019.03.077
•  Govender  N.,  Wilke  D.N.,  Wu  C.-Y.,  Rajamani  R.,  Khinast  J.,  Glasser  B.J., Large-scale GPU based DEM modeling of mixing using irregularly shaped particles, Advanced Powder Technology, 29 (2018c) 2476–2490. DOI: 10.1016/j.apt.2018.06.028
•  Hopkins  M.A.,  Hibler III  W.D.,  Flato  G.M., On the numerical simulation of the sea ice ridging process, Journal of Geophysical Research: Oceans, 96 (1991) 4809–4820. DOI: 10.1029/90jc02375
•  Kiangi  K.,  Potapov  A.,  Moys  M., DEM validation of media shape effects on the load behaviour and power in a dry pilot mill, Minerals Engineering, 46–47 (2013) 52–59. DOI: 10.1016/j.mineng.2013.03.025
•  Latham  J.-P.,  Munjiza  A.,  Garcia  X.,  Xiang  J.,  Guises  R., Three-dimensional particle shape acquisition and use of shape library for DEM and FEM/DEM simulation, Minerals Engineering, 21 (2008) 797–805. DOI: 10.1016/j.mineng.2008.05.015
•  Mishra  B.K.,  Rajamani  R.K., Motion analysis in tumbling mills by the discrete element method, KONA Powder and Particle Journal, 8 (1990) 92–98. DOI: 10.14356/kona.1990016
•  Mishra  B.K.,  Rajamani  R.K., The discrete element method for the simulation of ball mills, Applied Mathematical Modelling, 16 (1992) 598–604. DOI: 10.1016/0307-904x(92)90035-2
•  Morrison  R.D.,  Cleary  P.W., Using DEM to model ore breakage within a pilot scale SAG mill, Minerals Engineering, 17 (2004) 1117–1124. DOI: 10.1016/j.mineng.2004.06.016
•  Potapov  A.V.,  Campbell  C.S., A three-dimensional simulation of brittle solid fracture, International Journal of Modern Physics C, 07 (1996) 717–729. DOI: 10.1142/s0129183196000594
•  Rajamani  R.K.,  Mishra  B.K.,  Venugopal  R.,  Datta  A., Discrete element analysis of tumbling mills, Powder Technology, 109 (2000) 105–112. DOI: 10.1016/s0032-5910(99)00230-2
•  Rothenburg  L.,  Bathurst  R.J., Numerical simulation of idealized granular assemblies with plane elliptical particles, Computers and Geotechnics, 11 (1991) 315–329. DOI: 10.1016/0266-352x(91)90015-8
•  Tavares  L.M.,  André  F.P.,  Potapov  A.,  Maliska  C., Adapting a breakage model to discrete elements using polyhedral particles, Powder Technology, 362 (2020) 208–220. DOI: 10.1016/j.powtec.2019.12.007
•  Tavares  L.M., A review of advanced ball mill modelling, KONA Powder and Particle Journal, 34 (2017) 106–124. DOI: 10.14356/kona.2017015
•  Tavares  L.M.,  de Carvalho  R.M., Modeling breakage rates of coarse particles in ball mills, Minerals Engineering, 22 (2009) 650–659. DOI: 10.1016/j.mineng.2009.03.015
•  Venugopal  R.,  Rajamani  R.K., 3D simulation of charge motion in tumbling mills by the discrete element method, Powder Technology, 115 (2001) 157–166. DOI: 10.1016/s0032-5910(00)00333-8
•  Xu  L., Luo  K.,  Zhao  Y.,  Fan  J.,  Cen  K., Influence of particle shape on liner wear in tumbling mills: a DEM study, Powder Technology, 350 (2019) 26–35. DOI: 10.1016/j.powtec.2019.03.033

Authors’ Short Biographies

Alberto Puga

Alberto Puga completed his Bachelor’s degree in Mechanical Engineering from the University of Utah in Salt Lake City, in 2015. In 2016, he joined the Department of Metallurgical Engineering at the University of Utah in Salt Lake City as MS student. Later, he transferred to pursue a PhD in Metallurgical Engineering with an emphasis in mineral processing. As a student, he was a research assistant in comminution projects involving breakage mechanics in ball mills. He worked as a teaching assistant for the mineral processing and statistics courses. Because of his interaction with the students, he was awarded the “Teaching Assistant of the Year” in the department.

Nicolin Govender

Dr. Nicolin Govender is a principal scientist at Mondelez International in the position of global lead for granular material physics and GPU/HPC computing for the creating of digital twins. Prior to this, he held the prestigious EU Marie Curie Research Fellowship in the Department of Chemical and Process Engineering at the University of Surrey. Dr. Govender is the primary investigator (PI) of the NVIDIA GPU Research Centre for Discrete Element Particle Modeling, which he established in 2013 with the support for NVIDIA, specializing in GPU computing of particle systems. He has been working on GPU clusters since the advent of general-purpose GPU computing in 2008 and is an esteemed international expert in the field. He has solely developed the open source GPU-based DEM code “Blaze-DEM”, that is hosted by Elsevier, which is used by a growing number of institutes, nationally and internationally.

Raj K. Rajamani

He has been on the research and teaching faculty position at the University of Utah since 1979. Currently he is Professor of Material Science and Engineering. Dr. Rajamani’s research interests include: Breakage rate of mineral particle in ball mills and scale-up to plant scale mills, computational Fluid Dynamic study of hydrocyclones, Discrete Element Simulation of mineral ore grinding mills and electrodynamic sorting of auto-shredder scrap stream to recover aluminum and copper. Dr. Rajamani received the Antoine M. Gaudin Award, presented by the Society of Mining, Metallurgy and Exploration Engineers Inc. in the year 2009. The award plaque reads, “For his seminal work in the application of discrete element methods in the modeling of charge motion in semi-autogenous and ball mill grinding, and for his contribution to the basic science of comminution and classification.”