Modelling of the High-Energy Ball Milling Process

Mohsen Mhadhbi

Laboratory of Useful Materials, National Institute of Research and Physical-chemical Analysis, Ariana, Tunisia.

DOI: 10.4236/ampc.2021.111004                

Abstract

In this paper, the milling parameters of high energy ball mill (Fritsch Pulverisette 7) like vial geometry, number and size of balls and speed of the mill were modelled and discussed. Simulations through discrete element method (DEM) provide correlation between the milling parameters. A mathematical model is used to improve and develop this process. The results show that the loss of powder mass can remarkably improve the performance of milling. The balls made of stainless-steel have a positive effect on the milling efficiency. The simulation shows that the high ball milling velocities can contribute to faster particle size reduction.

Keywords

Ball Mill Modelling, Discrete Element Method, Planetary Ball Mill, High-Energy Ball Milling

1. Introduction

High-energy ball milling is a complicated process employed in solid reactions for obtaining nanostructured materials, in powder form, with an average particle size of less than 100 nm. The planetary mill is one of high-energy ball mills, which is used for efficient and precision milling. Thus, elastic/plastic deformation occur during the collisions of the powder particles between two colliding balls [1].

However, extensive research [1] [2] [3] [4] [5] has been done on the modelling of mechanical alloying, principally based on thermodynamic and mechanical approaches. Among them is discrete element modeling (DEM), which is a powerful process in predicting the interactions and the motion of balls, the power used, the energy transferred, the media and wear of the liner, and mill output [6] [7] - [12]. In this context, Lu et al. [13] proposed and validated kinetic models for simulating the trajectory energy transfer of single milling ball in a planetary ball mill. These models illustrate the relationships between kinetic and dynamic parameters. Additionally, Govender et al. [14] validated the predictions of a DEM model by comparison with the positron emission particle tracking (PEPT) realized in a wet tumbling ball mill. They revealed the capabilities of this model to generate the movement of the solid charge. Similarly, Zidane et al. [15] founded an equation explaining the contact temperature of the powders at the impact point of the balls as a function of some dynamic, geometrical, and thermo-physical parameters basing on a mathematical model. Hence, it has been reported that the impact angle and the impact velocity have a significant effect on the energy transferred to the powder during high-energy ball milling [16].

The present work aims to propose a numerical approach using discrete element method to analyze the dynamic aspects of a planetary ball milling process.

2. Method

Following, the model of the Fritsch Pulverisette 7 [17] planetary ball mill is presented. The mills consist of two vials rotating both around its own axis and the line of symmetry of the plate. The movement of the balls interior the vial is operated by the resultant of both Coriolis and centrifugal forces, which generate impacts that remove shear and compressive forces to a charge on powder particles. Thus, these forces generate variation in structure and microstructure and mechanochemical transformations [18] [19]. In fact, at optimum milling conditions the trajectory ideal for a single ball follows a point on the vial circumference through the resulting field of forces acts to carry the ball across the vial to the opposite point [20]. In practice, the trajectories, which depend on the interactions between several forces and balls, are exchanged in every possible direction with respect to the impact reference frame.

2.1. Numerical Model and Simulation Conditions

2.1.1. Discrete Element Method Representation

Discrete element method is a computational method that is used for predicting the flow of particulates in circumstances where their collisions are the dominant physical process. The ball motion in planetary ball mill is examined by DEM modelling during high-energy ball milling of a titanium carbide powder. The effect of velocities of vial and plate on milling efficiency has been studied.

Figure 1 shows the high-energy planetary ball mill (Fritsch Pulverisette 7 classic line, Idar-Oberstein, Germany) equipped with two vials used in this work. It has an external structure made of steel, which isolates the rotating parts from the external environment for protecting the operator from the risk of accidents from detached parts.

Moreover, computer simulations are a useful method to explore the milling process, providing information about kinematic and dynamic quantities [21]. In

Figure 1. Photograph of the high-energy planetary ball mill.

