A Novel Particle-Based Approach for Modeling a Wet Vertical Stirred Media Mill

Abstract

Modeling of wet stirred media mill processes is challenging since it requires the simultaneous modeling of the complex multiphysics in the interactions between grinding media, the moving internal agitator elements, and the grinding fluid. In the present study, a multiphysics model of an HIG5 pilot vertical stirred media mill with a nominal power of 7.5 kW is developed. The model is based on a particle-based coupled solver approach, where the grinding fluid is modeled with the particle finite element method (PFEM), the grinding media are modeled with the discrete element method (DEM), and the mill structure is modeled with the finite element method (FEM). The interactions between the different constituents are treated by loose (or weak) two-way couplings between the PFEM, DEM, and FEM models. Both water and a mineral slurry are used as grinding fluids, and they are modeled as Newtonian and non-Newtonian fluids, respectively. In the present work, a novel approach for transferring forces between grinding fluid and grinding media based on the Reynolds number is implemented. This force transfer is realized by specifying the drag coefficient as a function of the Reynolds number. The stirred media mill model is used to predict the mill power consumption, dynamics of both grinding fluid and grinding media, interparticle contacts of the grinding media, and the wear development on the mill structure. The numerical results obtained within the present study show good agreement with experimental measurements.

1. Introduction

In the mineral processing industry, comminution is the single most energy-intensive process and accounts for a large part of the capital and operating cost for mineral production [1,2]. It has been estimated that approximately 4% of the electric energy produced globally and around 50% of the energy at mine sites is consumed by comminution [3]. Tumbling mills are frequently used for comminution in the mineral processing industry. In a tumbling mill, much of the energy is absorbed in low-impact contacts that do not result in particle breakage [1]. The inefficient nature of the comminution process suggests that the mineral processing industry depends on improvement of the efficiency of this process. Compared to tumbling mills, stirred media mills are an attractive alternative since they have a significantly higher energy efficiency [4,5]. In a stirred media mill, the reduction of product particle size occurs inside a bed of grinding media. Comminution in a wet stirred media mill is obtained by pumping a slurry of product mixed with water through a grinding chamber. Inside the chamber, a shaft with attached agitator elements rotates, bringing the grinding media and slurry into motion. This motion induces many grinding media contacts, resulting in reduction of the product’s particle size. The interest in stirred media mills is reflected in a growing body of research, with recent contributions such as [6,7,8,9,10]. Despite the industrial importance, comminution in stirred media mills is a technology that is still poorly understood. To advance the knowledge and understanding of the operation of stirred media mills, modeling and simulation are powerful tools. A trustworthy computational model allows study of phenomena that are difficult or impossible to investigate experimentally. However, modeling of the dynamics of grinding media and slurry in a stirred media mill is a complex problem that requires robust and efficient numerical models.

The discrete element method (DEM) has been widely employed to model various granular materials, including grinding media in comminution. Early applications of the DEM to model dry stirred media mills, focusing on the dynamics of grinding media, are found in [11,12,13]. More recent works include [14], in which DEM was used for simulation of a lab scale dry ball mill, and [15], which combined DEM simulations and experiments to study liner design and liner wear in a semiautogenous (SAG) mill. A DEM study of the grinding media collision environment in a SAG mill under varying operating conditions is found in [16]. In [17,18], DEM is used for representation of realistic grinding media particle shapes in a SAG mill by using superquadratic discrete elements. Fukui et al. [19] presents models of horizontally and vertically oriented stirred mill designs, investigating the influence of agitator shaft direction on the grinding performance. Positron emission particle tracking (PEPT) is a technique that has been used with some success to experimentally study the dynamics of grinding media as an attempt to quantify the grinding performance [20,21]. In [22], a combination of PEPT and DEM was used to study the grinding media dynamics in a vertical attritor mill. Particle breakage in a ball mill is modeled using DEM and a bonded-cell method in [23]. In [24], a DEM model is used to study a batch vertical stirred mill. In [25], the computational aspects of DEM modeling of crushing and comminution of granular materials are addressed.

