Digital Twins with Distributed Particle Simulation for Mine-to-Mill Material Tracking

Abstract

Systems for transport and processing of granular media are challenging to analyse, operate and optimise. In the mining and mineral processing industries, these systems are chains of processes with a complex interplay among the equipment, control and processed material. The material properties have natural variations that are usually only known at certain locations. Therefore, we explored a material-oriented approach to digital twins with a particle representation of the granular media. In digital form, the material is treated as pseudo-particles, each representing a large collection of real particles of various sizes, shapes and mineral properties. Movements and changes in the state of the material are determined by the combined data from control systems, sensors, vehicle telematics and simulation models at locations where no real sensors could see. The particle-based representation enables material tracking along the chain of processes. Each digital particle can act as a carrier of observational data generated by the equipment as it interacts with the real material. This make it possible to better learn the material properties from process observations and to predict the effect on downstream processes. We tested the technique on a mining simulator and demonstrated the analysis that can be performed using data from cross-system material tracking.

1. Introduction

It is estimated that mining accounts for 6% of the global energy consumption. Grinding is the most energy-demanding process, accounting for 32% of the mine’s energy consumption, followed by haulage (24%), underground mine ventilation (9%) and digging (8%) [1]. The grinding efficiency is many times as low as 1% and depends strongly on the particle size distribution and ore hardness [2]. Therefore, there is a strong motivation for determining and controlling the properties of the ore fed into the mill. Mine-to-mill refers to a holistic approach to mine optimisation; see [3] for a comprehensive description. Instead of optimising each process individually, one takes the ore body characteristic into account and considers the throughput and energy consumption of the full chain of processes when deciding on target fragmentation and controlling the feed to the mill. The mine-to-mill approach has, however, proven difficult to implement in practice. One reason is that it requires an elaborate digital infrastructure for collecting, processing and communicating data between edge devices, distributed control systems and centralized data services. The technology, popularized as digital twins [4] or the Industrial Internet of Things (IIoT) [5], is currently under intensive development. A second challenge in implementing the mine-to-mill approach is the unknown variations in the properties of the ore body. The variations propagate through the process chain in a complicated way and add to the full-system dynamics.

A mine-to-mill process chain is illustrated in Figure 1. The performance of each individual process depends on its design and control and on the properties of the material feed. Each process also alters the state of the material. The complex interdependency among the process dynamics and variations in the material properties makes it hard to determine what is cause and what is effect. The mineral concentration, hardness, strength and toughness of the ore body is partially known in advance, from exploratory drilling, and stored in a block model [6]. These are input parameters to many models for blasting [7], crushing [8] or grinding [9], such as the Kuz–Ram blasting and JKMRC grinding models. Measurement-while-drilling (MWD) of blast holes produces data that contain information about the geo-mechanical properties at higher spatial resolution than the exploration data. This may be used for improving the blast plans and for optimizing the subsequent processes [10,11,12]. Drilling and blasting are conducted in cycles, each cycle producing a volume of fragmented rock that is loaded using bucket excavators or load-haul-dump trucks before the next blasting. The particle size and shape distribution are subsequently transformed by the crushers and grinders before entering the concentrator. Variations in the diggability of the blasted rock [13,14,15] are useful for improving blast plans or for predicting the throughput at the crusher and mill. Transport and storage systems alter the location of the material, but may also induce (unwanted) mixing and segregation that propagate to subsequent processes and affect the dynamics there. Mine-to-mill analysis and optimization require knowledge of these changes, but on-line measurements of granular media properties are inherently difficult. The material movements during vehicle and belt conveyor transport are easily tracked, and these systems may be instrumented with sensors for monitoring some of the material properties. However, in chutes, stockpiles and silos, little is known about the material beyond the knowledge of what is fed into them. If any sensors can be installed at all, the observations are either indirect or limited to the surface of the granular media, which constitutes a small and perhaps not representative fraction of the bulk. Conventional tracking systems, which rely on combinations of plug-flow models (or first-in/first-out) and perfect mixing models that are calibrated using tracer experiments [16], fall short on transport and storage systems with complex surface and discharge flow. This creates blind spots in the tracking of material movements and knowing the properties at any given location.

To remedy this, we explored digital twins with distributed particle simulation, of systems doing transport and processing of granular media. The idea was to represent the material and its movements in a structured data format based on a particle representation. In addition to the material’s identity and current position, the particle data structure supports the reading and writing of observations from sensors and equipment along the chain of processes. When connected equipment performs unit operations on the material, the digital copy is updated accordingly. When the material reaches blind spots in the system, for example a silo or a stockpile, the digital copy is driven by a simulation model fed with data from the control system and available sensors in real time. For the sake of computational performance and memory, the digital twin uses pseudo-particles that represent a large collection of real mineral particles with a distribution of grain size and shape. For fast simulation, we used the nonsmooth discrete element method, which allows large time-step integration [17]. To achieve real-time performance in large processes, a data-driven model order reduction technique [18] was employed. The combination of the two models in an executable application, which can be part of a distributed system like a digital twin, we call a granular surrogate.

