Grinding Kinetics Study of Tungsten Ore

: The European Commission (EC) maintains the consideration of tungsten as a critical raw material for the European industry, being the comminution stage of tungsten-bearing minerals an essential step in the tungsten concentration process. Comminution operations involve approximately 3–4% of worldwide energy consumption; therefore, grinding optimization should be a priority. In this study, the grinding behavior of tungsten ore from Barruecopardo Mine (Salamanca, Spain) is analyzed. A protocol based on Austin’s methodology and PBM is developed in order to study the inﬂuence of operational and geometallurgical variables on grinding kinetics. In addition to the kinetic parameters, the breakage probability ( S i ) and breakage function ( B ij ) is determined. The selection function was formulated for the Barruecopardo Mine with respect to the mill speed.


Introduction
The European Union has recently published the updated list of critical raw materials, in which tungsten (W) is included. This critical condition is defined by both the supply risk to the EU and the economic importance developed on the industrial value chains of the European Union [1].
Tungsten presents strategic applications on high strength alloys for machining tools, automotive and mobile phone sectors, among others [2]. Currently, there are several tungsten mines in Europe, some active and others on the exploration stage [3]. This is the case of the Barruecopardo mine in Salamanca (Spain), owned by Ormonde Mining PLC and currently administered by Saloro S.L., which is estimated to provide 11% of the non-Chinese global supply of tungsten [4]. The main minerals present in Barruecopardo are scheelite (CaWO 4 ) and wolframite ((Fe, Mn)WO 4 ), which constitute the ore. Arsenopyrite (FeAsS), pyrite (FeS 2 ), chalcopyrite (CuFeS2), and ilmenite (FeTiO 3 ) are also present as primary minerals of the gangue [5,6].
Regarding mineral benefit, the comminution stage represents 3-4% of the energy consumption worldwide and 40-70% of the energy consumed in a mineral processing plant [7,8]. In fact, ball-mill grinding is one of the most energy-consuming techniques. Therefore, setting the optimal values of the operational and mineralogical parameters both for the initial design and the process adaptation to ore variations [9].
Several researchers have investigated the influence exerted on kinetic conditions by operational parameters such as mill speed [10,11] and filling volume [12,13]. Other researchers devoted their work to study geometallurgical variables such as grain size, shape and roughness, specific surface area, orientation, hardness, fracture strength, feed particle size distribution, and mineralogy [14][15][16][17][18][19][20] using optical microscopy or more advanced techniques, such as Quantitative Microstructural Analysis (QMA) [21]. Consequently, a small improvement in machinery efficiency and an optimal design in the grinding system, taking into account the optimization of the above-mentioned parameters, would greatly cut down plant operational costs, impacting environmental issues and resource management [22,23].
This work aims to characterize the grinding kinetic behavior at a lab-scale of tungsten ore, due to its importance as a critical raw material to the EU, by determining the kinetic parameters following the Austin model. A second objective is to assess the influence of mill speed on the kinetic and geometallurgical parameters, being mill speed the more easily adjustable operational parameter at an industrial scale.

Theoretical Background
The population balance model (PBM) has been widely used in ball mills since its proposal by Austin [24]. This model is based on determining the particle size distribution grouped in size classes. A mass balance for the class i in a well-mixed grinding process is done by employing Equation (1), where a first-order kinetic fragmentation is assumed.
where w i (t) is the remnant mass fraction of particle size class i at grinding time t. The first term of the right-hand side is the mass fraction of the monosize i particles that break and thus no longer belong to that monosize, being S i the probability of fracture. The second term represents the contribution of all monosizes coarser than i that break to produce particles of monosize i. The fracture rate of a monosize material can be expressed by Equation (2): where S i is the probability of fracture or specific fracture rate, whose unit is t −1 . Assuming that S i does not change with time, the integral results in Equation (3).
where w i is the weight fraction of mineral feed into the mill having a size 1 for time t, and S i is the probability of fracture. According to the methodology proposed by Austin et al. [25], once S i values have been obtained through slope determination, they are plotted to the particle size, and Equation (4) is proposed to study the behavior of the probability of fracture S i .
where X i is the upper size limit of the interval in mm, and α T and α are model parameters that depend on the material properties and the grinding conditions. To find the second term of Equation (1), the fracture function b ii is defined. This function represents the particle fraction that belongs initially to interval j, after fracture falls in interval i. It is recommended to represent this value in cumulative form B ij , whose calculation is done with Equation (5).
That is, B ij is the sum of the mineral fractions finer than the upper limit of interval i as a result of the primary break of the size interval j. Austin et al. [25] showed that B ij values could be estimated from a size analysis of the products over short grinding times of an initial feed chiefly of size j through the method BII [25][26][27]. With the parameters of fracture function, B ij can be determined graphically through an empirical function like Equation (6).
where φ j , γ and β are parameters that depend on the material properties. The critical speed, Nc, is calculated using Equation (7).
where D is the mill diameter and d is the ball diameter [m]. Ball mill filling volume is calculated using Equation (8).
On the other side, Austin and Brame [28] calculated the sorting function α T in a general way through Equation (9).
where ν c is the mill speed expressed as a fraction of the critical speed. Finally, according to [29], the Froude number expresses the ratio of centrifugal acceleration to gravity acceleration at the perimeter of the mill chamber (Equation (10)). This number can be used to characterize the charge motion in the mill and the ball regime. Thus, in laboratory ball mills, it is recommended to define work conditions with centrifugal acceleration at the shell equalling 1 /2 of the acceleration due to gravity (F r = 0.5), corresponding to ν c = 70.7% [29].
A representative sample was characterized chemically through XRF, using a Bruker XRF S-4 Pioneer Advance, with sample preparation in a Claisse Perler, model M-4.

