Numerical Simulation of Mechanical Alloying in a Shaker Mill by Discrete Element Method

Modeling of Mechanical Alloying (MA), which is a solid-state powder processing technique, is carried out by examining one widely used laboratory scale milling device, the SPEX 8000 shaker mill. It is a vibratory mill; its vial is agitated at a high frequency in a complex cycle that involves motion in three orthogonal directions. In this work, a popular dynamic simulation technique, Discrete Element Modeling, is applied to examine dynamics of a SPEX 8000 shaker ball mill based on the movement of milling balls. The computational results for energy dissipation rate inside the mill are calculated for dif ferent ball sizes and varied total ball to powder mass ratios (charge ratios). The computational results are well correlated with the experimental results tracking milling dose (used to define the degree of milling) as a function of ball sizes and charge ratios. Moreover, the numerical (theoretical) milling dose that correlates well with its experimental analog was found to depend on the energy dissipation rate of the head-on ball collisions. The numerical simulations also indicated that the milling progress is most significantly affected by milling media collisions with the energy within a specific threshold, while the collisions with smaller and greater energies are less ef fective. Finally, discussion shows how this novel approach of correlating specific scaling terms between experiments and simulations can be applied to other powder processing equipment.


Introduction
Mechanical Alloying (MA) is a solid-state, highenergy ball milling technique used to produce powders with unique microstructures.Mechanical Alloying process has been widely employed in industry to synthesize a variety of commercially useful and scientifically interesting materials 1) .Generally, in MA process, a mixture of powder is loaded into a high-energy mill along with a suitable grinding medium (milling balls).Powder particles trapped between colliding media are constantly subjected to deformation, resulting particle f lattening, coalescence (cold welding) and/or fragmentation.
MA is a complex process.Aspects of the events that occur in MA process have been described qualitatively 2) .However, modeling the MA process is a difficult task, since a large number of variables need to be considered.Despite its complexity, modeling of the process has been actively pursued during the last 10 Ҁ 15 years, e.g., 3) .In general, modeling of MA approaches can be classified in either local or global terms.Local modeling describes effects (thermal and mechanical) and events (deformation, fracture and welding) when powder particles are entrapped between colliding or sliding media surfaces [4][5][6] , while global modeling considers aspects such as distribution of impact velocity, impact angles as well as heterogeneity of powder within the mill 7) .Attempt has been made to incorporate local and global approaches to examine mechanics of MA process 8) .While mod-W. Chn, M. Schoenitz, T.S. Ward

Abstract
Modeling of Mechanical Alloying (MA), which is a solid-state powder processing technique, is carried out by examining one widely used laboratory scale milling device, the SPEX 8000 shaker mill.It is a vibratory mill; its vial is agitated at a high frequency in a complex cycle that involves motion in three orthogonal directions.In this work, a popular dynamic simulation technique, Discrete Element Modeling, is applied to examine dynamics of a SPEX 8000 shaker ball mill based on the movement of milling balls.The computational results for energy dissipation rate inside the mill are calculated for different ball sizes and varied total ball to powder mass ratios (charge ratios).The computational results are well correlated with the experimental results tracking milling dose (used to define the degree of milling) as a function of ball sizes and charge ratios.Moreover, the numerical (theoretical) milling dose that correlates well with its experimental analog was found to depend on the energy dissipation rate of the head-on ball collisions.The numerical simulations also indicated that the milling progress is most significantly affected by milling media collisions with the energy within a specific threshold, while the collisions with smaller and greater energies are less effective.Finally, discussion shows how this novel approach of correlating specific scaling terms between experiments and simulations can be applied to other powder processing equipment.
erate success has been achieved in modeling the mechanics of the MA process based on a simplified mill 7) , the mechanics of MA remain incompletely understood.In the present work, a novel approach is considered where a parametric description of the milling progress of mechanical alloying in SPEX 8000 shaker mill is investigated in a manner such that experimental and theoretical results can be correlated.This paper describes the numerical modeling method used to develop theoretical description of the milling progress.
The SPEX 8000 series shaker mill is widely used in laboratories throughout the world.It allows relatively rapid milling of small (20 g) quantities of material, and has become a de-facto standard piece of equipment.Focusing on this type of mill will enable subsequent comparison of the computations with experimental results published by other researchers.Validation of the numerical approach for this type of mill is the first step required for transition to largerscale industrial equipment.
The SPEX 8000 shaker mill is a vibratory mill, where its vial is agitated at a high frequency in a complex cycle that involves motion in three orthogonal directions.Fig. 1 shows the schematic of a SPEX 8000 mill device that is used in our lab.The center of vial vibrates in two-dimensional mode with the same frequency and different amplitude, and its slanted axis rotates around the third direction.The movement of rotation and vibration has the same frequency.In this paper, a popular dynamic simulation technique, Discrete Element Modeling, is applied to 3D simulation of the SPEX 8000 shaker ball mill based on the interactions of balls with other balls as well as the vial boundaries.The purpose of the study is to examine dynamic impact inside this milling device and develop a correlation between the numerical results and the experimental data.As it is known, the driving force for the milling process is the countless impacts between ball and ball and ball and vial.Energy dissipation during the impact directly contributes to any changes in the milled powder, as powder is trapped between impacting surfaces.The cumulative work performed on the powder is described by the milling dose D m , quantified later in the paper.Therefore, through detailed simulations, energy dissipation rate inside the system is calculated and based upon that a numerical (i.e., theoretical) milling dose is defined.Calculations are carried out for different ball sizes and total ball-to-powder mass ratios (ball mass, m b , to powder mass, m p , the ratio is also called the charge ratio, C R ).This numerical milling dose is correlated with a suitably determined experimental milling dose to evaluate the milling progress.