We tested the feasibility of doing material tracking with this technique on an idealized mine simulator with both discrete and continuous transport processes from two different sources. The tests included the analysis of how crusher power draw and milling performance depend on variations in ore hardness and fragmentation at the blast location. From crusher time-series, it was possible to identify the relation between the ore hardness and measurement-while-drilling signal despite having multiple ore sources being mixed. The sensitivity to tracking errors was analysed by perturbing the pseudo-particle origins. The response in mill performance to a change in ore hardness at the source was tested with different states of an intermediate stockpile. The effect was significant, but it was still possible to reconcile the ore’s grindability, in terms of the bond ball mill index, to the locations of origin in the ore body thanks to the particle-based representation.

Related work includes the data-driven mine-to-mill framework proposed and tested in [19,20]. It was concluded that material tracking needs to be performed in a more comprehensive way to better connect mining and mineral processing data, as well as stockpile management, and the required infrastructure was found to be the main limitation for implementing the framework. A closed-loop framework for real-time reserve management was introduced in [21]. It incorporated sensor-based material characterisation, geostatistical modelling under uncertainty, data assimilation for a sequential model updating and mining system simulation and optimisation. Using a modified Kalman filter technique, online sensor and process information was back-propagated into the resource block model, which in turn was used for production planning and other operative decisions. Innes et al. [22,23] developed a method for representing, tracking and fusing the information of excavated material. The bulk material was represented as lumped masses, although not as particles, and no other properties than mass and location. Tracking along the transport and storage process was achieved using an augmented state Kalman filter with a mass preservation constraint. Reconciliation of the bond mill work index to the ore block model using mine-to-mill tracking was demonstrated in [24], using a first-in-first-out model for the stockpile. A 3D real-time stockpile model and mapping technique of product quality was developed in [25] for the purpose of planning and control in stacking, reclaiming and blending. The pile surface was reconstructed from the point-cloud of a mobile scanning device. A voxelized model was then associated with a reclaiming machine to achieve autonomous operation and predictive quality. The simulator in [26] demonstrated the combination of multiple, relatively simple, models for capturing the ore flow through a complex storage facility. No previous work was found where material tracking was performed using a particle representation.

2. Digital Twin as a Distributed Particle Simulation

From the material point of view, the digital twin of a mine may be considered a system of distributed particle simulation. For computational efficiency, we used pseudo-particles, each representing a large collection of real rock particles of various sizes, shapes and mineral compositions. Labelled n, a pseudo-particle has position ${\mathit{x}}_n\left(t\right)\in {\mathbb{R}}^3$, velocity ${\mathit{v}}_n\left(t\right)={\dot{\mathit{x}}}_n$, mass $m_n\in {\mathbb{R}}^1$, size distribution function ${\mathit{f}}_n\in {\mathbb{R}}^{N_{\mathrm{f}}}$, mineral concentration ${\mathit{c}}_n\in {\mathbb{R}}^{N_{\mathrm{c}}}$ and mechanical properties ${\mathit{h}}_n\in {\mathbb{R}}^{N_{\mathrm{h}}}$ that characterize the real particles it represents, each discretized in classes of dimensionality $N_{\mathrm{f}}$, $N_{\mathrm{c}}$ and $N_{\mathrm{h}}$, respectively. The particle state vector is represented as ${\mathit{X}}_n=[m_n,{\mathit{x}}_n,{\mathit{v}}_n,{\mathit{f}}_n,{\mathit{c}}_n,{\mathit{h}}_n]$. Each particle has a position of origin ${\mathit{x}}_n\left(t_0\right)$. The material properties may be known to some degree of certainty at that location from prior exploration and represented in a block model of grid cells indexed i, with centre point ${\mathit{x}}_i$ and volume $V_i$, or they are the subject for determination by measurements in the process downstream. The fragmentation is set by the blasting model and/or measurements after blasting.

The fixed and mobile mining equipment constitute a discrete set of assets that transform the material by a set of unit operations, the primary two being move and fragment. Each pseudo-particle was tracked, or simulated, by a virtual asset, which is a real-time model that mirrors a physical asset and its unit operations on the material. We distinguished between high-frequency micro-tacking inside each asset and macro-tracking of special events, e.g., when particles enter or exit an asset or pass a point of special interest. This is illustrated in Figure 2.

For a bucket excavator or truck asset, the micro-tracking can be accomplished simply by inheriting the position of the bucket or the truck during the transport from the loading to the dumping area. Particles transported on a belt conveyor may similarly be tracked between the loading and discharge points by simply tracking the movement of the belt using the conveyor control system. In the transporter’s virtual asset, the belt is discretized longitudinally in bins of known location and particle content. Chutes, silos and stockpiles are different. In these assets, the particle flow is passive, driven by gravity and feeders at the outlet for a controlled discharge rate. Micro-tracking requires a model that may range from simple first-in-first-out or perfect mixing models to a physics-based flow model.

