Pseudo-dynamic simulations applied to ball mill grinding circuit using population balance model and Monte Carlo Method

Abstract

Process simulations can be used to improve grinding circuit performance, which efficiently reduces operating costs. The population balance model (PBM) is widely accepted for grinding modeling because it can reproduce breakage events in tumbling mills, as described by Austin et al. (1984). In this study, a pseudo-dynamic model is introduced, integrating the PBM with the Monte Carlo Method to stochastically simulate variables in an industrial grinding circuit. This integrated approach enabled circuit simulations over a period of 2 hours, representing the operational variables as seen in historical data. Model validation showed a correlation of 0.74 in the product size distribution when comparing simulated outcomes with the original population.

Keywords:
comminution; modeling; stochastic process; PBM; MCM

1. Introduction

Grinding circuits are essential in mineral processing to increase the mineral surface or recover valuable metals from gangue. They have a direct influence on the efficiency of downstream metallurgical processes. However, grinding demands high energy, liners, and grinding media consumption (Herbst et al., 2003). Therefore, small improvements to this process can result in considerable economic returns.

Grinding simulations have been the subject of intensive research over the past decades, and several models have been developed to provide information on circuit design, optimization, and control (Fuerstenau et al., 2004; Shi & Xie, 2015; Le Roux et al., 2020). PBM is widely accepted for simulating the grinding process, where particles are transformed in terms of size and composition (King, 2012), and satisfactorily replicating particle breakage events in tumbling mills through selection and breakage functions (Austin et al., 1984).

Simulation is a tool that combines mathematical models to virtually reproduce the behavior of a system, phenomenon, or process. Moreover, it can be used to provide an “artificial” experience and to evaluate different operational conditions. Simulations have been used successfully in the mining sector at various stages of the production chain due to their ability to improve the performance of unit operations and circuits at a low investment cost. The term “process simulation,” either steady-state, dynamic, or pseudo-dynamic, in mineral processing implies predicting plant or circuit performance in terms of mass flow, solid concentration, particle size distribution, and other variables as functions of the intrinsic properties of a specific ore, equipment, and operational conditions (Napier-Munn et al., 1999). Unlike the steady-state model, which assumes the immutability of the process variables through a time interval, dynamic or pseudo-dynamic simulation is especially useful for predicting the behavior of a stochastic process, such as grinding.

A stochastic process is a family of random variables that represents the evolution of a system over time. This evolution can be modeled by mathematical models, such as the Monte Carlo Method (MCM). Moreover, the ability to analyze the unpredictable behavior of process variables is especially important in the context of modeling and simulation (Hodouin et al., 1988; Mishra, 2007; Ghaffari et al., 2012).

Recent advances in dynamic modeling for mineral processing simulations highlight the use of models created using Simulink in MATLAB (Liu & Spencer, 2004). These models have been employed to predict the dynamic behavior of plants (Khoshnam et al., 2015). Moreover, hybrid models, which blend dynamic and steady-state approaches, have been effectively employed for simulating mill-flotation circuits (Le Roux et al., 2020; Karelovic et al., 2016).

1.1 Population Balance Model

The PBM solution proposed by Reid (1965) for the size discretization of breakage events occurring in batch grinding processes is shown in Equation 1 (Fuerstenau et al., 2004):

(1) $\frac{d{m}_i\left(t\right)}{dt}=-S_i{m}_i\left(t\right)+{\displaystyle {\sum }_{j=1}^{i-1}{b}_{ij}{S}_j{m}_j\left(t\right),i=1,2,3\dots n}$

where m_i (t) represents the mass fraction of particles at size interval i and time t. S_i represents the breakage rate of particles at size interval i, and b_ij represents the fragment distribution after breakage, that is, the larger j material fraction that changes to size i after breaking. The breakage function (Equation 2) can be defined as the distribution of fragments appearing after the breakage of a single particle within a specific size interval (Austin et al., 1984):

(2) $B_{i,j}={\Phi }_j{\left(\frac{x_{i-1}}{x_j}\right)}^{\gamma }+\left(1-{\Phi }_j\right){\left(\frac{x_{i-1}}{x_j}\right)}^{\beta },i>j$

where B_ij is the breakage function and ɸ, γ, and ß are parameters to be adjusted to the experimental data. The selection function, or the specific rate of breakage, describes the rate at which particles within a particle size range are reduced and cross their lower size limit. Equation 3 presents a model that describes the selection function (Austin et al., 1984):

