**1. Introduction**

The increased infrastructure development of society has been accompanied by an increase in the use of natural materials for construction [1]. One resource that has seen massive extraction is natural sand deposits. The usefulness of sand in creating quality concrete products has led to the erosion of numerous coastlines caused by aggressive sand dredging [2]. Some governments such as Sweden are trying to incentivize new approaches with increased taxes to minimize the use of natural sand along with limiting permits for existing deposits.

Given the finite nature of natural sand and the increased demand on the market, replacement resources are becoming more common. One approach is replacing natural sand with manufactured sand. The process of creating manufactured sand for concrete includes vertical shaft impactors (VSIs), which are common comminution machines for crushing particles in the tertiary or quaternary stage of aggregate crushing circuits. Autogenous impact crushing results in particles that, compared to compressive crushing, have a smoother surface and a more spherical shape but with an increased amount of fine particulate matter in the process. Manufactured sand has successfully been used to replace natural aggregates in several applications, such as concrete production [3,4]. The main attributes sought after in natural sand particles are sphericity and high volume-to-surface ratio. These attributes minimize the use of cement in the concrete mixes and lead to desirable rheology for pouring [5].

Natural sand requires fewer processing steps than the manufacturing of artificial sand. To allow for a transition to more sustainable products, the process needs to become more

**Citation:** Grunditz, S.; Asbjörnsson, G.; Hulthén, E.; Evertsson, M. Fit-for-Purpose VSI Modelling Framework for Process Simulation. *Minerals* **2021**, *11*, 40. https:// doi.org/10.3390/min11010040

Received: 31 October 2020 Accepted: 28 December 2020 Published: 31 December 2020

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

efficient. By establishing a framework to adopt crusher models to existing production plants, simulations can be used to optimize settings at real plants, and therefore, a more efficient process can be achieved.

The dynamics of the particles has been shown to determine the product size distribution. A detailed description of the dynamics of the particles inside and exiting a rotor was provided by Rychel [6]. However, Kojovic [7] formulated a predictive model for the VSI crusher that correlated product size with power draw based on breakage of the ore in a drop-weight test. Attempts to model the VSI output by adapting the Whiten model as a basis for modeling have also been proven to be a viable approach [8].

According to Nikolov [9], the essential part of the performance of the VSI crusher is the relationships between the radius of the rotor, angular velocity and the rate of material being fed into the VSI crusher when given a specific particle size distribution for the feed material. Bengtson and Evertsson [10] then proposed a model that also considered the variations in rock properties and predicted capacity and power draw as well as the particle shape distribution. Lindqvist [11] also suggests that impact crushing can significantly reduce the energy required in grinding.

The attempts to model the performance of VSI crushers have often been limited by the lack of information about the interaction between particles inside the crusher. An approach to overcome this issue has been through the use of the discrete element method (DEM). With the DEM, every particle and its dynamics are simulated, making it possible to estimate the particle trajectories, residence times and collision energies when particles travel through the crusher. DEMs have been the main method of choice for other researchers such as Cunha, Cleary and Sinnott [12–14].

In Grunditz et al. [14], the data from DEM simulations of a VSI crusher were analyzed to gain an understanding of how the energy in a VSI crusher is distributed in relation to the different sizes of the particles fed into it and the number of impacts that occur. One of these energy collision spectra can be seen in Figure 1. In addition, an exploration of how the particles travel in the crushing chamber after leaving the rotor was conducted in conjunction with evaluating the residence time of each particle [15,16].

**Figure 1.** The distribution of particle collision energies inside a discrete element method (DEM) simulation of a vertical shaft impact (VSI). Two different rotor tip speeds were simulated.

To create a robust model, material from numerous sites was collected and used for the calibration and validation of the model. The model accurately predicts the product size distribution and is useful to give operators the opportunity to simulate changes to their feed size, material attributes and machine parameters before committing to them.

The use of data from several different sites of varying geological types enables the mapping of material properties to certain model parameters, achieved by deploying the same crushing process and equipment in the same configuration at multiple different crushing plants. These physical material properties were attributed to *fmat* and *E0 min* established in previous modelling efforts by Nikolov [9,17].

A selection–breakage approach is used to simulate the effects of impact energy present in VSI crushers. Vogel and Peukert established a master curve for breakage that considers the accumulation of impacts when observing particle breakage [18,19].

The purpose of this paper is to use previously established models and appropriate optimization routines for calibration of a fit-for-purpose VSI crusher model. A framework was proposed to enable crushing operators, plant owners and product managers to explore and test differing machine configurations to reach alternate product goals more systematically.

