A discrete elements simulation and analysis of a high energy stirred milling process

– Stirred bead mills used in industry allow to split up particules in suspension by agitating a milling medium. The multiple impacts may damage beads and thus reduce the stirred milling process eﬃciency. A discrete numerical approach including simple ﬂuid eﬀects is proposed and carried out in this paper to model and study ﬁnely the damage phenomenon. The simulation data allow to identify internal variables of the milling medium and to localize high energy zones. The order of magnitude of contact forces including ﬂuid contribution will enhance the bead wear law

bad-understood for beads.In fact, bead mass loss occurring by surface damages, lead to a progressive diminution of the bead size.These milling debris may contaminate the comminution particles.Moreover, experimental observations highlight that the milling efficiency can be optimized by modifying the bead size for a given operational parameter set [1].In this case, the milling efficiency is defined by the effective diminution of comminution particles median diameter.
To improve the milling performance, it is crucial to better understand the wear mechanisms and to take this phenomenon into account.Consequently, in the background of the identification of a wear law of milling medium, different types of ceramic beads of identical size have been used to mill two powders with different granulometry, shape and hardness in a laboratory stirred bead mill with agitating fingers.These experimental tests have shown different wear behaviors and resistances.Moreover, specific damages depending on comminution particles have been observed; this may characterize different wear mechanisms.Some scratching tests and Vickers indentation tests made on polished beads highlight two type of surface damages: padded stripping induced by plastic deformation and shelling by crack growth.These observations lead us to propose a wear law of milling beads based on the hypothesis that the comminution particles, which are caught in a bead vs. bead contact, are considered as sharp micro-indentor cf. Figure 1.This approach requires an accurate characterization of the beads vs. bead contacts during the stirred milling process and especially qualification and quantification of the contact reaction network.These forces in complement of geometrical and Article published by EDP Sciences material data allow to estimate the indentation contact surface, the indentation pad volume and to evaluate the statistics of communition particles caught between two beads.Indeed, the bead material leads to two limits of indentation loading; a plastic limit and a shell limit.To improve the contact study and therefore, the induced surface damages, we propose an analysis of the contact distribution in the milling chamber and the impact frequency above a given limit of force.The numerical approach of this milling problem have often been treated by fluid mechanics.In this case the mixed granular material and fluid is seen as an equivalent fluid with specific intrinsic properties [2][3][4][5].Nevertheless, these homogeneous simulations can not report the localized problem of contact which is a natural discrete problem.Only a discrete element method coupled with fluid actions is able to determine all internal variables of the milling granular medium.As the stress field is an internal variable of a continuum, the contact force network is that of a discrete medium.This network determination would considerably enhance the bead wear law.