(3) $S_i=\left[a{\left(\frac{x_i}{x_0}\right)}^{\alpha }\right]/\left[1+{\left(\frac{x_i}{\mu }\right)}^{\Lambda }\right],\Lambda \ge 0$

where x_i is the particle size; x_o is the standardization size, which is equal to 1 mm; ɑ is a constant that depends on the grinding conditions; α is a constant that depends on the characteristics of the material (positive number normally between 0.5 and 1.5); µ is the critical size at which the selection function achieves maximum value; and Λ is a positive number indicating the speed of decrease of the S_i with increase in size i. The S_i and B_ij parameters can be obtained through bench tests and/or industrial sampling. Herbst & Fuerstenau (1980) proposed that S depends on the consumed power and material mass in the mill. Thus, it is possible to obtain a specific selection function (S_i^E) using Equation 4:

(4) $S_i=S_i^E\left(\frac{P_{net}}{H}\right)$

where H is the material mass in the mill (Holdup), and P_net is the net consumed power.

1.2 Classification model

Sepúlveda & Gutiérrez (1986) proposed a modified version of Plitt's (1976) hydrocyclones classification model as described below. The cyclone feed pressure (C_p) is represented by Equation 5:

(5) $C_p=a_1\left(\frac{Q^{1.46}exp\left(-7.63\varphi +10.79{\varphi }^2\right)}{{\left(Dc\right)}^{0.20}{C}_h^{0.15}{\left(Di\right)}^{0.51}{\left(Do\right)}^{1.65}{\left(Du\right)}^{0.63}}\right)$

where a₁ is a constant dependent on the application; Dc is the cyclone diameter (in); Di is the cyclone inlet diameter (in); Do is the vortex (in); Du is the apex (in); C_h is the cyclone height (in); Q is the feed flow ( m³/h ) ; ɸ is the feed volume solids fraction. The adjusted cutting size, also known as d_50c, refers to the particle size that is evenly divided between the underflow and overflow streams, based on the adjusted efficiency curve. This parameter can be calculated using the following Equation 6:

(6) $d_{50c}=a_2\left(\frac{{\left(Dc\right)}^{0.44}{\left(Di\right)}^{0.58}{\left(Do\right)}^{1.98}exp\left(11.12\varphi \right)}{{\left(Du\right)}^{0.80}{C}_h^{0.37}{Q}^{0.44}{\left({\rho }_s-1\right)}^{0.5}}\right)$

where a₂ is a constant dependent on the application. The Equation 7 represents the flow partition, denoted by F_p, which is defined as the ratio of the underflow and overflow in a circulating load:

(7) $F_p=a_3\frac{C_h^{1.19}{\left(Du/Do\right)}^{2.64}exp\left(-4.33\varphi +8.77{\varphi }^2\right)}{C_p^{0.54}{\left(Dc\right)}^{0.38}}$

where a₃ is a constant dependent on the application. The efficiency of adjusted classification (E) and Plitt's (1976) m parameter can be calculated using the following Equations 8 and 9, respectively:

(8) $E_i^c=1-exp\left[-0.693{\left(\frac{d_i}{d_{50c}}\right)}^m\right]$

(9) $m=exp\left[a_4-1.58\frac{F_p}{\left(F_p+1\right)}\right]{\left[\frac{{\left(Dc\right)}^2{C}_h}{Q}\right]}^{0.15}$

where a₄ is a constant dependent on the application. The Equation 10 can be used to express the slurry short-circuit:

(10) $B_{pf}=\lambda B_{pw}$

where B_pf is the underflow short-circuit; B_pw is the water short-circuit; λ is constant depending on the application. The water short-circuit B_pw can be expressed by Equation 11:

(11) $B_{pw}=\left(\frac{F_p}{\left(F_p+1\right)}-\phi R_s^c\right)/\left(1-\phi \left[1-\lambda \left(1-R_s^c\right)\right]\right)$

where R_s^c is the total solids recovery (hypothetical) expressed by weight.

1.3 Monte Carlo Method

A stochastic process can be defined as the unpredictable evolution of a random time phenomenon. Despite the known initial conditions, there are many possible directions for the process evolution. In contrast, the deterministic process has only one solution. The stochastic process is represented by a set of time-indexed random variables. Stochastic processes are simulated using MCM, which are based on the generation of pseudorandom numbers (Gentle, 2003).