Modeling of wet comminution introduces the challenge of representing the slurry, namely the mix of product and water. The product particle size is typically small, and thus it is not feasible to model individual product particles with DEM. Instead, the slurry is commonly represented as a continuous fluid. For wet comminution, the smoothed particle hydrodynamics (SPH) method, developed by [26,27], is frequently used to model the slurry. Coupling numerical methods, such as the finite element method (FEM), DEM, and SPH are candidates for modeling of a wet comminution circuit. Jonsén et al. [28] employed a coupled SPH–FEM model for a tumbling mill process, which was later extended to include both slurry and grinding media using a coupled SPH–DEM–FEM approach [29,30]. In [31], a coupled DEM–SPH was applied on an industrial scale stirred media detritor. Computational fluid dynamics (CFD) is an alternative to SPH and coupled CFD–DEM models of stirred media mills can be found in [32,33,34]. Eulerian CFD simulations of solid–liquid flow in a horizontal stirred media mill are presented in [35].

Early applications of the particle finite element method (PFEM) combined with DEM for modeling particle-laden flows can be found in [36,37]. More recently, a PFEM–DEM approach was used by [38] to model particle-laden flows with free surface. Coupled particle finite element method (PFEM), FEM, and DEM models were used in a recent work by [39] for a tumbling mill and in [40] for a wet stirred media mill.

The present work presents numerical simulations of the interactions between grinding media, slurry, and mill structure in a pilot wet stirred media mill (Outotec HIG5) with 7.5 kW of installed power. The mill structure is modeled with FEM, grinding media with DEM, and slurry with PFEM. The PFEM employed in this work is based on an arbitrary Lagrangian–Eulerian (ALE) approach and a robust and automatic volume mesh generation based on Delaunay triangulation, originally developed by [41]. The simulations are realized with a two-way loosely (or weakly) coupled transient DEM–PFEM–FEM model. A novel approach for transferring forces between grinding fluid and grinding media is implemented. The force transfer is realized by specifying the drag coefficient as a function of the Reynolds number. The aim of the present work is to establish a numerical tool able to provide information on grinding media and slurry dynamics that are difficult or impossible to obtain experimentally. Furthermore, the aim is also to investigate wear development on shaft and agitators in the HIG5 mill. The results of the motion of grinding media and slurry, wear development, power draw, and grinding media contact energy spectra are presented. The proposed numerical approach is validated against experimental measurements of power draw during operation of the HIG5 mill.

2. Materials

In the present study, the grinding media are beads of a composite alumina and zirconia ceramic material. The composition is 75% Al₂O₃ and 25% ZrO₂ and SiO₂. The bead shape is spherical, and the particle size is in the range 2.0–4.0 mm. The solid and bulk densities are 3.9 g/cm³ and 2.4 g/cm², respectively, and the Vickers hardness is 1150 HV1. Two different grinding fluids are investigated, water and a mineral slurry. The slurry has a concentration of 46% solids by weight, a density of 1.48 g/cm³. P₈₀ is 81 μm, and 46% of the feed material has a particle size smaller than 20 μm. P₈₀ is the screen size through which 80% of the particles will pass.

3. Modeling and Simulation

In this section, the modeling approach of the present study is briefly reviewed. A more detailed presentation can be found in [40]. The grinding media is modeled with DEM, the grinding fluid is modeled with PFEM, and the mill structure is modeled with FEM. The grinding media is modeled without considering breakage or attrition, which is justified by the short duration of the simulations. The grinding fluid is represented as a continuous fluid governed by the Navier–Stokes equations. For the interaction between grinding media, fluid, and mill structure, a coupled DEM–PFEM–FEM model is used. The coupled model is implemented and solved in the nonlinear multiphysics code LS-DYNA, version R11.1 [42], using 32–64 cores with 2.6 GHz CPUs with $2\times 128$ GB local memory. The computational time was between 35 and 42 h, depending on the choice of fluid and agitator rotational speed.