this work, the high-energy planetary ball mill process is modeled using the software Automated Dynamic Analysis of Mechanical System MSC Adams [22]. Figure 2 illustrates a schematic of a high-energy planetary ball mill equipped with two cylindrical vials. The mill was powered by a 0.37 kW electric motor. The main geometric and operation parameters include: R_p is the distance between the global rotation axis and the rotation axis of the vial; R_v is the radius of the vial cylinder; ω and Ω are the angular velocities of the vial and plate, respectively. The direction of the vials revolution was set counter to that of the plate revolution. The geometrical and physical properties of the vial and ball used for mill are given in Table 1.

Figure 2. (a) Photograph of the two vials on the supporting disc and (b) schematic representation of the experimental planetary ball mill.

Table 1. Geometrical and physical properties of the vial and ball used for mill.

The velocities used in the force-displacement law are obtained based on Newton’s second law of motion applied to disc x:

$m_{\left(x\right)}\stackrel{••}{x_i}=m\frac{\text{d}{\nu }_{\left(x\right)}}{\text{d}t}=\sum F_{\left(x\right)i}$ (1)

$I_{\left(x\right)}\stackrel{••}{{\theta }_i}=\frac{\text{d}{\omega }_{\left(x\right)}}{\text{d}t}=\sum M_{\left(x\right)}$ (2)

where $m_{\left(x\right)}$ represents the mass of a particle, ${\nu }_{\left(x\right)}$ is the translational velocity, ${\omega }_{\left(x\right)}$ is the angular velocity, $I_{\left(x\right)}$ is the moment of inertia of disc x and ( $•$ ) denotes a time derivative.

In this work, the Hertz-Mindlin force model was used to resolve particles interactions [23] [24] [25]. We used the commercially available software Experts in Discrete Element Modeling (EDEM 2.4) [26].

Therefore, the total force acting on particle i is obtained by taking the summation of the above forces:

$F=F_N+F_T$ (3)

where $F_N$ is the normal force component and $F_T$ is the tangential force.

$F_N=\frac{4}{3}{E}^{*}\sqrt{R^{*}}{\delta }_N^{3/2}-2\sqrt{\frac{5}{6}}\frac{\ln e}{{\ln }^{2}e+{\pi }^2}\sqrt{2E^{*}{m}^{*}}{\left(R^{*}{\delta }_N\right)}^{1/4}{\nu }_N^{rel}$ (4)

where ${\delta }_N$ is the normal particle overlap, e is the restitution coefficient, ${\nu }_N^{rel}$ is the normal component of the relative velocity between two particles, $E^{*}$, $R^{*}$ and $m^{*}$ are the equivalent Young’s modulus, equivalent radius and equivalent mass of a single grinding ball, respectively, which are defined as following:

$E^{*}={\left(\frac{1-{\nu }_a^2}{E_a}+\frac{1-{\nu }_b^2}{E_b}\right)}^{-1}$ (5)

$R^{*}={\left(\frac{1}{R_a}+\frac{1}{R_b}\right)}^{-1}$ (6)

and

$m^{*}={\left(\frac{1}{m_a}+\frac{1}{m_b}\right)}^{-1}$ (7)

where $\nu$, E, R, m are the Poisson’s ratio, Young’s modulus, radius, mass of a particle a in contact with particle b, respectively.

The tangential force is a function of the elastic force and damping force.

$F_T=-8G^{*}\sqrt{R^{*}{\delta }_N}{\delta }_T-2\sqrt{\frac{5}{6}}\frac{\ln e}{{\ln }^{2}e+{\pi }^2}\sqrt{8G^{*}{m}^{*}}{\left(R^{*}{\delta }_N\right)}^{1/4}{\nu }_T^{rel}$ (8)

where $G^{*}$ is the equivalent shear modulus calculated from the following Equations:

$G^{*}={\left(\frac{1-{\nu }_a^2}{G_a}+\frac{1-{\nu }_b^2}{G_b}\right)}^{-1}$ (9)