Description of the Numerical Model
Most particle dynamic simulations use the Discrete Element Method (DEM).The term "discrete element" refers to the fact that the simulation models the particles as a system of individual element (called discrete elements) 9) .In other words, this technique simulates systems consisting of discrete particles in which motion of each individual particle is controlled by its interactions with other particles and with the system boundaries, which are also treated as discrete elements.Therefore, in a DEM simulation, discrete elements include particles and boundaries.For numerical simulation of the SPEX 8000 system by DEM, the boundary of the system is disassembled into three elements; one is represented as cylindrical boundary and the other two are represented as bounded top and bottom flat surfaces.In the modeling, the movements of boundary elements are described according to the actual motion of the SPEX 8000 vibrating mill during the experiments.Referring to Fig. 1, these include rotation of the vial, which is itself tilted around z-axis, and linear vibratory motion of vial center in x and y directions.As per its design and subsequent operating conditions, the value of rotation frequency is the same as that of vibration.Hence, the vibration movement in x and y directions can be expressed by the following simple equations, Ά x҃A x sin ω t (1) y҃A y cos ω t where, A x and A y are the amplitudes of motion in x and y directions respectively, and ω is the frequency of vibratory motion.
For dynamic simulation of a particulate system by DEM technique, the computation time is a critical issue, which is governed by the number of particles and size of the smallest particle in the system.In a given experiment, the physical size of powder to be mechanically alloyed is very small, and the size varies significantly, while the number of particles is very large.This poses a formidable computational challenge, and hence it is impractical to include individual powder particles in the modeling.The main purpose of this work is to compute useful information based on ball-ball and ball-boundary interactions in order to understand the mechanics of the SPEX 8000 mill.Therefore, it is not necessary to exactly model all the powder particles which are being milled, and only the milling media balls along with the vial motion are modeled.] , there is always some effect of powder on ball-ball and ball-boundary interactions, and therefore it is taken into consideration by allowing a change in certain properties of the ball particles.This is done based on the strategy suggested by Kano et al., through proper selection of restitution coefficient 10) .Moreover, previous systematic study has shown that the impact velocity, thickness and strength of the powder layer coating the balls and the ball size contribute to the variation in restitution coefficient 11) .The trends reported earlier were supported by our experimental evaluation of restitution coefficient between the ball and vial surface for a clean ball versus a milled powder coated ball 12) .Therefore, in our study (here as well as in Ref. 12) ), it is assumed that under the same amount of powder loading, restitution coefficient is inversely proportional to surface coverage of milling balls and changes within the range of 0.5 Ҁ 0.8 under different operating conditions.In our numerical simulation, Walton-Braun soft sphere model [13][14][15] is applied to ball-ball and ballboundary interactions. Te detailed description of the force model and the numerical implementation can be found in 16) .
The time step, ∆t, is an important quantity that affects the total computational time required for the simulation, and is calculated by the following expression [13][14][15] , ∆t҃ where e is the restitution coefficient of the milling balls, m i is the mass of one ball, K 1 is the spring stiffness during the loading (i.e., approach during the contact) and n is the desired number of time steps for one contact.This number is typically chosen to be between 20 and 60, and for this simulation n҃40.Equation (2) shows that the time step of the calculation is proportional to the ball mass.
The normal stiffness of the milling balls K 1 has been estimated to fall within the recommended guidelines [13][14][15] , allowing for a maximum deformation of 1 % of the ball diameter during a collision. Te normal stiffness has thus been calculated using the following equation, where E is Young's modulus and r i is the ball radius.