Calculation of Critical Speed and Initial Conditions for the Grinding Kinetic Tests
Mill critical speed was calculated using Equation (7). Table 1 shows the three milling speeds used in the tests. Table 1. Mill rotation speeds used in the grinding kinetic tests.

Mill Speed n [rpm] Froude Number
Grinding kinetic tests were run in a laboratory mill, 17.8 cm in diameter and 4500 cm 3 in capacity. Feed was of 900 cm 3 , and milling load consisted of 6.6 kg of steel balls with the following ball size distribution: 45 balls 19 mm in diameter, 23 balls 29.7 mm in diameter, and 17 balls 36.8 mm in diameter. Fill fraction was calculated using Equation (8). The feed consisted of samples of the 9 monosizes selected (5000/4000, 4000/3000, 3000/2000, 2000/1000, 1000/500, 500/250, 250/125, 125/75, 75/45 µm). Grinding times were 0.5, 1, 1.5, 3.5, 6.5, and 10.5 min. Each sample was dumped and, after performing a grain size analysis, W content was measured to assess the evolution of the W grade with respect to the kinetic parameters.
3.3. Determination of Fracture Probability (S i ), Fracture Function (B ij ), and the Kinetic Parameters (α T, α, φ j , γ, and β) 3.3.1. Fracture Probability (S i ) and Kinetic Parameters (α T, and α) After obtaining the oversize weights for each grinding time, and plotting the time function log (w i (t)/w i (0)) for each monosize, the equation for each curve, and consequently, the S i value, were calculated through the linear fitting using Equation (3). Then, S i , values for each monosize were plotted, and the parameters (α T and α) were calculated using Equation (4) for each mill speed condition to study the influence of this operational variable on the probability of fracture and on the kinetic parameters α T and α. The selection function α T , obtained through Equation (9) was calculated using an equation designed for this particular ore, as detailed in Section 4.2.

Determination of the Fracture Function (B ij ) and the Kinetic Parameters (φ j , γ and β)
These calculations were done for each feed monosize, i, at each mill speed condition and after the grinding time. Product particle size analysis was done by sieving with mesh sizes i to j−n. Weight of oversizes i and S j−n was determined, and fracture function values b ij and cumulative value B ij , were obtained using Equation (5). Finally, B ij data were plotted to relative size j/i for each monosize and mill speed conditions. The rest of the kinetic parameters were calculated through Equation (6).

Study of P 80 and the Ratio of Reduction R r
The evolution of some relevant parameters of the product was represented. These parameters were P 80 and the ratio of reduction, R r , after a 0.5 min grinding time at several mill speeds for each feed monosize. Likewise, the evolution of the probability of fracture, S i , and R r with mill speed was studied, for each feed monosize.

Chemical Characterisation of the Product
Milling products were chemically characterized with the equipment specified in Section 3.1.