Simulation of the stirred milling process 2.1 A discrete element method with fluid contribution
Simulation methods of fluid vs. particles interaction problems have been proposed in [6,7].Nevertheless, the first method is computational costly and both of them are very little adapted to determine the entire contact network of large multi-contact problems.In this study background we propose to use the non smooth contact dynamics theory a.k.a nscd [8,9].This method is able to solve multi-contact problems with the expression of dynamics equations in the α-contact local frame.Such transformations are achieved by the H α operator and H T α operator which transport velocities and forces from the particle frame to the α-contact frame and vice versa.For an integration time step i + 1 a linear relation is found between the relative velocity u α i+1 and the average impulsion p α i+1 on the interval [t i , t i+1 ].This one is associated with an interaction law such as, External forces F ext and dynamics effects are taken into account by the relative velocity free of contacts u α free,i and the Delassus operator W αβ = H T α M −1 H β which appears naturally.Several types of interaction law can be used between two contactors.The most used interaction law for granular materials is a Signorini contact condition on the normal direction 0 ≤ u n,i+1 ⊥ p n,i+1 ≥ 0 associated with a Coulomb friction law on the tangential direction.The multi-contact problem is solved implicitly by a Gauss-Seidel loop which goes through contacts, and stops when energetic and interpenetration criteria are satisfied.Despite the highly dynamic character of this problem, the energetic restitution of contacts is vanished.Indeed the energy of an impact between two particles in a viscous fluid is mainly dissipated by the generated flows moving away from the contact area.Therefore, the nscd method is well adapted to model contacts between rigid bodies.The beads vs. fluid interaction problem is simplified as a determination problem of fluid effects on beads.In fact, in a first approximation, we consider a stationary Couette flow between the enclosure and the shaft cf. Figure 2.This induces hydrodynamic forces on the granular medium.This simple fluid velocity profile is chosen to compare in a first step, contact forces orders of magnitude of dry milling simulations [10] and milling with fluid simulations.Moreover, effects of fluid disturbances behind agitating fingers on beads are neglected because of the few density of beads on this area.Thanks to this approach, position of a rigid body j can be determined as a function, among other things, of external forces F ext,j of j body including hydrodynamic contribution defined by Equation (2), with u f the velocity field of fluid, x j , u j , r j and m j position, velocity, radius and mass of the particle j, g the gravity acceleration, ρ the water specific mass and C 0. of fluid are stronger than viscous effects of fluid a priori.Indeed, if orders of magnitude of the drag efforts are compared: T μ from the Stokes theory of viscous fluid with low Re < 1 and T Γ from the Euler theory of perfect fluid, we found in our case, with μ the dynamic viscosity of water, r avg the average radius of beads and Δu avg the average difference between bead velocity and fluid velocity.Consequently, it seems reasonable to privilege the hydrodynamic contributions.

The stirred milling model
This numerical study simulates an industrial stirred bead mill available in laboratory.Its ∅97 mm cylindrical enclosure can hold two different granulometry of milling beads made in zircone ceramics, either of diameter close to 3 mm or to 1 mm.The milling granular medium and the suspension are moved by a ∅57.5 mm agitating shaft with fingers turning at 3000 tr.min −1 .Six ∅8 mm agitating fingers are uniformly distributed along the shaft circumference and repeat themselves each axial segment of 19 mm.Using symmetrical considerations, the simulated domain is limited between the median plane of fingers and the median plane between fingers.Axial limit conditions of this "layer" are given by two fictive planes which impede the axial egress of beads using frictionless contact conditions.This interaction does not influence radial and tangential displacements.Periodic conditions could be used but they require the simulation of a larger layer.Both population of beads have been tested and lead to simulations of 1400 and 38 500 bodies respectively.Nevertheless, the computation power available in the laboratory does not allow to finalize the simulation of milling beads of ∅1 mm in reasonable time.Consequently, the framework is restrained to milling beads of ∅3 mm.In contradistinction to a stirred ball mill, the stirred beads mill works with filling rates ranging from 80% to 90%.
Moreover the granular specificity of this problem requires the following simulation steps, 1.Generation of a granular sample -this must be included inside the restrained zone of the simulation and must correspond roughly to the filling rate.2. Deposit under gravity -the sample is submitted to a gravity field to check the filling rate.3. Milling simulation -the agitating shaft is rotated to simulate the stirred milling process.
The generation of the granular sample may be carried out with two methods.Pluviation under gravity allows to fill an enclosure bead by bead, preliminary knowing the granulometry.This method allows to generate directly a deposit sample under gravity and to check its filling rate in the same time.Nevertheless, this method is computational costly.The method used in this study is based on a geometric deposit of spherical beads [11].To verify an acceptable filling, the deposit volume is artificially increased by decreasing the agitating shaft volume.In a first step, particles are deposited by a geometric method which closes in on the real compacity under gravity.Then in a second step, they are submitted to a centrifugal gravity field.Once the beads deposited on the external cylinder of the stirred mill, the shaft recovers its initial dimension and the sample is submitted to a classic gravity field.This process allows to generate granular samples in a reasonable time, cf. Figure 4, filling until 85% of the mill enclosure.Once the filling is checked, the agitating shaft reaches a rotation speed of 314.16 rad.s −1 via a Heaviside function.The presence of fluid may dissipate a part of the impact energy but this first approach will consider that the contact laws between different elements correspond to Signorini-Coulomb contact condition with a friction coefficient of 0.1, excepted for contacts with the fictive planes.The milling simulation is made using 50 000 time steps of 2 μs, corresponding to 5 rotations of the agitating shaft.
The resulting data indicate that the first half of the simulation is transitory and depends on non realistic initiation of movement.However, the second half of the simulation emphasizes a stationary phase which testifies the elementary mechanisms occurring in the stirred milling process.To check effects of the velocity field of fluid, each simulation of "milling with fluid" is systematically compared with a similar simulation without fluid contribution, so called "dry milling".
3 Stirred bead milling analysis

