Heterarchical modelling of comminution for rotary mills: Part I — Particle crushing along streamlines

Rotary mills aim to eﬀectively reduce the size of particles through a process called comminution. Modelling comminution in rotary mills is a challenging task due to substantial material deformation and the intricate interplay of particle kinematics of segregation, mixing, crushing, and abrasion. Existing particle-based simulations tend to provide predictions that cannot cope with the large number of particles within rotary mills, their wide range of sizes, and the physics dictating the crushing of individual particles. Similarly, there is currently no deterministic modelling means to determine the evolving population of particle sizes at any point in time and space within the mill. The aim of this two-part contribution is to address these gaps by advancing a framework for a novel stochastic comminution model for rotary mills, which has a well-deﬁned deterministic continuum limit and can cope with arbitrarily large numbers of particles. This work describes the basic physics and structure of the new model within a heterarchical framework for ball and autogenous grinding mills. The primary focus of this Part I paper is to develop a computational model for the integration of motion of material along streamlines inside a mill. Coupled to this process is the kinetic physics dictating particle crushing. In a subsequent work, Part II, segregation and mixing will be added to this model such that realistic behaviour from the mill can be observed.


Introduction
Comminution is the process of reducing the size of particles, which in the mineral processing is achieved using different machinery such as crushers, rotary mills, vertical mills, and high pressure grinding rolls.These crucial unit operations are overwhelmingly energy intensive.The mineral processing industry is estimated to consume 3-4% of global energy [27,36], where nearly half of that energy is spent in the milling operation [58,26,57,6].Here, our focus is on the rotary mills.Cur-11 rently, there is no deterministic model that can 12 explain the particle size dynamics within these 13 mills.In this study, we introduce a novel approach 14 based on the multiscale heterarchical modelling this end, we follow the original ideas advanced by [40] to simultaneously describe the mechanisms of crushing, mixing, and segregation of particles.
Just like PBM models, the heterarchical model is capable of accommodating an arbitrarily large number of particles.However, unlike PBM models, the heterarchical approach can further provide information regarding the particle size distributions at any point in time and space within the rotary mill domains.This new approach follows a stochastic structure with a well-defined continuum limit for open system particle size dynamics such as segregation and mixing [40] and closed system ones such as crushing [33].This multiscale approach is an alternative to the other hierarchical multiscale approaches as here there is no scale separation between the representative volume and the continuum scale, both scales coexist within a single framework.This facilitates the transfer of information (i.e., particle size) from one scale to another as well as within the representative volume scale.Nevertheless, so far heterarchical models were limited to simple flow systems that may be idealised using only one spatial coordinate, such as those occurring during the granular avalanches, landslides, and debris flows [7,50,23,51].Furthermore, thus far the description of crushing in the framework ignored the role of kinetic particle collisions, dealing solely with quasi-static comminution.
In this study, we extend the heterarchical framework to track comminution in higher spatial dimensions.For the first time, the approach is used to model comminution in rotary mills.To deliver these goals the heterarchical approach is here coupled with the streamline method, originally called the strain path method in [5], which has been used in the past to integrate constitutive laws for predicting penetration resistance [24] in geomechanics.Here, rather than integrating constitutive models, the streamline method would be used to integrate the heterarchical method.
In summary, the current two-part study offers three significant novel contributions.Firstly, we extend the crushing model within the heterarchical approach to account for kinetic crushing, which goes beyond quasi-static conditions in previous publications.This factor is shown to be particularly important for rotary mills, where the particles are often advected through highly turbulent regions.Secondly, we integrate the heterarchy Fig. 1: Illustrative example for the structure of a heterarchical model.On the Left, a physical onedimensional polydisperse granular system.On the Right, an equivalent heterarchical model.The physical granular system is discretized into a number of representative volume elements (RVEs) along the physical coordinate.The particles within each RVE are then further mapped along the micro-structural coordinate into an array of particle sizes s t i = (s t i,1 , s t i,2 , ..., s t i,M ) at time t with periodic boundary (PB) conditions (colour figure online).To model comminution we need to define the conditions for the onset of particle crushing and the fragment size distribution for crushed particles.The onset of crushing occurs at any time t when the externally applied (bulk) stress (σ a ) is greater than or equal to the crushing strength (σ s i,z (t)) of the particles.If this criteria is met, at the next time step t + ∆t, the particle size at cell location (i, z) is reduced by some factor.Otherwise, the particle size in that cell is kept unchanged (see Fig. 2).The above comminution rules are applied to each cell of the RVE, and each cell of the RVE is subjected to the same external applied stress (σ a ).Accordingly, the particle size of the (i, z) th cell evolves as where X t i,z is the factor determining how much where x m = s min /s.The parameters k and λ can 268 be determined experimentally for the fragment 269 size distribution of any given mineral particle [64].
In the original work by Marks & Einav [40], for the z th cell of the i th RVE, the crushing strength of particles under quasi-static conditions was expressed as where ⃗ s = {s i,z−1 , s which can happen even when the external pressure 295 acting on them is relatively small.Fig. 3 depicts 296 such a collision between two particles, which 297 results in their crushing into smaller fragments.

298
Considering an assembly of jiggling particles, the effect of the collisions on the crushing strength under dynamic conditions (σ d i,z ) could be related to the (measurable [39]) mean fluctuating velocity (v ′ ), which represents the typical difference between particle velocities and their ensemble average at any particular point in space and time: where, σ s i,z is the crushing strength of particles  while K is introduced as a kinetic strength reduc-301 tion factor that is here taken to depend on the 302 fluctuating velocity v ′ and the particle size s i,z .

303
To determine this kinetic strength reduction factor one may consider the vast experimental observations [30,53,54], which show that the survival probability of the particles against crushing depends on their impact velocity.With the Weibull distribution [62] being a popular and successful representation of this survival probability, the kinetic strength reduction factor is here taken to follow the same statistics where, v ′ 50 = v ′ 50 (s i,z ) is the fluctuating velocity of a particle of size s i,z corresponding to 50% survival probability against crushing, and w v is the Weibull modulus for velocity.Both of these parameters (v ′ 50 and w v ) can be determined experimentally for a given mineral.Importantly, the work in [54] showed that v ′ 50 is also a function of the particle size (here s i,z ) and can be calculated as where a v and b v are correlation parameters and s r 304 is a reference particle size which in this case has Finally, notice that v ′ represents the magnitude of the velocity fluctuations, and thus can only be positive.Therefore, it is important to furthermore consider the typical direction of the fluctuations, since during expansion particles tend to move apart, so we would not expect collisions and associated crushing.Provided with the field of velocities (v) in the rotary mill, one can calculate the volumetric strain rate from εv = −∇ ⋅ ⃗ v to determine whether the medium expands or contracts.With that in mind, the particle size of a heterarchical cell in the general heterarchical comminution model for rate dependent, dynamic loading evolves according to the modified crushing criteria below where, H is the Heaviside step function, (i.e.dur-308 ing volumetric expansion εv < 0 and H( εv ) = 0, so 309 the particles do not crush, otherwise H( εv ) = 1, so 310 they may crush).

319
The methodology in this paper will follow a 320 one-way coupling scheme.Specifically, the use of 321 the kinematic and stress fields will be based on 322 either analytical or numerical proxies for the flow 323 of non-crushable particles.This has certain lim-324 itations since the fields retrieved ignore the fact 325 that the particles actually crush in those mills.

326
The tacit assumption here is that crushing would  386 Furthermore, in the case of the ball mill, the fluctuating velocity (v ′ ) of the particles (see Fig. 4e), as required from Section 2, can be derived as where P k is the kinetic pressure and ρ b is the bulk 387 density.

388
The kinetic pressure can be estimated for shear-induced flow in terms of the dimensionless kinetic number (I k ), as introduced by [2], which can be expressed as P k = P s I k , with P s being the static pressure.The calculation of bulk stresses acting on particles is discussed in the following subsection.Furthermore, the kinetic number can be approximated from [2] as I k = I, where I = εs s √ ρ m /P is the inertial number [20], εs the shear strain rate, s the mean particle size, ρ m the material density, and P the total pressure (sum of static pressure and kinetic pressure).Given these relations, the kinetic pressure may be estimated analytically from: Finally, the bulk static pressure (P s ) is assumed to be a linear function of depth (y cos ϕ) measured from the free surface, i.e.
where ρ b is the bulk density, g is the acceleration

AG mill 394
The numerical simulations from which we retrieve  have been established thanks to the coarse-grained velocity field derived from the DEM simulation.
Given this velocity field, the streamlines for the granular flow can be directly generated using MATLAB.However, given the inevitable numerical noise, the streamlines that are generated from MATLAB would not necessarily form the required set of closed streamlines.
To obtain closed streamlines, we adopt the following alternative approach.We start by selecting a subset of streamlines generated from MAT-LAB, ensuring they cover the entire flow region.
In order to close these streamlines, we strategi-

Material points on streamlines
Once the closed streamlines for the bulk material flow are obtained, the subsequent step is to allocate discrete set of material points on the streamlines.Each material point is defined as a

533
The position of a material point along the streamline at any time-step can be calculated using the Verlet integration scheme [34] as where ∆t is the time-step, v and a are the bulk The positions of the material points along the streamlines remain as such and do not need to be updated further during the simulation.

Tessellation into subvolumes
Given a set of material points along the streamlines, the material domain within the mill is tessellated into subvolumes.To compute the subvolume occupied by a material point in space we have used the Voronoi neighborhood method [1] as shown in Fig. 9.The volume will therefore be equal to the area of the Voronoi polygon assuming unit length along the mill axis.These Voronoi subvolumes are here required for the purpose of mass calculation, while in the following paper (Part II) they will be further required for implementing the segregation and mixing mechanisms.
Each of the Voronoi subvolumes are then further discretised along the micro-structural coordinate (z) into M cells as shown in Fig. 8.
Therefore, the whole volume of a material point is divided into M subvolumes, each containing particles of different sizes which evolve as the system undergoes crushing.

607
For the ball mill, the bulk density field is obtained by integrating the bulk density along the material points of a streamline using a suitable integrating scheme.Correspondingly, the material points are spaced using the analytical velocity field.As such, this method does in fact ensure the constancy of the material point masses.Here, using the forward Euler method, the bulk density field can be calculated as: where the indices p ∈ ∥1, N s ∥, and q ∈ ∥1, N p ∥ 608 represent the streamline and the corresponding 609 material point along it.N s is the total number of 610 Weibull modulus for velocity w v 6 6 Correlation parameters for impact velocity streamlines and N p is the total number of material points along a given streamline.εv is the volumetric strain rate and ∆t is the time step that we chose to discretise the streamlines into discrete material points.

Simulations and results
It is now possible to demonstrate how the newly defined 'heterarchical streamline method' may be used to analyse the process of comminution within due to impacts between the particles and therefore, the scaled mineral strength parameters in the simulation represent a weak crushable mineral.
For the AG mill simulation, we considered a mill with a diameter (d) of 6 m rotating at 13.5 rpm, and the feed consisting of particles uniformly distributed between 20 mm -100 mm in size.These operating parameters lie within the range that is typically associated with actual AG and SAG mills [43,10] of [30] for cement-mortar balls.Here, our assumption is that these values would not differ much for the minerals.In order to check the sensitivity of these impact velocity parameters on crushing, we do a parametric study towards the end of this paper.
For the conceptual ball mill simulation, we utilised a mill of diameter(d) of 30 cm rotating at 6.3 rpm speed, and the feed consisted of particles with a uniform size distribution ranging between 2 mm -3 mm.The value of the crushing strength parameter (σ m ) was chosen to represent a weak, crushable mineral and its value is taken two orders of magnitude lower than the crushing strength of the mineral used in the AG mill simulations.Similarly, we also scale the impact velocity parameters for the ball mill based on the corresponding value obtained for the AG mill.Given the parameters above, and the heterarchical streamline method, we can track the evolution of particle size at any point in space and time inside the ball mill.

693
The CDF plots provide clear insights into the influence of the fluctuating velocities on the crushing of particles in the granular flow.Notably, streamlines passing through regions of higher fluctuating velocity experience more crushing.In addition to generating CDF plots for individual streamlines, it is instructive to plot the overall or ultimate CDF that represents the entire material within each given mill, after one mill revolution.This product CDF is the weighted average by mass of the CDFs for each streamline, which is specified by where F p,q (s) is the cumulative mass distribution function, M p,q is the mass of particles in a material point for a given streamline and M total is the total mass of particles in the system which can be calculated as The results of the ultimate CDF of the par- The bulk density ρ b and velocity v fields are obtained using:  The fluctuating velocity v ′ is obtained using where ⟨v⟩ is linearly interpolated at the location

15 paradigm [ 40
] to unveil the particle size dynamics 16 within these mills.
and the streamline methods to investigate the pro-

Fig. 2 :
Fig. 2:The stochastic comminution rule of the heterarchical model for a cell (i, z) under externally applied stress σ a .At any given time t, if the external stress acting on the particles (σ a ) is greater than or equal to the crushing strength (σ s i,z ) of particles, the particle size in that cell is reduced by a factor X t i,z in the next time step t + ∆t.Otherwise, the particle size remains unchanged (colour figure online).

Fig. 3 :
Fig. 3: Schematic representation of kinetic crushing of particles during high deformations.Each particle experiences an independent velocity that fluctuates from the local mean velocity field, as depicted by arrows over the particles.The particles inside the highlighted box (left) have a higher relative velocity at time t, which may result in crushing.The figure on the right shows the crushed state of the particles at time t + ∆t.The length of the arrow (not to scale) represents the magnitude of the velocity of the particles.

305
been taken as equal to the maximum particle size, 306 i.e. s r = s m .

312
As shown in the previous section, the crushing 313 of particles within rotary mills depends on local 314 continuum information regarding kinematic enti-315 ties such as the velocities and strain rates, as well 316 as stress variables.Therefore, in the following, we 317 describe how these fields could be obtained for the 318 granular material within rotary mills.

327 not substantially affect the measured continuum 328 velocityFig. 4 :
Fig. 4: Kinematics for a ball mill using the analytical solution by Ding et al. [22].(a) Flow profile for the rolling regime in a ball mill.Contour plots: (b) Velocity magnitude (streamlines superimposed in white colour), (c) Volumetric strain rate (positive in compression), (d) Shear strain rate, and (e) Fluctuating velocity (colour figure online).

Fig. 5 :
Fig. 5: Coarse-grained kinematics from DEM for an autogenous mill (AG).(a) Snapshot of particle distribution and their sizes for the cataracting regime in an AG mill.Contour plots: (b) Velocity magnitude (streamlines superimposed in white colour), (c) Volumetric strain rate (positive in compression), (d) Shear strain rate, and (e) Fluctuating velocity (colour figure online).

395Fig. 7 :
Fig. 7: Streamlines in: (a) the ball mill obtained using the analytical velocity solution by Ding et al. [22], and (b) the AG mill obtained using Spline interpolation from the coarse-grained DEM velocity field.The streamlines are superimposed over the velocity vector data (shown in red colour) from DEM.The streamlines represented by dotted lines are the master streamlines, while the rest are interpolated using them.The number of streamlines in panels (a) and (b) are just for visualisation, while in the actual simulations the number of streamlines is much higher.

Fig. 8 :
Fig. 8: Schematic allocation of material points along streamlines for the flow within an example rotary mill.The material points are distributed by marching along the length of the streamlines in discrete timestep ∆t.The material points shown here are just for visualisation, while in the actual simulations they are distributed by adopting a very small time step.The granular mass within a material point is then further discretized along the micro-structural coordinate.Such a scheme is applied to both the ball and the AG mills (colour figure online).
cally add additional data points where needed and then generate a Spline function [25].These carefully selected streamlines are referred to as master streamlines for the bulk flow within the AG mill.By employing the Spline function, we can effectively interpolate any number of streamlines between any of these master streamlines, as illustrated in Fig. 7b.This technique allows us to obtain a set of closed streamlines that accurately describe the granular flow inside the AG mill.

534
velocity and acceleration, respectively, obtained 535 from the velocity field.The term1  2 a(t)∆t 2 ≈ 0 536 in Eq. 11 for very small accelerations, as long as 537 the time-step (∆t) chosen is very small.Therefore, 538 X(t + ∆t) = X(t) + v(t)∆t is also a good approx-539 imation to determine the position of a material 540 point in subsequent time-steps.Considering the 541 above, it is important to make a couple of observations with respect to the distribution of the material points along the streamlines: 1.For the ball mill case, near the vicinity of the zero-velocity line (α-line) the accelerations are actually significant since the velocities transform rapidly from v = 0 on the α-line to some non-zero value on the δ-line, given by the analytic velocity solution.Therefore, in Eq. 11 the term 1 2 a(t)∆t 2 ≠ 0, and does actually have a significant role.In fact, it is the acceleration that drives the material when it lies on the α-line, since v = 0 at that location.2. For the AG mill case, the streamlines are obtained using Spline interpolation.It was observed that the gradient of a streamline may not always be precisely oriented along the velocity vector obtained from the DEM simulation, as depicted in Fig. 7b, where the interpolated streamlines are superimposed over the velocity vectors obtained from DEM simulation.However, while distributing the material points along the streamlines, we neglect this second-order effect and assume that the velocity vectors are parallel to the streamlines.

Fig. 9 :
Fig.9: Volume of material points for the granular flow in AG mill calculated using the area obtained from Voronoi polygons[1].The length along the axis of the mill is assumed to be unity.The polygons are coloured randomly (colour figure online).

Fig. 10 :Fig. 11 :
Fig. 10: Cumulative distribution function (CDF) of product particle sizes along different streamlines after one mill revolution.(a) Ball mill, and (b) AG mill.The insets show the actual location of the corresponding streamlines in the flow region.An ultimate product size distribution for the entire minerals within the mills is also plotted.The number of cells used along the micro-structural coordinate for each simulation are M = 10 3 (colour figure online).

Fig. 12 :
Fig. 12: Evolution of mean particle size (d 50 ) normalized with respect to the feed mean particle size (d 0 50 ) for variation in values of impact velocity parameters a v , b v and w v .The base values of the parameters are listed in Table 1.Only one of the parameters is varied at a time keeping the other two parameters equal to their base values.[Top -ball mill] The horizontal axis represents the depth of the flow measured from the free surface and normalised with respect to the maximum depth of zero velocity line (α 0 ).[Bottom -AG mill] The horizontal axis corresponds to the number of streamlines denoted from 1 to N s as shown in the inset at the bottom figure e (colour figure online).

Fig. 10a and 10b 689 illustrate the cumulative distribution functions 690 (
CDFs) of particle size for different streamlines in 691 the flow region after one revolution of the ball mill 692 and the AG mill respectively.

Fig. 13 : 7 Conclusions
Figure12, we illustrate the variation of the mean particle size (d 50 ) with changes to these parameters for both the ball and the AG mills.It can be seen from the plots that for both the ball mill and the AG mill, the mean particle size varies

922 √ ( ε − εv 3 1 )
3) where x i , v i , m i are the position, velocity, and 915 mass of particle i, respectively, and x is the posi-916 tion at which the field is estimated.The velocity 917 field is then time-averaged to obtain ⟨v⟩, and the 918 spatial gradients are taken to obtain the strain-919 rate tensor ε = 1 2 (∇⟨v⟩ + ∇ T ⟨v⟩).The volumetric 920 strain rate is then calculated with εv = Tr ε, 921 and we define the shear strain rate as εs = ∶ ( ε − εv 3 1), with 1 the identity tensor. 923

924 of particle i from the field calculated at discrete 925 grid points. 926 Finally 1 0
, the static pressure P s = Tr σ is calculated from the contact forces f ij between particles i and j byσ = ∑ contacts f ij ⊗ (x j − x i ) ∫ w(|x i + s(x j − x i )|)ds,(B.5)with the integration performed numerically.927 Acknowledgments.The authors would like to 928 thank Molycop and The University of Sydney for 929 financially supporting this project.930 Declaration of interests.The authors did not 931 receive support from any organization for the sub-932 mitted work.No funding was received to assist 933 with the preparation of this manuscript.No fund-934 ing was received for conducting this study.No 935 funds, grants, or other support was received. 936 i,z , s i,z+1 } represents the parti-