3.1. Discrete Element Method

Originally formulated by [43], the DEM is based on a representation of the motion of individual granular particles using Newton’s second law of motion. The translation and rotation of a particle i with mass $m_i$ and inertia ${\mathit{I}}_i$ is modeled by

$$m_i\frac{d{\mathit{v}}_i}{dt}=\sum _j{\mathit{F}}_{ij}^c+\sum _k{\mathit{F}}_{ik}^{nc}+{\mathit{F}}_i^g$$

(1)

and

$${\mathit{I}}_i\frac{d{\mathbf{\omega }}_i}{dt}=\sum _j{\mathit{M}}_{ij},$$

(2)

where ${\mathit{v}}_i$ and ${\mathbf{\omega }}_i$ are the translational and angular velocities of the particle. Contact forces and torque acting on particle i by particle j are given by ${\mathit{F}}_{ij}^c$ and ${\mathit{M}}_{ij}$. Noncontact forces (e.g., capillary forces) acting on particle i by particle k are given by ${\mathit{F}}_{ik}^{nc}$ and the gravitational force is given by ${\mathit{F}}_i^g$. The evolution of particle velocities, locations, and contact forces is obtained by explicit time integration of Equations (1) and (2). Interparticle contact forces are modeled using a linear spring and damper approach. The contact force is given by $f^c=\delta k$, where $\delta$ is the overlap and k is a linear spring stiffness. Contact energy dissipation is governed by a damping model where the damping force is given by $f^d=vc$, v is the relative velocity of the contacting particles, and c is a damping coefficient. Furthermore, the interparticle friction is governed by a sliding friction parameter $\mu$. In the DEM, stability of the explicit time integration scheme is met by selecting the time step size as a fraction of the critical time step. In the present study, the critical time step is given by $\Delta t_c=c_t\sqrt{m/k}$, where $c_t$ is a user-defined constant. If differently size DEM particles are used, the critical time step is the minimum value calculated over all i particles.

3.2. Particle Finite Element Method

The particle finite element method (PFEM) was originally developed for free-surface fluid flow problems [44] and fluid–structure interaction [45]. Since its original development, the PFEM has been applied to a wide range of physical problems where large deformations are present. The main idea of the PFEM is a combination of a Lagrangian FEM with a powerful remeshing strategy. The computational domain is defined via a set of points/particles coinciding with the mesh nodes. The motion of the particles is treated in a Lagrangian manner. The particles contain both nodal variables, such as displacements or velocity, and physical properties, such as density and viscosity. The interactions between the particles are calculated using a finite element mesh.

3.2.1. Balance Equations in an ALE Framework

The continuity equation is given by

$$\frac{\partial \rho }{\partial t}+\frac{\partial }{\partial x_i}\left(\rho u_i\right)=0,$$

(3)

where $\rho$ is the fluid density, t is time, $u_i$ are the velocity components, and the index i refer to the space coordinates $x_i$. The flow is approximated as incompressible where $\rho$ is not a function of time or space. Thus, $\partial \rho /\partial t\cong 0$ and Equation (3) is reduced to a volume continuity equation

$$\frac{\partial u_i}{\partial x_i}=0.$$

(4)

The Navier–Stokes equations can be written as

$$\rho \left(\frac{\partial u_i}{\partial t}+u_j\frac{\partial u_i}{\partial x_j}\right)=\frac{\partial {\sigma }_{ij}}{\partial x_j}+\rho f_i,$$

(5)

where ${\sigma }_{ij}$ is the total stress tensor and $f_i$ are the components of external volume forces (including the external forces exerted by DEM particles on the fluid). The total stress tensor is given by