Feed Characterzsation
The sample mechanically prepared displays a particle size distribution (PSD), as shown in Figure 1.
F 80 for the Barruecopardo sample is 1690 µm. Table 2 shows W contents obtained from a mineralogical study for each size interval.  F80 for the Barruecopardo sample is 1690 µ m. Table 2 shows W contents obtained from a mineralogical study for each size interval.  Table 2 reflect that the sample comes from a waste dump of the Barruecopardo mine, and therefore from an area with low W contents.
The probability of fracture (Si) and the kinetic parameters (α, αT), a size of 80% of product undersize (P80), and the reduction ratio were obtained from this procedure (Rr): Figure 2 shows the values (Si) plotted to the particle size according to the mill speed. W contents shown in Table 2 reflect that the sample comes from a waste dump of the Barruecopardo mine, and therefore from an area with low W contents.
The probability of fracture (S i ) and the kinetic parameters (α, α T ), a size of 80% of product undersize (P 80 ), and the reduction ratio were obtained from this procedure (R r ): Figure 2 shows the values (S i ) plotted to the particle size according to the mill speed.  Figure 2 shows that comparing the three speeds tested, the highest probability of fracture occurred at 70% critical speed, agreeing with Gupta amd Sharma [9], Gupta [10], and Herbst and Fuerstenau [30], who did their experiments at lab scale, with several materials, simulations, and finally upscaling. This also verified the Steiner [29] recommendation of a configuration with Fr value close to 0.5. For intermediate and coarse grain sizes, Figure 2. Evolution of the probability of fracture (S i ) with particle size, in relation to mill speed. Figure 2 shows that comparing the three speeds tested, the highest probability of fracture occurred at 70% critical speed, agreeing with Gupta amd Sharma [9], Gupta [10], and Herbst and Fuerstenau [30], who did their experiments at lab scale, with several materials, simulations, and finally upscaling. This also verified the Steiner [29] recommendation of a configuration with F r value close to 0.5. For intermediate and coarse grain sizes, a rather linear trend was observed using a logarithmic scale. This coincided with Deniz [22] for this grain size interval.
Noteworthy was a sharp break in slope that occurred at around 250 µm for the three tested speeds, coinciding with the size fractions more enriched in W, as shown in Table 2. Moreover, the results backed Gupta and Sharma's [9] statements, pointing to the probability of fracture S i being one of the mill operational conditions more influenced by the mineralogical variability among monosize fractions.
On the other side, the evolution of S i , P 80 , and R r for each feed monosize was studied at different mill speeds as summarized in Table 3. The reduction ratio was the result of dividing the d 80 in the feed size (F 80 ) by the d 80 in the product size (P 80 ). It must be highlighted in Table 3 that once the total grinding time was reached, a finer P 80 and a higher R r, were obtained at 70% working speed for most of the monosizes. This confirmed that the mill speed affected the grinding product [11,31].
The kinetic parameters (α, α T ) obtained after linearization of Equation (4) are summarized in Table 4 and the selection function is represented in Figure 3.  Table 4 shows that the parameter α is coherent with what Austin et al. [25] reported. These authors pointed out that it usually ranges between 0.5 and 1.5 and that it depends only on the mineral. On the other side, Table 4 and Figure 3 show that the value of the selection function α T does not vary significantly despite the increasing speed, because the mill geometry remains unchanged. Figure 4 displays α T values calculated from Equation (9), as proposed by Austin and Brame [28]. It can be seen that this expression does not fit this case. This led to a polynomial adjustment using the experiment values, as it was shown in Equation (11) that it fitted better with the studied sample.
Metals 2021, 11, 71 7 of 13 marized in Table 4 and the selection function is represented in Figure 3.   Table 4 shows that the parameter α is coherent with what Austin et al. [25] reported. These authors pointed out that it usually ranges between 0.5 and 1.5 and that it depends only on the mineral. On the other side, Table 4 and Figure 3 show that the value of the selection function αT does not vary significantly despite the increasing speed, because the mill geometry remains unchanged. Figure 4 displays αT values calculated from Equation (9), as proposed by Austin and Brame [28]. It can be seen that this expression does not fit this case. This led to a polynomial adjustment using the experiment values, as it was shown in Equation (11)

Fracture Function, Kinetic Parameters ( , γ and β)
The values of fracture function Bij in relation to the particle size for each monosize as the mill speed varied were determined through Equation (5)