Mathematically, we represent a unit operation by an asset a as an operator ${\widehat{u}}_a(t,\Delta t)$ that propagates a particle state ${\mathit{X}}_n\left(t\right)$ forward in time into the new state ${\mathit{X}}_n(t+\Delta t)={\widehat{u}}_a(t,\Delta t){\mathit{X}}_n\left(t\right)$, with propagation time-step $\Delta t$. If the particle n enters the asset at time $t_{\mathrm{in}}$ in a state ${\mathit{X}}_n\left(t_{\mathrm{in}}\right)$, the exit state at time $t_{\mathrm{out}}$ is given by:

$${\mathit{X}}_n\left(t_{\mathrm{out}}\right)=\prod _{j=0}^{J_n^a-1}{\widehat{u}}_a(t_{\mathrm{in}}+j\Delta t){\mathit{X}}_n\left(t_{\mathrm{in}}\right),$$

(1)

where $J_n^a$ is the number of unit operations by the asset. Micro-tracking resolves the evolution of the particles this way. In general, the step size $\Delta t$ may vary, and the unit operation is a function of both time and the state ${\mathit{X}}_a$ of the particles in the asset, i.e., ${\widehat{u}}_a(t,{\mathit{X}}_a\left(t\right))$. Note that this representation supports both kinematic and dynamic models for the particle motion, for move operations. Discrete element simulations, taking real-time asset data as the input, are an example of the latter. Letting the particles co-move with a vehicle on the conveyor is an example of the former. The evolution of the particle size distribution may also be represented by this, for fragment operations.

Macro-tracking is represented as message passing, between virtual assets or a centralized tracking service, about events of special interest and the involved particles. The message should have a structured data format, such as:

$$\begin{array}{cc}\hfill \mathrm{event}=& \left\{\mathrm{id},\mathrm{type},\mathrm{location}\mathit{x},\mathrm{time}t,\mathrm{particles}=\left\{\left[n_1,{\mathit{X}}_{n_1}\left(t\right)\right],\left[n_2,{\mathit{X}}_{n_2}\left(t\right)\right],\cdots \right\}\right\}\hfill \end{array}$$

(2)

including an identifier for the event, its type, location and time, plus a set of particle identifiers and state tuples. The macro-tracking data may then be stored in a database for later analytics or other services based on tracking data. An illustration of the architecture of the system of assets and services is found in Figure 1.

Most assets are equipped with sensors that produce observations while handling the material. An observation is either a direct measurement of an intrinsic attribute or an indirect measurement. The latter may be registered as an extrinsic attribute that is assigned to the pseudo-particles at the time and location for the observation. Examples of extrinsic attributes are measurement-while-drilling (MWD) data, diggability, crushability and grindability. In complement real sensors, it is also possible to produce virtual sensors that enable monitoring or controlling a process at higher precision. A virtual measurement signal of an attribute p is accomplished by extracting the pseudo-particles $\mathcal{N}(\mathit{x},t)$ in a selected volume $V(\mathit{x},t)$, centred around $\mathit{x}$, at time t, and computing the mass-weighted attribute:

$$p(\mathit{x},t)={〈\frac{{\sum }_{n\in \mathcal{N}(\mathit{x},t)}{m}_n{p}_n}{{\sum }_{n\in \mathcal{N}(\mathit{x},t)}{m}_n}〉}_{t_{\mathrm{w}}},$$

(3)

with a moving average ${〈\cdots 〉}_{t_{\mathrm{w}}}$ of time window $t_{\mathrm{w}}$ for filtering the effects of coarse particles.

Each pseudo-particle has an individual size distribution, ${\mathit{f}}_n=[f_n^1,f_n^2,\cdots ,f_n^{N_{\mathrm{d}}}]$, discretized in $N_{\mathrm{d}}$ size classes, where $f_n^k$ denotes the mass percentage of particles in the size span $[d_{k-1},d_k]$. The averaged size distribution in a volume $V(\mathit{x},t)$ is computed using Equation (3) to $\mathit{f}(\mathit{x},t)$ and the mean particle diameter by $〈d〉(\mathit{x},t)=\mathit{d}·\mathit{f}$, where $\mathit{d}$ are the diameters of the particle bins. The cumulative size distribution $\mathit{F}$ is the mass fraction of particles equal to or larger than each size class $d_k$, such that $F^k={\sum }_{k^{\prime }=k}^{N_{\mathrm{s}}}{f}^{k^{\prime }}$. The $D_F$ value is the diameter d for which the cumulative size distribution equals F. It is easily determined from the cumulative size distribution F by direct inspection or interpolation of the nearest discrete size classes.