$${\sigma }_{ij}=-p{\delta }_{ij}+\mu \left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}-\frac{2}{3}\frac{\partial u_l}{\partial x_l}{\delta }_{ij}\right),$$

(6)

where p is the hydrostatic pressure, $\mu$ is the fluid dynamic viscosity, and ${\delta }_{ij}$ is the Kronecker delta. In a nearly incompressible flow, the last term inside the parenthesis in Equation (6) can be neglected and Equation (6) can be written as

$${\sigma }_{ij}\approx -p{\delta }_{ij}+\mu \left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right).$$

(7)

Thus, the viscous term in Equation (5) can be simplified, resulting in the following Eulerian form of the governing system of equations

$$\rho \left(\frac{\partial u_i}{\partial t}+u_j\frac{\partial u_i}{\partial x_j}\right)=-\frac{\partial p}{\partial x_i}+\mu \frac{{\partial }^2{u}_i}{\partial x_j\partial x_j}+\rho f_i,$$

(8)

$$\frac{\partial u_i}{\partial x_i}=0.$$

(9)

In the present study, an arbitrary Lagrangian–Eulerian (ALE) formulation is used for the fluid domain. The ALE description of motion requires some reformulation of Equation (8), resulting in

$$\rho \left(\frac{\partial u_i}{\partial t}+(u_j-v_j)\frac{\partial u_i}{\partial x_j}\right)=-\frac{\partial p}{\partial x_i}+\mu \frac{{\partial }^2{u}_i}{\partial x_j\partial x_j}+\rho f_i,$$

(10)

where $v_j$ is the velocity of the moving frame of reference [46]. The objective of using an ALE formulation is that the PFEM particles/nodes can move to a different position to maximize the mesh quality.

3.2.2. Meshing Procedure in the Particle Finite Element Method

The essence of the PFEM is the powerful mesh regeneration algorithm. The motion of the particles will eventually cause the mesh to distort and degenerate; this mesh is deleted, and a new mesh is constructed using the same set of nodes. For the mesh regeneration, a Delaunay triangulation algorithm is applied, and an alpha shape method is used to identify the internal and external boundaries. The general solution scheme of the PFEM is summarized below.

• The computational domain is filled with a set of points/particles.
• The particles are used as nodes to generate a finite element mesh using Delaunay triangulation [47,48].
• Internal and external boundaries are identified using an alpha shape scheme [49].
• The FEM is used to solve the governing equations on the mesh.
• Nodal positions are updated.
• Return to step 4. If a remesh is required, return to step 2.

The mesh regeneration is a key ingredient in the PFEM and requires a fast and robust algorithm. The PFEM uses a Delaunay triangulation approach for the mesh regeneration, and it is important to note that in this step, it is only the nodal connectivities that are updated; the nodes from the previous mesh are kept at the same location. For an extended overview of the theory and applications of the PFEM, the reader is referred to [50].

3.2.3. Time Integration

In the present study, the time integration of the Navier–Stokes equations is performed using the fractional step method, originally outlined independently by [51,52]. In the fractional step method, pressure and velocity are uncoupled, resulting in four linear systems of equations, three for the momentum equations (Equation (10)), and one for the continuity equation (Equation (9)). The fractional step method computes the velocity in three main steps. In the first step, a predictor velocity that does not satisfy the incompressibility condition is computed. In the second step, the predictor velocity is projected onto a space of divergence free vector fields, resulting in a Poisson equation of pressure. The computed pressure is then used to correct the velocity, resulting in a divergence free velocity. In the third step, the corrected velocity is used to move the particles to a new position. Convergence is obtained when the final position of the particles is stationary, within some error margin. More details of the fractional step method applied to incompressible flow problems can be found in [53,54,55].

3.2.4. Spatial Discretization by the Finite Element Method