Owing to various factors, such as circuit characteristics, control network, instrument accuracy, and response time, the mill process-controlled variable path oscillates unpredictably over time. The logic controller operates in the system with the aim of maintaining the mean of the control variable in the vicinity of a set point. Therefore, the Mean Regression Method (MRM) was chosen for simulation. This model is included in the MCM and can be used to randomly describe a stochastic process (Dias, 2005). The following formalization (Equation 12) of the MRM, that is,

the arithmetic Ornstein–Uhlenbeck process for stochastic variable X(t), is in accordance with the sequence described by Dias (2005):

(12) $dX=\eta \left(\overline{X}-X\right)dt+\sigma dz,$

where η is the regression parameter, X̅ is the mean, and σ is volatility. Equation 12 shows that a regression force acts on variable X, pulling it toward an equilibrium level (mean). The parameter η represents the speed of the regression process. The stochastic differential equation has the following solution (Equation 13).

(13) $X\left(T\right)=X\left(0\right)e^{-{\eta }^T}+\left(1-e^{-{\eta }^T}\right)\overline{X}+\sigma e^{-{\eta }^T}{\displaystyle {\int }_0^T{e}^{\eta t}dz\left(t\right)}$

With the functions of mathematical expectation (Equation 14) and variance (Equation 15), the variable X(T) has a normal distribution. Equation 16 shows the relationship between the mean reverting speed and half-life (h) of the regression process.

(14) $E\left[X\left(T\right)\right]=X\left(0\right)e^{-{\eta }^T}-\overline{x}\left(1-e^{-{\eta }^T}\right)$