#### **2. Methodology**

#### *2.1. Modeling*

Breakage modeling relies on descriptive functions of particle behavior when subjected to high energy impacts. The goal of these studies was to assess how much energy a particle requires to break relative to its initial size. In addition, the size of the daughter fragments following such an event was of interest.

Several previous studies and experiments [17] provide the basis to understand the role of particle size and impact energy on the breakage of the particular materials studied. Each material tested in this way provides properties that apply to only one specific quarry site. Having to run these exhaustive tests for each modelling attempt for a site is costly.

The different sizes into which a particle will statistically split are referred to as an appearance function. An appearance function was included in the optimization routine and assumed that the behavior of the breakage follows a polynomial function, as seen in Equation (1), where *ai,j* is the fraction of particles of one size being reduced into a smaller subset; *β<sup>i</sup>* values are coefficients, and *xi* is the size of the originating particle. The index *i* represents the number of particle sizes available in the PSD dataset starting from 1, the largest size, to N being the smallest size. This creates an appearance matrix for the selected particles sizes. An example of such an appearance matrix can be seen in Table 1 where each column for progeny Equation (1) is shown.

$$\boldsymbol{a}\_{i,j} = \beta\_1 + \beta\_2 \mathbf{x}\_i^1 + \beta\_3 \mathbf{x}\_i^2 + \beta\_4 \mathbf{x}\_i^3 \tag{1}$$

**Table 1.** The appearance matrix for a PSD with 5 different size bins. Each column for progeny Equation (1).


The most common approach to estimate the appearance function is with a drop weight test or a Julius Kruttschnitt rotary breakage tester (JKRBT) to generate a family of *t10* curves [20]. In [21], three different ore types were tested with the JKRBT [22], and their response was calculated with a modified equation based on the work of Vogel and Peukert [19] with satisfying results. A similar approach was implemented in this work by integrating the theories established by Vogel and Peukert and further expanded on by Bonfils into the calibration of an appearance function.