A simple four-noded linear tetrahedral finite element is used to discretize the velocity and the pressure fields. The Navier–Stokes equations in an ALE framework can suffer from two numerical instabilities. The first is associated with the incompressibility (the tetrahedral does not fulfill the inf-suf condition) and the second one with the dominant convection (usually at high Reynolds numbers). The inf-suf condition is stabilized using the finite calculus presented in [46,56,57] and the convection terms are stabilized using orthogonal subscale stabilization [46,58]. The mesh distortion is minimized through a nonphysical motion of the mesh nodes using the Laplacian operation [46,59], and in the case that large boundary deformation takes place, the alpha shape and Delaunay triangulation are used as was explained in Section 3.2.2.

3.3. Grinding Media, Slurry, and Mill Interaction

A partitioned (or staggered) approach is used for fluid structure interaction (FSI), in which the fluid and solid equations are uncoupled. The uncoupling allows using specifically written codes on the different domains, which is beneficial in terms of efficiency. FSI is realized by a loosely (or weakly) coupled scheme, in which only one solution of either field is required each time step. This is advantageous in terms of computational efficiency since it avoids the iterative step required for convergence. For additional details on the FSI formulation and its implementation in LS-DYNA, the reader is referred to [46].

The fluid and grinding media interaction is modeled by a two-way coupling between PFEM and DEM models. The grinding media affects the fluid by adding mass and velocity of the particles to the volume forces in Navier–Stokes equations. The fluid flow over the particles result in a drag force computed by a potential flow around a sphere

$$f_d=\frac{v_f^{2}A{\rho }_f{C}_d}{2},$$

(11)

where $v_f$ is the fluid velocity, $A=\pi r^2$ is the projected area of a sphere with radius r, ${\rho }_f$ is the fluid density, and $C_d$ is a drag coefficient. The drag coefficient is a dimensionless quantity that depends on the shape of the object and the flow conditions around it. The drag coefficient is typically dependent on the Reynolds number. In the present work, a relationship between the drag coefficient and the Reynolds number based on the work by [60] is used. The drag coefficient is related to the Reynolds number by the following expression

$$C_d=\frac{24}{\mathrm{Re}}+\frac{2.6\left(\frac{\mathrm{Re}}{5.0}\right)}{1+{\left(\frac{\mathrm{Re}}{5.0}\right)}^{1.52}}+\frac{0.411{\left(\frac{\mathrm{Re}}{2.63\times {10}^5}\right)}^{-7.94}}{1+{\left(\frac{\mathrm{Re}}{2.63\times {10}^5}\right)}^{-8.00}}+\frac{0.25\left(\frac{\mathrm{Re}}{{10}^6}\right)}{1+\left(\frac{\mathrm{Re}}{{10}^6}\right)}$$

(12)

where Re is the Reynolds number. A plot of the drag coefficient as a function of Reynolds number is shown in Figure 1.

The mineral slurry is modeled as a shear-thinning (pseudoplastic) fluid. A viscosity model originally proposed by [61] is used for the relationship between dynamic viscosity $\mu$ and shear rate $\dot{\gamma }$

$$\mu ={\mu }_{\infty }+\frac{{\mu }_0-{\mu }_{\infty }}{1+{\left(\lambda \dot{\gamma }\right)}^n}$$

(13)

where ${\mu }_0$ is the zero shear-rate viscosity, ${\mu }_{\infty }$ is the infinite shear-rate viscosity, n is a dimensionless constant governing the deviation from a Newtonian fluid, and $\lambda$ is a time constant.

The parameters for the viscosity model are obtained from a measured viscosity and shear-rate relationship (Figure 2). The viscosity and shear-rate relationship is measured at room temperature with an Anton Paar Rheolab QC rotational rheometer using a CC27 bob cup setup. A linear least squares method is used to fit the viscosity model to the experimental points, as shown in Figure 2. The model parameters obtained from the fit are given in

Table 1. Fluid model parameters used for water and dolomite slurry.