Simulation Approach
As mentioned before, because of the large difference between the particle size of the powder being milled and the size of the milling balls, only milling balls were considered.This simplification is reasonable; it allows focusing on the high-energy interactions between milling balls as well as milling ball and milling vial walls, which are responsible for mechanical alloying in a ball mill.During the simulation, energy dissipation in the system, which arises from ball-ball and ball-wall collisions, is calculated in the term of the energy dissipation rate by the following equation, where E 1 and E 2 are the energies of a binary impacting system before and after a collision, respectively; k is the collision index, and N c represents the total number of collisions during the time interval t s .For a ballball collision, the energy of the impacting system is calculated by: where subscripts 1, 2 refer to the impact energy before (E 1 ) and after the collision (E 2 ), respectively; i, j are indices of two colliding balls, v is the translational velocity, and w is the rotational speed of each impacting ball; m is the mass and I is the moment of inertia of the ball; and E p is the potential energy of the ball under gravity (E p ҃mgh, where h is height of mass center, and g is the gravitational acceleration); and for a ball-boundary interaction, where Dv i is the relative impact velocity between the ball i and the boundary it is colliding with.Table 1 lists the parameters and their values or variation ranges used in the simulation and in the corresponding experimental study 12) .The parameters related to the milling vial geometry and movement come from empirical measurements for a specific SPEX 8000 series shaker mill used in experiments 12,17) .Milling ball sizes and loadings are also the same as for the corresponding experimental investigation.Milling balls used in the study are made of hardened steel (density, r҃7.86g/cm 3 ) with four different sizes.

Numerical Results
Numerical simulations are carried out to represent closely the experimental studies of the effect of milling ball size and ball loading, represented by C R (total ball mass to powder mass, termed as the charge ratio), on the milling time.According to the experimental study, 1 2 1 2 three values of C R , i.e., C R ҃2.5, 5 and 10 are used for each ball size.As stated in the previous sections, because of the computational limitation, powder particles are not considered directly in the numerical analysis.However, the effect of powder to ball-ball and ball-wall interactions is taken into consideration by the appropriate selection of the restitution coefficient.It is assumed that under the same powder loading, the restitution coefficient is inversely proportional to the surface coverage of the milling balls, varying in the range of 0.5 Ҁ 0.8.For the milling process, the required milling time is correlated with the impact energy consumption rate, i.e., the energy dissipation rate.The larger the energy consumption rate, the shorter will be the milling time required to achieve the same state in the milled powder.The dynamics of the SPEX 8000 shaker mill is studied in detail in terms of the distribution of energy dissipation per impact and the impact angle.Fig. 2 represents the histograms of impact energy dissipation inside the system corresponding to four different ball sizes with C R ҃5.0, and m p ҃2 g.The vertical scale is the fraction of all impacts occurring at a given energy dissipation level.For the system loaded with the smallest balls considered (d҃2.36 mm), the value of energy dissipation for most impacts is less than 10e-5 J.For the two ball diameters, d҃3.16 mm and d҃4.76 m, most impacts involve energy dissipation at the similar energy levels that do not exceed 10e-2 J.For the system loaded with the largest ball size (d҃9.52mm), some impacts involve the energy levels exceeding 10e-2 J. Fig. 3 illustrates distribution of the dissipated energy as a function of collision angles for the systems with different milling ball sizes.Generally, for all ball sizes, most collisions occur at an angle that is in range 80 Ҁ 90°.Such impacts can be classified as oblique or glancing impacts.However, a trend can be noticed that with the increased ball size, the fraction of the energy dissipated in the collisions occurring at smaller angles increases.For the largest milling balls (d҃9.52 mm), the fraction of the energy dissipation from the collisions occurring at the angles less than 30°(which can be called head-on collisions) approaches 10%.These observations made through numerical simulations show that the processing of powder in a SPEX 8000 mill involves a large fraction of glancing impacts.Relatively few impacts are accompanied with high levels of energy dissipation.These conclusions are in good agreement with previously observed experimental results 2) , suggesting that this numerical simulation strategy can capture the experimental phenomena reasonably well.

