Motion Analysis in Tumbling Mills by the Discrete Element Methodt

The media motion in a ball mill is simulated using a numerical algorithm known as the Discrete Element Method. The motion of the charge is modelling by considering the forces acting at each contact point during a collision and following the movement of individual balls as per Newton's law. First, experimental verification of the model is shown. Then a few simulated results are shown for some interesting cases. For the first time it is possible to get information about the distribution of energy in the col lision between balls. Also the cataracting and cascading motion of charge in large mills can be simulated accurately with this simulation program.


Introduction
Grinding is an important unit operation in metal production from mined ores. Crushed ore in the size range of I~ 3 centimeters is ground to a size finer than 72 microns. In almost all plants a tumbling ball mill is used for grinding. Importantly, this operation consumes a large amount of energy; the power cost may represent more than half the total processing cost. Therefore considerable attention has been focused on the breakage rates and energy consumption in tumbling mills. Our knowledge of comminution (the effect of operating variables on mill performance, modelling, simulation and control) has greatly widened, yet we face fresh problems that require both practical as well as theoretical solutions. One of the current problems is predicting the power draft of the mill. This would entail knowing precisely the distribution of forces during the tumbling motion of the balls. A good approach to the above problems clearly resides in the understanding of the dynamic motion of the balls within the mill, which is the genesis of the forces experienced by the particles to be ground.
Motion analysis of the charge of particulate material in a tumbling mill was reported, by White 1 ) and Davis 2 ), among others and has been reviewed by Rose and Sullivan 3 >. Although simplistic in nature, these works established the motion of the charge inside a ball mill approximately.
Once the charge motion was established, researchers took different approaches to predict the power draft of the mill. In one such approach, certain assumptions were made with regard to the dynamic charge: the entire mill charge composed of a single mass whose profile is determined by a line joining the toe and shoulder of the ball mass and the mechanical energy sustaining the mass at this position continuously. While the mill is at rest, the center of mass of the charge is on the vertical line passing through the center of the mill axis, but as the mill rotates the center of mass of the charge assumes a new position to the right of the vertical axis (for counterclockwise mill rotation). The torque required to maintain the charge in its offset position is then computed from the load and the offset distance from the mill center. The power draft is then related to the torque through the rotational speed of the mill.
The power-draft calculation described is quite simple, and hence the resulting analytical equations are widely used in the mineral processing industries. However, a host of factors like the profile of the charge, the natural angle of response of the charge, frictional forces acting at the mill wall, the shape of the lifter bars, the size distribution of the ball charge, etc., influence power draft so much that the simple analysis begins to fail when higher accuracy is demanded. Liddell and Moys 4 ) have proved this point: none of the torque formulae could predict their carefully measured power draft of the 0.545-m X 0.308-m mill. In their recent works Mishra and Rajamani 5 ) have shown that the center of gravity of the charge is not a fixed point in space; rather, it changes its position as the mill rotates. They have also shown 6 ) that the profile of the charge changes with the shape of the lifter bars. Therefore, any torque analysis treating the ball charge as a single mass is unrealistic.
The breakage rate of particles, wear rate of balls and liners, and power draft are known to be intimately related to the motion of the charge. Hence to accurately predict these quantities, motion of individual balls must be modelled. An analytical solution within the realm of classical mechanics is impossible. Therefore, the intention here is to explore some other avenues of fundamental inquiry, by means of which motion and collision of each and every ball in the ball charge can be tracked. Such an analysis is possible with the Discrete Element Method, which is a numerical tool for simulating multi-body collisions. In this analysis, the motion of the charge is treated as a multi-body collision problem, for which only the material properties with regard to impact must be known. The purpose of this paper is to present the results obtained through numerical simulation and discuss the trends in the results in comparison to physical experiments.