Fluid	Viscosity Model	$\mathit{\rho }$ (g/cm3)	$\mathit{\mu }$ (Pa·s)	${\mathit{\mu }}_0$ (Pa·s)	${\mathit{\mu }}_{\mathit{\infty }}$ (Pa·s)	$\mathit{\lambda }$ (s)	n (-)
Water	Newtonian	1.00	$1.052\times {10}^{-3}$	-	-	-	-
Slurry	Non-Newtonian *	1.48	-	1.498	0.006	0.579	1.331

* The non-Newtonian viscosity model is a model for shear-thinning fluids formulated by [61].

.

3.4. Simulation Procedures

In the present study, the stirred media mill is a 7.5 kW Outotec HIG5. The simulations are performed on a geometrically simplified version of the mill (see the CAD model in Figure 3). The specifications of the HIG5 mill are presented in

Table 2. Outotec HIG5 mill specifications.

Property	Value
Mill volume (10${}^{-3}$ m3)	6.2
No. of discs	9
Installed power (kW)	7.5

. The volume enclosed by the mill outer structure is 6.2 × 10${}^{-3}$ m³. The centrally located shaft is fitted with eight identical agitator discs and one bottom disc with the same diameter but of slightly different design. The shaft is connected to an electric motor. The simulations are run for a total of 5.0 s and the rotational speed is ramped up from 0 to the final rotational speed during the first 0.5 s. In the present study, the final rotational speed is varied between 300–600 rpm. The mill structure is discretized using triangular FE shell elements and the casing and agitator materials are modeled with a rigid material description.

The amount of grinding media 8.93 and 7.66 kg for the water and mineral slurry simulations, respectively, correspond to fill levels of approximately 60% and 50% of the grinding chamber volume. The mill is filled with approximately 270,000 discrete elements for a grinding media mass of 8.93 kg and 120,000 discrete elements for a grinding media mass of 7.66 kg. The diameter of the discrete elements is between 2.0–3.0 mm, following a Gaussian distribution for the water simulations. For the mineral slurry simulation, 50% of the grinding media has a diameter between 2.0–3.0 mm and 50% between 3.0–4.0 mm. The sliding friction coefficients are obtained experimentally, and since the grinding media is spherical, a low value of the rolling friction coefficient is assumed (see

Table 3. Physical properties for the discrete element method (DEM) model of the grinding media common for all simulations.

Property	Value
Particle density (g/cm3)	3.9
Young’s modulus (MPa)	300
Poisson’s ratio	0.21
DE–DE frict. coeff., sliding	0.11
DE–DE frict. coeff., rolling	0.01
DE–FE frict. coeff., sliding	0.5
DE–FE frict. coeff., rolling	0.01
Damping coefficient	0.5

). The sliding friction coefficient for DE–DE contact is determined using a CSM Instruments micro scratch tester in combination with a CSM Instruments micro indenter. The DE–FE sliding fiction coefficient is determined with a UMT-2 tribometer using a ball on disc configuration. Young’s modulus of grinding media is approximately 300 GPa. Since the DEM in the present study uses explicit time integration, the critical time step size is proportional to the Young’s modulus of the grinding media. In DEM simulations, an increase of the critical time step size is usually realized by decreasing the contact stiffness between particles, which is proportional to the Young’s modulus. To increase the critical time step size, a value of Young’s modulus three orders of magnitude smaller than the real value is used in the present study. It has been shown that values of Young’s modulus in the range 10${}^7$–10${}^{11}$ Pa has a negligible effect on the granular material bulk flow properties [62,63].

PFEM in the ICFD module in LS-DYNA is used to model the grinding fluids. When the grinding fluid is water, it is modeled as a Newtonian fluid with a constant dynamic viscosity of 1.052 × 10${}^{-3}$ Pa·s and density of 1.0 g/cm³. The mineral slurry density is 1.48 g/cm³. The slurry is modeled as a non-Newtonian fluid with a dynamic viscosity dependent on the shear rate. Fluid boundary conditions are applied according to the settings used in the HIG5 mill. The inlet flow rate is 0.24 m³/h for both water and mineral slurry. A prescribed pressure equal to zero is used for the outlet. Both inlet and outlet are circular. A no-slip condition is applied between mill structure and fluid.