Fracture Function, Kinetic Parameters (φ j , γ and β)
The values of fracture function B ij in relation to the particle size for each monosize as the mill speed varied were determined through Equation (5) and are shown in Figures 4-6 for 60%, 70%, and 80% mill speed, respectively.    Figure 5. Behavior of the fracture function (B ij ), with respect to particle size. 70% working speed Nc (Barruecopardo ore).
Metals 2021, 11, x FOR PEER REVIEW 9 of 14 Figure 6. Behavior of the fracture function (Bij), with respect to particle size. 80% working speed Nc (Barruecopardo ore).
The kinetic parameters of fracture function ( , γ, and β) are shown in Table 5. According to Austin et al. (1984), and β are parameters that depend on the material. Regarding γ and β, these authors propose that their values were usually in the range of 0.5-1.5 and 2.5-5.0, respectively. Figures 4-6 and Table 5 show that Bij depends on the feed grain size for parameters of 60, 70, and 80% of critical speed. The influence that mill speed exerts on Bij, can also be noticed by comparing the different monosizes: a greater difference existed for coarser sizes, whereas it was lesser for finer sizes. This was due to the fact that coarse sizes not only possessed a higher Si, but also were more prone to yield new finer particles (progeny). That meant that Bij depended on the feed particle size, as Ipek and Goktepe [32] observed, which was also influenced by the mill speed, and concurred with results by Deniz [22]. Nevertheless, this variation was not as significant as reported by Austin et al. The kinetic parameters of fracture function (φ j , γ, and β) are shown in Table 5. According to Austin et al. (1984), φ j and β are parameters that depend on the material. Regarding γ and β, these authors propose that their values were usually in the range of 0.5-1.5 and 2.5-5.0, respectively.  Table 5 show that B ij depends on the feed grain size for parameters of 60%, 70%, and 80% of critical speed. The influence that mill speed exerts on B ij, can also be noticed by comparing the different monosizes: a greater difference existed for coarser sizes, whereas it was lesser for finer sizes. This was due to the fact that coarse sizes not only possessed a higher S i , but also were more prone to yield new finer particles (progeny). That meant that B ij depended on the feed particle size, as Ipek and Goktepe [32] observed, which was also influenced by the mill speed, and concurred with results by Deniz [22]. Nevertheless, this variation was not as significant as reported by Austin et al. [25]. Table 5 shows parameter γ, which represents the fineness factor. In Figure 7, the γ values are plotted against mill speed for two feed particle sizes (4000 and 500 µm, respectively). Figure 7 depicts that γ values are influenced by both mill speed and feed particle size. Smaller γ values were related to coarse particles (4000 µm), which meant that more fine particles were generated. Conversely, finer particles (500 µm) generated a lesser proportion of fine particles, agreeing well with results by Ipek and Goktepe [32] and Austin et al. [25].  Table 5 shows parameter γ, which represents the fineness factor. In Figure 7, the γ values are plotted against mill speed for two feed particle sizes (4000 and 500 µ m, respectively).    particles were generated. Conversely, finer particles (500 µ m) generated a lesser proportion of fine particles, agreeing well with results by Ipek and Goktepe [32] and Austin et al. [25].     As can be seen in Figures 8-10, for most monosizes, W content increased with decreasing particle size. This could be partly explained because feed grain sizes around 250 µm already yielded higher W content, as shown in Table 2. Nevertheless, maximum values after grinding were 6-7 times higher than original values pointing undoubtedly to a differential grinding effect leading to W mineral grains mainly falling in the 250-125 µm interval. The accumulation of these W-enriched particles, which were more difficult to grind, supports the mineralogical explanation of the aforementioned breakage probability reduction at 250 µm size. In addition, this interval always presented S i values higher than the coarser and finer intervals. This would suggest that particles of this size would have a higher probability of fracture compared to the adjacent size intervals. Indeed, in all cases, a decrease in W content could be observed down to 75 µm size, followed by an increase of further finer particles. It must be highlighted that the highest W content was yielded for grinding at 70% Nc. This could be explained because, under these grinding conditions, the mill performance was more efficient due to a more adequate charge regime (F r close to 0.5), leading to better grinding kinetics.

Conclusions
The experimental work here presented and its further analysis permits to draw the following conclusions:

•
Austin's methodology has allowed us to characterize the kinetic behavior of tungsten ore by determining the kinetic parameters (α, α T , ϕ j , γ, β), concluding that the values S i and B ij do not vary significantly with time.

•
The probability of fracture, S i , is highest at 70% critical speed. Fracture function B ij , does not vary significantly with mill speed. Nevertheless, it is affected by the feed particle size, becoming higher for coarser sizes. • Equation (10) is proposed as the best α T fit, specifically for the studied ore.

•
Values of parameter γ are influenced by both mill speed and feed particle size. Coarser particles yield smaller γ values, i.e., they produce more fines, whereas finer particles produce lesser quantities. The parameters φ j and β depend on the features of the ore.

•
The chemical characterization and the study of the evolution of the relevant grinding parameters, such as P 80 and R r , to the grinding time have demonstrated that, first, the highest probability of fracture occurs at 70% of the critical speed and; second, the effect of differential grinding is evidenced between the W-bearing species and the gangue. The latter results in an increase of the W grade in the monosize 250/125 µm.