The Model
The Discrete Element Method (DEM) is a numerical scheme pioneered by Cundall and Strack 7) for simulating the behavior of systems of discrete interacting bodies. It has gained popularity over the past decade in several areas of science where behavior of particulate material had to be studied taking an alternate route to the continuum approach. This scheme allows finite displacements and rotation of discrete bodies, where complete loss of contact and formation of new contacts take place as the calculation progresses. The program DISC KONA No. 8 (1990) developed by Corkum and Ting 8 ) implements the DEM algorithm for two-dimensional discshaped bodies; it has been used to study basic soil behavior and geotechnical problems. This code has been modified by Mishra 9 ) to simulate media motion in tumbling mills.
The numerical scheme adopted in the DEM formulation applies Newton's second law to the discs and force displacement law at the contacts. Newton's second law gives the motion of the discs resulting from the forces acting on it. The forces developing at the contacts are modelled by a pair of normal and tangential spring-dashpots at every contact point. The forces are found from the overlaps of the discs at the contact points by using a linear force displacement law; the contact forces are proportional to the amount and the rate of overlap between the discs. The provision of viscouscontact damping models the system realistically, because the coefficients of restitution for different contacts are used to compute the damping parameters. Also, friction damping occurs during sliding in any time step, when the absolute value of the shear force at any contact exceeds a maximum value found by the product of the normal force and the coefficient of friction. As shown in Fig. 1, every disc is identified separately; they are allowed to overlap at the contact points; and the calculation is done for each disc in turn. The relative velocity of the disc i with respect to other discs k, l, m, and n is determined for each contact. These relative velocities for every contact of disci are resolved in the normal (along the line drawn through the disc centers) and shear direction, and the force calculation is then done for each contact in the following way: where 6Fn and 6Fs are the incremental forces due to the springs, 6dn and 6ds are the incremental forces due to the dashpots, dvn and dvs are incremental velocities, K and C are the spring and dash pot constants, respectively, and M is the incremental time step. Then the contact forces on a disc are added to get the net out-of-balance force acting on it. The acceleration of disc i of mass m; is given by where .X, y, 8 are the accelerations in the x, y and e directions, respectively, 1 0 ; is the moment of inertia of disc i, and M 0 is the total moment acting on the disc. Then by integrating acceleration with a central difference scheme, the velocity and displacement are determined. Thus, by repeating this calculation for each disc, the new position of all the discs is determined, and, by repeating the entire set of calculation at successive time steps, the movement of the charge in the x, ycoordinate space is given.

Experimental Work
The majority of the experimental work was carried out in a small ball mill, and specialized equipment was used to determine material parameters. The principal aim of the experimental work is to verify the numerical simulations. The parameters needed for the numerical calculation are the coefficient of restitution, materi-94 al stiffness and coefficient of friction. The material stiffness is determined by dropping the ball over a layer of particles situated on an anvil. A strain gauge embedded in the anvil records the strain waves. The analysis of the data gives stiffness.
Similarly in a pendulum experiment, the coefficient of restitution for the ball-ball collision is determined. Drop ball experiments are done to find the same for ball-wall collisions. These values of the coefficient of restitution are used to estimate damping parameters.
The coefficient-of-friction data are obtained from the results of the experiments done by Rose and Sullivan 3 ). The data used in the simulations are given in Table 1.  IO) for studying the movement of ball media in vibration mills.
For recording the motion of ball charge experimentally, a ball-mill/video-camera system was built. The video camera records events on a video tape, and a frame-grabber hardware transmits the image to the computer memory. Then a software program tracks the x, y-coordinates of the moving balls. In this manner the motion of the entire ball mass is compared with the simulation results. A schematic of this equip-

Ball M1ll
Frame grabber Fig. 2 Experimental apparatus ment is shown in Fig. 2. The ball mill is a thin cylinder I centimeter thick and 30 centimeters in diameter, fitted with six equally spaced lifter bars. The discs are 3.4 centimeters in diameter and 0.8 centimeters thick. The discs are positioned in the mill, and a plexiglas end plate is bolted to the lifter bars. The plexiglas plate allows video recording of the charge motion. Thus, the experimental mill is made to correspond to the simulated mill.