4. Results and Discussion

A number of simulations were set up to investigate the ability of the coupled PFEM–DEM–FEM model to reproduce the operation of the HIG5 stirred media mill. In the subsequent simulations, the dynamics of the fluid and grinding media and the performance of the stirred media mill were investigated for different operational conditions.

4.1. Steady Flow Past a Static Sphere

To validate the implementation of Equation (12), the flow over a stationary sphere located within a rectangular channel was simulated. The $C_d$ values from the simulations were obtained from Equation (11). As shown in Figure 4, the simulated results agree well with results predicted by Equation (12).

4.2. Power Consumption

During the operation of the mill, a driving torque is required to maintain the rotational velocity of the shaft. To keep the rotational velocity constant, the torque will vary due to the loads on the shaft and agitator from the interactions between grinding media, fluid, and the agitator elements. The grinding mill power output can be calculated as the product of the driving torque and the rotational velocity of the shaft. The power consumption can be studied to give an estimate on the performance of the grinding mill under varying operational conditions as well as for different fluids and grinding media. It is important to note that the power draw obtained from the numerical simulations of the present study is a theoretical (pure) power draw for just rotating the agitator. Thus, no mechanical or electrical losses are accounted for, nor any product breakage.

In Figure 5, the simulated theoretical (pure) power draw is presented for three rotational velocities with water as grinding fluid. Experimentally measured power draw for the 600 rpm case is shown in the same figure. A clear correlation between rotational velocity and power draw is observed; increasing rotational velocity results in a significantly increased power draw. Simulated and experimentally measured power draw agrees well for the 600 rpm case. The pure power draw for the mineral slurry was simulated using the non-Newtonian fluid model; the result for a rotational speed of 550 rpm is shown in Figure 6 together with the experimentally measured power draw. Simulated and experimental power draw shows excellent agreement. When comparing Figure 5 and Figure 6, it is observed that a the power draw is slightly underpredicted for the water case and somewhat overpredicted for the mineral slurry case. This might be partly explained by the fact that the rheology of the slurry was determined experimentally at room temperature. The temperature of the grinding fluid will likely increase during the operation of the mill. However, this was not investigated experimentally, and in the present study, no temperature dependency was included in the rheology model.

The power draw is an important aspect in wet comminution, and the proposed modeling approach provides opportunity to determine the power draw for different grinding fluids and operational conditions. Thus, it constitutes a valuable tool for stirred media mill operators. It is important to note that the simulated power draw in the present study is a relative calculation that only considers the resistance of the grinding media and the fluid.

4.3. Grinding Media and Fluid Dynamics and Wear Prediction

With the coupled model developed in the present study, detailed information on the fluid and grinding media dynamics can be obtained for varying mill operating conditions. In Figure 7, the grinding media distribution and velocity are shown for the case where water is used as grinding fluid and the rotational speed is 300 rpm. The results are shown at steady state. From Figure 7, it is observed that the highest grinding media velocities are located close to the perimeter of the agitator discs and that the velocity decreases rapidly with increasing radial distance from the agitator discs. This result indicates that most of the feed material particle breakage occurs between the discs. In Figure 8, the fluid velocity is illustrated for a case where the fluid is water and the rotational speed is 300 rpm. The result is presented as multiple semitransparent isosurfaces of fluid flow velocity. The fluid has its highest velocity close to the agitator discs and the velocity decreases with the radial distance from the discs, similar to the result for the grinding media. To further visualize the grinding fluid motion, a line integral convolution (LIC) was generated for the case where the fluid is water and the rotational speed is 300 rpm (see Figure 9). The LIC technique convolves noise with a vector field, producing streaking patterns that follow vector field tangents. Thus, in a LIC visualization, the flow of the velocity vector field is visualized showing the direction of the velocity field vectors. In Figure 9, the strength of the field (velocity magnitude) is shown by color where blue represents stationary fluid and red represents fast-moving fluid. From the LIC, some recurring patterns are observed, such as the formation of vortices above and below the agitator elements and irregular turbulent flow close to the agitators.