When particle simulation was needed for micro-tracking, we used the discrete element method (DEM) with spherical pseudo-particles. The simulation attributes include particle diameter $d_{\mathrm{p}}$, coefficient of restitution $e_{\mathrm{p}}$, friction coefficient ${\mu }_{\mathrm{p}-\mathrm{t}}$, rolling resistance coefficient ${\mu }_{\mathrm{p}-\mathrm{r}}$ and particle cohesion $c_{\mathrm{p}}$. The parameters were calibrated to values that approximated the mechanics and flow dynamics of the real material on the bulk level, as measured by the bulk mass density ${\rho }_{\mathrm{b}}$, internal friction ${\mu }_{\mathrm{b}}$, cohesion $c_{\mathrm{b}}$ and angle of repose. The particle diameter was chosen as large as possible, but small enough for resolving the important flow features. Even with pseudo-particles as coarse as a loader bucket, the number of particles in a silo or stockpile asset may exceed what can be simulated in real time with present hardware. To guarantee the real-time performance of the particle dynamics, we employed data-driven model order reduction [18]. The idea was to run numerous DEM simulations in advance, covering the relevant state-space, as well as possibly using authentic CAD models and control signals for in- and out-flow to the asset. From the simulation data, a low-dimensional model was trained to predict the flow field that occurred in a stockpile or silo as a function of the control signal and mass distribution. The computationally intense process of computing all contact forces was substituted by a simple evaluation of the predicted velocity field at the positions of the particles, which were simply advected using this velocity. This was applied to the particles in the bulk. To accurately capture the dispersion of particles impacting the pile surfaces, resolved DEM was used. The combination of the two models in an executable application, part of a distributed system, is referred to as granular surrogates.

3. Integration Test

