Heterarchical modelling of comminution for rotary mills: Part II --- Particle crushing with segregation and mixing

.


Introduction
Granular materials often segregate by size or density when flowing.This segregation is typically not perfect, and some level of mixing is present counteracting this segregation.These two competing mechanisms create patterns in terms of the spatial distributions of size or density of particles.When particles also can crush, these patterns can become  In natural granular flows such as fault gauges [1,2], snow avalanches [3], debris flows [4,5], and rock avalanches [6,7], particle segregation is a key mechanism driving the dynamics of particle sizes.Segregation is also responsible for stratification in several rock formations [8].The process of segregation in granular flows has been extensively studied, and various theories have been developed to model the behaviour of mixtures of two [9][10][11][12] or more [13][14][15][16][17][18][19][20] types of particles.
These models describe segregation in simplified flow geometries, such as inclined planes [15,16,18], which is somewhat representative of granular avalanches, landslides, and debris flows where segregation occurs in the direction normal to the flow [21].Other models have focused on the fronts of these avalanches [22].Recently, models have been proposed for bounded heap flow, rotating tumblers, and hopper discharge [17,19,20].
In addition to segregation, another significant mechanism that occurs during granular flows is the mixing of particles.Indeed, complete segregation where particle species are fully separated from one another is rarely observed due to the random fluctuations of particles during the flow, which act as a de-segregation mechanism.In the past, the phenomenon of mixing has been studied analytically [12,23,24] and explored using experiments and numerical simulations [19,25,26].
In this work, we will focus on the modelling of the interior of autogenous grinding (AG) mills, where the particles undergo substantial deformations.AG mills are specifically designed to produce severe particle crushing, as first described in Part I of this contribution.In this Part II, we additionally explain the physics of the segregation and mixing mechanisms within the heterarchical stochastic formalism and study the effects of segregation and mixing on the system.
In Part I of this contribution, we delved into the dominant mechanism of kinetic crushing, which is prevalent in industrial mills.However, that paper highlighted the clear necessity to further treat the segregation and mixing mechanisms as part of the modelled process of comminution.
By solely relying on kinetic crushing, the extent of crushing within the system diminished over time, as individual particles do not explore the whole system, but instead follow closed streamlines through the mill and reach a state where they are cushioned against further crushing.In practice, however, this cushioning is disturbed by

Swapping of particles with frequency 𝒇 𝒎𝒊𝒙𝒊𝒏𝒈
Swapping of particles with frequency

Crushing Mixing Segregation
Fig. 1: The stochastic rules for crushing, segregation, and mixing mechanisms in the general heterarchical framework.The crushing rule is applied along the internal micro-structural coordinate while the rules for segregation and mixing are applied along the physical external coordinate(s).This figure shows the evolution of the initial state over time when the rules for crushing, segregation, and mixing are applied (colour figure online).
framework [18] outlines the rules for the segregation and mixing of particles in space, which are briefly discussed here.
As described in [18] and in Part I of this contribution, in the heterarchical framework a polydisperse granular system (i.e. one composed of many particle sizes) is replaced by an equivalent stochastic lattice.Each cell in the lattice represents particles of a specific size and the position of a cell in the lattice is defined by both external physical coordinates and an internal micro-structural coordinate.The rule for crushing is applied along the internal micro-structural coordinate and the advective rules of segregation and mixing are applied along the external physical coordinates.The rules for segregation and mixing are such that the particle size in a given cell of the lattice is swapped with the particle size of its neighbouring cells with a certain frequency (see Fig. 1).While the direction of the swapping of particle sizes is random for the mixing process, it is dependent on the neighbouring particle sizes for the segregation process.
In this section, we present the mathematical formulation for implementing the particle size dynamics of segregation and mixing in rotary mills using the heterarchical model.The formulation is presented for a polydisperse particle size distribution.