Results
It is easy to calculate the trajectory of a  (Fig. 3a) while the other rests on the outermost layer of balls (Fig. 3b). The mill is filled with 24 balls of 3.4 centimeters diameter resulting in a 40 percent mill filling. Experimentally determined positions are compared against numerical simulation results for two revolutions of the mill. Even though collision events are stochastic in nature, the comparison is surprisingly close. It is difficult to say, however, that predictions would be equally good in the third revolution.
Since the disc surface is rough, the actual trajectories deviate slightly from the theoretical, and as time passes these deviations build up to a large extent. It is believed that by correctly choosing the parameters the collision events can be simulated within the limits of acceptable error. Even though the trajectories of individual balls become stochastic at longer periods of time, the mean tendency of the charge motion can be predicted with DEM. After a number of raised. Therefore, to produce mostly cascading motion and very little cataracting of balls, the mill should be operated at an optimum speed.
It is possible to ascertain this speed for any size mill with this simulation procedure. After having built some confidence in the computer code, a few numerical experiments were carried out to compare the trends in the simulation results with actual observations. Lifter-bar shape tremendously influences the profile of the charge, but due to wear lifter-bar shape changes after prolonged hours of milling, causing a change in the ball trajectories; these changes in turn affect the power drawn. As an illustration, Fig. 5 shows snapshots of the charge profile as it evolves, for two different face angles of lifter bars -134 degrees and 90 degrees. A comparison of power draft indicates that a 30-centimeter diameter mill, operating at 60 percent of its critical .>peed, 40 percent filling, with 90-degree face-angle lifter bars, draws twice as much power when compared to the same mill fitted with 130-degree face-angle lifter bars. How much of this extra power is utilized in particle breakage is not known.
It is possible to calculate the energy distribution of collisions in the ball charge. What makes this possible is the contact model using normal and tangential spring-dashpots. The frequency distribution of energy expended in various collisions reveals more diverse information about the behavior of the mechanical process. To il---::::: IMPACT ENERGY (mJ) Fig. 6 Frequency distribution of energy between collisions for different mill speeds lustrate this, a 55-centimeter diameter mill with 40 percent ball filling was simulated. Figure 6 shows the frequency of impacts for a single ball during one revolution, as a function of impact energy, for different rotational speeds.
It is seen that the number of impacts corresponding to small energy values is more than the number of impacts corresponding to that of large energy values. This happens when the cascading motion predominates the overall motion of the charge; this fact was readily noticed in the simulation-animation results which are not shown here.
It is important to know how the impact frequency corresponding to a given energy value changes with the rotational speed; with this information the particle breakage process can be controlled. Figure 7 shows the cumulative number of impacts greater than a given energy plotted against the rotational speed of the mill. The number of low-energy impacts ( ~ 10 mj) decreases with increasing rotational speed. The frequency of medium-(~50 mj) and high-energy ( ~ 100 mj) impacts increases with increasing rotational speed and subsequently decreases with speed. Also, low-energy impacts are greater in number than the high-energy impacts at any rotational speed of the mill. As the balls begin to centrifuge, the total number of impacts decreases. A combination of these impacts influ-KONA No. 8 (1990) ences the mode of breakage of the material and the power draft. Also from the figure, it is safe to predict that the maximum in power draft of the mill would occur at about 70 percent of the critical speed. It is known that mill capacity is at its highest when the power draw is highest. For a given liner configuration and ball size distribution, it is possible to predict power draft with this DEM simulator. The critical speed of the mill applies only to the outer layer of the balls, i.e., the ones which are in contact with the mill shell. This outer layer of balls cataracts, and the inner layers of balls then tend to cascade, as they would in a mill operated at a low speed. So, at any given time, some of the balls would be in flight (cataracting) while others would not be. The proportion of these balls with respect to one another determines the frequency distribution of impacts. However, these energies of impacts cannot be determined experimentally, and one has to rely on numerical schemes as presented here. Rolf and Simonis 11 > measured the energies of impact in a carefully planned experiment. The balls were fitted with strain gauges and memory circuits and after going through collisions the information was transmitted to a computer for analysis. The simulated trends of energy frequency agree with Rolf's experimental findings.

Conclusions
An analytical description of multi-body collision is difficult, but the Discrete Element Method enables numerical simulation of ball charge motion. Currently, the simulation requires hours of supercomputer time, depending on the number of balls to be simulated. Published models describing the motion of the charge in a mill are invariably based on empiricism, and hence their use is limited. DEM simulation can give collision frequency and impact energy frequency. The collision frequency can be subdivided into ball-ball collisions and ballliner (shell) collisions. Hence, with these frequencies, rates of both liner and ball wear and particle breakage can be modelled accurately.
The simulation results pertaining to two important aspects of ball motion are checked against laboratory experiments: First, the position of two arbitrary balls are tracked and found to match with the simulated trajectories, and second, the toe and shoulder positions defining the profile of the charge as a whole as predicted by simulation agree with the experiments. Also the charge motion was simulated for two different shapes of the lifter bars. The snapshots of a charge motion depict the differences in the overall motion of the charge.
The key to any plant operation lies in a good understanding of the individual unit operations; ball mills in a concentrator are no exception. Very little is known about the processes taking place inside the mill due to the lack of sensors which can function in the harsh environment. Moreover, for process simulation, an adequate theory capable of describing the individual motion of balls is yet to emerge. As described here, it is possible to simulate multi-body collisions. Some simulations were carried out to study the impact process occurring within the mill. The low-energy impacts are far larger in number than high-energy impacts for all speeds. Furthermore, the number of impacts for a specified energy value reaches a maximum and then drops off as speed increases.
The results shown here, though in no way exhaustive, are quite fundamental in nature, and unveil the modes of charge motion and hence the mechanism of breakage. For the first time, the frequency distribution of energy of collisions is numerically calculated. This work is currently being pursued at greater lengths, to establish a mechanism of particle breakage with respect to the modes of charge motion. It will then be possible to predict the size distribution of the particles being milled. We believe that this work should open a new dimension of further investigation into comminution in ball mills.