The selection–breakage method was implemented to determine whether a particle sustained enough impact energy to suffer a breakage event [10]. By using the kinetic energy and size of particles, a probability distribution of particle breakage can be predicted (see Equations (2) and (3) where B, the breakage, is the appearance matrix. The index *i* represents the number of particle sizes available in the PSD dataset starting from 1, the largest size, to N being the smallest size. The model essentially treats particles that enter the crusher as a binary option: either they are crushed to breakage or the particle withstands the impact and retains its size. The particles that are chosen for breakage are then reduced in size according to the appearance function. The product is then added to the unbroken particles that pass through. This process can be seen in the flowchart in Figure 2. The number of impacts a particle, following its exit from the rotor, is assumed to be 1 based on the results of DEM simulations analyzing the collision energy spectrum. The majority of the collisions are under 10 J/kg and are assumed not to influence the impact breakage. These lower energy collisions are assumed to be after the initial collision with the rockbed in the crushing chamber.

$$P = \begin{bmatrix} BSf + (I - S)f \end{bmatrix} \tag{2}$$

$$S = \left[1 - \exp\left(f\_{mat}\mathbf{x}\_i^a k \left(E\_{CS} - \frac{E\_{min}^0}{\mathbf{x}\_i^a}\right)\right)\right] \tag{3}$$

**Figure 2.** The selection–breakage process was used to determine the product generated from the feed [10].

After the model had been calibrated for different sites, an effort to validate it was performed. The feed data and machine parameters for the VSI crushers were then put into the model, and the product was predicted. The speed of the particles at the rotor tip is transformed to kinetic energy, *ECS*, as seen in Equation (4).

$$E\_{CS} = \frac{v\_{rotor}^2}{2} \tag{4}$$

The predicted particle size distributions (PSDs) for the product and the experimental product were then compared to determine how much the model varied from the experiments. The model was based on crushers that employ a rockbox configuration. The input parameters and the overall process of running the model assumed that no prior material characterization had been performed. The approach can be seen in Figure 3.

**Figure 3.** The overall modelling process for a new site where the material has not been calibrated and then adequately mapped to the material database.

Once the material parameters have been obtained, the model can be run with alternative machine parameters to achieve different rotor tip speeds, as seen in Figure 4.

**Figure 4.** A model overview showing the parameters required to run the VSI crusher model at alternate speeds to predict the product PSD. The model requires a material database with *fmat*, *E0 min* and α values obtained from previous calibration efforts.

#### *2.2. Sampling Procedures*

Due to the enclosed design and destructive behavior of comminution machines, it is difficult to directly observe what occurs inside with changes in operational parameters. Experimental data that are measurable for a VSI crusher include the incoming material feed rate and particle size distribution, outgoing mass flow and particle size distribution, rotor tip speed, rotor radius and material properties. These parameters are often quantified on a site-by-site basis to establish the range of products that can be expected with the current setup of machines and the specific ore characteristics. For this model, datasets were collected from different sites throughout Sweden to obtain a comprehensive overview of the performance of the VSI crusher. The material from these sites was crushed with a Metso BARMAC 5100 SE VSI crusher (Metso, Helsinki, Finland), which has a diameter of 510 mm and speeds between 69 and 79 m/s. The conveyer before and after the VSI crusher was stopped to allow material sampling. A conveyor length of 1 m was isolated, and all the material was taken according to the SS-EN 933-1 standard. The intricacies of variation between the sites in terms of materials and machine parameters can create uncertainties in the analysis of its effects. Each source of variation, originating from either equipment, operator, material or site production level, contributes to a potential confounding of the results. The locations of the different sites can be seen in Figure 5 and the overall plant layout used in Figure 6.

Previous studies [13,15,17] have investigated the phenomena, conditions and effects of crushing particles at single-energy levels. Several of these studies have been conducted to analyze the inputs of specific energy levels on particles of different sizes. The results have often been reported in useful units such as J/kg to make comparisons between different particle sizes easier. More recent studies have shown simulations and observations of the different energies that the particles are exposed to while being expelled from the rotor and into the crushing chamber [12].

**Figure 5.** Location of the quarries tested in Sweden. Circles denote where the mobile VSI plant was set up, and the colored dots show where the material originated from.

**Figure 6.** The mobile VSI crusher plant layout used at the different quarries tested.

#### *2.3. Optimization Framework*

The purpose of this framework is to enable the rapid deployment of an accurate model without the need for exhaustive levels of sampling and ore characterization. The initial design space for the framework was intentionally kept small to reduce the computational load to run the optimization. Each round of sampling takes considerable time and effort on site for collection. With the data at hand and prior models to simulate certain aspects of impact crushing, a surrogate model was established for the appearance function and the selection breakage. Both of these sets of values were assessed for accuracy, and an error function was formulated (see Equation (5)), where *Pexp,High* and *Pexp,Low* are PSD data from experiments, *Psim,High* and *Psim,Low* are PSD data from simulations using the model and the *i* index is for the different particle sizes. The two datasets were run at two different rotor tip speeds, 69 and 79 m/s, which is indicated with the Low and High indexes, respectively. The framework for running the optimization utilizes the genetic algorithm for minimizing the error function.

$$\min\_{\{\beta\_{i,f\_{\text{mid}},F\_{\text{mid}}^{0},a}\}} E = \sum\_{i=1}^{n} \left( \left( P\_{i\_{Exp,High}} - P\_{i\_{Sim,High}} \right)^2 - \left( P\_{i\_{Exp,Low}} - P\_{i\_{Sim,Low}} \right)^2 \right) \tag{5}$$

The different bounds used in conjunction with the Genetic Algorithm can be seen in Table 2. A population size of 8000 was initiated with a set of 800 generations. The limits for the coefficients for the appearance generation polynomial were reached by initial testing to end up with a feasible yet large solution space. The values for *fmat*, *E0 min* and α were derived through the implications of their properties and then expanded on the spectrum at which these types of materials have been documented.

**Table 2.** Listing the lower and upper bounds of the 6 different variables used during the optimization routine.


The error function takes the result of the product estimation function using the same optimizer for two different speeds and compares each of them to their experimental data according to Equation (5). The genetic algorithm will then continue to the next generation until it reaches the tolerance limit and finds an optimum.

#### **3. Results**

Running calibrations on multiple datasets from multiple sites enabled the quantification of the material properties for the different geological rock types. The property values obtained from a wide array of sites all converge on values within the range that is to be expected of rock materials. In Table 3, the summary of the calibrated material properties per site is listed, where "Default Feed" indicates that only one sample was taken for the feed into the crusher, and the variation was assumed to be negligible, while the term "Individual Feeds" denotes that each rotor tip speed configuration run has matching feed sample data.

**Table 3.** Material data from all the sites used in the study, along with their locations and calibrated values.


The results for each dataset originating from one site can be visualized, as in Figure 7. The dashed lines represent the estimation of the model for the product, and the solid lines are the tested product curves for the corresponding speeds. Lastly, the black solid line indicates the feed or feeds in the dataset. Of the three rotor tip speeds, only two were used to train

the model, the highest and lowest speeds, while the middle speed was used to verify the response of the model.

**Figure 7.** Particle size distributions from one of the sites and the corresponding simulation results. Case data are from Tierp, as seen in Table 3.

By running all the datasets through the same optimization routine, similar results were obtained for each site, like the results achieved for Figure 7 but with different levels of magnitude relating to the feed used and material properties. While these properties include a multitude of factors that affect particle breakage and size progeny, some parameters are more prevalent than others and will vary less between sites. For instance, the feed size PSD has a large impact on the product, but the majority of sites surveyed used an 11.2-mm top size feed. Three factors were calibrated from each site with the use of experimental data.

*E0 min* is a characterization of the energy level required to achieve particle breakage, while the *fmat* property encompasses not only the material properties but also the properties arising from different particle shapes and *α* is a size dependency parameter. These parameters are dependent on many variables within the material but generate a perception of the strength of the different material types.

The models created and calibrated with these methods have created significantly accurate models within normal operative speeds. In each dataset, a speed, intermediate, was excluded from the training data to be used as verification of the model. Values of P20, P50 and P80 were extracted from the simulated product data and were compared to similar points interpolated from the experimental data. The values can be seen in Figure 8. These sets were compared to each other to assess the accuracy and range of the model. In large part, due to the optimization method, the values for P50 were the least differing, while the tail ends were often more spread away from the experimental data. Some variation in the experimental data is also assumed to be linked to the belt cut method of samples taken from the sites but is generally considered to be negligible.

To test the robustness of the model, PSD predictions of speeds exceeding the rotor tip speeds the model was trained with were made. In Figure 9, predictions were made for a tip speed of 85 m/s (15% higher kinetic energy than training data) and 65 m/s (12% less kinetic energy than training data), respectively.

**Figure 8.** A parity plot of the P20, P50 and P80 values from the experimental data compared to the simulated data using the developed models.

**Figure 9.** Particle size distributions from one of the sites and the corresponding simulation results for speeds exceeding the training data of the model. Note that the training speeds were still 69 and 79 m/s, not shown in the figure.

#### **4. Discussion**

The distribution of the datasets for the different sites and materials in Figure 10 displays an interesting pattern. Materials on the large scale tend to either end up with higher values for *E0 min* and lower values for *fmat* or, on the opposite side, with high values for *fmat* and low values for *E0 min*. Since *fmat* can be regarded as the resistance of the material against fracture, materials with lower values will be more pliable, and a higher level of particle size reduction can be expected when breakage does occur. *E0 min* is regarded as the energy threshold of a material, the level it can withstand before succumbing to comminution breakage. Materials with higher values tend to require higher levels of energy to break, ultimately requiring higher rotor tip speeds relative to materials with lower values.

The model that has been presented can be used to optimize existing plants. The machine parameters, the particle size distribution for the feed and product and the rock type are entered into the model, and the site operator can then change the flow rates, rotor speed or sizes to reach their desired production goals. The practical application of the model is to enable plant owners to estimate the feasibility of operating a VSI crusher on site and to see what different types of products they can expect to achieve and in what tonnages.

**Figure 10.** A scatter plot of the *fmat* and *E<sup>0</sup> min* values for different rock types and sites analyzed.

One possible shortcoming of the model and the theory it uses as a basis is that the properties of materials are summarized into two parameters and the size dependency into one parameter. While useful to make the optimization and helpful to characterize a material, it is possible that certain details or complexities are overlooked or disregarded due to this approach. For instance, the sphericity and flakiness of material going into the crusher compared to what goes out are often very different due to the abrasive nature of the crusher. This factor is included in the material parameters but not isolated as a separate measurement.

**Author Contributions:** Conceptualization, S.G. and E.H.; methodology, S.G.; software, S.G. and G.A.; writing—original draft preparation, S.G.; writing—review and editing, S.G., E.H. and M.E.; supervision, G.A., E.H. and M.E. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Swerock, Chalmers University of Technology, VINNOVA, SBUF and MinFo/Hesselman.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Acknowledgments:** The authors would also like to acknowledge the financial support from Swerock through the 3B project, Vinnova, SBUF and MinFo/Hesselman. The personnel at Swerock are also gratefully acknowledged for their support and insights. This work was carried out within the Sustainable Production Initiative and the Production Area of Advance at Chalmers University of Technology.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **Abbreviations**



#### **References**