Particle size distribution 142
It is useful to represent a polydisperse granular medium using a continuum particle size coordinate s [15].The fraction of particle sizes in the range (s a < s < s b ) can then be calculated as where φ is the probability density function that characterises the particle size distribution, such that while s is the average particle size.When dealing with polydisperse media, it is also useful to define three different forms of density other than the density of the material (ρ m ).The intrinsic density ρ * (s) represents the mass of a unit volume of particle size s.On the other hand, the partial density ρ(s) signifies the contribution of the intrinsic density of a given particle size s towards the average (or bulk) density ρ b .The relationships between those densities are given by: where ν is the solid fraction.The governing equation for the conservation of mass [27] of a given class of particle sizes, in the reference frame of the barycentric motion (i.e.having already accounted for the bulk motion of the material) can be expressed as where ⃗ F (ρ) is the corresponding mass flux, which is generally given by ⃗ F (ρ) = ρû + D∇ρ. ( Notice that the mass flux consists of two com- 151 The segregation velocity of particles [9,11,15,17,19] has been shown to be a function of at least the shear strain rate, particle size ratio, concentration of the particle sizes, particle shape and density.Existing segregation models typically consider particle segregation driven by either the contrast in their size or density.Most of these models cannot capture the segregation behaviour when both the particle size and the particle density vary.As an alternative, [28] proposed a framework to explain the segregation behaviour in granular flows using the granular temperature and the kinetic pressure, which are functions of the measurable velocity fluctuation of the particles.The gradient of the granular temperature, or of the associated kinetic pressure, can be used to explain the direction and the rate of segregation in systems where both particle size and particle density vary.Building upon this concept [29] have expressed the segregation velocity as a function of the kinetic pressure where c 0 is a parameter controlling the rate of segregation, and s −1 is the local harmonic average particle size.This harmonic average particle size is a weighted average with respect to the partial density (ρ) of each particle size (s) and is given as The expression in Eq. 6 allows us to calculate the segregation velocity for any arbitrary particle size in a polydisperse granular system.For the current case of an AG mill (where only mineral particles exist), by assuming that all the particles have the same intrinsic density, the above expression can be simplified to Following the work of [30,31], the diffusivity (D) can be calculated as where l is a non-dimensional parameter controlling 152 the rate of diffusivity, εs is the shear strain rate, 153 I = εs s √ ρ m /P is the inertial number [32] previ-   For a given material point (i, j) and particle size (s), Eq. 4 can be written as where ρ i,j,s is the discretised array of partial densities.Upon integrating the above equation over the volume of a material point, we get Using the Gauss divergence theorem, the second term in the Eq.11 can be written to give where n is the surface unit normal.Assuming that the density within a given cell is homogeneous, we thus find We approximate the surface integral ( ∮ S ) by a summation (Σ) over the sides of the Voronoi polygons.Accordingly, Eq. 13 becomes where the index k denotes the edges of a Voronoi polygon with its neighbouring cells, l ijk is the length of the corresponding sides of the polygon, and V i,j is the volume of the polygon.
The frequency of swapping particle sizes s along the k th edge of the Voronoi polygon with its neighbouring cell is then defined as Notice that the mass flux ( ⃗ F i,j,s ) in the above Once the frequency of swapping is obtained, we need to determine the probability with which the particle sizes are swapped in a given time increment ∆t.The probability of swapping a particle size in a Voronoi volume along a neighbouring cell from one of its edges is obtained as where ⟨⟩ is the Macaulay function (defined as ⟨x⟩ = x for x ≥ 0, and ⟨x⟩ = 0 for x < 0).Note that the expression for swapping frequency in Eq. 15 also carries a direction that is embedded in its sign.We choose to only propagate material in the positive direction and use the Macaulay function to satisfy this requirement.This is justified as the frequencies that have a negative direction for the current cell, will be taken care of in the adjacent cells where they will have a positive direc-  that is given by Eq. 15.
Fig. 3 shows the initial and final states after the segregation is completed.After some time, as the particles begin to segregate, they gradually develop into a state where all the large particles are found at the top of the lattice, while the smaller ones reside at the bottom.
A spatiotemporal plot of the average particle size over time along the physical vertical direction is also shown in Fig. 4 for varying values of M .As the averaging is done over an increasing number of cells along the heterarchical coordinate, the spatiotemporal plot converges to the analytical solution given in [9,11], as also shown in [18,34].The concentration shocks as seen in the bottom right of Fig. 4 represent the analytic solution (the superimposed black lines) of the following partial differential equation

279
To validate the mixing mechanism, we take the

293
A spatiotemporal plot of the average particle size with time along the physical coordinate in the vertical direction is shown in Fig. 6.Similar to the case of segregation, as the averaging is done over an increasing number of cells along the heterarchical coordinate, the plot converges to the analytical solution given by [12] as also shown by [18], which is found by solving the diffusion equation

Unstructured lattice 294
After validating the heterarchical approach using corresponding results for the structured lattice and are in agreement with the analytical solution of the corresponding continuum equations for the advection and diffusion problems (Eq.17 and Eq.18), respectively.This validates our approach to model segregation and mixing using the heterarchical method.The next task is thus to adopt this method on the unstructured lattice used for rotary mills.

Rotary drums without crushing: lab-scale
The segregation and mixing models contain two free parameters, one controlling the rate of segregation (c o ) in Eq. 6 and one controlling the rate of diffusivity (l) in Eq. 9. To calibrate values for these models, a qualitative comparison is made here with experimental results for mixing and segregation in the literature [25,35].Fig. 7: Bidisperse segregation and mixing in an unstructured lattice.The cells are randomly distributed in the lattice and are filled with small (yellow) and large (blue) particles in equal proportion.A number of such lattices are considered along the heterarchical coordinate.The particle size is averaged along the heterarchical coordinate for each cell.For the segregation, we start with the initial condition of a uniform distribution of particle size in the lattice and for the mixing, we start with an initial condition of a perfectly segregated state.The rules for segregation and mixing are then applied (colour figure online).

336
The mixing of particles in rotary drums is sim-

368
To simulate the mixing mechanism, the frequency of swapping cells (here with the same particle size, though different colours) can be calculated using Eq. 15.To obtain the swapping probability, the critical time step can be calculated as done before using max k (f swap i,j,s,k ) ∆t c < 1/2 for each material point, and considering the minimum value.However, in doing so, there are very few material points in the system for which the critical time step is substantially small (about ≈ 10 4 order smaller than the rest of the points).Adopting such time step will cause the simulations to run very slowly.In order to overcome this, the critical time step for mixing (∆t m ) is taken as the volume average of the critical time steps of all the material points in the system and is given as where ∆t i,j c is the critical time step for the mate- The diffusivity can be calculated using Eq. 9.

374
The free parameter l in the diffusivity equation  whose exact form depends on the degree of filling and the particle size ratio.However, at higher rotational speeds, these streaks get suppressed, while instead the system quickly develops a steady state motion with radially segregated pattern, where the smaller particles tend to occupy the inner core of the bulk.This work focuses on the high rotational speeds associated with rotary mills.Therefore, here the heterarchical approach will be adopted, to see if it can recover the radially segregated pattern in binary systems of particle sizes.In so doing, our primary objective is to calibrate the parameter (c o ), which controls the rate of segregation (see Eq. 6).The formation of radial streak formations using heterarchical approach is beyond the scope of the present study and would require a different approach than integration along 416 steady state streamlines.

417
The work of [35] provides a particularly valu-418 able observation for the current study.Specifically, 419 that work illustrated the impact of the degree of 420 filling and particle size ratio on the number of revolutions required to achieve close to complete segregation of bidisperse particles in rotary drums.
The simulation parameters based on [37] using glass beads for segregation in a rotating drum are listed in Table .2. The critical simulation time is obtained by volume averaging.
Based on the degree of filling (f = 40%) and the particle size ratio (= 0.5) the number of revolutions required for achieving complete segregation is about ≈ 4 revolutions.Fig. 9 shows the evolu- tion of the segregated state in rotary drum.The parameter controlling the rate of segregation (c o ) as defined in Eq. 6 is varied until complete segregation is achieved in the aforementioned number of revolutions.The calibrated value c 0 = 0.05 s −1 is obtained from the simulations.
As done previously for the mixing of particles, It is important to note that in this simulation, no mass flux is explicitly modelled inside the mill.
The material is fed into the mill at the start of the simulation, and the crushed finer particles are removed only at the end of the simulation which in this case happens after 13 revolutions of the mill.The exit of fine material through the discharge end of the mill is facilitated by using a grated screen [42].This screen consists of a series of openings distributed over its surface, and covers 2-12% of the mill cross-sectional area.We consider one design of such grates from [42], which is shown in Fig. 10.The size of the grates on the screen controls the size of the fine particles filtering out of it.Once this material comes out through the grated screen, it is then passed into the discharge trunnion with the help of pulp lifters [43].Details regarding how we measure the product particle size distribution exiting through this grated screen will be discussed later in this section.
The simulation parameters for the AG mill are the same as those used for the simulations in part I of this contribution, as recapped in Table .3.
The segregation velocity of the particles is still calculated using Eq.6 where the adopted value of constant (c o ) is taken to be the same as calibrated in Sec 4. The segregation velocity is also dependent on the gradient of the kinetic pressure field.The kinetic pressure field for the AG mill is obtained from the discrete element method (DEM) simulation using coarse-graining, as shown in Fig. 11a and explained in Part I of this contribution.The gradient of the kinetic pressure can thus be obtained.For implementing the mixing mechanism, the diffusivity field (see Fig. 11b) is calculated using Eq. 9 where the value of the constant (l) is used as calibrated in Sec 4.
For the case of the current simulation, the particles residing on a material point of a given streamline not only undergo crushing but also advect and diffuse in space due to segregation and mixing events.This is illustrated for the two stochastic simulations shown in Fig. 12, where the trajectory of a tracer parcel of particles is plotted over time inside the mill.As done in Part I of this contribution, when only the crushing mechanism    particles, the fluctuating velocity of particles, and the kinetic pressure field.At this stage, we do not delve deep into these physics and rather adopt a simple probabilistic approach that allows us to estimate the particle size distribution of the fine material that exits through the grate.
Specifically, we consider a grated screen that has radially distributed holes covering the region from the periphery to a certain distance from the centre, which in this case is assumed to be a distance of two-thirds from the centre of the screen.
These apertures are assumed to be circular and have a diameter s g for simplicity.The effect of the shape of apertures, their distribution, and the distance from centre of the screen needs to be further investigated in the future and is beyond the present scope of this study.
Therefore, the available area (a) for a particle of size s ≤ s g to pass through this aperture is a = π(s g − s) 2 /4.The probability (p) of passing a particle of size s through the aperture of size s g is then, p = available area aperture area where ⟨⟩ is the Macaulay function.Note that

721
The probability density function P (s) for the above function can be obtained by normalising the function with the total area under the curve and can be expressed as

8
more pronounced and complex.The mechanisms 9 of particle segregation, mixing, and crushing have 10 many important implications for different natu-11 ral and industrial processes.The purpose of this 12 paper is to address the coupling between these 13 mechanisms within the context of particle milling. 14

154
ously defined in Part I of the contribution, and P 155 is the total pressure.

156 2 . 4
Stochastic swapping rules 157 As detailed in Part I of this two-part contribu-158 tion, the comminution problem within rotary mills 159 is tackled by integrating the physics dictating the 160 particle size dynamics along a set of closed set of 161 streamlines.This set is taken consistently from a 162 known velocity field, obtained either analytically if 163 such an equation is available or by coarse-graining 164 particle-based simulations.These streamlines are 165 then discretised into equal mass material points 166 along the length of the streamlines, consistent 167 with discrete time steps.The volume (assuming a 168 unit length into the page) corresponding to these 169 material points is obtained using a Voronoi tessel-170 lation [33] as shown in Fig. 2. The mass balance 171 equation (Eq.4) is solved for each of these material 172 points at each time step. 173

Fig. 2 :
Fig. 2: Tessellation of space and flux calculation of material points for the segregation and mixing mechanisms in the mill.Left: The space around each material point is discretised by a Voronoi tesselation.The areas are coloured randomly.Right: The flux of material moves in and out of a given cell to its neighbouring cells through the sides of the Voronoi cell.The highlighted cell in this figure has six neighbouring cells numbered from k = 1 to 6, where nijk is the unit normal vector on each side of the polygon and l ijk is the corresponding length of the side of the polygon (colour figure online).

178Fig. 3 :
Fig.3: Bidisperse segregation in a structured lattice.An regular lattice is randomly filled with small (yellow) and large (blue) particles in equal proportion.A number of such lattices are considered along the heterarchical coordinate as shown in the left.The particle size is averaged along the heterarchical coordinate for each cell.The rule for segregation is applied to this lattice and as segregation takes place we achieve a completely segregated state as shown in right.The animation showing the evolution of average particle size from initial uniform to a perfectly segregated state can be found in SI Video 1 (colour figure online).

229
scheme for unstructured grids to account for the 231

3. 1
Structured lattice 3.1.1Segregation We begin with a regular square lattice, with each cell randomly filled with either small particles (s a = 0.5 mm) or large particles (s b = 1 mm) in equal proportion.A number of such lattices are considered along the heterarchical coordinate, as shown in Fig. 3.The average particle size at a given cell location in physical space is obtained by averaging all of the particle sizes along the heterarchical coordinate, as defined by s i,j = 1 M ∑ M z=1 s i,j,z .We set an initial condition such that the average particle size at each cell location in the physical space is uniform and the mass of particles in each cell is equal.Another condition on the segregation rule is implemented for a simple case where the gradient of kinetic pressure in space is constant and is given as ∇p k = 1 ĵ N/m 3 , where ĵ is the unit vector along the vertical spatial coordinate j.The parameter c 0 in Eq. 6 is here taken as 1 s −1 for simplicity.The boundary condition is set such that the advective flux ( ⃗ F s = ρû) is zero for the top and bottom cells of the lattice.Neighbouring cells swap particle sizes with a frequency

Fig. 4 :
Fig. 4: Time evolution of average particle size (s) in a structured lattice for bidisperse segregation with varying the number of cells (M ) along the heterarchical coordinate.For each case, the gradient of kinetic pressure field ∇p k = 1 ĵ N/m 3 and the lattice consists of 100 × 100 cells along the two spatial coordinates.From (a) to (d) the number of cells along the heterarchical coordinate are varied as M = 1, 10, 10 2 and 10 3 .As the number of cells along the heterarchical coordinate increases the solution converges to the equivalent analytical solution of the 1D advection equation (Eq.17) (colour figure online).

280Fig. 5 :
Fig.5: Bidisperse mixing in a structured lattice.An initial uniform lattice is filled in equal proportion with small (yellow) and large (blue) particles.A number of such lattices are considered along the heterarchical coordinate as shown in the left.The particle size is averaged along the heterarchical coordinate for each cell.The rule for mixing is applied for this lattice and as the particles mix we achieve a complete uniform state from an initial perfectly segregated state as shown in right.The animation showing the evolution of average particle size from the initial perfectly segregated state to a uniform state can be found in SI Video 2 (colour figure online).

295 7 310Fig. 6 :
Fig. 6: Time evolution of average particle size (s) in a structured lattice for bidisperse mixing with varying number of cells along the heterarchical coordinate (M ).For all cases the diffusivity D = 0.01 m/s 2 and the lattice consists of 100 × 100 cells along the two spatial coordinates.From (a) to (d) the number of cells along the heterarchical coordinate are varied as M = 1, 10, 10 2 and 10 3 .As the number of cells along the heterarchical coordinate increases the solution converges to the equivalent analytical solution of the 1D diffusion equation.Bottom right shows the superimposed contour lines in gray obtained from the analytical solution of the 1D diffusion equation (Eq.18) (colour figure online).

337Fig. 8 :
Fig. 8: Spatiotemporal plot of average particle size for bidisperse mixture in an unstructured lattice.Left Segregation, Right Mixing.For segregation, the gradient of kinetic pressure field ∇p k = 1 ĵ N/m 3 and for mixing, diffusivity D = 0.01 m/s 2 .For both cases M = 10 3 .The superimposed lines show the corresponding analytical solution of the advection (Eq.17) and the diffusion equation (Eq.18) (colour figure online).

375
Fig.9i, where the particle size is averaged along 391

3924. 2 Fig. 9 :
Fig.9: Comparison of mixing and segregation for a lab scale rotating drum: (a-d) The evolution of mixing for a rotating drum in an experiment[25].For the experimental rotating drum, the drum is initially filled in two halves using identical particle sizes coloured white and black.The particles are initially placed horizontally and as the drum rotates, the particles get inclined at the angle of repose of the material.In all simulations M = 1 unless specified otherwise.(e-i) Evolution of mixing in a simulated rotating drum.The drum is initially filled in two halves using identical particle sizes coloured yellow and blue.The experiment starts horizontally whereas the simulated sample begins at the initial avalanche condition.(jn) Evolution of segregation in a simulated rotating drum.The drum is initially filled with small particles (yellow) and large particles (blue) in equal proportion and distributed uniformly in the drum.(i,n) The continuum limit of the heterarchical model is illustrated with M = 10 3 (colour figure online).

Fig. 10 :
Fig.10: Schematic of an Autogenous grinding mill and a typical grate design for the discharge of finely crushed material through the mill.Note that in the current 2-dimensional (2D) simulation for the AG mill, we use time (or mill revolution) as a proxy for the processing time of the mineral particles inside the mill.

Rate of diffusivity l 10 isFig. 12 :
Fig. 12: Trajectory of a tracer parcel of particles for the simulation of comminution in an AG mill.The segregation and mixing events cause the tracer to move from one streamline to another.Otherwise, in the absence of segregation and mixing the tracer traverses only along its original given streamline.Panels (a) and (b) show the trajectory for the same tracer parcel of particles for two simulations.The stochastic nature of the segregation and mixing mechanisms cause the tracer to follow different trajectories (colour figure online).

Fig. 13 : 10 !Fig. 16 :
Fig. 13: Cumulative distribution function (CDF) of the product particle sizes for the whole material inside a simulated AG mill (colour figure online).

713
the expression for probability in Eq. 20 is an 714 approximation where we consider the probabil-715 ity of passing finer particles through the grate 716 as a function of only particle size and the size 717 of the grate.The actual physics governing this 718 mechanism which might be additionally related 719 to the kinetic pressure or fluctuating velocities of 720 particles is not accounted for here.

2 .Fig. 17 :
Fig. 17: Cumulative distribution function (CDF) of the product particle size inside the mill, as well as outside the mill after screening through the grates.A typical grate size of diameter s g = 50 mm is considered here (colour figure online).

Table 1 :
Simulation parameters for mixing[25] 369 rial point and V i,j is the volume of the material 370 point.This choice is validated by reducing the 371 time step further and verifying that the results do 372 not practically vary.373

Table 3 :
Simulation parameters for an AG mill