Size of vial (mm) f38҂57
Angle of vial axis with rotating axis (°) Discussing further the results of the numerical simulations, one can consider the frequency of collisions as a function of the number of balls loaded in the milling vial for different ball sizes.The results are plotted in Fig. 4. It is observed that the collision frequency inside the SPEX 8000 mill is approximately proportional to the number of balls.This trend remains valid for different ball sizes, allowing some generalization of this result followed from the numerical modeling.Fig. 4 also shows that the number of collisions per second is quite large, and that the collision frequency remains linearly related to the number of balls while the number of balls increases in a wide range, as illustrated by the log-log plot.Although it is not shown explicitly, the number of collisions per ball increases dramatically as the total number of balls within a system increases.While this is a useful result, it is unclear whether an increase in the number of balls within the same size milling vial would be advantageous for the milling progress.Hence, there is a need to develop parametric information that can allow generalization as well as consolidation of the results for different operating parameters such as ball size and number of balls, and can allow for an easy comparison with the experimental results.
One such function that can be formulated and subsequently studied through simulations is proposed in terms of the energy dissipation rate, E d , and the charge ratio, C R , of the system, and is expressed as 12) .The powder mass, m p is incorporated in the function as a scaling factor to correlate with the experimental study and characterize the energy dissi-

Parametric Investigation of Mechanical Alloying
Recent experimental studies 12,17) on reactive milling using a SPEX 8000 mill suggested that the milling progress can be described using the specific milling dose, D m given by: where W is the work performed on the powder, and m p is the powder mass.
The studies involved milling of reactive materials such as thermites, where a self-sustaining reaction is triggered after a certain period of milling.For the purpose of the comparison of experimental and numerical results, it was assumed that the self-sustaining reaction is triggered when a specific degree of refinement of the powder components was achieved.The work performed on the powder in the milling system can be expressed through the energy dissipation rate, E d , i.e., W҃E d t init , where t init is the milling time required to initiate (or trigger) the reaction 12,18) .This leads to, From DEM simulations, it was observed that E d is proportional to the mass of the milling balls m b , hence, This equation provides an easy method for experimental measurements of the milling dose for a specific milled material.Based on the experimental results for milling of 2AlѿFe 2 O 3 thermite 12,17) , the product of the measured milling time (leading to initiation of the selfsustained reaction) and the charge ratio, D m ҃C R t init , is plotted in Fig. 7 as a function of the ball size for different charge ratios.The plots of C R t init for different C R values are similar to one another, especially for the intermediate ball sizes, and the overall value of D m does not change significantly.Therefore, the experimentally determined milling dose shows the behavior that is qualitatively similar to its numerically calculated equivalent, cf.Fig. 5.
However, the experimentally determined milling dose is not always constant, particularly for smaller and larger ball sizes.While the deviations may be explained qualitatively by collective ball motions or vial size limitations, further details could be determined from the numerical simulations.For instance, in the numerical simulation, the effect of the collision impact angle on the milling progress can be readily examined.Following the initial analysis presented in Fig. 3, the collisions occurring at different impact angles can be classified, in the first approximation, as head-on impacts (the impact angle is less than 30°) and glancing impacts (the impact angle is greater than 30°).Therefore, the energy dissipation rates for each type of collisions can be computed and respective milling progress functions can be considered and  compared to the experimental data shown in Fig. 7.
The sum of the head-on and glancing energy dissipation rates is the total energy dissipation rate, shown in Figs. 5 and 6.In contrast, Fig. 8 illustrates numerical milling dose, D' m based on the energy dissipation rate from the head-on collisions only as a function of ball sizes for a series of charge ratios corresponding to the experimental data in Fig. 7.For comparison, a similar plot for the numerical milling dose based on energy dissipation rate in glancing impacts only is shown in Fig. 9.
The numerical milling corresponding to the total E d , shown in Fig. 5, reveals a relatively constant milling dose as a function of ball size, while the milling dose based on the head-on energy dissipation rate shows a different trend which correlates better with the experimental milling dose as shown in Fig. 7. Conversely, the calculation based on the glancing energy dissipation reveals different value of milling dose upon different powder loading, where larger loading results in smaller milling dose, which is not supported by the experimental results.The numerical milling dose based on the head-on dissipation energy rate has a divergence similar to the experimental milling dose for small diameter balls, but demonstrates a less efficient milling (shown as larger average milling dose value) for C R ҃10 rather than C R ҃2.5 as the experimental results.For the intermediate ball sizes (d҃3.16mm and d҃4.76 mm), which generate collisions with very similar levels of energy dissipa-tion (c.f.discussion in section 4 for Fig. 2), the numerical milling dose agrees well with the experimental milling dose.This may imply that the specific range of dissipation energies is critical for reactive milling and, by implication, for the progress in mechanical alloying.
The validity of experimental milling dose and numerical milling dose based on the energy dissipation rate for the head-on collisions is also well represented in the DEM results with a near constant milling dose for 4.76 mm balls among all charge ratios.This intersection point represents the most efficient milling condition for the examined charge ratios.Therefore, head-on collisions are found to be more significant in defining the milling progress.This finding agrees with reports stating that the glancing impacts do not contribute significantly to the deformation, coalescence, and fragmentation 4,7) .The results shown so far indicate that this methodology of combined numerical and experimental study can be used for the analysis of mechanical alloying suitable for the further process development and scale-up.The numerical (or theoretical) milling dose, although expressed in terms of energy, behaves similarly to the experimental milling dose which is expressed in terms of time.This demonstrates the DEM methodology's ability to capture the physical interactions during milling.In a scaled up operation, the optimal conditions, observed here for 4.76 mm balls, may actually occur at a different ball size, but in general, numerical modeling can be expected to identify the optimal ball size and help identify the additional parameters required for successful production.The same idea could be applied to the analysis of other dynamic systems, such as mechanofusion, where a parameter similar to milling dose can be used to evaluate the coating degree.That topic will be discussed in a future paper.