(15) $Var[X\left(T\right)=\left(1-e^{-{\eta }^T}\right)\frac{{\sigma }^2}{2\eta }$

(16) $h=\left(ln\left(2\right)\right)/\eta$

The parameter h can be defined as the expected time for the stochastic variable X to reach half way toward the equilibrium level X̅. Equation 17, discretized in time, enables the simulation.

(17) $X_t=X_{t-1}{e}^{-\eta \Delta t}+\overline{X}\left(1-e^{-\eta \Delta t}\right)+\sigma \sqrt{\left(1-exp\left(-2\eta \Delta t\right)/\left(2\eta \right)\right)}N\left(0,1\right)$

With a normal distribution and mean equal to zero, Equation 17 is treated as the sum of mathematical expectation and a random term. The model identifies the long-term equilibrium level for the random variable as a very important reference, such that the variable P follows a mean regression towards an equilibrium level given by Equation 18:

(18) $\overline{X}=\ln \left(\overline{P}\right)$

where P = e(X) represents the long-term equilibrium constant of a random variable. The mean regression model defines the simulated mean as E [P(T)] = e{x[X(T)]} , de monstrating that the ratio between variables X and P it is possible to obtain the following (Equations 19 and 20) expected value for the random variable at time T during the simulation:

(19) $E\left[P\left(T\right)\right]=exp\left\{X\left(0\right)e^{-\eta T}+\overline{X}\left(1-e^{-\eta T}\right)\right\}$

(20) $P\left(T\right)=exp\left\{X\left(t\right)-0.5Var\left[X\left(t\right)\right]\right\}$

where Var[X(T)] is defined in Equation 15. By combining Equations 15, 19, and 20, it is possible to simulate the paths of the random variable P, following a mean regression process. This allows the direct simulation of the real process for P (t) using Equation 21 (Dias, 2005).

(21) $P\left(T\right)=exp\left\{\left[ln\left[P\left(t-1\right)\right]exp\left[-\eta \Delta t\right]\right]+\left[ln\left(P\right)\left(1-\text{exp}\left[-\eta \Delta t\right]\right)\right]-\left[\left(1-\exp \left[-2\eta \text{t}\right]\right){\sigma }^2/4\eta \right]+\sigma \sqrt{(\left(\left(1-\exp \left[-2\eta \Delta \text{t}\right]/2\eta \right)\right)}N\left(0,1\right)\right\}$

1.4 Objective

The purpose of this study is to develop a tool that simulates continuous grinding processes using well-known mathematical models. Incorporating both the PBM and the MCM, the pseudo-dynamic model aims to reproduce operational variables similar to the original data. To validate and demonstrate the model's efficacy, a ball mill circuit from a gold ore plant was selected for modeling and simulation.

2. Materials and methods

2.1 Grinding data and ore sampling

The gold plant process can be summarized into three stages of crushing and one stage of grinding to obtain the required product sizes for the downstream leaching process. The mill has an internal grate discharge and is in a direct closed circuit classified by cyclones (Table 1).

Table 1
Grinding circuit parameters.

Ball mill / Ore
Effective Diameter/Length (feet)	12/17	Interstitial slurry (%)	100
Balls charge filling (%)	40	Lift angle (°)	31.9
Critical speed (%)	73	Installed power (kW)	1150
Make-up ball size (mm)	63	Balls density (t/m³)	7.75
Ore density (t/m³)	2.75	New feed setpoint (t/h)	115
Cyclone cluster
Vortex (mm)	140	Number in operation/total	3/5
Apex (mm)	100	Diameter/Height/Inlet (mm)	500/1090/195

A sampling campaign was performed to obtain data for the mass balance and ore for laboratory tests. New feed, mill discharge, cyclone overflow and underflow were collected every 5 min for 2 h. Throughout this time, the plant maintained a stable operation with a consistent new feed rate of 115t/h, and control setpoints remained unchanged. At the end of the sampling period, a large sample of ore feed (≈ 500kg) was collected directly from the feed belt. Samples were weighed, dried, and divided into aliquots for screening and grinding batch tests in mono size fractions (Austin et al., 1984). All data for the grinding process variables were collected during the sampling time from the plant PLC (programming logic controller).

2.2 Grinding tests in mono-size fractions

Batch grinding tests with mono size fractions, are, in general, short-time tests with the purpose of preventing regrind. The material to be tested has a pre-classified predominant size, is sieved between mesh ranges, and has a √2 ratio. The material was ground at different time intervals, and a particle size analysis was performed at each step, including the original “zero time” sample. The tests were conducted in a mill of 254 mm diameter and length, with 70% critical speed, 20% ball filling, 25.4 mm single size balls and 50% material filling (Austin et al., 1984).

2.3 Griding Modeling

For the pseudo-dynamic simulation, an algorithm was developed in Visual Basic for Application (VBA) - Excel to calculate the model of N individual steady-state simulations between short fractions of the total simulated time and to obtain the expected outputs. To identify the optimal values for the circuit's operating parameters, a search algorithm was employed during the simulation.

To estimate the circuit's mass balance, the simulation's time interval was set equal to the residence time required for a specific fraction of the total mass introduced into the grinding circuit.

This fraction ensures the complete filling of the effective ball circuit volume. The developed model considers an accumulation or a decrease in the mass in the circuit over time, in contrast to one of the basic premises of steady-state simulation. The mass balance can be summarized using Equation 22.

(22) $\text{INPUT-OUTPUT}=\text{ACCUMULATION}$

The control variables for the circuit, specifically new feed and cyclone feed density, were selected for stochastic simulation over time using MRM and subsequently utilized as inputs to the grinding model. To test the applicability of MRM, it was essential to ensure that these variables followed a normal density function. Therefore, an Anderson-

Darling test for normality was conducted, with a significance level set at p > 0.05 (Umaporn, 2011).

Following is the description of the algorithm to perform pseudo-dynamic model calculations for the grinding circuit:

(a) The MRM statistical parameters are entered to calculate the throughput rate and cyclone feed density variables: mean to be maintained (setpoint), volatility, and half-life time (Dias, 2005).

(b) The simulation period is set to approximately 2h (the same as the sampling campaign).

(c) N interactions loops (N = 1, 2 ,3,4, ...) are initiated: I. First, the circuit residence and simulation times are calculated:

a. Simulation time evolution is calculated using Equation 23:

(23) $\{\begin{array}{c}\begin{array}{cc}t\left(N\right)=\text{Initial}\text{start}\text{time,}& N=1\end{array}\\ \begin{array}{cc}t\left(N\right)=\left[{\displaystyle {\sum }_1^N\tau \left(N\right)}\right]+t\left(N=1\right),& N>1\end{array}\end{array}$

b. The residence time in the circuit can be calculated as a function of the total mill-fed mass rate (new through-put + underflow) using Equation 24: where τ(N) represents the residence time (in minutes) at cycle N, F(N) denotes the total feed mass rate during cycle N (in t/h), P(N) is the total slurry rate in cycle N (in m³/h), VC is the occupied volume of the mill discharge box (m³), and H signifies the mill hold-up (in tons).

(24) $\tau \left(N\right)=\left(\frac{H}{F\left(N\right)}+\frac{VC}{P\left(N\right)}\right)60$

II. The MRM is used for a random simulation of the selected process variables (new feed and cyclone feed density), or the values are searched in the database, depending on the simulations described in section 3.2.

III. The minimization routine of the objective function (Equation 26) at time t(N) complies with the following mass balance Equation 25:

(25) $\frac{NF\left(N\right)+UF\left(N\right)}{\tau \left(N\right)}-\frac{OF\left(N\right)+UF\left(N\right)}{\tau \left(N\right)}=\pm CLV$

(26) $F_{obj}=\sqrt{{\left(\frac{NF\left(N\right)+UF\left(N\right)}{\tau \left(N\right)}-\frac{OF\left(N\right)+UF\left(N\right)}{\tau \left(N\right)}\pm CLV\right)}^2}$

where NF(N) is the dry mass of the circuit feed in the mill in cycle N, OF(N) is the dry mass in the cyclone overflow in cycle N, UF(N) is the dry mass in the cyclone underflow in cycle N. CLV represents the difference in dry mass of the circulating load between consecutive cycles (N and N-1), indicating either an accumulation or a decrease.

A search algorithm was applied to minimize the objective function (Equation 26). This algorithm adjusts input model variables, including the percentage of solids in mill discharge, cyclone underflow and overflow, and water addition, aiming to identify the system's optimal solution.

To calculate the net consumed power of the model was used the Equation 27 proposed by Hogg & Fuerstenau (1972):

(27) $P_{net}=0.238D^{3.3}\left(L/D\right)N_c{\rho }_{ap}\left(J-1.065J^2\sin \delta \right)$

where D and L are the mill diameter (ft) and length (ft), respectively, N_c is the tumbling speed, expressed as a fraction of the critical speed, ρ_ap is apparent load density, (t /m³) is the apparent volumetric fractional mill filling, and δ represents charge lifting angle.

3. Results and discussion

3.1 Model parameters estimation

The breakage function was estimated using mono-size grinding tests on the sampled ore, employing the BII method developed by Austin et al. (1984). This technique directly calculates the breakage function values and includes corrections for fractures caused by undesirable material regrinding. Figure 1 shows the distribution of fragments from the normalized breakage of the 3 mono size tests. Values are plotted in cumulative form, B_ij, versus size as a fraction of the breaking size, d_i /d_j (Austin et al., 1984).

Figure 1
Cumulative progeny fragment distributions from breakage of 3 mono size grinding tests and estimating the breakage function estimation.

Parameter calculations of the classification and selection functions were performed using Moly-Cop Tools™ . This simulator is provided with parameter estimation routines for analyzing data produced in the laboratory and/or industrial sampling. The grinding model parameters are listed in Table 2 and the sample’s particle size distribution and the fitted results are plotted in Figure 2.

Table 2
Classification, selection, and breakage parameters.

Classification		Selection		Breakage
a 1	6.54	a	8.07×10–5	Φ	0.62
a 2	1.31	α	1.5	γ	0.59
a 3	12.93	Λ	2.98	β	5.13
a4	0.72	μ	1331	@
λ	1.31	@

Figure 2
Particle size distribution: sampling and fitting results.

The MRM fitting (constants presented in Table 3) was obtained by statistical analysis of the original population, empirical modeling through a set of inputs, and observation of random outputs.

Table 3
Regression model constants.

Random variables	Mean	Volatility	Regression parameter
Mill feed rate (t/h)	115	1.90	0.9
Cyclones feed density (t/m³)	1.581	3.0 × 10–3	1.9

3.2 Pseudo-dynamic model validation

The pseudo-dynamic model and MRM were validated by comparing the behavior of the simulated results with the data measured by instruments/sensors installed in the processing plant, at the same time of the sampling. An evaluation was performed to determine whether the proposed models could produce data similar to the original data. First, to validate the pseudo-dynamic grinding model, the mill control variables of the cyclone feed density and mill feed rate in the S1 simulation (Table 4) were set to the same value by the algorithm as those measured by online instruments. The variation in size distribution results sampled from the mill feed was also used as input for the model (Figure 5-a).

Table 4
Summary of the pseudo-dynamic simulations (mean results).

Unit		Measured(A)	Pseudo-dynamic simulations
			S1	MCM model
				S2
Process variables
Throughput	t/h	114.95	115.0	115.1
Standard deviation	@	1.88	1.93	2.02
Coefficient of variation	%	1.63	1.68	1.76
Cyclones feed density	t/m³	1.5806	1.581	1.582
Standard deviation (10–3)	@	2.91	2.94	3.21
Coefficient of variation	%	0.190	0.186	0,20
Mill product (passing 106 µm)	%	82.04(B)	81.87	81.48
Coefficient of variation	%	1.61	1.59	1.42
Circulating load	(%)	3.55(C)	3.27	3.21
Coefficient of variation	%	-	5.51	4.42
Nominal power	kW	1057	1056	1057
Coefficient of variation	%	0.23	0.45	0.43
Specific energy consumption	kWh/t	9.18	9.19	9.19

The behavior replication test was applied to one output variable of the model - product size distribution - comparing simulated results with measured values (Figure 3). Disregarding the discrepant values of the data, good correlation and R-squared results were observed for the product size distribution, with values of 0.7376 and 0.5441, respectively. To identify outliers, the data points that significantly influenced the correlation coefficient were ranked, and the IQR (Interquartile Range) method was employed to calculate the lower and upper limits.

Figure 3
Correlation between measured and simulated (S1) mill product size (passing 106 µm).

In the second part of the validation, MRM was used in the S2 simulation (Table 4) to randomly calculate the evolution of the selected process variables (cyclone feed density and mill feed rate) over time, using the statistical model parameters (Table 3): mean to be maintained (setpoint), volatility, and regression parameter.

Statistical analysis was performed using histograms and normal distributions of the database and randomly estimated populations for mill throughput (Figure 4 – a and b) and cyclone feed density (Figure 4 – c and d). The analysis was performed to verify the MRM capacity to reproduce the behavior of the original data. The S2 simulation provided the mean and standard deviation results that were very close to the real data for the analyzed variables, and both populations presented similar normal distribution curves. The population of these measured and random variables passed the Anderson–Darling normality test (p-value > 0.05), which is a requirement for MRM (Figure 4).

Figure 4
Histogram and normal distribution for the measured data (a and c) and values estimated using MRM (b and d) in the S1 simulation for mill throughput and cyclone feed density. Normality test results indicated (p-value > 0.05).

The output behavior of the pseudodynamic model, as represented by the S1 and S2 simulations, was plotted alongside the curves of the measured data from the grinding circuit (Figure 5). The results demonstrated that the model could identify the time-evolution behavior of process variables. Specifically, the product size of the S1 simulation exhibited similar trends, periods, and amplitudes when compared to the data (Figure 5 - b). The S2 simulation results indicate that the MRM successfully generated a population mirroring the behavior of the original data. This capability paves the way for testing alternative setpoints for cyclone feed density and new mill feeds in subsequent studies, all while preserving the intricate variability of the circuit.

Figure 5
Time evolution comparison of grinding circuit variables: S1 simulation vs. S2 (MCM) simulation vs. instrument measurements.

(A) Plant database. (B) The mean of the online measurement of particle size in individual cyclone overflow. (C) Results from mass balance derived from the sampling data.

4. Conclusion

In short, the study developed and tested a pseudo-dynamic model to accurately simulate the grinding process. Highlights include the fact that the proposed models underwent behavior replication tests, which proved their efficacy in simulating the grinding process. Additionally, the research presented a pseudo-dynamic model that combined PBM and MCM. This approach enabled the simulation of time-indexed random variables, resulting in outcomes with dispersion parameters that were comparable to the measured data. For future studies, the recommendation is to employ the developed tool for optimizing operational parameters and formulating effective control strategies for the grinding circuit. This approach allows the evaluation of the effects of deviations in grinding variables on downstream processes.

References

• AUSTIN, L. G.; KLIMPEL, R. R.; LUCKIE, P. T. Process engineering of size reduction: ball milling. New York: Society of Mining Engineers of the AIME, 1984. 561 p.
• DIAS, M. A. G. Hybrid real options with petroleum applications. 2005. Tese (Doutorado em Engenharia Industrial) - Pontifícia Universidade Católica do Rio de Janeiro, Rio de Janeiro, 2005.
• FUERSTENAU, D. W.; DE, A.; KAPUR, P. C. Linear and nonlinear particle breakage processes in comminution systems. International Journal of Mineral Processing, v. 74, p. 317-327, 2004.
• GENTLE, J. E. Random number generation and Monte Carlo methods. 2. ed. New York: Springer, 2003. 377 p.
• GHAFFARI, A.; HAYATI, M.; SHEKHOLESLAMI, A. Probability and sensitivity analysis in flotation circuit of bama lead and zinc processing plant using Monte Carlo simulation method. Mineral Processing and Extractive Metallurgy Review, v. 33, n. 6, p. 416-426, 2012.
• HERBST, J. A.; FUERSTENAU, D. W. Scale-up procedure for continuous grinding mill design using population balance models. International Journal of Mineral Processing, v. 7, p. 1-31, 1980.
• HERBST, J. A.; LO, Y. C.; FLINTOFF, B. Size reduction and liberation. In: FUERSTENAU, M. C.; HAN, K. N. Principles of mineral processing. Littleton, Colorado: Society for Mining, Metallurgy, and Exploration, Inc., 2003. cap. 3, p. 61-118.
• HODOUIN, D.; DUBÉ, Y.; LANTHIER, R. Stochastic simulation of fltering and control strategies for grinding circuits. International Journal of Mineral Processing, v. 22, p. 261-274, 1988.
• HOGG, R.; FUERSTENAU, D. W. Power relations for tumbling mills. Trans. SME-AIME, v. 252, p. 418-432, 1972.
• KARELOVIC, P.; PUTZ, E.; CIPRIANO, A. Dynamic hybrid modeling and simulation of grinding–flotation circuits for the development of control strategies. Minerals Engineering, v. 93, p. 65-79, 2016.
• KHOSHNAM, F.; KHALESI, M. R.; DARBAN, A. K.; ZAREI , M. J. Development of a dynamic population balance plant simulator for mineral processing circuits. International Journal of Mining and Geo-Engineering, v. 49, p. 143-153, 2015.
• KING, R. P. Modeling and simulation of mineral processing systems. 2. ed. Englewood, USA: Society for Mining, Metallurgy and Exploration, 2012. 466 p.
• LE ROUX, J. D.; STEINBOECK, A.; KUGI, I.; CRAIG, I. A. Steady-state and dynamic simulation of a grinding mill using grind curves. Minerals Engineering, v. 152, n. 106208, 2020.
• LIU, Y.; SPENCER, S. Dynamic simulation of grinding circuits. Minerals Engineering, v. 17, p. 1189-1198, 2004.
• MISHRA, B. K. Monte Carlo method for the analysis of particle breakage. In: SALMAN, A. D. et al. Handbook of powder technology: particle breakage. Amsterdam, The Netherlands: Elsevier, 2007. v. 12, cap. 15, p. 637-660.
• NAPIER-MUNN, T. J.; MORRELL, S.; MORRISION, R. D.; KOJOVIC, T. Mineral comminution circuits: their operation and optimization Indooroopilly, Queensland, Australia: Julius Kruttschnitt Mineral Research Center, 1999. 413 p.
• PLITT, L. R. A mathematical model for the hydrocyclone classifer. CIM Bulletin, v. 69, p. 114-123, 1976.
• REID, K. J. A solution to the batch grinding equation. Chemical Engineering Science, v. 20, p. 953-963, 1965.
• ROUX, J. D. L.; STEINBOECK, A.; KUGI, A.; CRAIG, I. K. Steady-state and dynamic simulation of a grinding mill using grind curves. Minerals Engineering, v. 152, p. 106-208, 2020.
• SEPÚLVEDA, J. E.; GUTIÉRREZ, L. R. Dimensionamiento y optimización de plantas concentradoras mediante técnicas de modelación matemática. Santiago: Centro de Investigación Minera y Metalúrgica, 1986.
• SHI, F.; XIE, W. A. Specific energy-based size reduction model for batch grinding ball mill. Minerals Engineering, v. 70, p. 130-140, 2015.
• UMAPORN, C. Efficiency comparisons of normality test using statistical packages. Thammasat International Journal of Science and Technology, Bangkok, v. 16, n. 3, p. 9-25, 2011.

Publication Dates

• Publication in this collection
  19 Apr 2024
• Date of issue
  Apr-Jun 2024

History

• Received
  29 Mar 2023
• Accepted
  29 Nov 2023