$G=\frac{E}{2\left(1+\nu \right)}$ (10)

The torque applying on a particle due to contact with a surface is given by:

${\tau }_a=-{\mu }_R{F}_N{R}_a{\omega }_a$ (11)

where ${\mu }_R$ is the coefficient of rolling and ${\omega }_a$ is the angular velocity.

Table 2 summarizes the simulation parameters and their values used in this study.

Table 2. Simulation parameters used in this study [4].

2.1.2. Ball-Particle Collisions

The model proposed by Barrios et al. [27] assumes that the particles can be modeled as a monolayer bed. If the particles have spherical shapes and arranged according to a dense hexagonal packing, the number of particles trapped in each collision as a function of radius could be estimated by:

$N_{cap,i}=\frac{1}{4}+\frac{3}{4}{\left(\frac{2r_c}{d_i}\right)}^2$ for $r_c\ge \frac{d_i}{2}$ (12)

$N_{cap,i}=1$ for $r_c<\frac{d_i}{2}$

where r_c is the radius of the bed and d_i is the mean size of the particles caught by the colliding balls.

They also proposed that the radius of capture can be obtained by adding the radius of contact due to the elastic deformation during the contact, and the radius of contact caused by the geometry of media (balls) and particle powders.

$r_c={\left(\frac{15E}{8K_e}\right)}^{1/5}{K}_g^{2/5}+{\left[K_g^2-{\left(K_g-h\frac{\Delta }{h}\right)}^2\right]}^{1/2}$ (13)

where K_e is the elastic constant of the contact is given by the following equation [28]:

$K_e=\frac{K_1{K}_2}{K_1+K_2}$ (14)

and K_g is the geometric constant of the contact given by the equation:

$K_g=\frac{r_{c1}{r}_{c2}}{r_{c1}+r_{c2}}$ (15)

K_i is elastic stiffness of each of the bodies in contact, r_ci is the radius of curvature and Δ is the maximum deformation of the bed during collision.

Therefore, in the case of monolayer bed the initial bed height may be estimated equal to the mean size of the particles caught by the colliding balls, d_i.

Thus, the ratio Δ/h represents the maximum relative deformation of the bed [27], given by:

$\frac{\Delta }{h}=0.362\left[1-\exp \left(-0.152{\left(\frac{E}{E_{50b}}\right)}^{0.413}\right)\right]$ (16)

This model includes grinding media (two balls) and particles (TiC powder) since the DEM software records only collisions and contact involving two bodies. In this approach, only ball-particle collisions are extracted (see Figure 3). Thus, the presence of the particles represents the powder captured between grinding media. When the powder particles (TiC) were milled by large balls (the

Figure 3. Schematic illustration of ball to powder particles collision during milling.

case of our study), only some particles are collided with the balls because of the difference in size between balls and powder particles.

The most important ingredient in the model of a ball mill is the contact law. In this context, Gilardi et al. [29] have reported the existence of both models (discrete and continuous).

The contact force is defined by the hard-coded impact function [22], based on the non-linear spring (F_k) and linear damper scheme (F_d) proposed by Dubowsky and Freudenstein [30] from following relation:

$F_c=F_k+F_d=k{u}^n+c\stackrel{\dot{}}{u}$ (17)

where u and $\stackrel{\dot{}}{u}$ the relative displacement and velocity of the colliding bodies, k and c are the spring generalized stiffness and the damping coefficient. With respect to the above formulation, to prevent discontinuities, the impact function, defined within MSC Adams, implements a damping coefficient which depends on the relative displacement of colliding bodies, given by:

$c=-c_m{\left(\frac{u}{d}\right)}^2\left(3-2\frac{u}{d}\right)$ (18)

c varies from zero to cm, assuming the latter value when the relative displacement is greater than or equal to d. With the exponent n = 3/2, an estimate of k can be derived from the Hertzian theory of contact [31], as given in the following relation:

$k=\frac{4}{3\pi }\left(\frac{1}{{\eta }_i+{\eta }_j}\right){\left(\frac{R_i{R}_j}{R_i+R_j}\right)}^{1/2}$ (19)

for the case of contact between sphere i and j, with radius R, and

$k=\frac{4}{3\pi }\left(\frac{1}{{\eta }_i+{\eta }_j}\right){\left(R_i\right)}^{1/2}$ (20)

for the case of contact between sphere i and a plane.

In fact, the parameter, accounts for elastic properties of materials, is given by the following equation:

${\eta }_i=\frac{1-{\nu }_i^2}{\pi E_i}$ (21)

where n is the Poisson ratio and E represents the Young modulus. The solution of the classical problem of damped vibration of a mass spring system provides a guess of maximum the damping parameter [30],

$c_m=-2\ln \left(e\right){\left[\frac{1}{{\ln }^2\left(e\right)+{\pi }^2}\right]}^{1/2}{\left(m_e{k}^{\ast }\right)}^{1/2}$ (22)

being e = (h/h₀)^1/2 the restitution coefficient, accounting for the difference from initial (h₀) and final (h) height in a drop test, 1/m_e = 1/m_i + 1/m_j and k^* [32]:

$k^{*}={\left(\frac{320}{81}{m}_e{\stackrel{\dot{}}{u}}_{n0}{E}_e^4{R}_e^2\right)}^{1/5}$ (23)

Geometrical and elastic properties are expressed, respectively, as 1/R_e = 1/R_i + 1/R_j and 1/E_e = 1/E_i + 1/E_j and ${\stackrel{\dot{}}{u}}_{n0}$ is the relative velocity during the collision.

2.2. Parameters of Planetary Ball Mill

The planetary ball mill process involves loading of the powders with the milling medium in a vial with severe deformation. The main parameters of planetary ball mill are [33]:

- Type of mill (mixer or planetary);

- Height and volume of vial;

- Operation: Milling duration, rotation speed, kinetic shock energy, shock frequency;

- Number, size and type of milling balls, ball to powder mass ratio;

- Other parameters: temperature of milling, milling atmosphere (e.g. argon) and process control agent.

Moreover, all these variables are not completely independent. For example, the optimum milling time depends on the type of ball mill, temperature of milling, ball-to-powder ratio, etc. However, other parameters have no significant effect on the final product obtained after milling [33].

In this work, the rotation speed, the number of the balls, the diameter of the balls, and the time of milling were varied.

2.3. Experimental

Elementary Ti (<40 μm, 99.9%) and C (5 μm, 99.9%) powder mixture was sealed into a stainless-steel vial with 5 stainless-steel balls (15 mm in diameter) in a glove box filled with purified argon to avoid oxidation. The ball to powder weight ratio was 70:1. The milling process was performed at room temperature using a high-energy ball mill (Fritsch Pulverisette-7 planetary mill).

The crystalline properties of the powders were identified through X-ray diffraction technique using a Panalytical XPERT PRO MPD diffractometer with Cu K_α radiation. The existing phases were resolved by the HighScore Plus program based on the ICDD PDF2. The morphology of the milled powder was carried out using an FEI Quanta 200 environmental scanning electron microscope.

3. Results and Discussion

3.1. XRD and SEM Study

Figure 4 provides the X-ray diffraction patterns of the elemental powders of Ti and C and the elemental mixture of Ti and C powders. All the expected peaks for all the as-received powders are present.

Figure 5 illustrates SEM micrograph of TiC powder milled for 20 h. The powder is characterized by fine agglomerated particles having homogenous shape and with average particle size less than 1 µm. The state is marked by particle welding predominance.

Figure 4. XRD pattern of the elemental powder of Ti, the elemental powder of C and the elemental mixture of Ti and C powders.

Figure 5. SEM micrograph of TiC powder milled for 20 h.

3.2. The Influence of Friction on Milling Process