For the case where water is used as fluid, the abrasive wear distribution was calculated on the mill casing, agitator discs, and shaft. The results are shown in Figure 10. The wear was predicted by an implementation of Archard’s wear model [64], in which the wear is proportional to the normal contact pressure and sliding velocity. In Figure 10, it is observed that the wear increases close to the bottom of the HIG5 mill. This effect is due to gravity since the highest contact pressures occur at the bottom of the mill. Furthermore, it is observed that the outer perimeter of the agitator discs is heavily worn. This wear pattern is in line with observations from industrial use of the HIG5 mill.

Using the proposed coupled model, it is thus possible to predict the complex multiphysics of the stirred media mill in terms of grinding media and fluid dynamics for a wide range of operating conditions and materials. The inclusion of an efficient wear prediction model makes the model an efficient tool useful for design and optimization of stirred media mills.

4.4. Contact Energy Spectra

The contact energy dissipation associated with the collisions between grinding media was calculated for the HIG5 stirred media mill model. Interparticle collisions are associated with energy dissipation, and to quantify these collisions during operation of the HIG5 mill, a contact energy spectrum was calculated. The spectrum gives the frequency distribution of the normal component of the contact energies (see Figure 11). To obtain this contact energy spectrum, a methodology outlined by [34] was used. The normal contact energy distribution in Figure 11 is shown for all cases simulated in the present study. From this, it is observed that increasing the rotational speed results in more energy intensive interparticle collisions. The model with the mineral slurry as grinding fluid gives a completely different contact energy spectrum, with a significant shift to more energy-intensive contacts. It should be noted that the rotational speed for the mineral slurry case was 550 rpm and the grinding media size distribution is different from that of when water was used as grinding fluid. In Figure 11, a lower limit of contact energy of ${10}^{-15}$ J was used, as suggested by [34]. Since feed material particle breakage occurs mainly due to grinding media interparticle collisions, the contact energy distribution is an interesting parameter to study, and the efficiency of a stirred media mill can be improved by shifting the contact energy spectra to more useful contact events. Modeling of breakage of feed material particles was not included in the present study.

5. Conclusions

In the present study, a particle-based modeling approach based on fundamental, measurable physical relations is used to model the interactions between fluid, grinding media and mill structure in a HIG5 pilot stirred media mill. The modeling and simulations are based on a loose (or weak) two-way coupling PFEM–DEM–FEM approach where PFEM is used for the grinding fluid, DEM for the grinding media, and FEM for the mill structure. Furthermore, a novel coupling between the grinding media and fluid models based on the Reynolds number was utilized. With the proposed modeling strategy, the dynamics of fluid and grinding media, power draw, wear distribution, and interparticle collisions were investigated. Two grinding fluids were considered, water and a mineral slurry. Water was modeled as a Newtonian fluid and the mineral slurry was modeled as a non-Newtonian fluid with properties obtained from rheology measurements. With the proposed modeling approach, the complex multiphysics of the HIG5 mill could be accurately predicted and the model was validated against experimentally measured power consumption for both water and mineral slurry. The model predicts high velocity of the grinding media and fluid close to the agitator discs as is expected from experimental observations. Furthermore, the wear prediction agrees with experimental observations of wear patterns after prolonged use of the HIG5 mill. The model developed in the present study is a powerful tool for increasing the knowledge and understanding of the operation of stirred media mills. Thus, it should be valuable for industrial users of stirred media mills, allowing inexpensive study of how different materials and process parameters affect the grinding efficiency.