Kinetic analysis
After two simulations of a milling with fluid and a dry milling of a granular sample of ∅3 mm, the velocity distribution and the contact forces are given by Figure 5 for a time step included in the stabilized phase.In the dry milling simulation, the granular medium behaves like rigid body turning around the mill axis.Frictions between beads and the cylindrical enclosure and the shaft, decrease the rotation of granular set in comparison with the rotation imposed by the agitating shaft.However, beads are in quite stable configurations in zone far away from the agitating fingers.It is not the same for the milling with fluid simulation.The fluid field presence, and particularly the fact that it vanishes on the cylindrical enclosure, slows down external beads and privileges the slippage of radial layers.This phenomenon is highly influenced by the fluid velocity profile.On the other hand, we observe, in specific time steps, some locally more stable granular configurations, i.e. blockage configurations.These structures better resist to loadings and generate strong contact reaction chains as shown in Figure 5.The force chains are privileged by the fluid that decreases the granular mobility on the cylindrical enclosure.To compare the kinetic of the granular set, an average rotation rate Ω a , around the agitating shaft axis e a , is calculated post mortem for each saved time step.For a n s particle set of total volume V and neglecting the self-rotation effect, Ω a is defined as, with ω j a the rotation rate of one bead j around the axis e a and χ j a ratio between the volume of a bead j and the average volume of beads.In our case, the bead granulometry is quasi constant and the ratio χ j is considered equal to 1.The rotation rate ω j a depends on the radial excentration and the tangential velocity of a particle j.
In our case we write, t )e j t = ω j a e a ∧ (x j • e j r )e j r e a = (1, 0, 0), e j r = (0, x j 2 , x j 3 ) , e j t = (0, −x j 3 , x j 2 ) with e j r the direction between e a and the particle j and e j t the tangential direction defining (e a , e j r , e j t ) as an orthonormal frame, cf. Figure 6.The evolution of this rotation rate during the test is given in Figure 7 for three different simulations: a milling with fluid, a dry milling with dense beads and a dry milling with light particles.This last simulation is carried out to check the global kinetic behavior to the specific mass of beads.Rough kinetic data converge around a rotation rate constant value.An asymptotic smooth function is used to compare results and is defined as )), with t the time, c 0 the stabilized rotation rate, c 1 the characteristic time and c 2 an updating time.The milling with fluid plot tends to stabilize sooner with a constant value of rotation rate lower than the dry mills.This emphasizes the breaking effect of the fluid on the granular medium.In every cases, the milling regime is considered stable after five characteristic times, i.e. in the second half of the simulation.Furthermore, we notice that for the industrial stirred mill studied, the granular kinetic behavior is not influenced by the specific mass of beads.

Contact forces analysis
To establish a descriptive bead wear law, it is of primary importance to characterize contact forces in the stationary milling regime with regards to their nature, their intensity and their frequency.The second half of the simulation is considered stable, so the investigations on contact forces are carried out only on this restrained domain.Data are averaged on this stable phase to smooth discrete effects.In Table 1, we give especially the maximal force that loads a bead during the stirred milling process and the average force associated to the entire granular set.Comparing either normal or tangential part of these forces, we notice that the stress level inside the medium is stronger in the simulation with fluid.These values correspond to observations made in Section 3.1 which associate a kinematic decrease to a contact forces increase.We notice that the contact forces inside the stirred mill with fluid are distributed on a lower contact number.This privileges the bead wear.We propose to study the impact frequency distribution of beads which are in a contact loading state larger than a critical reference force F c , for different F c values. Figure 8

