- Email: [email protected]
Brownian-like kinematics of ball milling for particulate structural modeling
Powder Technology 301 (2016) 1077–1084
Contents lists available at ScienceDirect
Powder Technology journal homepage: www.elsevier.com/locate/powtec
Brownian-like kinematics of ball milling for particulate structural modeling Constantine C. Doumanidis a, Hamda A. Al Kaabi b, Abdelaziz S.M. Alzaabi b, Ibrahim E. Gunduz c,⁎, Claus Rebholz d, Charalabos C. Doumanidis b a
Pankyprion Gymnasion, Nicosia, Cyprus Department of Mechanical Engineering, Khalifa University, Abu Dhabi, UAE c School of Mechanical Engineering, Purdue University, W. Lafayette, Indiana, United States d Department of Mechanical & Manufacturing Engineering, University of Cyprus, Nicosia, Cyprus b
a r t i c l e
i n f o
Article history: Received 21 March 2016 Received in revised form 20 May 2016 Accepted 9 July 2016 Available online 19 July 2016 Keywords: Ball milling Brownian motion Modeling Microstructure NiAl
a b s t r a c t Ball milling motion has been previously studied through computationally expensive, off-line experimental video processing and numerical simulations by the discrete element method. This research establishes a more efficient formulation of the ball energetics and kinetics similar to the Brownian kinetic theory of statistical mechanics. Based on assumptions of thermomechanical equilibrium, negligible gravitational, aerodynamic and surface condition effects, and decoupled impact interaction among balls and with milled particulates, this model proposes mono-parametric spectral energy and velocity probability density functions akin to Maxwell-Boltzmann statistics, along with uniformly distributed impact directionality. The model predictions are calibrated and validated by comparison with published experimental measurements and computationally derived spectra. This descriptive Brownian-like motion model enables effective simulation of contact and impact, material deformation and micro-joining of ball milled bimetallic powders. A comprehensive simulation of the evolving internal fractal microstructure of the processed particulates is implemented at real-time computation speed, and its predictions are compared with experimental micrographs of ball milled Ni\\Al particulates. © 2016 Elsevier B.V. All rights reserved.
1. Introduction Milling of particulates, originally powered by humans, animals, water and wind, has been one of the prehistoric fundamental manufacturing methods for processing of rock and soil, construction materials, foodstuff and medicines [1]. Through compressive impact and/or shearing between hard surfaces, resulting in plastic deformation and/or crushing fracture, milling has been used for comminution, i.e. reduction of volumetric size and increase in active surface, coagulation or sharpening of particle shape, as well as mixing, agglomeration and compounding of diversified solid matter. Material synthesis through rolling, tumbling, rod and ball milling resulting in mechanical alloying of metals has been of recent interest, particularly in fabrication of reactive particulates out of bi-material powder mixtures such as Al\\Ni [2– 5]. The resulting multi-scale Apollonian packing or lamellar micro/ nano-structure of the milled particulates, upon their thermal or electrical ignition, determines the kinetics of self-propagating exothermic
⁎ Corresponding author. E-mail addresses: [email protected] (C.C. Doumanidis), [email protected] (H.A. Al Kaabi), [email protected] (A.S.M. Alzaabi), [email protected] (I.E. Gunduz), [email protected] (C. Rebholz), [email protected] (C.C. Doumanidis).
http://dx.doi.org/10.1016/j.powtec.2016.07.033 0032-5910/© 2016 Elsevier B.V. All rights reserved.
reactions (SPER) [6,7] with highly-localized and time-controlled enthalpy release and high-temperature synthesis (SHS) of refractory compounds. Application of such nanoheater particulates [8] in rapid thermal processing, micro-joining and micro-soldering of electronics, self-sintering powders and self-healing composites, all the way to biomedical thermal ablation of tumors, has kindled research in understanding the process-structure relationships of milling energetics with the resulting self-similar (fractal) bimetallic microstructures in mechanically alloyed particulates. Ball milling in rotating uniaxial cylindrical vial or planetary settings (Fig. 1a) has been implemented for this purpose in both high-energy (interrupted adiabatic or actively cooled continuous) and low-energy (continuous or interrupted isothermal) configurations, and has been extensively investigated by laboratory studies and computational simulation. Experimental observation of the ball motion in the milling vial through a transparent lid by high-speed (HS) optical recording [9–13] has elucidated and empirically described their kinematics in stagnant, cascading, cataracting and centrifuging modes of ball trajectories [5]. However, laboratory process monitoring via HS video is computationally costly in its requisite image post-processing. Nevertheless this provided the basis for identification of (visco)elastic hysteresis/plastic collision and frictional parameters of material models in the theoretical linear or non-linear dynamic equations of the ball motion [14–17], such
1078
C.C. Doumanidis et al. / Powder Technology 301 (2016) 1077–1084
Fig. 1. a. Planetary ball mill setup, b. Milling ball-particulate impact.
as coefficients of restitution, rebound angles, contact duration and forces. Such collision formulations enabled calibration and validation of multi-body interaction simulations in the milling vial by the discrete element method (DEM) [18–24]. Numerical DEM models have provided valuable information on impact frequencies, ball wear rates, heat generation and power requirements for the ball milling process. However, complex Lagrangian simulations of multiple moving balls in a rotating vial by DEM are computationally inefficient, typically requiring parallel or super-computing resources for prolonged ball milling process times, while often missing description of mechanical interactions with and structural effects to the milled particulates. Therefore current numerical models, along with laboratory sensing techniques, are suitable for offline estimation and statistical process control (SPC) only, requiring counter-productive and invasive process interruption for sample removal and characterization to correct ball milling conditions [25]. Recently, real-time modeling of the particulate microstructure for in-process feedback control of ball milling was introduced, based on a computationally cost-effective semi-analytical predictive simulation of the dynamic product geometry and material composition [26]. The latter hinges on an Eulerian formulation, following and describing a typical particulate as it develops during the process, through its collisions with the balls and vial. These result in impact loading, contact micro-bonding with other particulates, material deformation, heat dissipation, gradual diversification of particulate species and their population proliferation in the final milled product. The sequence of these phenomena commences with and critically depends on the energetics and kinetics of the milling balls and vial. For compatibility with this real-time framework and instead of a full DEM simulation, a comprehensive yet parsimonious model of ball milling motion is required. This paper elaborates on such an ab-initio analytical, theoretical physics-based, experimentally calibrated and computationally validated descriptive model of the milling ball energy and velocity stochastic distributions, along with its integration into the particulate structure simulation. 2. Theory 2.1. Ball milling motion During steady-state ball milling, the previous experimental and numerical literature elucidates a continuous, time-invariant, chaotic motion pattern of the ball population in a confined vial, consisting of deterministic free trajectory segments between random collisions with other balls and vial walls, under predictable impact conditions. Upon distant zoomed-out observation from a larger dimensional (macro) scale, this Wiener process/Lévy flight pattern is reminiscent of the familiar nanoscale Brownian motion of ideal gas atoms/molecules at thermal equilibrium, with well-established descriptions by kinetic theory [27–29]. Despite similarities across scales, however, a number
of salient differences in milling ball motion behooves certain assumptions before proceeding to modeling analogies. 2.1.1. Steady state thermomechanical equilibrium The Brownian model involves fully reversible, perfectly elastic collisions between spherical particles and with stationary, adiabatic container walls. In ball milling, irreversible mechanical energy losses into heat during impacts, as denoted by the coefficient of restitution [22–24], are compensated during collisions with the mechanically moving (inducting work), thermally abducting vial walls. Thermomechanical equilibrium assumed at steady state implies constancy of the total mechanical energy of the balls, as well as the enthalpy of the system, i.e. isothermal ball milling conditions analogous to those in an ideal gas. This also ensures constant material properties of the balls and vial. 2.1.2. Negligible gravitational, aerodynamic and surface condition effects In a fixed-elevation ball milling setup, gravity introduces no mechanical work to the average ball energy, while limited vial height (Fig. 1a) has negligible influence to the standard deviation, especially in high-energy systems. In low-energy systems gravity bias may cause visibly parabolic trajectory segments of the balls instead of the straight ones in Brownian motion. Aerodynamic drag and lift (through the vorticity of spinning balls) interactions with the inert gas in the setup are also negligible, since they scale as the mass of gas volume swept between collisions over that of the balls (i.e. by about 3 orders of magnitude lower). Dry ball milling without liquid process control agents is assumed, avoiding complex non-Newtonian viscous and surface tension/capillary flow effects. 2.1.3. Decoupled impact interaction among balls and with milled particulates As per kinetic theory, decoupled ball motion requires that their mean free path between successive collisions as analyzed below, should be larger than the characteristic geometric size, i.e. the harmonic mean of the ball and vial dimensions. While temporarily adopted for simplicity (i.e. balls only in vial), the assumption on ball-particulate interactions will be explicitly relaxed in the sequel to reflect transport of mechanical energy to the processed particulates during collisions. The (normal and shear) impact energy transported to the particulate contact area, and the corresponding restitution of the impacting ball or vial, are determined by the relative linear and spin velocity components of the impactor surface at the contact location. Such an effective composite velocity v at the surface of a vial wall, e.g. in a planetary mill with angular velocity Ωο and radius Ro of the sun, and respectively Ω and R on the planet surface as in Fig. 1a, is: V 2 ¼ ðΩRÞ2 þ ðΩo Ro Þ2
ð1Þ
C.C. Doumanidis et al. / Powder Technology 301 (2016) 1077–1084
1079
Similarly, at the surface of a ball of radius r and mass m the effective velocity v corresponds to its total kinetic energy E (translational and rotational), i.e. (Fig. 1b): E¼
1 1 1 2 1 2 mυ2 þ Jω2 ¼ m υ2 þ r 2 ω2 ¼ mv2 →v2 ¼ υ2 þ ðωrÞ2 ð2Þ 2 2 2 5 2 5
where J = mρ2 is the ball moment of inertia around its rotation axis, pffiffiffiffiffiffiffiffiffiffi and ρ ¼ 2=5 r its gyration radius. 2.2. Ball energy and velocity distribution Besides the deterministic motion of the vial, the experimental evidence [5,12] points to a statistical distribution of the ball effective velocity v and energy E similar to that of Brownian motion. According to statistical mechanics [29,30], in an ideal gas or solution, the number N of atoms/molecules with kinetic energy E is proportional to its inverse exponential, i.e. N ∝ exp (− E/kT), where k is the Boltzmann constant and T the gas temperature. Integration over the entire population of atoms/molecules yields their average kinetic energy hEi ¼ 12 mhv2 i ¼ 32 kT for 3-D motion. For each 1-D motion component (k = x,y,z), the average energy hEk i ¼ 12 mhv2k i ¼ 12 kT ¼ 12 ma2 , where the parameter a pffiffiffiffiffiffiffiffiffiffiffiffi ¼ kT=m. By analogy to statistical mechanics, the energy E of the milling balls of mass m follows a probability density function (pdf) f(E) (Fig. 2a): 2 f ðEÞ ¼ pffiffiffi 2 πma
rffiffiffiffiffiffiffiffiffi E E exp − ma2 ma2
ð3Þ
In this distribution, the most probable energy is Ep = ma2, and the average energy hEi ¼ 32 ma2 = 1.5Ep. qffiffiffiffiffiffiffiffiffi The parameter a ¼ hv2k i of this distribution is the root mean square (RMS) velocity in one direction of motion k, and can be conveniently measured experimentally from a video sequence of the ball motion if available. The measured a = γV can be related to the vial cylinder wall velocity V (Eq. (1)) via an empirical efficiency factor γ. Alternatively a can be determined in the laboratory through the steady-state average energy 〈E〉 of the balls in the vial, from power measurements during the initial transient to steady operation or the final transient to rest of the milling system (Fig. 3). The steady state mechanical energy stored in N balls N〈E〉 , the rotating vials (Eq. (1a) and other parts EV = Σ1/ 2JjΩ2j of the milling system, plus its enthalpy H = CT at thermal equilibrium, originates from the power P = MΩm supplied by the motors minus the thermal power Q = T / Rq removed from the setup (shaded
Fig. 3. Initial and final transient of ball milling power.
area in Fig. 3): Z τ X 1 3 NhEi þ Ev þ H ¼ J ω2 þ CT ðP−Q Þdt→N ma2 þ j2 j j 2 0 Z τ ¼ MΩm −T=Rq dt 0
ð4Þ
where Jj and Ωj are the moments of inertia and angular velocities of rotating elements, M and Ωm the motor torque and rotation speed (measured via the electrical current and voltage), C the heat capacity of the system at uniform equilibrium temperature T (measured e.g. via IR thermocouples), and Rq the total heat loss resistance of the system, with C and Rq being experimentally measurable parameters. pffiffiffiffiffiffiffiffiffiffiffiffi Through its relation with energy E, the effective velocity v ¼ 2E=m at the ball surface (Eq. (2)) is also statistically distributed, similar to a Maxwell-Boltzmann pdf [30,31] derived from Eq. (3) (Fig. 2b): pffiffiffi 2 2 v2 v2 ϕðvÞ ¼ pffiffiffi exp − 2 2 2a πa 2a
ð5Þ
pffiffiffi The most probable velocity in this pdf is vp ¼ 2a, the average velocqffiffi qffiffi pffiffiffiffiffiffiffiffiffi pffiffiffi ity is hvi ¼ 2 π2a ¼ p2ffiffiπ vp, and the RMS velocity hv2 i ¼ 3a ¼ 32vp in 3-D ball motion. 2.3. Impact with particulates In the presence of milled particulates in the vial along with the balls, it is now assumed additionally that both these elements follow kinetic energy E and effective velocity v distributions with pdfs similar to Eqs. (3) and (5) respectively, with separate parameters. Moreover, for simplicity the Eulerian collision simulation described below [26] follows
Fig. 2. a. Milling ball energy pdf (m = 0.03). b. velocity pdf.
1080
C.C. Doumanidis et al. / Powder Technology 301 (2016) 1077–1084
the motion of one evolving particulate, assumed stationary and impacted by impinging vial walls and ball surfaces (Fig. 1) moving at relative velocities vi also according to a Maxwell-Boltzmann-like pdf. The equivalent (or reduced) mass m* of the colliding two-body system is the harmonic sum of the particulate mass μ and the impactor (ball or vial) m [32]: m ¼
μm μ þm
ð6Þ
typically resulting to m* ≈ μ since μ ≪ m in ball milling. Therefore the kinetic energy of the impactor before collision Ei ¼ 12 m v2i and after it Eo ¼ 12 m v2o depends on the relative input vi and output vo velocities, with the coefficient of restitution defined as e = vo/vi. During each such collision, from the impinging vial or ball energy Ei the particulate receives an amount of elastic work w = Eo which is reversibly restored as kinetic energy to the impactor upon separation, as well as an inelastic work W retained irreversibly by the particulate and eventually converted to heat. The latter consists of energy Wr stored in hysteretic elastic and residual stress fields (which is subsequently thermalized by stress relieving, annealing, etc.), plastic deformation work Wp to the materials, as well as frictional work Wf at the contact surfaces. This irreversible (positive) work W is therefore: Ei −Eo ¼
1 2 1 2 m vi − m vo ¼ W ¼ W r þ W p þ W f N0 2 2
Fig. 5. Particulate surrounded by inner cluster and outer medium layers while impacted by milling ball for its structural simulation.
ð7Þ
implying that vi N vo and therefore e b 1. Therefore for a certain impact event, while the input velocity vi can be randomly selected from the full range of pdf ϕ(vi) (Eq. (5)), the resulting output velocity is statistically limited in the range [0...vi] corresponding to the already selected vi. Therefore vo must follow a finite-range pdf ϕ(vo) / Φ(vi), normalized by the value of the cumulative distribution function (cdf) Φ(v) at vi: Z ΦðvÞ ¼
v 2 v v2 ϕðυÞdυ ¼ erf pffiffiffi − pffiffiffi pffiffiffi exp − 2 2a π 2a 2a 0 v
ð8Þ
This also implies that, while the corresponding incoming energy Ei to the collision is randomly selected from the full-range pdf f(Ei) (Eq. (3)), the concomitant out-coming energy Eo is limited in the range [0…Ei] of the pdf f(Eo) / F(Ei), i.e. normalized by the value of cdf F(E) at Ei: Z F ðEÞ ¼
rffiffiffiffiffiffiffiffiffiffiffi! rffiffiffiffiffiffiffiffiffiffiffi E 2 E E p ffiffiffi f ðε Þdε ¼ erf − exp − m a2 m a2 π m a2 0 E
ð9Þ
The kinetic energy Ei of the impactor due to its incoming velocity vi is imparted to the particulate via a normal vn and a tangential vt component (Fig. 1b), corresponding to a normal En ¼ 12 m v2n at and a tangential
Fig. 4. 2-D warped ellipsoidal primitive of particle domain for particulate structure simulation.
Fig. 6. a. DEM simulation [22] vs kinetic spectra (continuous line) of ball energy during ball milling. b. Experimental (1st column) and DEM simulation (2nd column) [19] vs kinetic energy spectra.
C.C. Doumanidis et al. / Powder Technology 301 (2016) 1077–1084
Et ¼ 12 m v2t impact energy: 1 2 1 2 1 2 m vi ¼ m vn þ m vt ¼ En þ Et ¼ W þ Eo and v2i 2 2 2 ¼ v2n þ v2t
Ei ¼
ð10Þ
The random incidence angle θ = tan−1(vn / vt) of the impactor on the local contact surface of the particulate is assumed to follow a uniform pdf within the range [−π/2…π/2]. Finally the particulate surface is hit by an impactor which is either one of N balls (of curvature 1/r), a vial top or bottom cap (curvature 0) or a vial cylinder (curvature 1/R). Because of the uniform pdf of the incidence angle θ, the respective probabilities are proportional to the solid angles subtended by each impactor type to the particulate, i.e.
1081
~ Nr2 / 4 bD N2 for a ball, ~ 2(R/h)2 for a vial cap and ~ 1/2(h/R) for a vial cylinder. The collision frequency f = N/b tN of N balls at all energy levels in the vial of volume O = πR2h depends on the average flight pffiffiffiffiffiffiffiffiffi time b tN = b DN/ hv2 i of a single impactor between successive collipffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffi sions across a mean free path hDi≈ 3 O=N , with hv2 i the RMS ball velocity. 2.4. Particulate contact, deformation and joining Statistical selection of impact energy Ei or velocity vi and incidence angle θ according to the previous experimentally validated analytical pdfs is instrumental in the implementation of a real-time simulation of particulate microstructural evolution in ball milling [26]. Consistent
Fig. 7. Experimental [5] vs Maxwell-Boltzmann (dashed line) velocity spectra for various values of speed ratio K = Ω / Ωo of planetary ball milling.
1082
C.C. Doumanidis et al. / Powder Technology 301 (2016) 1077–1084
with experimental planar micrographs of such self-similar Apollonian and lamellar structure formations in bimetallic ball-milled materials [2–8], this model is a 2-dimensional, dynamic simulation of the crosssectional fractal geometry and material distribution in one typical particulate as it develops. The model assumes a versatile warped ellipsoidal primitive geometry for each monometallic particle domain, adequately describing its deformation from initially spheroidal powder to the final lamellar layer form, through variations of its position, orientation and dimensional descriptive parameters (Fig. 4). It also adopts a linear elasto-plastic strain-hardening behavior of each material volume, with dry Coulomb friction on its surface, and temperature-invariant (steady-state isothermal), impact high-rate specific parameters [24]. Before each impact event, this original primitive or developed composite particulate is surrounded by statistically selected additional clusters, similar to the current particulate and/or its prior development generations as explained below, and brought into random but stable contact by relative surface rolling of each cluster (Fig. 5). During each collision, this adjacent inner layer is considered for potential integration of each neighboring cluster to the particulate by micro-welding at their contact locations. The entire assembly is encircled by an outer boundary layer of particles, with contact loading insufficient for joining to the particulate during the impact event, but providing mechanical motion confinement and energy dissipation. A ball or vial impactor subsequently collides with the stationary particulate at a random contact site of its surface according to the previously established kinetics and energetics. During the progression of each impact event, an initially elastic Hertzian contact at the particulate-impactor ellipsoidal interface patch with respective equivalent contact curvatures and stiffnesses, absorbs the impact energy components (Eq. (10)) over the particulate bulk material [33,34]. This provisional elasticity assumption yields the normal and shear load distributions over the contact patch, along with its progressing footprint ellipsoid radii, penetration approach and displaced interface geometry. This contact loading of the particulate and surrounding inner cluster layer determines the normal and shear reaction distributions over their boundary surface with the outer restraint layer of materials, considered as an equivalent elastic circumscribing continuum with an ellipsoidal inclusion [35,36]. This outer boundary reaction loading is used for computation of the inner contact loads between the particulate and adjacent clusters (Fig. 5) under mechanical equilibrium conditions, along with all requisite internal elastic stress fields. Since such redundant contact conditions lead to indeterminacy of their solutions, Castigliano's theorem is employed to determine elastic deformation displacements at the contact locations and their resulting Hertzian conditions, and hence the tentative deformations of individual particulate and cluster domains and their internal elastic strain fields [37]. Next, the equivalent elasticity assumption above is relaxed and the total work w + W absorbed by each domain of the particulate assembly is differentiated into a reversible elastic component w and an irreversible inelastic one W as explained above. This is done according to the stress-strain diagram of its bulk material, therefore defining plastic work Wp and eventual deformation of the domain, as well as according to the load-slip conditions at its contact surfaces, defining frictional work Wf. When the specific thermalized energy distribution Wp + Wf over a certain contact surface area A locally exceeds a critical activation energy density Ew/A for joining of the bi-material system, a micro-weld is established over A. The load-carrying ability of this micro-joint is determined by the specific residual strain energy distribution, i.e. inelastic work Wr = W-Wp-Wf over A. Such surface joining of domain interfaces in the particulate also results in its integration with additional surrounding clusters; while separation of particulate clusters also occurs when joint strength is exceeded upon subsequent impact loading. Therefore, sequential collisions with the impactors cause a progressive evolution of the external geometric size and shape along with internal material structure of the particulate under study, yielding successively new particulate generations along with previous ones in
the setup, all the way back to original mono-material powders. Since the other balls and vial impactors, besides the one involved in the collision with the studied particulate, also act simultaneously during ball milling, the populations of particulate generations grow in parallel according to a logistic parabola law. The pdf of each particulate generation is composed of products of the population sizes for combinations of its prior generations. These pdfs are used in assembling the adjacent cluster layer around the particulate before each collision; and they reflect the eventual demographics of particulate diversity in the ball milled product.
3. Results and discussion Fig. 6a compares the energy spectrum f(E) of the physics-based pdf in Eq. (3) (a = 0.36 m/s, γ = 0.43) with those derived in [22] by experimentally calibrated DEM simulation for a ∅90 × 14 cm cylindrical vial with 32 ∅5.08 cm stainless steel balls rotating at 18 rpm. The statistical model satisfying the decoupled ball interaction condition fits rather well, but somewhat over-predicts impact frequencies at lower ball energies, while it under-predicts especially at the very high energy end. Fig. 6b also shows a comparison of the energy distribution predictions (a = 0.99 m/s, γ = 0.71) with those in [19] by experiment and DEM simulation, for a ∅90 × 4 cm vial with 150 ∅2.9 cm alumina balls, rotating at 30 rpm for 30 s. The statistical model of decoupled impacts fits the measured impact frequencies within ~10%, again with over-prediction at low ball energies (possibly due to missed impacts by the image processing [19]), while it somewhat under-predicts the DEM results. Fig. 7 illustrates a comparison of the velocity predictions ϕ(v) of Eq. (5) with experimental measurements in [5] for planetary ball milling (Fig. 1), with Ro = 10.4 cm, R = 4 cm, h = 5 cm, ball r = 3 mm and volume fraction 0.35 in the vials, and sun Ωo = 694 rpm, for various planet rotation speeds Ω = KΩo within 0–1388 rpm. The distribution parameters fitting the laboratory data are: Speed ratio K = Ω/Ωo Pdf parameter a (m/s) Efficiency parameter γ
1.0 1.0 0.12
1.5 0.92 0.094
1.8 0.65 0.071
1.9 1.98 0.21
2.0 2.41 0.25
It can be seen that the statistical distribution slightly over-predicts the ball fraction at very low velocities as explained before, while it somewhat under-predicts at high velocities. A plausible reason for these deviations, also manifested in the variance of the efficiency parameter γ, is the high bulk charge challenging the decoupledness assumption for ball impacts. The variation in pdf parameter a and efficiency parameter γ can also be related to the nature of ball motion observed in [5]. The first two speed ratios correspond to mostly rolling
Fig. 8. Experimental (x) and DEM simulated (solid line) [24] vs kinetic theory predictions (dashed line) of coefficient of restitution e = vo/vi during stainless steel ball milling.
C.C. Doumanidis et al. / Powder Technology 301 (2016) 1077–1084
1083
Table 1 3D Fractal dimensions of experimental and simulated sections at 2, 4, 6 and 8 ball milling hours. BM time/fractal dimension Capacity Dimension Information Dimension Correlation Dimension Probability Dimension
Fig. 9. DEM simulated [23] vs kinetic theory predicted (dashed lines) of normal impact velocity vn. during planetary ball milling with 10 and 20 balls.
motion (cascading) of balls on balls with similar parameters a and γ. The last two ratios correspond to a falling motion (cataracting) where the balls are rolling mostly on vial surfaces and crossing across the vial walls during a part of the rotational cycle. These motions can be expected to produce both higher average ball speeds and high efficiency parameters due to longer and more intimate contact with the vial walls. The comparison of coefficient of restitution e = vo/vi determined by sequential selections of V = vi and vo according to the pdfs in Section 2.3 (with a = 0.31 m/s, γ = 0.042), to measured values and DEM predictions in [24] for stainless steel ball milling is shown in Fig. 8. The decoupled statistical formulation again results to some over-prediction of e for low velocities V, and a slight under-prediction for higher values. Fig. 9 shows a comparison of the statistical distribution of the normal impact velocity vn = vi sinθ (with a = 0.34 or 0.26 m/s, γ = 0.096 or 0.074), with the DEM predictions in [23] for planetary milling with Ro = 12.5 cm, R = 4 cm, h = 6 cm, Ωo = 300 rpm and Ω = −750 rpm, using N = 10 or 20 stainless steel balls of r = 5 mm. This confirms that components of the impact velocity under decoupled ball interaction conditions are also distributed randomly according to the Maxwell-Boltzmann pdf of Eq. (5), therefore justifying the assumption of uniform distribution of θ over the contact half-space. The calculated impact frequencies f ≈ 120 Hz for the conditions of [22] and f ≈ 25 Hz for those of [19], according to the Section 2.3, are consistent with the values of Fig. 6a and b, respectively. To validate the ball milling energetics formulation, the previously described simulation model [26] with decoupled impact kinetics following a Maxwell-Boltzmann pdf (Eq. (5)) was computationally
2 Hours Exp
4 h
6 h
8 h
Model Exp
Model Exp
Model Exp
Model
2.54 2.36 2.01 2.04
2.74 2.56 2.59 2.45
2.77 2.57 2.59 2.45
2.79 2.69 2.73 2.61
1.69 1.78 2.26 2.23
2.07 2.19 2.45 2.41
2.14 2.15 2.44 2.47
2.46 2.31 2.50 2.49
implemented for the experimental conditions of [4]. In [4], a low-energy planetary ball milling system with cylindrical stainless steel vial of volume 80 ml (Ro = 121 mm, R = 20 mm) containing 5 stainless steel balls of 2r = ∅10 mm in inert gas (N2) was employed, at room temperature and at rotation speed of 300 rpm (K = 1.82). This was loaded with a total mass 32 g of Al (99.5%) and Ni (99.95%) powders of 325 mesh size at a molar ratio of 1:3. The powder load was processed by continuous ball milling, briefly interrupted after 2, 4, 6 and 8 h for partial sample removal and analysis. These samples were sectioned, polished and characterized by scanning electron microscopy (SEM, back scattered electron BSE-first row) in Fig. 10. The micrographs illustrate a typical particulate structure evolution, from initial coarse globular Apollonian assemblies of harder Ni particles surrounded by softer Al phase, to planarized finer lamellar multi-layers of the two metals in the subsequent and final sections. The second row of Fig. 10 comparatively shows the respective particulate structure sections as predicted by the ball milling simulation model for the same processing times as before, with its kinetic parameters (a = 0.87 m/s, γ = 0.22) calibrated from the laboratory conditions. Similar to the experimental BSE micrographs above, initially larger warped ellipsoids of Ni (bright domains) appear as spheronized by the omni-directional random impactor collisions, directionally bent and joined with smaller Ni particles through intermediate highly deformed Al (dark domains) in an Apollonian-packed fractal architecture. With progression of the process, further impact to the particulate surface leads to more bonding and unification of similar metal domains into continuous layers, as well as dispersion of dissimilar layers with clearly discernible bilayer interfaces in lamellar formations. These multi-layers tend to planarize due to the higher probability of collisions normally to the flatter, larger radius surfaces of an ellipsoidal particulate (Fig. 4) compared to its curved, small radius sides. This is because of their respective difference in solid angles subtended to the directionally uniform incidence of impactors. Further deformation and refinement of the lamellae is evident in the finally produced microstructures by ball milling. For a quantitative comparison of the fractal dimension metrics in the experimental and simulated particulate sections of Fig. 10 (first and
Fig. 10. Particulate microstructure after 2, 4, 6 and 8 h of planetary ball milling Experimental SEM micrographs [4]: BSE (1st row) versus simulated internal structure of impacted particulate (2nd row).
1084
C.C. Doumanidis et al. / Powder Technology 301 (2016) 1077–1084
second row) after 2, 4, 6 and 8 h of ball milling, a box-counting algorithm [38] was implemented for their off-line image analysis with different measures. The respective attributes in Table 1 show comparative agreement of the model predictions to the laboratory data within ± 7% of the parameter values. This indirectly confirms sufficient calibration and suitability of the Brownian impactor kinetics for the particulate structure simulation. Therefore, unlike prior DEM simulations of ball milling motion, the computational efficiency of this kinetic theory formulation enables real-time implementation of the particulate structural model on a standard desktop computer running at 130 GIPS. 4. Conclusions The previous Brownian-like model of statistical ball milling motion is described by cost-efficient analytical probability density functions for the ball energy and velocity, and by spherically uniform directionality of impactor incidence. A single calibration parameter in a MaxwellBoltzmann expression, i.e. the RMS velocity of the balls in a single direction, can be either calibrated from video-measured energy or velocity laboratory spectra, or computed from the initial or final transients of the process conditions. The estimates of this descriptive model directly fit published energetic and kinetic spectra for various decoupled ball milling process conditions, while they also enable simulation of impacted particulate structures in low-energy planetary processing of bimetallic powders at real-time computation speeds on ordinary hardware. The convenience of this probabilistic kinetic theory-based model of ball milling is attractive for description of other milling processes, potentially challenging its original assumptions and calling for appropriate adjustment of its formulation and/or parameters. High-energy uncooled processing, for example, may result in dominant thermomechanical transients departing from steady-state equilibrium and isothermal assumptions. Temperature variation of the material parameters may also alter impact restitution and affect spectral distribution parameters. Use of liquid process control agents may also affect such parameters through concomitant viscosity and capillary surface tension effects during contact. Finally coupled interaction of ball impactors at high bulk volume charges or with larger milled particulates of commensurate mass and energy, as well as complex surface geometries and impact property distributions, could invalidate the decoupledness assumption of ball motion made above. Such conditions may necessitate coupledspecies, modified energetic and kinetic spectral expressions with dependencies on non-uniform impact directionality distributions, to be established through subsequent investigation. Acknowledgments This research was supported in part by a Khalifa University Internal Research Fund award (Level 1) and by the US Department of Energy, National Nuclear Security Administration under Award Number DENA0002377. Appendix A. Supplementary data Supplementary data to this article can be found online at http://dx. doi.org/10.1016/j.powtec.2016.07.033. References [1] A.J. Lynch, C.A. Rowland, The History of Grinding, Society for Mining, Metallurgy and Exploration, 2005. [2] A. Hadjiafxenti, I.E. Gunduz, C. Tsotsos, T. Kyratsi, S.M. Aouadi, C.C. Doumanidis, C. Rebholz, Synthesis of reactive Al/Ni structures by ball milling, Intermetallics 18 (2010) 2219–2223.
[3] K.V. Manukyan, B.A. Mason, L.J. Groven, Y.C. Lin, M. Cherukara, S.F. Son, A. Strachan, A.S. Mukasyan, Tailored reactivity of Ni + Al nanocomposites: microstructural correlations, J. Phys. Chem. C 116 (2012) 21027–21038. [4] A. Hadjiafxenti, I.E. Gunduz, C.C. Doumanidis, C. Rebholz, Exothermic reaction characteristics of continuously ball-milled Al/Ni powder compacts, Vacuum 96 (2013) 73–78. [5] A.S. Rogachev, D.O. Moskovskikh, A.A. Nepapushev, T.A. Sviridova, S.G. Vadchenko, S.A. Rogachev, A.S. Mukasyan, Experimental investigation of milling regimes in planetary ball mill and their influence on structure and reactivity of gasless powder exothermic mixtures, Powder Technol. 274 (2015) 44–52. [6] I.E. Gunduz, K. Fadenberger, M. Kokonou, C. Rebholz, C.C. Doumanidis, T. Ando, Modeling of the self-propagating reactions of nickel and aluminum multilayered foils, J. Appl. Phys. 105 (2009) 074903. [7] K. Fadenberger, I.E. Gunduz, C. Tsotsos, M. Kokonou, S. Gravani, S. Brandstetter, A. Bergamaschi, B. Schmitt, P.H. Mayrhofer, C.C. Doumanidis, C. Rebholz, In-situ observation of rapid reactions in nanoscale Ni-Al multilayer foils using synchrotron radiation, Appl. Phys. Lett. 97 (2010) 144101. [8] Z. Gu, Q. Cui, J. Chen, J. Buckley, T. Ando, D. Erdeniz, P. Wong, A. Hadjiafxenti, P. Epaminonda, I.E. Gunduz, C. Rebholz, C.C. Doumanidis, Fabrication, characterization and applications of novel nanoheater structures, Surf. Coat. Technol. 215 (2013) 493–502. [9] P. LeBrun, I.L. Froyen, L. Delaey, The modeling of the mechanical alloying process in a planetary ball mill: comparison between theory and in-situ observations, Mater. Sci. Eng. A161 (1993) 75–82. [10] R.K. Rajamani, B.K. Mishra, R. Venugopal, A. Datta, Discrete element analysis of tumbling mills, Powder Technol. 109 (2000) 105–112. [11] H. Mori, H. Mio, J. Kano, F. Saito, Ball mill simulation in wet grinding using a tumbling mill and its correlation to grinding rate, Powder Technol. 143-144 (2004) 230–239. [12] S. Rosenkranz, S. Breitung-Faes, A. Kwade, Experimental investigation and modeling of the ball motion in planetary ball mills, Powder Technol. 212 (2011) 224–230. [13] H. Shin, S. Lee, H.S. Jung, J.-B. Kim, Effect of ball size and powder loading on the milling efficiency of a laboratory scale wet ball mill, Ceram. Int. 39 (2013) 8963–8968. [14] M. Abdellaoui, E. Gaffet, The physics of mechanical alloying in a modified horizontal rod mill: mathematical treatment, Acta Mater. 44 (1996) 725–734. [15] C. Thornton, Z. Ning, A theoretical model for the stick/bounce behavior of adhesive, elastic-plastic spheres, Powder Technol. 99 (1998) 154–162. [16] D.A. Gorham, A.H. Kharaz, The measurement of particle rebound characteristics, Powder Technol. 112 (2000) 193–202. [17] P.P. Chattopadhyay, I. Manna, S. Talapatra, S.K. Pabi, A mathematical analysis of milling mechanisms in a planetary ball mill, Mater. Chem. Phys. 68 (2001) 85–94. [18] T. Inoue, K. Okaya, Grinding mechanism of centrifugal mills – a simulation study based on the discrete element method, Int. J. Miner. Process. 44-45 (1996) 425–435. [19] S. Agrawala, R.K. Rajamani, P. Songfack, B.K. Mishra, Mechanics of media motion in tumbling mills with 3D discrete element method, Miner. Eng. 10 (1997) 215–227. [20] P.W. Cleary, Prediction charge motion, power draw, segregation and wear in ball mills using discrete element method, Miner. Eng. 11 (1998) 1061–1080. [21] H. Mio, J. Kano, F. Saito, K. Kaneko, Optimum revolution and rotational directions and their speeds in planetary ball milling, Int. J. Miner. Process. 74S (2004) S85–S92. [22] B.K. Mishra, A review of computer simulation of tumbling mills by the discrete element method: part II – practical applications, Int. J. Miner. Process. 71 (2003) 95–112. [23] Y.T. Feng, K. Han, D.R.J. Owen, Discrete element simulation of the dynamics of high energy planetary ball milling processes, Mater. Sci. Eng. A 375-377 (2004) 815–819. [24] H. Kruggel-Emden, E. Simsek, S. Rickelt, S. Wirtz, V. Scherer, Review and extension of normal force models for the discrete element method, Powder Technol. 171 (2007) 157–173. [25] S. Lu, P. Zhou, T. Chai, W. Dai, Modeling and simulation of whole ball mill grinding plant for integrated control, IEEE Trans. Autom. Sci. Eng. 11 (2014) 1004–1019. [26] C.C. Doumanidis, I.E. Gunduz, C. Rebholz, C.C. Doumanidis, Real-time computational model of ball-milled fractal structures, ASME J. Nanotechnol. Eng. Med. 6 (2015) 031001-031001-6. [27] J. Perrin, Les Atomes (Trans. D.L. Hammick), Constable & Co Ltd, London, 1916. [28] A. Einstein, Investigations of the Theory of Brownian Movement, Dover, 1956. [29] D. Revuz, M. Yor, Continuous Martingales and Brownian Motion, second ed. Springer-Verlag, 1994. [30] F. Mandl, Statistical Physics, Manchester Physics, second ed. J. Willey & Sons, 2008. [31] J.C. Maxwell, Illustrations of the dynamic theory of gases part I: on the motions and collisions of perfectly elastic spheres, Phil. Mag. 19 (1860) 19–32. [32] J.R. Forshaw, A.G. Smith, Dynamics and Relativity, Wiley, 2009. [33] V.L. Popov, Rigorous Treatment of Contact Problems-Hertzian Contact, Chapter 5 in Contact Mechanics and Friction, Springer-Verlag, 2010. [34] J.A. Williams, R.S. Dwyer-Joyce, Contact between Solid Surfaces(Chapter 3) CRC Press, 2001. [35] J.D. Eshelby, The elastic field outside an ellipsoidal inclusion, Proc. R. Soc. Lond. A252 (1959) 561–569. [36] Y.N. Podil'chuk, Stress state of a transversely-isotropic body with elliptical inclusion, Int. Appl. Mech. 33 (1997) 881–887. [37] S.P. Timoshenko, J.N. Goodier, Theory of Elasticity, McGraw-Hill, 1971. [38] J. Li, Q. Du, C. Sun, An improved box-counting method for image fractal dimension estimation, Pattern Recogn. 42 (2009) 2460–2469.
Contents lists available at ScienceDirect
Powder Technology journal homepage: www.elsevier.com/locate/powtec
Brownian-like kinematics of ball milling for particulate structural modeling Constantine C. Doumanidis a, Hamda A. Al Kaabi b, Abdelaziz S.M. Alzaabi b, Ibrahim E. Gunduz c,⁎, Claus Rebholz d, Charalabos C. Doumanidis b a
Pankyprion Gymnasion, Nicosia, Cyprus Department of Mechanical Engineering, Khalifa University, Abu Dhabi, UAE c School of Mechanical Engineering, Purdue University, W. Lafayette, Indiana, United States d Department of Mechanical & Manufacturing Engineering, University of Cyprus, Nicosia, Cyprus b
a r t i c l e
i n f o
Article history: Received 21 March 2016 Received in revised form 20 May 2016 Accepted 9 July 2016 Available online 19 July 2016 Keywords: Ball milling Brownian motion Modeling Microstructure NiAl
a b s t r a c t Ball milling motion has been previously studied through computationally expensive, off-line experimental video processing and numerical simulations by the discrete element method. This research establishes a more efficient formulation of the ball energetics and kinetics similar to the Brownian kinetic theory of statistical mechanics. Based on assumptions of thermomechanical equilibrium, negligible gravitational, aerodynamic and surface condition effects, and decoupled impact interaction among balls and with milled particulates, this model proposes mono-parametric spectral energy and velocity probability density functions akin to Maxwell-Boltzmann statistics, along with uniformly distributed impact directionality. The model predictions are calibrated and validated by comparison with published experimental measurements and computationally derived spectra. This descriptive Brownian-like motion model enables effective simulation of contact and impact, material deformation and micro-joining of ball milled bimetallic powders. A comprehensive simulation of the evolving internal fractal microstructure of the processed particulates is implemented at real-time computation speed, and its predictions are compared with experimental micrographs of ball milled Ni\\Al particulates. © 2016 Elsevier B.V. All rights reserved.
1. Introduction Milling of particulates, originally powered by humans, animals, water and wind, has been one of the prehistoric fundamental manufacturing methods for processing of rock and soil, construction materials, foodstuff and medicines [1]. Through compressive impact and/or shearing between hard surfaces, resulting in plastic deformation and/or crushing fracture, milling has been used for comminution, i.e. reduction of volumetric size and increase in active surface, coagulation or sharpening of particle shape, as well as mixing, agglomeration and compounding of diversified solid matter. Material synthesis through rolling, tumbling, rod and ball milling resulting in mechanical alloying of metals has been of recent interest, particularly in fabrication of reactive particulates out of bi-material powder mixtures such as Al\\Ni [2– 5]. The resulting multi-scale Apollonian packing or lamellar micro/ nano-structure of the milled particulates, upon their thermal or electrical ignition, determines the kinetics of self-propagating exothermic
⁎ Corresponding author. E-mail addresses: [email protected] (C.C. Doumanidis), [email protected] (H.A. Al Kaabi), [email protected] (A.S.M. Alzaabi), [email protected] (I.E. Gunduz), [email protected] (C. Rebholz), [email protected] (C.C. Doumanidis).
http://dx.doi.org/10.1016/j.powtec.2016.07.033 0032-5910/© 2016 Elsevier B.V. All rights reserved.
reactions (SPER) [6,7] with highly-localized and time-controlled enthalpy release and high-temperature synthesis (SHS) of refractory compounds. Application of such nanoheater particulates [8] in rapid thermal processing, micro-joining and micro-soldering of electronics, self-sintering powders and self-healing composites, all the way to biomedical thermal ablation of tumors, has kindled research in understanding the process-structure relationships of milling energetics with the resulting self-similar (fractal) bimetallic microstructures in mechanically alloyed particulates. Ball milling in rotating uniaxial cylindrical vial or planetary settings (Fig. 1a) has been implemented for this purpose in both high-energy (interrupted adiabatic or actively cooled continuous) and low-energy (continuous or interrupted isothermal) configurations, and has been extensively investigated by laboratory studies and computational simulation. Experimental observation of the ball motion in the milling vial through a transparent lid by high-speed (HS) optical recording [9–13] has elucidated and empirically described their kinematics in stagnant, cascading, cataracting and centrifuging modes of ball trajectories [5]. However, laboratory process monitoring via HS video is computationally costly in its requisite image post-processing. Nevertheless this provided the basis for identification of (visco)elastic hysteresis/plastic collision and frictional parameters of material models in the theoretical linear or non-linear dynamic equations of the ball motion [14–17], such
1078
C.C. Doumanidis et al. / Powder Technology 301 (2016) 1077–1084
Fig. 1. a. Planetary ball mill setup, b. Milling ball-particulate impact.
as coefficients of restitution, rebound angles, contact duration and forces. Such collision formulations enabled calibration and validation of multi-body interaction simulations in the milling vial by the discrete element method (DEM) [18–24]. Numerical DEM models have provided valuable information on impact frequencies, ball wear rates, heat generation and power requirements for the ball milling process. However, complex Lagrangian simulations of multiple moving balls in a rotating vial by DEM are computationally inefficient, typically requiring parallel or super-computing resources for prolonged ball milling process times, while often missing description of mechanical interactions with and structural effects to the milled particulates. Therefore current numerical models, along with laboratory sensing techniques, are suitable for offline estimation and statistical process control (SPC) only, requiring counter-productive and invasive process interruption for sample removal and characterization to correct ball milling conditions [25]. Recently, real-time modeling of the particulate microstructure for in-process feedback control of ball milling was introduced, based on a computationally cost-effective semi-analytical predictive simulation of the dynamic product geometry and material composition [26]. The latter hinges on an Eulerian formulation, following and describing a typical particulate as it develops during the process, through its collisions with the balls and vial. These result in impact loading, contact micro-bonding with other particulates, material deformation, heat dissipation, gradual diversification of particulate species and their population proliferation in the final milled product. The sequence of these phenomena commences with and critically depends on the energetics and kinetics of the milling balls and vial. For compatibility with this real-time framework and instead of a full DEM simulation, a comprehensive yet parsimonious model of ball milling motion is required. This paper elaborates on such an ab-initio analytical, theoretical physics-based, experimentally calibrated and computationally validated descriptive model of the milling ball energy and velocity stochastic distributions, along with its integration into the particulate structure simulation. 2. Theory 2.1. Ball milling motion During steady-state ball milling, the previous experimental and numerical literature elucidates a continuous, time-invariant, chaotic motion pattern of the ball population in a confined vial, consisting of deterministic free trajectory segments between random collisions with other balls and vial walls, under predictable impact conditions. Upon distant zoomed-out observation from a larger dimensional (macro) scale, this Wiener process/Lévy flight pattern is reminiscent of the familiar nanoscale Brownian motion of ideal gas atoms/molecules at thermal equilibrium, with well-established descriptions by kinetic theory [27–29]. Despite similarities across scales, however, a number
of salient differences in milling ball motion behooves certain assumptions before proceeding to modeling analogies. 2.1.1. Steady state thermomechanical equilibrium The Brownian model involves fully reversible, perfectly elastic collisions between spherical particles and with stationary, adiabatic container walls. In ball milling, irreversible mechanical energy losses into heat during impacts, as denoted by the coefficient of restitution [22–24], are compensated during collisions with the mechanically moving (inducting work), thermally abducting vial walls. Thermomechanical equilibrium assumed at steady state implies constancy of the total mechanical energy of the balls, as well as the enthalpy of the system, i.e. isothermal ball milling conditions analogous to those in an ideal gas. This also ensures constant material properties of the balls and vial. 2.1.2. Negligible gravitational, aerodynamic and surface condition effects In a fixed-elevation ball milling setup, gravity introduces no mechanical work to the average ball energy, while limited vial height (Fig. 1a) has negligible influence to the standard deviation, especially in high-energy systems. In low-energy systems gravity bias may cause visibly parabolic trajectory segments of the balls instead of the straight ones in Brownian motion. Aerodynamic drag and lift (through the vorticity of spinning balls) interactions with the inert gas in the setup are also negligible, since they scale as the mass of gas volume swept between collisions over that of the balls (i.e. by about 3 orders of magnitude lower). Dry ball milling without liquid process control agents is assumed, avoiding complex non-Newtonian viscous and surface tension/capillary flow effects. 2.1.3. Decoupled impact interaction among balls and with milled particulates As per kinetic theory, decoupled ball motion requires that their mean free path between successive collisions as analyzed below, should be larger than the characteristic geometric size, i.e. the harmonic mean of the ball and vial dimensions. While temporarily adopted for simplicity (i.e. balls only in vial), the assumption on ball-particulate interactions will be explicitly relaxed in the sequel to reflect transport of mechanical energy to the processed particulates during collisions. The (normal and shear) impact energy transported to the particulate contact area, and the corresponding restitution of the impacting ball or vial, are determined by the relative linear and spin velocity components of the impactor surface at the contact location. Such an effective composite velocity v at the surface of a vial wall, e.g. in a planetary mill with angular velocity Ωο and radius Ro of the sun, and respectively Ω and R on the planet surface as in Fig. 1a, is: V 2 ¼ ðΩRÞ2 þ ðΩo Ro Þ2
ð1Þ
C.C. Doumanidis et al. / Powder Technology 301 (2016) 1077–1084
1079
Similarly, at the surface of a ball of radius r and mass m the effective velocity v corresponds to its total kinetic energy E (translational and rotational), i.e. (Fig. 1b): E¼
1 1 1 2 1 2 mυ2 þ Jω2 ¼ m υ2 þ r 2 ω2 ¼ mv2 →v2 ¼ υ2 þ ðωrÞ2 ð2Þ 2 2 2 5 2 5
where J = mρ2 is the ball moment of inertia around its rotation axis, pffiffiffiffiffiffiffiffiffiffi and ρ ¼ 2=5 r its gyration radius. 2.2. Ball energy and velocity distribution Besides the deterministic motion of the vial, the experimental evidence [5,12] points to a statistical distribution of the ball effective velocity v and energy E similar to that of Brownian motion. According to statistical mechanics [29,30], in an ideal gas or solution, the number N of atoms/molecules with kinetic energy E is proportional to its inverse exponential, i.e. N ∝ exp (− E/kT), where k is the Boltzmann constant and T the gas temperature. Integration over the entire population of atoms/molecules yields their average kinetic energy hEi ¼ 12 mhv2 i ¼ 32 kT for 3-D motion. For each 1-D motion component (k = x,y,z), the average energy hEk i ¼ 12 mhv2k i ¼ 12 kT ¼ 12 ma2 , where the parameter a pffiffiffiffiffiffiffiffiffiffiffiffi ¼ kT=m. By analogy to statistical mechanics, the energy E of the milling balls of mass m follows a probability density function (pdf) f(E) (Fig. 2a): 2 f ðEÞ ¼ pffiffiffi 2 πma
rffiffiffiffiffiffiffiffiffi E E exp − ma2 ma2
ð3Þ
In this distribution, the most probable energy is Ep = ma2, and the average energy hEi ¼ 32 ma2 = 1.5Ep. qffiffiffiffiffiffiffiffiffi The parameter a ¼ hv2k i of this distribution is the root mean square (RMS) velocity in one direction of motion k, and can be conveniently measured experimentally from a video sequence of the ball motion if available. The measured a = γV can be related to the vial cylinder wall velocity V (Eq. (1)) via an empirical efficiency factor γ. Alternatively a can be determined in the laboratory through the steady-state average energy 〈E〉 of the balls in the vial, from power measurements during the initial transient to steady operation or the final transient to rest of the milling system (Fig. 3). The steady state mechanical energy stored in N balls N〈E〉 , the rotating vials (Eq. (1a) and other parts EV = Σ1/ 2JjΩ2j of the milling system, plus its enthalpy H = CT at thermal equilibrium, originates from the power P = MΩm supplied by the motors minus the thermal power Q = T / Rq removed from the setup (shaded
Fig. 3. Initial and final transient of ball milling power.
area in Fig. 3): Z τ X 1 3 NhEi þ Ev þ H ¼ J ω2 þ CT ðP−Q Þdt→N ma2 þ j2 j j 2 0 Z τ ¼ MΩm −T=Rq dt 0
ð4Þ
where Jj and Ωj are the moments of inertia and angular velocities of rotating elements, M and Ωm the motor torque and rotation speed (measured via the electrical current and voltage), C the heat capacity of the system at uniform equilibrium temperature T (measured e.g. via IR thermocouples), and Rq the total heat loss resistance of the system, with C and Rq being experimentally measurable parameters. pffiffiffiffiffiffiffiffiffiffiffiffi Through its relation with energy E, the effective velocity v ¼ 2E=m at the ball surface (Eq. (2)) is also statistically distributed, similar to a Maxwell-Boltzmann pdf [30,31] derived from Eq. (3) (Fig. 2b): pffiffiffi 2 2 v2 v2 ϕðvÞ ¼ pffiffiffi exp − 2 2 2a πa 2a
ð5Þ
pffiffiffi The most probable velocity in this pdf is vp ¼ 2a, the average velocqffiffi qffiffi pffiffiffiffiffiffiffiffiffi pffiffiffi ity is hvi ¼ 2 π2a ¼ p2ffiffiπ vp, and the RMS velocity hv2 i ¼ 3a ¼ 32vp in 3-D ball motion. 2.3. Impact with particulates In the presence of milled particulates in the vial along with the balls, it is now assumed additionally that both these elements follow kinetic energy E and effective velocity v distributions with pdfs similar to Eqs. (3) and (5) respectively, with separate parameters. Moreover, for simplicity the Eulerian collision simulation described below [26] follows
Fig. 2. a. Milling ball energy pdf (m = 0.03). b. velocity pdf.
1080
C.C. Doumanidis et al. / Powder Technology 301 (2016) 1077–1084
the motion of one evolving particulate, assumed stationary and impacted by impinging vial walls and ball surfaces (Fig. 1) moving at relative velocities vi also according to a Maxwell-Boltzmann-like pdf. The equivalent (or reduced) mass m* of the colliding two-body system is the harmonic sum of the particulate mass μ and the impactor (ball or vial) m [32]: m ¼
μm μ þm
ð6Þ
typically resulting to m* ≈ μ since μ ≪ m in ball milling. Therefore the kinetic energy of the impactor before collision Ei ¼ 12 m v2i and after it Eo ¼ 12 m v2o depends on the relative input vi and output vo velocities, with the coefficient of restitution defined as e = vo/vi. During each such collision, from the impinging vial or ball energy Ei the particulate receives an amount of elastic work w = Eo which is reversibly restored as kinetic energy to the impactor upon separation, as well as an inelastic work W retained irreversibly by the particulate and eventually converted to heat. The latter consists of energy Wr stored in hysteretic elastic and residual stress fields (which is subsequently thermalized by stress relieving, annealing, etc.), plastic deformation work Wp to the materials, as well as frictional work Wf at the contact surfaces. This irreversible (positive) work W is therefore: Ei −Eo ¼
1 2 1 2 m vi − m vo ¼ W ¼ W r þ W p þ W f N0 2 2
Fig. 5. Particulate surrounded by inner cluster and outer medium layers while impacted by milling ball for its structural simulation.
ð7Þ
implying that vi N vo and therefore e b 1. Therefore for a certain impact event, while the input velocity vi can be randomly selected from the full range of pdf ϕ(vi) (Eq. (5)), the resulting output velocity is statistically limited in the range [0...vi] corresponding to the already selected vi. Therefore vo must follow a finite-range pdf ϕ(vo) / Φ(vi), normalized by the value of the cumulative distribution function (cdf) Φ(v) at vi: Z ΦðvÞ ¼
v 2 v v2 ϕðυÞdυ ¼ erf pffiffiffi − pffiffiffi pffiffiffi exp − 2 2a π 2a 2a 0 v
ð8Þ
This also implies that, while the corresponding incoming energy Ei to the collision is randomly selected from the full-range pdf f(Ei) (Eq. (3)), the concomitant out-coming energy Eo is limited in the range [0…Ei] of the pdf f(Eo) / F(Ei), i.e. normalized by the value of cdf F(E) at Ei: Z F ðEÞ ¼
rffiffiffiffiffiffiffiffiffiffiffi! rffiffiffiffiffiffiffiffiffiffiffi E 2 E E p ffiffiffi f ðε Þdε ¼ erf − exp − m a2 m a2 π m a2 0 E
ð9Þ
The kinetic energy Ei of the impactor due to its incoming velocity vi is imparted to the particulate via a normal vn and a tangential vt component (Fig. 1b), corresponding to a normal En ¼ 12 m v2n at and a tangential
Fig. 4. 2-D warped ellipsoidal primitive of particle domain for particulate structure simulation.
Fig. 6. a. DEM simulation [22] vs kinetic spectra (continuous line) of ball energy during ball milling. b. Experimental (1st column) and DEM simulation (2nd column) [19] vs kinetic energy spectra.
C.C. Doumanidis et al. / Powder Technology 301 (2016) 1077–1084
Et ¼ 12 m v2t impact energy: 1 2 1 2 1 2 m vi ¼ m vn þ m vt ¼ En þ Et ¼ W þ Eo and v2i 2 2 2 ¼ v2n þ v2t
Ei ¼
ð10Þ
The random incidence angle θ = tan−1(vn / vt) of the impactor on the local contact surface of the particulate is assumed to follow a uniform pdf within the range [−π/2…π/2]. Finally the particulate surface is hit by an impactor which is either one of N balls (of curvature 1/r), a vial top or bottom cap (curvature 0) or a vial cylinder (curvature 1/R). Because of the uniform pdf of the incidence angle θ, the respective probabilities are proportional to the solid angles subtended by each impactor type to the particulate, i.e.
1081
~ Nr2 / 4 bD N2 for a ball, ~ 2(R/h)2 for a vial cap and ~ 1/2(h/R) for a vial cylinder. The collision frequency f = N/b tN of N balls at all energy levels in the vial of volume O = πR2h depends on the average flight pffiffiffiffiffiffiffiffiffi time b tN = b DN/ hv2 i of a single impactor between successive collipffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffi sions across a mean free path hDi≈ 3 O=N , with hv2 i the RMS ball velocity. 2.4. Particulate contact, deformation and joining Statistical selection of impact energy Ei or velocity vi and incidence angle θ according to the previous experimentally validated analytical pdfs is instrumental in the implementation of a real-time simulation of particulate microstructural evolution in ball milling [26]. Consistent
Fig. 7. Experimental [5] vs Maxwell-Boltzmann (dashed line) velocity spectra for various values of speed ratio K = Ω / Ωo of planetary ball milling.
1082
C.C. Doumanidis et al. / Powder Technology 301 (2016) 1077–1084
with experimental planar micrographs of such self-similar Apollonian and lamellar structure formations in bimetallic ball-milled materials [2–8], this model is a 2-dimensional, dynamic simulation of the crosssectional fractal geometry and material distribution in one typical particulate as it develops. The model assumes a versatile warped ellipsoidal primitive geometry for each monometallic particle domain, adequately describing its deformation from initially spheroidal powder to the final lamellar layer form, through variations of its position, orientation and dimensional descriptive parameters (Fig. 4). It also adopts a linear elasto-plastic strain-hardening behavior of each material volume, with dry Coulomb friction on its surface, and temperature-invariant (steady-state isothermal), impact high-rate specific parameters [24]. Before each impact event, this original primitive or developed composite particulate is surrounded by statistically selected additional clusters, similar to the current particulate and/or its prior development generations as explained below, and brought into random but stable contact by relative surface rolling of each cluster (Fig. 5). During each collision, this adjacent inner layer is considered for potential integration of each neighboring cluster to the particulate by micro-welding at their contact locations. The entire assembly is encircled by an outer boundary layer of particles, with contact loading insufficient for joining to the particulate during the impact event, but providing mechanical motion confinement and energy dissipation. A ball or vial impactor subsequently collides with the stationary particulate at a random contact site of its surface according to the previously established kinetics and energetics. During the progression of each impact event, an initially elastic Hertzian contact at the particulate-impactor ellipsoidal interface patch with respective equivalent contact curvatures and stiffnesses, absorbs the impact energy components (Eq. (10)) over the particulate bulk material [33,34]. This provisional elasticity assumption yields the normal and shear load distributions over the contact patch, along with its progressing footprint ellipsoid radii, penetration approach and displaced interface geometry. This contact loading of the particulate and surrounding inner cluster layer determines the normal and shear reaction distributions over their boundary surface with the outer restraint layer of materials, considered as an equivalent elastic circumscribing continuum with an ellipsoidal inclusion [35,36]. This outer boundary reaction loading is used for computation of the inner contact loads between the particulate and adjacent clusters (Fig. 5) under mechanical equilibrium conditions, along with all requisite internal elastic stress fields. Since such redundant contact conditions lead to indeterminacy of their solutions, Castigliano's theorem is employed to determine elastic deformation displacements at the contact locations and their resulting Hertzian conditions, and hence the tentative deformations of individual particulate and cluster domains and their internal elastic strain fields [37]. Next, the equivalent elasticity assumption above is relaxed and the total work w + W absorbed by each domain of the particulate assembly is differentiated into a reversible elastic component w and an irreversible inelastic one W as explained above. This is done according to the stress-strain diagram of its bulk material, therefore defining plastic work Wp and eventual deformation of the domain, as well as according to the load-slip conditions at its contact surfaces, defining frictional work Wf. When the specific thermalized energy distribution Wp + Wf over a certain contact surface area A locally exceeds a critical activation energy density Ew/A for joining of the bi-material system, a micro-weld is established over A. The load-carrying ability of this micro-joint is determined by the specific residual strain energy distribution, i.e. inelastic work Wr = W-Wp-Wf over A. Such surface joining of domain interfaces in the particulate also results in its integration with additional surrounding clusters; while separation of particulate clusters also occurs when joint strength is exceeded upon subsequent impact loading. Therefore, sequential collisions with the impactors cause a progressive evolution of the external geometric size and shape along with internal material structure of the particulate under study, yielding successively new particulate generations along with previous ones in
the setup, all the way back to original mono-material powders. Since the other balls and vial impactors, besides the one involved in the collision with the studied particulate, also act simultaneously during ball milling, the populations of particulate generations grow in parallel according to a logistic parabola law. The pdf of each particulate generation is composed of products of the population sizes for combinations of its prior generations. These pdfs are used in assembling the adjacent cluster layer around the particulate before each collision; and they reflect the eventual demographics of particulate diversity in the ball milled product.
3. Results and discussion Fig. 6a compares the energy spectrum f(E) of the physics-based pdf in Eq. (3) (a = 0.36 m/s, γ = 0.43) with those derived in [22] by experimentally calibrated DEM simulation for a ∅90 × 14 cm cylindrical vial with 32 ∅5.08 cm stainless steel balls rotating at 18 rpm. The statistical model satisfying the decoupled ball interaction condition fits rather well, but somewhat over-predicts impact frequencies at lower ball energies, while it under-predicts especially at the very high energy end. Fig. 6b also shows a comparison of the energy distribution predictions (a = 0.99 m/s, γ = 0.71) with those in [19] by experiment and DEM simulation, for a ∅90 × 4 cm vial with 150 ∅2.9 cm alumina balls, rotating at 30 rpm for 30 s. The statistical model of decoupled impacts fits the measured impact frequencies within ~10%, again with over-prediction at low ball energies (possibly due to missed impacts by the image processing [19]), while it somewhat under-predicts the DEM results. Fig. 7 illustrates a comparison of the velocity predictions ϕ(v) of Eq. (5) with experimental measurements in [5] for planetary ball milling (Fig. 1), with Ro = 10.4 cm, R = 4 cm, h = 5 cm, ball r = 3 mm and volume fraction 0.35 in the vials, and sun Ωo = 694 rpm, for various planet rotation speeds Ω = KΩo within 0–1388 rpm. The distribution parameters fitting the laboratory data are: Speed ratio K = Ω/Ωo Pdf parameter a (m/s) Efficiency parameter γ
1.0 1.0 0.12
1.5 0.92 0.094
1.8 0.65 0.071
1.9 1.98 0.21
2.0 2.41 0.25
It can be seen that the statistical distribution slightly over-predicts the ball fraction at very low velocities as explained before, while it somewhat under-predicts at high velocities. A plausible reason for these deviations, also manifested in the variance of the efficiency parameter γ, is the high bulk charge challenging the decoupledness assumption for ball impacts. The variation in pdf parameter a and efficiency parameter γ can also be related to the nature of ball motion observed in [5]. The first two speed ratios correspond to mostly rolling
Fig. 8. Experimental (x) and DEM simulated (solid line) [24] vs kinetic theory predictions (dashed line) of coefficient of restitution e = vo/vi during stainless steel ball milling.
C.C. Doumanidis et al. / Powder Technology 301 (2016) 1077–1084
1083
Table 1 3D Fractal dimensions of experimental and simulated sections at 2, 4, 6 and 8 ball milling hours. BM time/fractal dimension Capacity Dimension Information Dimension Correlation Dimension Probability Dimension
Fig. 9. DEM simulated [23] vs kinetic theory predicted (dashed lines) of normal impact velocity vn. during planetary ball milling with 10 and 20 balls.
motion (cascading) of balls on balls with similar parameters a and γ. The last two ratios correspond to a falling motion (cataracting) where the balls are rolling mostly on vial surfaces and crossing across the vial walls during a part of the rotational cycle. These motions can be expected to produce both higher average ball speeds and high efficiency parameters due to longer and more intimate contact with the vial walls. The comparison of coefficient of restitution e = vo/vi determined by sequential selections of V = vi and vo according to the pdfs in Section 2.3 (with a = 0.31 m/s, γ = 0.042), to measured values and DEM predictions in [24] for stainless steel ball milling is shown in Fig. 8. The decoupled statistical formulation again results to some over-prediction of e for low velocities V, and a slight under-prediction for higher values. Fig. 9 shows a comparison of the statistical distribution of the normal impact velocity vn = vi sinθ (with a = 0.34 or 0.26 m/s, γ = 0.096 or 0.074), with the DEM predictions in [23] for planetary milling with Ro = 12.5 cm, R = 4 cm, h = 6 cm, Ωo = 300 rpm and Ω = −750 rpm, using N = 10 or 20 stainless steel balls of r = 5 mm. This confirms that components of the impact velocity under decoupled ball interaction conditions are also distributed randomly according to the Maxwell-Boltzmann pdf of Eq. (5), therefore justifying the assumption of uniform distribution of θ over the contact half-space. The calculated impact frequencies f ≈ 120 Hz for the conditions of [22] and f ≈ 25 Hz for those of [19], according to the Section 2.3, are consistent with the values of Fig. 6a and b, respectively. To validate the ball milling energetics formulation, the previously described simulation model [26] with decoupled impact kinetics following a Maxwell-Boltzmann pdf (Eq. (5)) was computationally
2 Hours Exp
4 h
6 h
8 h
Model Exp
Model Exp
Model Exp
Model
2.54 2.36 2.01 2.04
2.74 2.56 2.59 2.45
2.77 2.57 2.59 2.45
2.79 2.69 2.73 2.61
1.69 1.78 2.26 2.23
2.07 2.19 2.45 2.41
2.14 2.15 2.44 2.47
2.46 2.31 2.50 2.49
implemented for the experimental conditions of [4]. In [4], a low-energy planetary ball milling system with cylindrical stainless steel vial of volume 80 ml (Ro = 121 mm, R = 20 mm) containing 5 stainless steel balls of 2r = ∅10 mm in inert gas (N2) was employed, at room temperature and at rotation speed of 300 rpm (K = 1.82). This was loaded with a total mass 32 g of Al (99.5%) and Ni (99.95%) powders of 325 mesh size at a molar ratio of 1:3. The powder load was processed by continuous ball milling, briefly interrupted after 2, 4, 6 and 8 h for partial sample removal and analysis. These samples were sectioned, polished and characterized by scanning electron microscopy (SEM, back scattered electron BSE-first row) in Fig. 10. The micrographs illustrate a typical particulate structure evolution, from initial coarse globular Apollonian assemblies of harder Ni particles surrounded by softer Al phase, to planarized finer lamellar multi-layers of the two metals in the subsequent and final sections. The second row of Fig. 10 comparatively shows the respective particulate structure sections as predicted by the ball milling simulation model for the same processing times as before, with its kinetic parameters (a = 0.87 m/s, γ = 0.22) calibrated from the laboratory conditions. Similar to the experimental BSE micrographs above, initially larger warped ellipsoids of Ni (bright domains) appear as spheronized by the omni-directional random impactor collisions, directionally bent and joined with smaller Ni particles through intermediate highly deformed Al (dark domains) in an Apollonian-packed fractal architecture. With progression of the process, further impact to the particulate surface leads to more bonding and unification of similar metal domains into continuous layers, as well as dispersion of dissimilar layers with clearly discernible bilayer interfaces in lamellar formations. These multi-layers tend to planarize due to the higher probability of collisions normally to the flatter, larger radius surfaces of an ellipsoidal particulate (Fig. 4) compared to its curved, small radius sides. This is because of their respective difference in solid angles subtended to the directionally uniform incidence of impactors. Further deformation and refinement of the lamellae is evident in the finally produced microstructures by ball milling. For a quantitative comparison of the fractal dimension metrics in the experimental and simulated particulate sections of Fig. 10 (first and
Fig. 10. Particulate microstructure after 2, 4, 6 and 8 h of planetary ball milling Experimental SEM micrographs [4]: BSE (1st row) versus simulated internal structure of impacted particulate (2nd row).
1084
C.C. Doumanidis et al. / Powder Technology 301 (2016) 1077–1084
second row) after 2, 4, 6 and 8 h of ball milling, a box-counting algorithm [38] was implemented for their off-line image analysis with different measures. The respective attributes in Table 1 show comparative agreement of the model predictions to the laboratory data within ± 7% of the parameter values. This indirectly confirms sufficient calibration and suitability of the Brownian impactor kinetics for the particulate structure simulation. Therefore, unlike prior DEM simulations of ball milling motion, the computational efficiency of this kinetic theory formulation enables real-time implementation of the particulate structural model on a standard desktop computer running at 130 GIPS. 4. Conclusions The previous Brownian-like model of statistical ball milling motion is described by cost-efficient analytical probability density functions for the ball energy and velocity, and by spherically uniform directionality of impactor incidence. A single calibration parameter in a MaxwellBoltzmann expression, i.e. the RMS velocity of the balls in a single direction, can be either calibrated from video-measured energy or velocity laboratory spectra, or computed from the initial or final transients of the process conditions. The estimates of this descriptive model directly fit published energetic and kinetic spectra for various decoupled ball milling process conditions, while they also enable simulation of impacted particulate structures in low-energy planetary processing of bimetallic powders at real-time computation speeds on ordinary hardware. The convenience of this probabilistic kinetic theory-based model of ball milling is attractive for description of other milling processes, potentially challenging its original assumptions and calling for appropriate adjustment of its formulation and/or parameters. High-energy uncooled processing, for example, may result in dominant thermomechanical transients departing from steady-state equilibrium and isothermal assumptions. Temperature variation of the material parameters may also alter impact restitution and affect spectral distribution parameters. Use of liquid process control agents may also affect such parameters through concomitant viscosity and capillary surface tension effects during contact. Finally coupled interaction of ball impactors at high bulk volume charges or with larger milled particulates of commensurate mass and energy, as well as complex surface geometries and impact property distributions, could invalidate the decoupledness assumption of ball motion made above. Such conditions may necessitate coupledspecies, modified energetic and kinetic spectral expressions with dependencies on non-uniform impact directionality distributions, to be established through subsequent investigation. Acknowledgments This research was supported in part by a Khalifa University Internal Research Fund award (Level 1) and by the US Department of Energy, National Nuclear Security Administration under Award Number DENA0002377. Appendix A. Supplementary data Supplementary data to this article can be found online at http://dx. doi.org/10.1016/j.powtec.2016.07.033. References [1] A.J. Lynch, C.A. Rowland, The History of Grinding, Society for Mining, Metallurgy and Exploration, 2005. [2] A. Hadjiafxenti, I.E. Gunduz, C. Tsotsos, T. Kyratsi, S.M. Aouadi, C.C. Doumanidis, C. Rebholz, Synthesis of reactive Al/Ni structures by ball milling, Intermetallics 18 (2010) 2219–2223.
[3] K.V. Manukyan, B.A. Mason, L.J. Groven, Y.C. Lin, M. Cherukara, S.F. Son, A. Strachan, A.S. Mukasyan, Tailored reactivity of Ni + Al nanocomposites: microstructural correlations, J. Phys. Chem. C 116 (2012) 21027–21038. [4] A. Hadjiafxenti, I.E. Gunduz, C.C. Doumanidis, C. Rebholz, Exothermic reaction characteristics of continuously ball-milled Al/Ni powder compacts, Vacuum 96 (2013) 73–78. [5] A.S. Rogachev, D.O. Moskovskikh, A.A. Nepapushev, T.A. Sviridova, S.G. Vadchenko, S.A. Rogachev, A.S. Mukasyan, Experimental investigation of milling regimes in planetary ball mill and their influence on structure and reactivity of gasless powder exothermic mixtures, Powder Technol. 274 (2015) 44–52. [6] I.E. Gunduz, K. Fadenberger, M. Kokonou, C. Rebholz, C.C. Doumanidis, T. Ando, Modeling of the self-propagating reactions of nickel and aluminum multilayered foils, J. Appl. Phys. 105 (2009) 074903. [7] K. Fadenberger, I.E. Gunduz, C. Tsotsos, M. Kokonou, S. Gravani, S. Brandstetter, A. Bergamaschi, B. Schmitt, P.H. Mayrhofer, C.C. Doumanidis, C. Rebholz, In-situ observation of rapid reactions in nanoscale Ni-Al multilayer foils using synchrotron radiation, Appl. Phys. Lett. 97 (2010) 144101. [8] Z. Gu, Q. Cui, J. Chen, J. Buckley, T. Ando, D. Erdeniz, P. Wong, A. Hadjiafxenti, P. Epaminonda, I.E. Gunduz, C. Rebholz, C.C. Doumanidis, Fabrication, characterization and applications of novel nanoheater structures, Surf. Coat. Technol. 215 (2013) 493–502. [9] P. LeBrun, I.L. Froyen, L. Delaey, The modeling of the mechanical alloying process in a planetary ball mill: comparison between theory and in-situ observations, Mater. Sci. Eng. A161 (1993) 75–82. [10] R.K. Rajamani, B.K. Mishra, R. Venugopal, A. Datta, Discrete element analysis of tumbling mills, Powder Technol. 109 (2000) 105–112. [11] H. Mori, H. Mio, J. Kano, F. Saito, Ball mill simulation in wet grinding using a tumbling mill and its correlation to grinding rate, Powder Technol. 143-144 (2004) 230–239. [12] S. Rosenkranz, S. Breitung-Faes, A. Kwade, Experimental investigation and modeling of the ball motion in planetary ball mills, Powder Technol. 212 (2011) 224–230. [13] H. Shin, S. Lee, H.S. Jung, J.-B. Kim, Effect of ball size and powder loading on the milling efficiency of a laboratory scale wet ball mill, Ceram. Int. 39 (2013) 8963–8968. [14] M. Abdellaoui, E. Gaffet, The physics of mechanical alloying in a modified horizontal rod mill: mathematical treatment, Acta Mater. 44 (1996) 725–734. [15] C. Thornton, Z. Ning, A theoretical model for the stick/bounce behavior of adhesive, elastic-plastic spheres, Powder Technol. 99 (1998) 154–162. [16] D.A. Gorham, A.H. Kharaz, The measurement of particle rebound characteristics, Powder Technol. 112 (2000) 193–202. [17] P.P. Chattopadhyay, I. Manna, S. Talapatra, S.K. Pabi, A mathematical analysis of milling mechanisms in a planetary ball mill, Mater. Chem. Phys. 68 (2001) 85–94. [18] T. Inoue, K. Okaya, Grinding mechanism of centrifugal mills – a simulation study based on the discrete element method, Int. J. Miner. Process. 44-45 (1996) 425–435. [19] S. Agrawala, R.K. Rajamani, P. Songfack, B.K. Mishra, Mechanics of media motion in tumbling mills with 3D discrete element method, Miner. Eng. 10 (1997) 215–227. [20] P.W. Cleary, Prediction charge motion, power draw, segregation and wear in ball mills using discrete element method, Miner. Eng. 11 (1998) 1061–1080. [21] H. Mio, J. Kano, F. Saito, K. Kaneko, Optimum revolution and rotational directions and their speeds in planetary ball milling, Int. J. Miner. Process. 74S (2004) S85–S92. [22] B.K. Mishra, A review of computer simulation of tumbling mills by the discrete element method: part II – practical applications, Int. J. Miner. Process. 71 (2003) 95–112. [23] Y.T. Feng, K. Han, D.R.J. Owen, Discrete element simulation of the dynamics of high energy planetary ball milling processes, Mater. Sci. Eng. A 375-377 (2004) 815–819. [24] H. Kruggel-Emden, E. Simsek, S. Rickelt, S. Wirtz, V. Scherer, Review and extension of normal force models for the discrete element method, Powder Technol. 171 (2007) 157–173. [25] S. Lu, P. Zhou, T. Chai, W. Dai, Modeling and simulation of whole ball mill grinding plant for integrated control, IEEE Trans. Autom. Sci. Eng. 11 (2014) 1004–1019. [26] C.C. Doumanidis, I.E. Gunduz, C. Rebholz, C.C. Doumanidis, Real-time computational model of ball-milled fractal structures, ASME J. Nanotechnol. Eng. Med. 6 (2015) 031001-031001-6. [27] J. Perrin, Les Atomes (Trans. D.L. Hammick), Constable & Co Ltd, London, 1916. [28] A. Einstein, Investigations of the Theory of Brownian Movement, Dover, 1956. [29] D. Revuz, M. Yor, Continuous Martingales and Brownian Motion, second ed. Springer-Verlag, 1994. [30] F. Mandl, Statistical Physics, Manchester Physics, second ed. J. Willey & Sons, 2008. [31] J.C. Maxwell, Illustrations of the dynamic theory of gases part I: on the motions and collisions of perfectly elastic spheres, Phil. Mag. 19 (1860) 19–32. [32] J.R. Forshaw, A.G. Smith, Dynamics and Relativity, Wiley, 2009. [33] V.L. Popov, Rigorous Treatment of Contact Problems-Hertzian Contact, Chapter 5 in Contact Mechanics and Friction, Springer-Verlag, 2010. [34] J.A. Williams, R.S. Dwyer-Joyce, Contact between Solid Surfaces(Chapter 3) CRC Press, 2001. [35] J.D. Eshelby, The elastic field outside an ellipsoidal inclusion, Proc. R. Soc. Lond. A252 (1959) 561–569. [36] Y.N. Podil'chuk, Stress state of a transversely-isotropic body with elliptical inclusion, Int. Appl. Mech. 33 (1997) 881–887. [37] S.P. Timoshenko, J.N. Goodier, Theory of Elasticity, McGraw-Hill, 1971. [38] J. Li, Q. Du, C. Sun, An improved box-counting method for image fractal dimension estimation, Pattern Recogn. 42 (2009) 2460–2469.