The effect of the friction conditions for the ball motion on the planetary ball milling is very important. Thus, we study the effects of mill feed by altering the friction coefficients in the DEM simulation. As shown in Figure 6, the simulated milling ball motion as well as the location and shape of the bulk of balls change with increasing both friction coefficients. Therefore, the static and rolling friction coefficients have a remarkable effect on the ball motion. For the initial simulation (μ_S = 0.2 and μ_R = 0.05) the balls motion is very low and the powder (small sizes) and balls (big sizes) are distributed over a small fraction in inner of the vial. When we have increased the two coefficients (μ_S = 0.5 and μ_R = 0.5) the powder and the balls are distributed over a larger fraction in inner of the vial. In addition, the transfer of kinetic energy (in the tangential direction) during the collision of balls becomes more significant with increasing friction coefficients, which is similar to the experimental results. These results are also obtained by Zeng et al. [34] who studied the effects of static friction coefficients (between particle and cylinder sieve wall-and between particles) on uniformity of motion in a horizontal rice mill.

Figure 6. Effect of friction coefficients on ball motion in a planetary ball mill.

3.3. The Influence of Speed Ratio on Milling Process

The speed ratio, k, is an important factor that influences milling ball motion. In this work, the speed ratio was varied between −3 and −1. To study the influence of rotational speed, the filling ratio of powders and balls were kept constant at 20 and 15%, respectively. Planetary ball mill partially filled by spherical particles powder was simulated. Figure 7 shows the outcome of selected simulations, illustrating trajectories of the balls for different motion regimes. For the first case, motion is not random and space regions of higher occupancy are present. For

Figure 7. Schematic showing the typical motion regime of the balls in the vial: (a) cascading, (b) cataracting, and (c) rolling.

the second case, the motion described within the second column is more random, with trajectories across the whole vial. The last case, corresponding to high vial velocity, as clearly shown by trajectories, is characterized by balls chiefly sticking to the vial surfaces.

Figure 8 shows the influence of varied speed ratio on milling ball motion (k = −3, k = −2.5, k = −2, and k = −1). The results of the DEM model show that no visible difference in milling balls motion pattern when the speed ratio is varied. Thus, for both cases (k = −3 and k = −1) cascading of the balls (in read) was observed. In this study, the friction was not considered. DEM model is in agreement with the experimental observations.

Figure 8. Influence of varied speed ratio on milling ball motion (in red color) in the presence of TiC powder (in green and yellow colors).

4. Conclusion

Mechanical alloying process through high-energy planetary ball mill is complicated and reliant on the conducted production, the milled materials, and the desired final product. We validate that DEM is able to simulate the dynamics of planetary ball milling process. During this approach, we show that the effects of any change of operation conditions on the dynamics of the system can be also studied. DEM modelling clarifies the efficiency of mechanical alloying process for low velocities ratios. The milling includes more movement of the balls in the inner and possibly very powerful interactions. Finally, further study on several aspects requires to be followed.

Acknowledgements