Conclusions
The numerical study of the mechanical alloying in a SPEX 8000 shaker mill is carried out by DEM technique.The movement of the vial is implemented into the modeling according to the experimental measurements.During the simulation, the energy dissipation rate of the system is calculated for four different ball sizes and three charge ratios.The numerical milling dose is thereafter constructed based on the energy dissipation rate and correlated with the experimental milling dose, which is a useful indicator of the milling progress for reactive milling.There is a good agreement between numerical simulation and experimental data indicating a trend for a constant milling dose for a specific material.Furthermore, the numerical milling dose based on the energy dissipation rate from headon impacts demonstrates a better correlation with the experimental milling dose than the similarly defined milling dose calculated considering the energy dissipation rates for all or only glancing impacts.This classification of the energy dissipation rate according to head-on and glancing impacts can only be accomplished by numerical simulation, emphasizing its importance for the analysis of the actual milling mechanism.Particularly good description of the milling progress by both the experimental milling dose and its numerical analog is found for intermediate ball sizes, which result in collisions within a specific range of energies.The collisions with energies which are higher or lower than the observed threshold, and which occur when larger or smaller balls are used, appear to be less effective in achieving the milling progress.Because the numerical simulation successfully predicts the trends observed in the experiments, it is suggested that the model can accurately describe the progress of mechanical alloying and reactive milling, and is useful for further scale-up studies.

Fig. 6
Fig. 6 Simulation results showing the energy dissipation rate as a function of the ball mass for different ball sizes.

Fig. 9
Fig. 9 Theoretical/numerical milling dose based on the glancing energy dissipation rate as a function of ball size.

Table 1
Input parameters of simulation Computed changes of the function C R • m p /E d for different ball sizes and different charge ratios.
Experimental data on reactive milling of 2AlѿFe 2 O 3 composition presented as the product of the charge ratio and milling time plotted versus ball size for different charge ratios.
Fig. 8 Theoretical/numerical milling dose D m based on the headon energy dissipation rate as a function of ball size.