Rotation axis
e a e j r e j t x j (x j • e j r )e j r v j (v j • e j t )e j t Fig. 6.Projections of position x j and velocity v j on the rotation axis normal plane (e j r , e j t ).
frequencies distribution on the milling beads in two application cases.The impact frequencies of the entire granular set are given in Table 2. Weak impacts appear more frequently in the dry mill than in the mill with fluid.However, for strong impacts, their frequency is larger in the milling with fluid simulation.It is not obvious to consider a "structure" in such a dynamic simulation, nevertheless, we observe a bead segregation which has a sensitivity to impacts function of their position in the milling chamber.This leads to an identification of bead subset that have a larger risk of premature wear.To highlight this phenomenon, we propose to analyze the distribution of volumetric density power dissipated by friction according to the radial granular layers.Friction being the major mechanism of the stirred milling process, this mapping allows to localize the milling efficient zone.This power density is calculated on 12 radial layers of the stirred mill and results are given in Figure 9.For the total volume, the friction power is larger in the dry mill, cf.Table 1.This is explained by high relative velocities and in particularly near the static enclosure.The mill with fluid keeps a high friction power close to the cylindrical enclosure, but there is a non negligible activity near the end of the agitating fingers influence zone.With the lack of geometric elements on the cylindrical enclosure which would disturb the granular structure of the external ring, beads are driven away by the centrifugal effect and become trapped in a zone frequently loaded by important contact forces.This privileges their premature wear.

Conclusions and prospects
In the framework of the wear study of milling bodies inside high energy stirred bead mills, this work allowed to carry out some discrete simulations including hydrodynamic contribution of a fluid in movement.Comparisons have been established to identify tangible effects of a simple fluid contribution on a granular medium.This new numerical approach of the milling problem complements knowledge of stirred bead mills with some accurate investigations on the kinetic behavior and on the contact network of these granular sets.We showed that the simulated industrial stirred milling process generates preferential friction zones leading to a premature bead wear.Another numerical campaign including beads of ∅1 mm is being carried out to complete this study cf. Figure 10.Collaborations are made to enhance the fluid velocity field, taking into account the vortex zone below the agitating fingers.This will allow to simulate more realistic bead kinetics, to identify milling power and to localize their efficient zone.The fluid effect, while beads bring together, modifies the granular medium behavior and cannot be modeled by a classic Signorini contact condition.This phenomenon requires new interaction laws inside the discrete model including a dissipation coefficient lower than 1.Different geometries of agitating shafts could also be studied to optimize the stirred milling process.

Fig. 1 .
Fig. 1.Cut view of a milling granular medium including comminution particles.

Fig. 3 .
Fig. 3. Some details about the stirred beads mill and delimitation of the simulation zone.

Fig. 5 .
Fig. 5. Velocity fields and contact network after milling simulations with beads of ∅3 mm, with ||uj || the velocity norm of a particle j.

Fig. 7 .
Fig. 7. Evolution of the rotation rate of the granular medium according to time for several simulations.

Fig. 8 .
Fig. 8. Impact frequency[Hz]  for each bead with intensity larger than the critical reference force Fc.These frequencies are ordered from the most impacted to the least impacted.

Fig. 9 .
Fig. 9. Distribution of the volumetric density power dissipated by friction [mW.mm −3 ] of the granular medium for several simulations.The influence zone of agitating fingers is limited by the dashed line.

Fig. 10 .
Fig. 10.Velocity field of a stirred mill layer with a granular medium of ∅1 mm.

Table 1 .
Normal and tangential parts of contact forces [N], volumetric density of power [mW.mm−3] and contact number, averaged on stationary part of the numerical test.

Table 2 .
Total contact frequency on the simulated mill layer.