We would like to express our special appreciation and thanks to Ms. Melissa Luo, AMPC Editorial Office in Advances in Materials Physics and Chemistry, for her guidance and support.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Maurice, D. and Courtney, T.H. (1994) Modeling of Mechanical Alloying: Part I. Deformation, Coalescence, Band Fragmentation Mechanisms. Metallurgical and Materials Transactions A, 25, 147-158. https://doi.org/10.1007/BF02646683
[2]	Bhattacharya, A.K. and Arzt, E. (1992) Temperature Rise during Mechanical Alloying. Scripta Metallurgica et Materialia, 27, 749-754. https://doi.org/10.1016/0956-716X(92)90500-E
[3]	Murty, B.S., Mohan, M.D., Mohan Rao, M. and Ranganathan, S. (1995) Milling Maps and Amorphization during Mechanical Alloying. Acta Metallurgica et Materialia, 43, 2443-2450. https://doi.org/10.1016/0956-7151(94)00402-1
[4]	Rosenkranz, S., Breitung-Faes, S. and Kwade, A. (2011) Experimental Investigations and Modelling of the Ball Motion in Planetary Ball Mills. Powder Technology, 212, 224-230. https://doi.org/10.1016/j.powtec.2011.05.021
[5]	Zhao, X. and Shaw, L. (2017) Modeling and Analysis of High-Energy Ball Milling through Attritors. Metallurgical and Materials Transactions A, 48, 4324-4333. https://doi.org/10.1007/s11661-017-4195-6
[6]	Ward, T.S., Chen, W., Schoenitz, M., Dave, R.N. and Dreizin, E.L. (2005) A Study of Mechanical Alloying Processes Using Reactive Milling and Discrete Element Modeling. Acta Materialia, 53, 2909-2918. https://doi.org/10.1016/j.actamat.2005.03.006
[7]	Jiang, X., Trunov, M.A., Schoenitz, M., Dave, R.N. and Dreizin, E.L. (2009) Mechanical Alloying and Reactive Milling in a High Energy Planetary Mill. Journal of Alloys and Compounds, 478, 246-251. https://doi.org/10.1016/j.jallcom.2008.12.021
[8]	Santhanam, P.R. and Dreizin, E.L. (2012) Predicting Conditions for Scaled-Up Manufacturing of Materials Prepared by Ball Milling. Powder Technology, 221, 403-411. https://doi.org/10.1016/j.powtec.2012.01.037
[9]	Santhanam, P.R., Ermoine, A. and Dreizin, E.L. (2013) Discrete Eement Model for an Attritor Mill with Impeller Responding to Interactions with Milling Balls. Chemical Engineering Science, 101, 366-373. https://doi.org/10.1016/j.ces.2013.06.048
[10]	Feng, Y.T., Han, K. and Owen, D.R.J. (2004) Discrete Element Simulation of the Dynamics of High Energy Planetary Ball Milling Processes. Materials Science and Engineering A, 375-377, 815-819. https://doi.org/10.1016/j.msea.2003.10.162
[11]	Yang, R.Y., Jayasundara, C.T., Yu, A.B. and Curry, D. (2006) DEM Simulation of the Flow of Grinding Media in Isa Mill. Minerals Engineering, 19, 984-994. https://doi.org/10.1016/j.mineng.2006.05.002
[12]	Cleary, P.W. and Morrison, R.D. (2011) Understanding Fine Ore Breakage in a Laboratory Scale Ball Mill Using DEM. Minerals Engineering, 24, 352-366. https://doi.org/10.1016/j.mineng.2010.12.013
[13]	Lu, S.Y., Mao, Q.J., Peng, Z., Li, X.D. and Yan, J.H. (2012) Simulation of Ball Motion and Energy Transfer in a Planetary Ball Mill. Chinese Physics B, 21, 566-574. https://doi.org/10.1088/1674-1056/21/7/078201
[14]	Govender, I., Cleary, P.W. and Mainza, A.N. (2013) Comparisons of PEPT Derived Charge Features in Wet Milling Environments with a Friction-Adjusted DEM Model. Chemical Engineering Science, 97, 162-175. https://doi.org/10.1016/j.ces.2013.04.023
[15]	Zidane, D., Bergheul, S., Rezoug, T. and Hadji, M. (2011) Study of the Temperature Variations of Al and Ti in Mechanical Alloying. Journal of Materials Engineering and Performance, 20, 1109-1113. https://doi.org/10.1007/s11665-010-9721-2
[16]	Gaffet, E. (1991) Planetary Ball-Milling: An Experimental Parameter Phase Diagram. Materials Science and Engineering: A, 132, 181-193. https://doi.org/10.1016/0921-5093(91)90374-V
[17]	Fritsch Planetary Micro Mill PULVERISETTE 7. http://www.fritsch-milling.com/products/milling/planetary-mills/pulverisette-7-classic-line/ description
[18]	Harris, K.D.M. (2013) Mechanochemical Synthesis: How Grinding Evolves. Nature Chemistry, 5, 12-14. https://doi.org/10.1038/nchem.1539
[19]	Frišvcić, T., Halasz, I., Beldon, P.J., Belenguer, A.M., Adams, F., Kimber, S.A.J., Honkimäki, V. and Dinnebier, R.E. (2013) Real-Time and in Situ Monitoring of Mechanochemical Milling Reactions. Nature Chemistry, 5, 66-73. https://doi.org/10.1038/nchem.1505
[20]	Burgio, N., Iasonna, A., Magini, M., Martelli, S. and Padella, F. (1991) Mechanical Alloying of the Fe-Zr System. Correlation between Input Energy and End Products. IL Nuovo Cimento D, 13, 459-476. https://doi.org/10.1007/BF02452130
[21]	Broseghini, M., D’Incau, M., Gelisio, L., Pugno, N.M. and Scardi, P. (2016) Effect of Jar Shape on High-Energy Planetary Ball Milling Efficiency: Simulations and Experiments. Materials & Design, 110, 365-374. https://doi.org/10.1016/j.matdes.2016.06.118
[22]	MSC Adams. http://www.mscsoftware.com/it/product/adams
[23]	Hertz, H. (1882) On the Fixed Elastic Body Contact. Journal für die reine und angewandte Mathematik (Crelles Journal), 92, 156-171. https://doi.org/10.1515/crll.1882.92.156
[24]	Mindlin, R.D. (1949) Compliance of Elastic Bodies in Contact. Transactions of the ASME, Journal of Applied Mechanics, 16, 259-268.
[25]	Tsuji, Y., Tanaka, T. and Ishida, Y. (1992) Lagrangian Numerical Simulation of Plug Flow of Cohesion Less Particles in a Horizontal Pipe. Powder Technology, 71, 239-250. https://doi.org/10.1016/0032-5910(92)88030-L
[26]	http://tm.spbstu.ru/images/2/28/EDEM2.4_user_guide.pdf
[27]	Barrios, G.K.P., Carvalho, R.M. and Tavares, L.M. (2011) Modeling Breakage of Mono Dispersed Particles in Unconfined Beds. Minerals Engineering, 24, 308-318. https://doi.org/10.1016/j.mineng.2010.09.018
[28]	Tavares, L.M. and King, R.P. (2004) Measurement of the Load-Deformation Response from Impact-Breakage of Particles. International Journal of Mineral Processing, 74, S267-S277. https://doi.org/10.1016/j.minpro.2004.07.017
[29]	Gilardi, G. and Sharf, I. (2012) Literature Survey of Contact Dynamics Modelling. Mechanism and Machine Theory, 37, 1213-1239. https://doi.org/10.1016/S0094-114X(02)00045-9
[30]	Dubowsky, S. and Freudenstein, F. (1971) Dynamic Analysis of Mechanical Systems with Clearances Part 1: Formation of Dynamic Model. Journal of Engineering for Industry, 93, 305-309. https://doi.org/10.1115/1.3427895
[31]	Goldsmith, W. (2001) Book on the Impact: The Theory and Physical Behaviour of Colliding Solids. Dover Publications, Mineola, 416 p.
[32]	Di Maio, F.P. and Di Renzo, A. (2004) Analytical Solution for the Problem of Friction Al-Elastic Collisions of Spherical Particles Using the Linear Model. Chemical Engineering Science, 59, 3461-3475. https://doi.org/10.1016/j.ces.2004.05.014
[33]	Suryanarayana, C. (2001) Mechanical Alloying and Milling. Progress in Materials Science, 46, 1-184. https://doi.org/10.1016/S0079-6425(99)00010-9
[34]	Zeng, Y., Jia, F., Meng, X., Han, Y. and Xiao, Y. (2018) The Effects of Friction Characteristic of Particle on Milling Process in a Horizontal Rice Mill. Advanced Powder Technology, 29, 1280-1291. https://doi.org/10.1016/j.apt.2018.02.021