An integration test of using granular surrogates for micro- and macro-tracking was performed; see Figure 3. A virtual stockpile asset was created using AGX Dynamics for the particle simulation. To accelerate the simulations to real-time performance, a reduced-order model was trained from ground truth simulation data, as described in [18], with a tracking accuracy of 10%. The storage asset was connected to an ABB 800xA process simulator [27] using OPC [28] for receiving a discharge control signal and for exchanging macro-tracking data. The ABB simulator represents the assets of the incoming and outgoing conveyors, discretized in bins moving with the controlled conveyor belt speed. Each bin on the incoming conveyor contains material with a specific mass, ore grade and a tracking id. That information is exchanged with the stockpile asset when a bin reaches the discharge point. The storage asset responds by emitting particles with the corresponding total mass and mineral concentration. The particles fall, impact and disperse over the pile surface as simulated by the physics-based stockpile model. The storage asset receives a discharge rate control signal from the ABB simulator. This drives a flow field in the bulk, and particles are propagated in real time by the data-driven reduced-order model. The particles that reach the outlet are eliminated, and a macro-tracking message is passed so that the ABB simulator receives information of the material that enters the outgoing conveyor. That information includes tracking ids, mass and ore grade. The micro-tracking data reside in the storage asset. A monitoring service was created as a stand-alone web application. It receives virtual sensor information from the stockpile asset using a REST (representational state transfer (REST): https://en.wikipedia.org/wiki/Representational_state_transfer; date accessed on 1 April 2021) interface. The monitoring information includes the total mass in the storage, the average grade and the spatial distributions of these quantities, as well as the predicted velocity field in the bulk of the pile. A video that demonstrates the integration test is included as Supplementary Material Video S1. It should be pointed out that the ABB 800xA simulator can be switched to sharing real-time data from a mine with the same layout.

4. Simulator

Before implementation in a real mine, particle-based material tracking was tested in a simulator of a simplified 2D model of a mine shown in Figure 4. A video of the simulator is included as Supplementary Material Video S2. Each asset can easily be replaced with more elaborate 3D models or virtual assets mirroring real physical assets. There are two sources of ore from which haul trucks transport blasted material to a single crusher. The crushed material is fed to a stockpile for intermediate storage, from which material is drawn and fed to a mill. Belt conveyors transport the material among the crusher, storage and mill. The simulator was used for testing the analysis that was enabled by the particle-based tracking and for understanding the sensitivity to tracking errors.

Macro-tracking events include the haul trucks dispatching from the loading point, dumping at the crusher station, material passing the crusher zone, the inlets and outlets to the stockpile and the inlet to the grinder. Model order reduction was applied on the belt conveyors. The simulator was implemented using the physics engine AGX Dynamics [29], using a nonsmooth (or time-implicit) DEM [17] for the particles and rigid multibody dynamics for the asset’s moving parts using feedforward controllers. The pseudo-particle size ranged between $0.2$ and $0.48$ m, which corresponded to a mass range between 12 and 160 kg. The distance between the two loading points was roughly 150 m, and the stockpile was 60 m tall. The simulations were run with a time-step of 11 ms and 100 projected Gauss–Seidel solver iterations, which corresponded to an error tolerance of around 5% for the stockpile bulk behaviour [17]. At the sources, each particle was assigned an identity and a location of origin. The left source was assigned a mineral concentration of $0.2$, while the right source was assigned the value $0.8$. The right source was given a mineral hardness of $1.0$, while the hardness at the left source varied over the ore body according to some function $h\left(\mathit{x}\right)$. It was tested how a step in the hardness distribution function, from a value of $1.0$ to $2.0$, in the left source propagated through the system. It was also investigated how well a hardness profile in the ore body could be backtracked from observations at the crusher and the grinder. The (true) particle size distribution was discretized in $N_{\mathrm{d}}=16$ size classes with a logarithmic binning. The initial size distributions is shown in Figure 5, including also the average size distribution after the crusher. Note that the pseudo-particle diameter, used for visualization and contact dynamics, had no correspondence to its internal particle size distribution.

Each haul truck received a load of roughly 15 tons (roughly 150 pseudo-particles) at the loading point. For simplicity, the loading equipment was not modelled explicitly. The trucks dumped the load on the crusher inlet, and there, mixing of the material sources occurred from here on. The crusher was modelled as a jaw crusher, drawing material from above and crushing it with a jaw driven in an oscillatory motion. The pseudo-particles retained their original diameter, but when they passed the jaw, the internal particle size distribution was transformed. The fractions larger than the crusher gap width, $d_{\mathrm{cr}}=0.25$ m, broke into a uniform distribution of the smaller size classes. This unit operation was represented by ${\mathit{f}}_n:={\mathit{f}}_n-\Delta {\mathit{f}}_{n-}+\Delta {\mathit{f}}_{n+}$, with the change vector for crushed (oversize) fractions, $\Delta f_{n-}^k=f_n^k$ for $d_k\ge d_{\mathrm{cr}}$ (otherwise zero) and the change vector for the $n_{\mathrm{d}}$ undersized fractions $\Delta f_{n+}^k=\frac{1}{n_{\mathrm{d}}}{\sum }_{k^{\prime }}\Delta f_{n-}^{k^{\prime }}$ (otherwise zero). The resulting size distribution is shown in Figure 5. We assumed that the torque for driving the jaw crusher was in proportion to the amount of mass in the crusher zone $V_{\mathrm{cr}}\left(t\right)$, its hardness and size distribution, according to the following model:

$${\tau }_{\mathrm{cr}}\left(t\right)=\sum _{n\in {\mathcal{N}}_{V_{\mathrm{cr}}}}{m}_n{h}_n{\mathit{k}}_{\mathrm{cr}}·{\mathit{f}}_n,$$

(4)

where the crushing coefficient $k_{\mathrm{cr}}^n=1.0$ for the oversized and $k_{\mathrm{cr}}^n=0.1$ for the undersized. The momentaneous power draw of the crusher was computed by $P_{\mathrm{cr}}\left(t\right)={\tau }_{\mathrm{cr}}{\omega }_{\mathrm{cr}}$, where ${\omega }_{\mathrm{cr}}=10$ rad/s was the crusher speed.

The stockpile held up to 2000 tons (roughly 20,000 pseudo-particles). There were two entry discharge points and four exit discharge points, but only one was used for the presented tests. When material was discharged from the bottom outlet, at a controlled speed ranging up to 2 m/s (or two tons/s), a funnel flow field was activated in the pile, reaching the surface where a depletion was gradually formed, unless new material was fed into the storage at the same rate or faster. Two time instances of pile states are shown in Figure 6. For the pointy pile state, one can expect more mixing of the incoming material and longer residence time, contrary to the depleted pile state where material flowed more directly through the storage and with less mixing.

An idealized mill model was assumed, guided by a model that was used for mill control at Boliden Aitik [30]. The idealized model consisted of a single grinder, although most real mills consist of a circuit of grinders and classifiers. When a pseudo-particle n reached the the grinder, a size reduction model began to operate on its size distribution. Each size fraction k was assigned a probability per unit time, ${\mathcal{P}}_n^k$, for breaking into smaller fractions and a probability per unit time, ${\mathcal{Q}}_n^k$, for leaving the mill. At each time-step, the mass fractions carried by each pseudo-particle were redistributed according to ${\mathit{f}}_n:={\mathit{f}}_n-\Delta {\mathit{f}}_{n-}+\Delta {\mathit{f}}_{n+}$, where the broken and passing fractions were:

$$\Delta f_{n-}^k=({\mathcal{P}}_n^k+{\mathcal{Q}}_n^k)f_n^k\Delta t.$$

(5)

and the reduced fractions were:

$$\Delta f_{n+}^k=\frac{1}{N_{\mathrm{d}}-k}\sum _{k^{\prime }=k+1}^{N_{\mathrm{d}}}{\mathcal{P}}_n^{k^{\prime }}{f}_n^k\Delta t.$$

(6)

Modelling a size classifier, the probability per unit time for the finest fractions, smaller than $d_{\mathrm{cl}}$, to leave the mill was ${\mathcal{Q}}_n^k=k_{\mathrm{esc}}/(1+exp\left[k_{\mathrm{cl}}(d_k/d_{\mathrm{cl}}-1)\right])$. We used $d_{\mathrm{cl}}=0.025$ mm, $k_{\mathrm{cl}}=10$ and $k_{\mathrm{esc}}={10}^3$. The breakage probability depended on the (true) particle size, impact energy, hardness and the number of impacts [31], which in turn depended on the mill’s angular velocity, ${\omega }_{\mathrm{gr}}$, and the active torque, ${\tau }_{\mathrm{gr}-\mathrm{a}}$. We assumed the following model for breakage probability per unit time:

$${\mathcal{P}}_n^k=k_{\mathrm{gr}}\frac{E_{\mathrm{imp}}^k}{E_{n,\mathrm{brk}}^k}{\left(\frac{{\tau }_{\mathrm{gr}-\mathrm{a}}{\omega }_{\mathrm{gr}}}{{\tau }_{\mathrm{gr}-\mathrm{dr}}^0{\omega }_{\mathrm{crit}}}\right)}^2,$$

(7)

where $E_{\mathrm{imp}}^k=E_{\mathrm{imp}}^{\infty }\left(e-exp[1-D_{50}/d_k]\right)d_k/d_{\mathrm{gr}}$ is an impact energy distribution function, $E_{n,\mathrm{brk}}^k=E_{\mathrm{brk}}^0{h}_n/[1+{(d_k/d_{\mathrm{gr}})}^{1/4}]$ is a specific breakage energy distribution, $D_{50}$ is the (time-dependent) mean diameter of the material currently in the grinder, $d_{\mathrm{gr}}=1.0$ mm is a grinding length-scale parameter and $k_{\mathrm{gr}}=1.0$ a dimensionless breakage coefficient. The internal drive friction, ${\tau }_{\mathrm{gr}-\mathrm{dr}}^0$, of an empty mill and the critical mill speed, ${\omega }_{\mathrm{crit}}=1.5$ rad/s, were used for normalization. The mill breakage model can be seen in Figure 7, with $\tilde{\mathcal{P}}\left(d_k\right)\equiv \mathcal{P}\left(d_k\right)/\frac{{\tau }_{\mathrm{gr}-\mathrm{a}}{\omega }_{\mathrm{gr}}}{{\tau }_{\mathrm{mill}}{\omega }_{\mathrm{crit}}}$. The exit probability $\tilde{\mathcal{Q}}\left(d_k\right)\equiv \mathcal{Q}\left(d_k\right)/k_{\mathrm{esc}}$ is also included.

The total grinding torque was the sum of three contributions, ${\tau }_{\mathrm{gr}}={\tau }_{\mathrm{gr}-\mathrm{dr}}+{\tau }_{\mathrm{acc}}+{\tau }_{\mathrm{gr}-\mathrm{a}}$. The total mass, $m_{\mathrm{tot}}\left(t\right)=m_0+m_{\mathrm{gr}}\left(t\right)$, was the mass of the moving part of the mill, $m_0={10}^3$ kg, plus the mass of the material in the grinder, $m_{\mathrm{gr}}\left(t\right)={\sum }_n{m}_n$. The internal drive friction was modelled ${\tau }_{\mathrm{gr}-\mathrm{dr}}={\mu }_{\mathrm{dr}}{m}_{\mathrm{tot}}g{r}_{\mathrm{mill}}$, and the torque required for accelerating the grinder was ${\tau }_{\mathrm{acc}}=m_{\mathrm{tot}}{r}_{\mathrm{mill}}^2{\dot{\omega }}_{\mathrm{gr}}$, where we denote the internal friction coefficient ${\mu }_{\mathrm{dr}}=0.02$, gravity acceleration $g=9.8$ m/s${}^2$ and mill radius $r_{\mathrm{mill}}=5$ m. The active torque depended on the charge’s mass, $m_{\mathrm{gr}}={\sum }_n{m}_n$, size distribution and dynamic angle, which depended on the mill angular velocity. We assumed, again guided by [30], that the active torque could be modelled ${\tau }_{\mathrm{gr}-\mathrm{a}}={\mu }_{\mathrm{gr}-\mathrm{a}}{m}_{\mathrm{gr}}g{r}_{\mathrm{mill}}$ with effective internal friction of the charge ${\mu }_{\mathrm{gr}-\mathrm{a}}=k_{\mathrm{gr}-\mathrm{a}}\Theta ({\omega }_{\mathrm{gr}}/{\omega }_{\mathrm{crit}})\frac{m_{\mathrm{gr}-\mathrm{a}}{\omega }_{\mathrm{gr}}}{m_{\mathrm{gr}}{\omega }_{\mathrm{crit}}}$. Here, $m_{\mathrm{gr}-\mathrm{a}}={\sum }_{n,k\ge N_{\mathrm{gr}-\mathrm{a}}}{m}_n^k$ is the mass of the charge fractions larger than an active grinding size $d_{\mathrm{gr}-\mathrm{a}}$, with index $N_{\mathrm{gr}-\mathrm{a}}$. We used $d_{\mathrm{gr}-\mathrm{a}}=1.0$ mm, which was $N_{\mathrm{gr}-\mathrm{a}}=7$ in the 16 class size representation, and $k_{\mathrm{gr}-\mathrm{a}}=1$. With $\Theta \left(x\right)=Re\left[{(1-x^4)}^{1/2}\right]$, the effective internal charge friction and active torque dropped to zero at and above the critical velocity. The power draw of the grinder was computed as $P_{\mathrm{gr}}\left(t\right)={\tau }_{\mathrm{gr}}{\omega }_{\mathrm{gr}}$. When nothing else was stated, the mill was run at constant nominal speed ${\omega }_{\mathrm{gr}}=1.2$ rad/s.

5. Results

5.1. Step Response

The simulator was initialized with ore of constant hardness $h=1$ transported from both sources. A step change to hardness $h=2$ was introduced at the left source, and the response in the crusher may be seen in Figure 8. The haul truck reached this harder material at time $t=350$ s. About 60 s later (five seconds more than one truck cycle), the crusher power draw increased by 40%, from $1.5$ kW to $2.1$ kW. A virtual sensor in the crusher revealed that the mineral hardness varied between one and two, with a time-averaged hardness around $1.4$. The large oscillations in hardness and power draw, with a 55 s time period, was because of the trucks dumping interchangeably from the left and right source with hardness of one and two, respectively.

For the step response in the mill, we analysed the effect of the ore passing through the stockpile in Figure 6. The discharge of the piled material started at 500, and the feed reached the empty mill shortly thereafter. The material in the mill remained as hardness $h=1$ for 90 more seconds until it started increasing at time 590 s up to $h=1.5$ with only small oscillations (well mixed). At time 875 s, the storage level was down at zero (apart from the dead zones), and the rate of mass flowing into the mill dropped quickly to the rate of mass leaving the crusher and flowing into the storage. Therefore, after 875 s, the oscillations in hardness and power again reflected the truck cycles. For reference, we ran an identical simulation with constant hardness $h=1$ at both sources. As seen in Figure 9, the difference in the two power draw signals followed the increase in material hardness inside the mill. The increase was about 20% on average. Note that the mechanism for the higher power draw was that the harder ore took a longer time to grind, and the mill therefore accumulated more mass over time.

5.2. Identification of Hardness from Measurement-While-Drilling

To test the feasibility of identifying the rock hardness from the MWD signal, we assumed the quadratic hardness profile $h\left(x\right)$ in Figure 10 at the left source, ranging continuously between one and two. While drilling, a measurement profile $s\left(x\right)$ was measured. For the test, we assumed $s=h^{0.4}$ and introduced a measurement location error such that $s_n=h{(x_n+\mathcal{R}(-\Delta x,\Delta x))}^{0.4}$, where $\mathcal{R}(-\Delta x,\Delta x)$ is a uniform disturbance. We used $\Delta x=12.5$ m, which corresponded to six haul truck cycles back and forth to the crusher. The right source had constant hardness $h=1$. Ore from the two sources entered the crusher, partially mixed, and gave rise to the torque signal in Figure 11.

Assuming the crusher torque is related to rock hardness according to Equation (4), any ansatz $s(h,\beta )$ can be explored and optimized with respect to model parameter $\beta$. Different outcomes for the ansatz $s=h^{\beta }$ are shown in Figure 11. The best fit was found at $\beta =0.402$, which deviated by $0.5$% from the ground truth because of the tracking location error that was introduced. Note that empirical models between the ore tensile strength and adjusted drill penetration rate have this form in the literature [12]. Once a relation between the MWD signal and hardness was established, it was possible to predict the hardness of the material that was about to enter the crusher given the MWD signal measured at its location of origin. This is illustrated in Figure 12 over a sampling duration of 2500 s, or roughly 40 haul trucks from the left source.

5.3. Predicting the Properties of the Mill Feed

With particle-based tracking and Equation (3), it was straight-forward to monitor the properties of the ore kept in the storage as a function of time and space. Given a storage flow model, it was also possible to estimate which particles in the funnel zone would leave the storage in the near future, say the next 30 s. This way, the properties of the ore feed that would reach the mill could be predicted in advance, even before exiting the storage. Sample spatial distributions are shown in Figure 13. The time evolution of the properties in the pile and of the future feed are shown in Figure 14.

Backtracking the Grindability

The grindability of an ore may be measured by the amount of energy, per unit mass, that is required for reducing the fragmentation down to a specific particle size. For the sake of production planning, the grindability is often estimated in advance using the relatively sparse exploration data of the ore body. The Bond ball mill work index (BWI) is defined as:

$$W=\frac{\mathrm{BWI}}{3.6}\left[\frac{10}{\sqrt{〈D_{80}^{\mathrm{out}}〉}}-\frac{10}{\sqrt{〈D_{80}^{\mathrm{in}}〉}}\right],$$

(8)

where the work per unit mass (W/kg) is computed as $W\left(t\right)=〈P_{\mathrm{gr}}〉/〈{\dot{m}}_{\mathrm{gr}-\mathrm{out}}〉$ using a moving average with time-window $t_{\mathrm{av}}$. $D_{80}$ is the 80% passing diameter (m), computed from the cumulative size distribution function for mass that enters and leaves the mill. The bond grindability is then defined as the inverse $G\equiv {\mathrm{BWI}}^{-1}$. We demonstrated that the grindability measured in the mill can be backtracked to the ore’s location of origin for the purpose of improving future estimates for grindability. Assume there is an a priori, $G_i^0$, estimate of the grindability of the ore in block i. At time $t_k$, a grindability $G\left(t_k\right)$ was observed at the mill. The average mass in the mill originating from the block i was $〈m_{\mathrm{gr}}^i〉\left(t_k\right)$. We reconciled the grindability at block i and time $t_k$ by the weighted average:

$$G_i\left(t_k\right)=\frac{{\alpha }_i^0{G}_i^0+{\sum }_{k^{\prime }=1}^k{\alpha }_i\left(t_{k^{\prime }}\right)G\left(t_{k^{\prime }}\right)}{{\alpha }_i^0+{\sum }_{k^{\prime }=1}^k{\alpha }_i\left(t_{k^{\prime }}\right)},$$

(9)

where the weight factor is given by the ratio of the mass with origin i to the mass in the mill and relative to the block mass ${\alpha }_i\left(t_{k^{\prime }}\right)={\langle m_{\mathrm{gr}}^i\rangle }^2/m_i^0\langle m_{\mathrm{gr}}\rangle$. Here, $\langle m_{\mathrm{gr}}\rangle \left(t_k\right)$ is the amount of mass in the mill and $m_i^0$ is the total mass of the ore block i. Note that the weight factor is small (or zero) when there is little (or no) mass from the considered block. The weight factor for the original grindability reflects the confidence of that estimate, ${\alpha }_i^0\ll 1$ if the confidence is low and ${\alpha }_i^0\gtrsim 1$ if the confidence is high. We tested backtracking the grindability of the ore originating from the left source for the case of a quadratic hardness profile shown in Figure 10, but mixed in the mill with ore from the right source. We used a block size of 10 m and ${\alpha }^0=0.1$. The BWI signals are shown in Figure 15, and the grindability backtracked to the source using Equation (9) is found in Figure 16.

6. Discussion

A digital twin of a mine may be realised as a system of particle simulations, distributed over virtual assets that evolve according to their physical counterparts using real-time particle simulation where no real sensors or control signals are available. This enables material tracking of ore from multiple sources passing through complex storage systems, such as stockpiles with multiple inlets and discharge outlets. Material tracking is necessary for understanding how the mill performance depends on natural variations in the raw material at the different sources and on the actions taken at drilling, blasting and crushing.

A concrete application example is to optimise the mill’s throughput energy consumption by controlling the fragmentation of the ore in relation to its hardness, known as mine-to-mill optimisation. It is widely accepted that energy savings up to 30% can be obtained by doing so [32]. First, the grindability of the ore body must be understood, by measuring the energy consumption in the mill and back-tracking it to the source. The next step is to determine how the grindability depends on the fragmentation and hardness of the material. The latter can be predicted from MWD data in combination with exploration data [12]. An optimal target fragmentation can then be predicted given the local material properties and realised by a specific drill and blast plan. The next level of optimisation is to control the stockpile flow so that the mill feed maintains an even and high grindability or pre-emptive regulation of the grinder in cases where variations are unavoidable. In [33], a 60% increase in the plant throughput was achieved this way, by using physical tracers to track the material from different sources to the mill via crushers and storages. That is, however, not sustainable over time, and a digital tracking system must be used in place.

This application puts a number of requirements on the digital material model used for tracking the material through the stockpile. To resolve the flow from different sources, the pseudo-particles must at least be smaller than a truck load. When the material lands on the stockpile surface, it is dispersed, which requires yet smaller pseudo-particles. The dispersion depends on the current shape of the pile surface, which may be complex when there are multiple discharge outlets, and on the material properties. When material is drawn, a flow channel develops, forming an inverted cone at the surfaces. If the incoming flow lands directly on the inverted cone, there is little dispersion, but if it lands on a slope or a conical top, the dispersion may be large: this is why a first-in-first-out or perfect mixing models is not suitable. Consequently, a grinder may at any time contain ore from many different trucks, and its content needs to be discretized sufficiently. Assuming the mill has a significant contribution of mass from 100 different trucks, thanks to dispersion and mixing in the stockpile, the mill mass must be resolved with at least 1000 pseudo-particles. For a medium-sized SAG mill, that implies a particle size less than $0.5$ m. A medium-sized haul truck is then resolved by more than 100 particles that may disperse on the pile. This size is also small enough to pass through a stockpile outlet and resolve the funnel flow in the pile.

There are a number of uncertainties to consider. The shape of the stockpile surface needs to be captured sufficiently. This can be done using the pseudo-particle model alone or by combining it with online measurements of the surface height to prevent the model drifting away from the true state. The particle contact parameters should be calibrated such that the dispersion over the surface is captured. This can be performed in advance using small test piles of the real material. The model for the funnel flow needs to be verified and calibrated. A measurement campaign with physical traces is the recommended way, measuring the residence time for each tracer and the evolving shape of the pile surface. In a cold climate, it might be necessary to carry this out several times to capture the temperature effect on the size of the dead zones and funnel flow, especially for temperatures below the freezing point. One limitation with the pseudo-particle model is that it does not automatically support size segregation effects, which may be significant when material is dispersed on a pile. The sieving mechanism makes large particles slide and roll further down the material bed than small ones. This can be dealt with by giving the pseudo-particles different internal size distributions and masses (maintaining the net size distribution and mass) and giving coarse and fine pseudo-particles different sizes and contact parameters such that the segregation effect evolves naturally.