Analysis of the Dynamics of Rougher Cells on the Basis of Phenomenological Models and Discrete Event Simulation Framework

Abstract

Due to the increase in the amount of copper sulphide minerals processed through concentration processes and the need to improve the efficiency of these production processes, the development of theoretical models is making an important contribution to generating a better understanding of their dynamics, making it possible to identify the optimal conditions for the recovery of minerals, the impact of the independent variables in the responses, and the sensitivity of the recovery to variations in both the input variables and the operational parameters. This paper proposes a method for modeling, sensitizing, and optimizing the mineral recovery in rougher cells using a discrete event simulation (DES) framework and the fitting of analytical models on the basis of operational data from a concentration pilot plant. A sensitivity analysis was performed for low, medium, and high levels of the operative variables and/or parameters. The outcomes of the modeling indicate that the optimum mineral recovery is reached at medium levels of the flow rate of gas, bubble size, turbulence dissipation rate, surface tension, Reynolds number of bubble, bubble–particle contact angle, superficial gas velocity and gas hold-up in the froth zone. Additionally, the optimal response is reached at maximum levels of particle size and density and at minimum levels of bubble speed, fluid kinematic viscosity and fluid density in the sampled range. Finally, the recovery has an asymptotic behavior over time; however, the optimum recovery depends on an economic analysis, examining the marginalization of the response over time in an operational context.

1. Introduction

Globally, approximately 20 million tons of copper were produced from mines, along with 25 million tons from refinery production, in 2020 [1]. Chile is the largest copper producer worldwide, with a 27.9% share [2,3,4] of production and 29% of the reserves of this commodity [5]. Copper oxides processed by hydrometallurgy are becoming increasingly scarce in Chile (from 30.8% in 2015 to 12% in 2027), while copper sulphides are being found in more significant amounts [6,7,8]. The National Service of Geology and Mining [2] (SERNAGEOMIN, by his acronym in Spanish) has indicated that 39.2% of copper production takes place through hydrometallurgical processes, but most of the production (60.8%) takes place through flotation processes. A report by the Chilean Copper Commission [6] (COCHILCO, by his acronym in Spanish) proposed steadily increasing the stock of copper concentrates in Chile, and it is estimated that growth will double between 2014 and 2026, accounting for 88% of domestic mining production. This means going from 3.9 to 5.4 million tons of concentrate. Therefore, the Chilean copper industry has indicated the concentration process to be the future production method; however, this process involves the generation of environmental liabilities. It is estimated that for every ton of Cu produced by flotation processes, 151 tons of tailings are generated [9,10,11,12]; there are currently 92 mining worksites that have been defined as environmental liabilities [13]. As part of sulphide mineral processing, the flotation, smelting and electrorefining processes are disaggregated. Flotation is defined as the process of concentration of minerals in which useful ore particles are separated from sterile ones or gangues, wherein by means of physical–chemical treatment, it is possible to modify the surface tension of the particles, causing them to adhere to bubbles of air, to later be recovered in the froth phase. Flotation is a key process in copper recovery from sulphide minerals, and determines the final recovery of ore, meaning that study of the dynamics of this process has the potential to improve performance and contribute to the generation of value in the mineral processing industry.

Essential variables that influence the performance of flotation operations include the size and number of bubbles, the particle sizes, and hydrophobic conditions [14]. Copper flotation operations have attracted special interest in the scientific community. Some works include the adjustment of phenomenological models through mathematical programming [15], conceptual and analytical designs of concentration circuits [16,17,18], sensitivity analysis [19], the performance evaluation of flotation banks and circuit designs using machine learning techniques [20,21,22,23], studies of the dynamic behavior of particles and bubbles with respect to recovery efficiency [24,25,26,27,28,29,30], the determination of the influence of factors on recovery [31], the impact of gas dispersion measures [32], and effect of clay minerals [33]. Foremost among the simulation and optimization methodologies are the use of mixed-integer linear and nonlinear programming algorithms [34], the optimization of recovery performance on the basis of the depth of the froth [35], methodologies for optimizing grinding and flotation using automated mineralogy data [36], reagent dose control methods based on optimal bubble size distribution [37], and optimization methodologies applied to frothing recovery, superstructures, mathematical models, and algorithms. This also takes into consideration the epistemic and stochastic uncertainties in both modeling and process performance [38,39].

The optimization models mentioned above are successful at determining the optimal conditions of the analytical or kinetic models that represent the flotation process or some of the threads that comprise it. However, these do not offer a comprehensive approach for rougher cells, indicating both the impact of the variability of the feed variables and the process control parameters and their respective optimal values. Therefore, in this work, we present an analysis of the modeling, simulation, and optimization of rougher cells in a flotation circuit on the basis of a discrete event simulation (DES) framework, representing an interesting opportunity to study the dynamics of the process and to identify target values of operational parameters in order to optimize the mineral recovery. The simulated cell is integrated with analytical models that represent the dynamics of the flotation process with rougher cells, previously validated in the literature, fitted to operational data in a high-level design for rougher cells in a pilot plant. Then, the objective of this work is to evaluate the ability of the tool (analytical expressions and simulation using DES) to represent the dynamics of a single rougher cell, simulate its operation, and optimize the response.

This article describes the materials and methods, details the simulation framework that was used (discrete event simulation; DES), the high-level superstructure of a rougher cell, the analytical formalization of the process, and the superstructure adjustment in the Arena simulation software. The Results and Discussion establishes the domain of the operating variables, the comparative analysis of the simulations for low, medium, and high levels of the operative variables, and the results of the optimization of the simulated model. Finally, the Conclusions and Future Works section presents the main findings of the work development and possible future research lines.

2. Materials and Methods

There are several mineral flotation processes, mainly varying with respect to their physical and chemical aspects, such as metal solubility, solution kinetics, reagent consumption, etc. In this section, a description of the discrete event simulation framework (DES), the mathematical modeling of a rougher cell and the explanation of the superstructure adjustment in the Arena simulation software is provided [40].

2.1. Discrete Event Simulation

Discrete event simulation comprises a collection of techniques that, when applied to the study of a dynamic events, generates sequences called sampling routes that characterize their behavior [41]. One or more phenomena of interest change their values at discrete points over time, making it possible to model all those systems in which events change state at time-spaced moments [42]. Although seemingly simple, such systems have the potential to model complex phenomena [43]. The main feature of a DES is that the system consists of a sequence of events that occur in time and changes to the system state occur at times ${\mathrm{t}}_1,{\mathrm{t}}_1,\dots ,{\mathrm{t}}_{\mathrm{n}}$ [44]. When generating a DES program, an approximation of the operation can be developed, describing the sequences of events and activities performed by the entities and how they are exposed to modifications when completing each of the phases that make up the simulated model [43].

2.2. Superstructure

The superstructures proposed in the literature usually correspond to a set of connected equipment; when the streams of concentrate and tailings leave one, they enter another, either downstream or upstream. Superstructures are important because they define the alternatives and size of the problem [38]. According to Sepulveda et al. [45], circuits of the flotation process generally use three to five stages, i.e., the superstructure of a generic circuit generally considers a rougher stage (R), a cleaner stage (C), and a scavenger stage (S), as shown in Figure 1.

For the purpose of mathematical modeling, simulation and optimization using a DES framework, only the rougher cell is considered, as described in more detail in the following sections.

2.3. Mathematical Modeling of the Flotation Process in Rougher Cells

For modeling, scaling and analysis purposes, two main zones are distinguished in the flotation cell [15]:

• Collection Zone: where the bubble–particle aggregate (true flotation) is formed; referred to as the pulp interface. This area is described by collection recovery (R_C).
• Cleaning Zone: located between the pulp–froth interface and the overflow of the concentrate, where the dragged particles have the possibility of returning to the collection area. This area is described by froth recovery (R_F).

Overall cell retrieval $R_G$ (see Equation (1)) [15] is related to the recovery of the collection $R_C$ area and froth recovery $R_F$, considering the overall mass balance of each cell. Then, the generalization of the concentrate output ($C$) and the tails ($T$) of the $k$ species in each cell $i$ are presented in Equations (2) and (3), respectively.

$$R_G=\frac{R_C\times R_F}{1-R_C\left(1-R_F\right)}$$

(1)

$$C_{ik}=R_{ik}\times F_{ik}, C_i={\displaystyle \sum }_{\forall k}{C}_{ik}$$

(2)

$$T_{ik}=\left(1-R_{ik}\right)\times F_{ik},T_i={\displaystyle \sum }_{\forall k}{T}_{ik}$$

(3)

where $F_{ik}$ is the mass flow of the feed of species $k$ in the stage of flotation $i$, and $C_{ik}$ is the mass flow of the $k$ species in the stream of the concentrate at the $C_i$ in flotation stage $i$, while $T_{ik}$ is the mass flow of the $k$ species in the tail flows, $T_i$ at the flotation stage $i$, and $R_{ik}$ is the recovery of the $k$ species in flotation stage $i$ (equal to Equation (1)), where $k=1,2,\dots ,m$, $i\in I$ and $m$ is the number of species to recover those present in the feeding [15].

2.3.1. Pulp Flotation Models

The recovery of ore $R_C$ in the collection area can be modeled analytically by the general expression shown in Equation (4) [15,46,47].

$$R_C=R\underset{0}{\overset{\infty }{{\displaystyle \int }}}\underset{0}{\overset{\infty }{{\displaystyle \int }}}\left(1-e^{-k\times t}\right)\times E\left(t\right)\times F\left(k\right)dtdk$$

(4)

where E(t) is the residence time distribution (RTD) function for continuous processes with different mixing characteristics, F(k) is the distribution of the flotation rate for different species of minerals with different flotation ratios, and $R$ is a factor that represents the maximum recovery of ore in infinite time [48]. The term $\left(1-e^{-k\times t}\right)$ represents mineral retrieval as a first-order process with a kinetic constant $k$, such as a time function [17,49].

The structure of the cell residence time distribution model (Equation (5)) used to describe the mixing regime in the flotation cells used for modeling of this work is a float cell model reported in [46,47,48], considering a single flotation cell and a perfect mixer with dead time [50,51].

$$E\left(t\right)=\frac{1}{{\tau }_{Mc}}{e}^{-\frac{\left(t-{\tau }_{Dc}\right)}{{\tau }_{Mc}}}$$

(5)

where ${\tau }_{Dc}$ is the delay, i.e., the time at which E begins to change ($E\left({\tau }_D\right)=0$), and ${\tau }_{Mc}$ is the average residence time of feeding in the cell, assuming that the ore inputs are evenly distributed. The residence time is one of the factors affecting the grade of concentrate and mainly the recovery of minerals. The mean residence time of the pulp in the collection area can be estimated by the relationship between the effective volume of the collection area, $V_{eff}$, and the volumetric flow of the feed pulp, $Q_f$ [52].

At the same time, the distribution of the flotation rate F(k) depends mainly on the mineralogical characteristics of the feed, which float at different speeds. Therefore, several models have been proposed in the literature that include normal, gamma or rectangular distributions in order to make the description of the flotation phenomena more flexible. In this study, the model developed by Yianatos et al. [48] was used to model the kinetic distribution F(k), which is a model widely used in the scientific literature [49,53,54,55,56,57,58]. Then, the rectangular model for the ith mineral class and for the nth cell, given by Equation (6), was used to simulate the flotation rate distribution in a distributed manner due to its flexibility and low number of parameters (${\mathrm{k}}_{\max , \mathrm{i}}^{\mathrm{n}}$) [17,48].

$${\mathrm{F}}_{\mathrm{i}}^{\mathrm{n}}\left(k\right)=\left\{\begin{array}{ccc}\frac{1}{{\mathrm{k}}_{\max , \mathrm{i}}^{\mathrm{n}}}& , \mathrm{si}& 0<k<{\mathrm{k}}_{\max , \mathrm{i}}^{\mathrm{n}}\\ 0& , \mathrm{si}& k>{\mathrm{k}}_{\max , \mathrm{i}}^{\mathrm{n}}\end{array}\right.$$

(6)

The kinetic constant under feeding conditions and operational conditions is given by the analytical model developed by Pyke et al. [59] (see Equation (7)).

$$k=2.39\frac{G_{fr}}{d_b{V}_{cell}}\left[0.33\frac{{\epsilon }^{4/9}{d}_b^{7/9}}{{\nu }^{1/3}}{\left(\frac{{\rho }_p-{\rho }_f}{{\rho }_f}\right)}^{2/3}\frac{1}{v_b}\right]E_{coll}$$

(7)

The first term, $2.39\frac{G_{fr}}{d_b{V}_{cell}}$, considers the mechanical characteristics of the flotation cells as a function of the gas flow rate and the cell volume. $G_{fr}$ is the gas flow rate (${\mathrm{cm}}^3/\min$),$d_b$ is the bubble diameter (cm) and $V_{cell}$ is the reference volume of the flotation cells (${\mathrm{cm}}^3$). The second term, $\left[0.33\frac{{\epsilon }^{4/9}{d}_b^{7/9}}{{\nu }^{1/3}}{\left(\frac{{\rho }_p-{\rho }_f}{{\rho }_f}\right)}^{2/3}\frac{1}{{\nu }_b}\right]$, indicates the turbulent nature of the flow field within the tank, where $\epsilon$ represents the turbulent dissipation rate (${\mathrm{cm}}^2/{\mathrm{s}}^3$), reflecting the rate at which smaller turbulent eddies are converting kinetic energy into thermal energy, ν is the kinematic viscosity of the fluid (${\mathrm{cm}}^2/\mathrm{s}$), ${\rho }_p$ is the particle density ($\mathrm{g}/{\mathrm{cm}}^3$), and ${\rho }_f$ is the fluid density ($\mathrm{g}/{\mathrm{cm}}^3$). Furthermore, the term $1/{\nu }_b$ is the effect of bubble speed on the nature of the flow, for which Pyke et al. [59] used an individual value equal to the average velocity of the fluid. Finally, the last term of Equation (6), $E_{coll}=E_c·E_a·E_s$, consists of the micro processes that occur between solid particles and air bubbles in the cell.

As detailed in Equation (7), the most crucial term in the estimation is $E_{coll}=E_c·E_a·E_s$, where the following assumptions are normally considered [56]: particles and bubbles are spherical; particles are finer than bubbles; and the flow around the bubble is modeled as if the bubbles were stationary in a flow field, resulting in an equivalent rate of bubble rise. Then,$E_c$, $E_a$ and $E_s$ are consecutive sub-processes comprising particle–bubble collision, attachment and stability, respectively.

The model developed by Yoon and Luttrell [60] ($E_c^{YL}$) has been broadly accepted by researchers [59,60,61,62], who have applied the interception mechanism and an empirical flow function valid for intermediate flow conditions in which the Reynolds number of bubbles is between 1 and 100. However, it is only applicable to particles smaller than 100 μm and bubbles smaller than 1 mm with immobile surfaces, and does not cover ranges of particles and bubbles under flotation conditions [56]. The analytical model proposed by Yoon and Luttrell [60] is presented in Equation (8), for bubbles with intermediate Reynolds numbers in the range $1<R{e}_b<100$.

$$E_c^{YL}=\left(\frac{3}{2}+\frac{4R{e}_b^{0.72}}{15}\right){\left(\frac{d_p}{d_b}\right)}^2$$

(8)

On the other hand, $E_a$ is dependent on the characteristics of the mineral surface and the degree of adsorption of the collector on the mineral surface [63], which has been extensively studied in the literature [59,60,61,64]. In this case, the Yoon and Luttrell model [60] was used in the present work, with the $E_a$ factor defined as presented in Equation (9).

$$E_a^{YL}=si{n}^2\left(2arctan\left(exp\left[\frac{-v_b{t}_{ind}\left(45+8R{e}_b^{0.72}\right)}{15d_b\left(d_b/d_p+1\right)}\right]\right)\right)/si{n}^{2}90$$

(9)

where $t_{ind}$ represents the induction time, as defined by the expression $t_{ind}=\frac{75}{\theta }{d}_p^{0.6}$ ($\theta$ in degrees and $d_p$ in m) by Koh and Schwarz [65], where $\theta \left(°\right)$ is the bubble–particle contact angle. The model presented in Equation (9) is described in greater detail in the work reported by Dai et al. [24]. A modified stability model $E_s$ [66] for determining the efficiency in separating bubbles and particles is shown in Equation (10).

$$E_s=1-exp\left(A_s\left(1-\frac{F_{att}}{F_{det}}\right)\right)$$

(10)

The value $A_s$ is considered to be a constant of mathematical adjustment [67], while $F_{att}$ is defined as the adhesion force and $F_{det}$ is defined as the separation force. Then, the $F_{det}/F_{att}$ ratio is defined as presented in Equation (11) [56,68], where $\mathsf{\sigma }\left(\mathrm{N}/\mathrm{m}\right)$ is the surface tension and $\mathrm{g} \left(\mathrm{m}/{\mathrm{s}}^2\right)$ is the gravitational constant.

$$\frac{F_{det}}{F_{att}}=\frac{d_p^2\left[\left({\rho }_p-{\rho }_f\right)g+1.9 {\rho }_p{\epsilon }^{\frac{2}{3}}{\left(\frac{d_p}{2}+\frac{d_b}{2}\right)}^{-\frac{1}{3}}\right]+1.5d_p\left(\frac{4\sigma }{d_b}-d_b{p}_{f}g\right)si{n}^2\left(\pi -\frac{\theta }{2}\right)}{\left|6\sigma sin\left(\pi -\frac{\theta }{2}\right)sin\left(\pi +\frac{\theta }{2}\right)\right|}$$

(11)

2.3.2. Froth Flotation Models

The cleaning zone is characterized by the froth recovery, which can be estimated from the measurement of the bubble charge; from modeling point of view, froth recovery, taking into account the effect of froth transportation, can be expressed using Equation (12).

$${\mathrm{R}}_{\mathrm{F}}=\int R_f\left(\beta ,t_f\right)·E\left(t_f\right)d{t}_f$$

(12)

where $R_f\left(\beta ,t_f\right)$ describes the detachment mechanism in the froth zone as a function of froth residence time ($t_f$) and froth stability factor ($\beta$), while $E\left(t_f\right)$ is the distribution function of froth residence time. Additionally, neither $R_f\left(\beta ,t_f\right)$ nor $E\left(t_f\right)$ are independent of each other, and they are strongly influenced by operation cell conditions [69]. The froth recovery $R_f$ is given by Equation (13) [54,69,70,71].

$$R_f\left(\beta ,t_f\right)=\alpha e^{-\beta {\tau }_f}$$

(13)

where α is the maximum recovery in the froth zone, β is the froth stability factor, and ${\tau }_F$ is the gas mean residence time in the froth [15,54]. Equation (13) is a practical approach to representing froth recovery, as is demonstrated later on the basis of the experimental results. The ${\tau }_F$ factor is given by Equation (14) [15,72], where $H_F$ is the froth depth and $J_{Gf}$ is the superficial gas velocity in the froth zone. The β parameter of Equation (13) is modeled as shown in Equation (15).

$${\tau }_f=\frac{H_F{\epsilon }_{Gf}}{J_{Gf}}$$

(14)

$$\beta =\vartheta M_C{R}_C$$

(15)

where $M_C$ is the feeding of valuable minerals into the cell (that is, the product of the flow rate entering the cell and the concentration of valuable minerals, $M_C=Q\times C$) and $\vartheta$ is a mathematical fit parameter. Then, replacing Equations (14) and (15) in Equation (13), the separation mechanism in the froth phase is given by Equation (16).

$$R_f=\alpha e^{-\vartheta QC{R}_C\frac{H_F{\epsilon }_{Gf}}{J_G}}$$

(16)

On the other hand, the residence time distribution function in the froth phase maintains the dynamics of the collection phase, with different estimates of mean times (${\tau }_{Mf}$) and delay (${\tau }_{Df}$), giving as a result the mathematical model presented in Equation (17).

$${\mathrm{R}}_{\mathrm{F}}=\int \alpha e^{-\vartheta QC{R}_C\left(\frac{H_F{\epsilon }_{Gf}}{J_G}\right)}\left(\frac{1}{{\tau }_{Mf}}\right)e^{-\left(\frac{t_f-{\tau }_{Df}}{{\tau }_{Mf}}\right) }d{t}_f$$

(17)

2.4. Adjustment of Superstructure in the Arena Simulation Software

Once the flotation process workflow has been mathematically modeled, it is possible to model and simulate the flotation stage in the Arena simulation software. The updating of mineral recovery over time is simulated by parameterizing the analytical models (see Equations (1)–(17)) recovered from the literature, and incorporating them into the Arena simulation software [40]. A high-level schematic of the simulation model is presented in Figure 2. Updating of the copper concentration status is carried out continuously throughout the course of each production campaign. At the same time, the use of operational parameters is updated in the module “cell parameters update”. The concentration under operational parameters is obtained from the analytical models derived from Equation (1) (function of time, mineralogical conditions of the feed, and maximum concentration of each mineral under operational conditions, among others).

The flotation process of copper minerals is based on production campaigns, which are determined by the availability of inventories from the comminution phase, the development of the campaign corresponding to the production of each stack, the limited production capacity according to the flotation cells, the downstream flow, and the storage capacity. For each feed run, the expected concentration of ore is measured on the basis of the adjusted analytical models. The storage of minerals is carried out according to the logic of inventory theory [73], where the crushing product is kept on hold until the end of each campaign. The main operative variables were obtained from the operational data of pilot plant of a mining company in the Antofagasta region, Chile. These variables are used as input variables for the analytical models in order to estimate the expected copper recovery in rougher cells under operating conditions.

3. Results and Discussion

3.1. Practical Validity of Modeled System Equations

In evaluating the practical validity of the mathematical models of the studied system extracted from the literature (as presented in Equations (1)–(17)), the estimates of the expected mineral recovery and the real recovery were contrasted in the context of operational data, considering the sampled operational variables (gas flow rate, bubble size, particle size, bubble speed, particle density, fluid density and superficial gas velocity) and estimates of non-sampled variables based on a review of the literature (turbulence dissipation rate, surface tension, particle-bubble contact angle, kinematic fluid viscosity, Reynolds number of the bubble and gas hold-up), with the contrast between the expected recovery of minerals in the rougher cell ($R_G$) versus the predicted recovery (${\widehat{R}}_G$) (see Figure 3a) indicating that the model is effective for estimating the responses, despite the variability of the predictions with respect to the real recovery, which can be attributed to the assumptions regarding the unsampled variables (mean values). Then, the distribution of the standardized residuals (see Figure 3b) indicates that they exhibit a normal trend, while the normal probability distribution plot of the residuals (see Figure 3c) validates the goodness of fit of the analytical models of mineral recovery from the rougher cell, contrasting the real recovery with the expected recovery and verifying the assumption of normality of the residuals, given that the p-value is greater than the significance level (α = 0.05).

3.2. Establishing the Domain of the Variables of Interest

A rougher cell is simulated in order to evaluate variations flotation performance by incorporating analytical models within a DES framework, integrating the mineralogical characteristics under uncertain conditions.

Table 1. Schematic summary of the independent variables.

Parameter and Variables	Units	Lower Limit	Upper Limit
Gas flow rate ($G_{fr}$)	[${\mathrm{cm}}^3/\min$]	3000	4000
Bubble Size ($d_b$)	[cm]	0.05	0.10
Particle size ($d_p$)	[cm]	0.0020	0.0075
Bubble Speed ($v_b$)	[cm/s]	10	30
Turbulence dissipation rate ($\epsilon$)	[W/kg] [${\mathrm{m}}^2/{\mathrm{s}}^3$]	18	30
Kinematic fluid viscosity ($\nu$)	[${\mathrm{cm}}^2/\mathrm{s}$]	0.008	0.01
Particle Density (${\rho }_p$)	[$\mathrm{g}/\mathrm{c}{\mathrm{m}}^3$]	1.3	4.1
Fluid Density (${\rho }_f$)	[$\mathrm{g}/{\mathrm{cm}}^3$]	0.9	1.1
Surface tension ($\sigma$)	[mN/m]	35	75
Reynolds number of the bubble ($R{e}_b$)	[dimensionless]	0.2	100
Bubble–particle contact angle (θ)	[◦]	30	90
Superficial gas velocity ($J_G$)	[cm/s]	1.5	2.2
Gas hold-up (${\epsilon }_{Gf}$)	[dimensionless]	0.92	0.98

defines the theoretical domain of the variables involved in the process for modeling, simulation and optimization purposes based on operational data, the knowhow of the authors, and a literature review, while additionally including ranges of variables not measured in the field, such as the bubble–particle contact angle, surface tension, and turbulence dissipation rate. (The inclusion of variables not measured experimentally is only for analytical purposes.)

The limits of the operational parameters shown in Table 1 correspond to the domain of the flotation process in the rougher cell of a pilot plant in a mine in Chile’s Antofagasta region. Additionally, for modeling, simulation and optimization purposes, the parameters average residence time in collection and froth phase (${\tau }_{Mc},{\tau }_{Mf}$), delay time in the collection and froth phase (${\tau }_{Dc},{\tau }_{Df}$), gas concentration in the collection phase (${\epsilon }_{Gc}$), cell volume ($V_{cell}$), material flow rate at the entrance to the cell (Q), concentration of valuable mineral at the entrance to the rougher cell (C), maximum recovery in the froth zone ($\alpha$) And depth of the froth phase ($H_F$) are considered to be constant, with values of 540 and 24 [s], 20 and 2 [s], 0.08, 0.5616 [m²], 1.5 [ton/h], 0.6 [%], 95 [%] and 4 [cm], respectively.

3.3. Sensitization of Mineral Recovery

Analytical models (see Equation (1)) for recovery from grounded ore in a rougher cell were implemented using an optimization model. By this means, the copper recovery of the rougher cell can be modeled as an inverse exponential model, while the sensitivity analysis for low, medium and high levels of the variables of interest is presented in Figure 4.

The analysis of the parameters from the sensitization of the analytical model for the rougher cell (see Equation (1)) indicates that particle size (see Figure 4c), bubble velocity (see Figure 4d) and particle density (see Figure 4h) have a significant impact on the response. Meanwhile, it is also indicated that variables like gas flow rate (see Figure 4a), bubble size (see Figure 4b), turbulence dissipation rate (see Figure 4e), fluid kinematic viscosity (see Figure 4f), fluid density (see Figure 4h) and bubble Reynolds number (see Figure 4j) have a moderate impact in response. In contrast, the surface tension (see Figure 4i), bubble–particle contact angle (see Figure 4k), superficial gas velocity (see Figure 4l) and gas hold-up in froth zone (see Figure 4m) have a marginal impact on mineral recovery.

It is important to note that the sensitivity analysis developed in Figure 4 is for analytical purposes only, and the variations are not necessarily applicable to real flotation circuits, given the impossibility of modifying a given parameter without affecting the other operative variables in a mining worksite. Estimates for the recovery of the rougher cell with different levels of the independent variables are presented in

Table 2. Sensitivity analysis of mineral recovery for low, medium and high levels of the independent variables.

Variable	Units	Levels [Low, Medium, High]	Mineral Recovery of Rougher Cell [Time = 60 min]
			Low	Medium	High
$\mathrm{Gas} \mathrm{flow} \mathrm{rate} (G_{fr}$)	[${\mathrm{cm}}^3/\min$]	[3000, 3500, 4000]	0.411678	0.430877	0.447594
$\mathrm{Bubble} \mathrm{Size} (d_b$)	[cm]	[0.065, 0.075, 0.085]	0.417253	0.430877	0.426852
$\mathrm{Particle} \mathrm{size} (d_p$)	[cm]	[0.001, 0.005, 0.010]	0.144555	0.430877	0.437685
$\mathrm{Bubble} \mathrm{Speed} (V_b$)	[cm/s]	[10, 20, 30]	0.487739	0.430877	0.376153
Turbulence dissipation rate ($\epsilon$)	[W/kg] [${\mathrm{m}}^2/{\mathrm{s}}^3$]	[18, 25, 30]	0.412690	0.430877	0.441012
Kinematic fluid viscosity ($\nu$)	[${\mathrm{cm}}^2/\mathrm{s}$]	[0.008, 0.009, 0.010]	0.435843	0.430877	0.426446
$\mathrm{Particle} \mathrm{Density} ({\rho }_p$)	[$\mathrm{g}/{\mathrm{cm}}^3$]	[1.3, 2.5, 4.1]	0.300519	0.430877	0.491779
$\mathrm{Fluid} \mathrm{Density} ({\rho }_f$)	[$\mathrm{g}/{\mathrm{cm}}^3$]	[0.9, 1.0, 1.1]	0.445052	0.430877	0.417223
Surface tension ($\sigma$)	[mN/m]	[35, 40, 45]	0.430877	0.430877	0.430877
Reynolds number of the bubble	[dimensionless]	[10, 50, 100]	0.404272	0.430877	0.396849
Bubble–particle contact angle ($\theta$)	[◦]	[30, 60, 90]	0.430877	0.430877	0.430877
$\mathrm{Superficial} \mathrm{gas} \mathrm{velocity} (J_G$)	[cm/s]	[1.5, 1.8, 2.2]	0.430685	0.430685	0.430685
Gas hold-up in froth zone (${\epsilon }_{Gf}$)	[dimensionless]	[0.92, 0.95, 0.98]	0.430861	0.430861	0.430861

.

3.4. Comparation of Simulated Scenarios

Developing the model in the Arena simulation software and incorporating the inherent uncertainty of the input parameters, an analysis of variance was carried out to determine the impact of these on the response variable. The contrast of the different levels of each factor in the response was studied using the Games–Howell test, a non-parametric approach used to compare combinations of samples or groups, as well as the differences between them. Because the Games–Howell test is not based on equal sample sizes and variances, it is often recommended over other approaches for performing comparisons between sample data. The hypothesis test considered the following research hypotheses:

• $Ho: \Delta E\left[Mineral Recovery\right]\left(i,j\right)=0 \forall i\ne j\in \left\{low, medium, high\right\}$
• $Ha: \exists i\ne j\in \left\{low, medium, high\right\} | \Delta E\left[Mineral Recovery\right]\left(i,j\right)\ne 0$

Then, the recovery behavior with respect to changes of the independent variables is presented in

Table 3. ANOVA of recovery distributions of the rougher cell for low, medium and high levels of independent variables.

Indicator	ANOVA		Comparisons
Parameter/Level	F Value	p-value	Games-Howell
Gas flow rate ($G_{fr}$)	22.107	0.000	Differ
Bubble Size ($d_b$)	10.024	0.008	Differ
Particle size ($d_p$)	35.263	0.000	Differ
Bubble Speed ($v_b$)	27.334	0.000	Differ
Turbulence dissipation rate ($\epsilon$)	4.9163	0.003	Differ
Kinematic fluid viscosity ($\nu$)	3.8574	0.083	They do not differ
Particle Density (${\rho }_p$)	31.284	0.000	Differ
Fluid Density (${\rho }_f$)	3.6159	0.048	Differ
Surface tension ($\sigma$)	1.2438	0.092	They do not differ
Reynolds number of the bubble ($R{e}_b$)	13.894	0.003	Differ
Bubble–particle contact angle ($\theta$)	1.1742	0.115	They do not differ
Superficial gas velocity ($J_G$)	1.1836	0.107	They do not differ
Gas hold-up in froth zone (${\epsilon }_{Gf}$)	1.0564	0.145	They do not differ

.

The simulation results considering the distributions of response to low, medium, and high levels of the parameters of interest are presented in Table 3. It can be seen that with varying gas flow rate, bubble size, particle size, bubble speed, kinematic bubble viscosity and particle density, the mean recovery differs significantly at low, medium, and high levels of the parameters in question. In contrast, there is insufficient evidence that the average recovery differs with variations in dissipation rate, kinematic viscosity, fluid density and surface tension, as indicated by the Games–Howell comparison test [74].

3.5. Optimization

After modeling the rougher cell in Arena, the model is optimized using the Arena OptQuest module [40]. The integration of the simulation with the optimization makes it possible to evaluate the system performance, measuring and inferring correlations of interest between the explanatory variables. Opquest is an Arena module that expands the analysis capabilities, making it possible to automatically search for optimal solutions in a simulation model by searching for the values of the control variables that maximize the objective.

The model restrictions with respect to the operating variables and the optimal operating results obtained by Optquest are summarized in

Table 4. Optimal levels of operating parameters.

Parameters	Dimensionality	Optimal Value
Gas flow rate	[${\mathrm{cm}}^3/\min$]	3491.83
Bubble size	[cm]	0.0708
Particle size	[cm]	0.01
Bubble speed	[cm/s]	10.00
Turbulence dissipation rate	[W/kg] [${\mathrm{m}}^2/{\mathrm{s}}^3$]	24.99
Fluid kinematic viscosity	[${\mathrm{cm}}^2/\mathrm{s}$]	0.008
Particle density	[$\mathrm{g}/\mathrm{c}{\mathrm{m}}^3$]	4.1
Fluid density	[$\mathrm{g}/{\mathrm{cm}}^3$]	0.9
Surface tension	[mN/m]	39.863
Reynolds number of the bubble	[dimensionless]	41.885
Bubble–particle contact angle	[◦]	59.739
Superficial gas velocity	[cm/s]	1.8
Gas hold-up in froth zone	[dimensionless]	0.95
Time	[s]	3600

, while a discussion of the results is presented in the following section.

4. Discussion

On the basis of the mathematical modeling of a rougher cell, and its simulation and optimization, the following discussion is presented.

• The optimum operating levels are reached at a medium gas flow rate (within the domain described in Table 2), indicating that the gas flow velocity is a measure of the aeration capacity of the flotation equipment, and, at the same time, has a considerable impact on the efficiency of the circuit. In other words, the gas flow rate and the bubble velocity indicate the way in which the air entering the flotation cell is dispersed through it [32]. It should be noted that with higher values of gas flow rate, there are two effects with conflicting impacts on mineral recovery: on the one hand, the larger the bubble size, the lower the collision efficiency (see Equation (8)); and on the other hand, the increased availability of bubble area to capture minerals should result in an increase in recovery rate. Thus, the maximum mineral recovery rate is obtained at medium gas flow rates [75], as shown by the optimal configurations shown in Table 4.
• The efficiency of the gas flow rate, in addition to the bubble velocity, is dependent on the geometry of the flotation cell. For example, there is a high superficial gas velocity in the corners of the tank [76]. Therefore, the optimal bubble size, considering that the PDF of the distribution is non-Gaussian, together with the analysis of the sample data and the corresponding literature [25,37,77], is 0.0708 cm. The distribution characteristics of bubble size reveal a skew distribution to the left, which leads to the logical conclusion that the smaller the bubble size, the greater the stability. This means a lower probability of the existence of samples with larger surface bubbles, which is in line with the analysis of the analytical model for low, medium, and high levels presented in Figure 4b, where it is evidenced that at larger bubble sizes, the degree of mineral recovery is lower than at medium levels, which can be explained by the inverse relationship between the bubble size and their stability and collision efficiency [75]. Analyzing the gas flow rate, it can be seen that the flotation rate, and therefore the mineral concentration, is directly proportional to the gas flow rate for particle sizes > 10 µm (for finer particles, it is not possible to affirm the direct proportionality or the difference) [61]. Controlling the gas flow rate parameter can be achieved by generating a control mechanism for the air supply and the impeller rotation rate.
• Particle size is one of the most important parameters in flotation processes, and the recovery dynamics based on this parameter can be divided into three size classes: an intermediate size range where the recovery is high; a fine size range where the recovery is lower and increases over time; and a thick size range where recovery is again lower, but less affected by time. The size divisions are not precise, and will vary with the type of mineral [78]. Control of particle size is dependent on the comminution stages that are executed prior to the flotation stage [79,80]. Considering the above, particle size is an essential parameter in the efficiency of the flotation process [81], resulting in an optimal value of approximately 100 µm, which is within the optimal particle range, which is usually between 30 and 120 µm. Additionally, Wyslouzil et al. [82] indicated that the efficiency of the flotation process is negatively impacted when operating at the extremes, that is, when working with particles that are too fine (<10 µm) or too coarse (>200 µm). Meanwhile, Pyke et al. [59] indicated that the constant float velocity versus particle size curve increases almost linearly with particle size. Several methods have been proposed as measures for monitoring bubble size, among which is the one developed by Cheng et al. [83], consisting of inserting a viewing camera made of transparent acrylic into a PVC collection tube inserted diagonally into the flotation cell, while the size can be manipulated by means of the control of the added frothers and the dispersion mechanism [84].
• Flotation requires a certain degree of turbulence for several reasons, including the maintenance of solids suspended in the pulp phase, introducing air into the pulp so that it disperses into bubbles, and mixing the aerated pulp to achieve the distribution and conditioning of reagents, as well as providing opportunities for collisions between particles and bubbles [85]. This variable is directly related to the speed of the impeller, i.e., the higher the speed of the impeller, the greater the separation of solid particles from air bubbles, and therefore the less mineral recovered [61]. Therefore, the higher the turbulence dissipation rate, the lower the probability of maintaining solids in suspension, while the lower the dissipation rate, the greater the probability of the recovery of solid particles from the froth. The control of the turbulence dissipation rate parameter depends on several factors [86], with the main ones being fluid density, bulk flow velocity, and impeller dynamics, among others [87]; it is possible to control these variables to impact the response. Analysis of the turbulence dissipation rate shows that at low levels, the floating rate tends to remain constant, while at high turbulent dissipation energies, the floating rate constants increase in magnitude because of the increased frequency of bubble–particle collision, while they also decrease due to bubble–particle instability. Since the latter effect is of greater magnitude for coarser particles, these calculations produce the characteristic shape of the float rate constant against the particle size curve, with the maximum being obtained at intermediate particle sizes, as has been reported in practice [59]. Therefore, the optimal results are achieved at mean levels of the variable, with the local optimum being obtained when the value of the turbulence dissipation rate is approximately 25 ${\mathrm{m}}^2/{\mathrm{s}}^3$.
• Chemical modifiers are chemical reagents that alter the properties of fluid flow and in particular the suspension of solids in liquid media [88], and handling the response depends mainly on the type and amount of chemical additives used to modify the flow properties. Regarding the study of the densities of both the particles and the fluid, the values of particle density or the bubble–particle collision must be greater than that of the fluid, as confirmed in [56]. A study by Mehrotra and Kapur [89] indicated that the lower the pulp density, the higher the flotation ratio, and therefore, the higher the mineral concentration. Other studies show that the pulp density affects the key parameters of the pulp and froth phases, since the kinetic rates of the pulp phase decrease at high density, which is attributed to a reduction in turbulence within the cell [90]. Particle density was assumed to be ≤4.1 $\mathrm{g}/{\mathrm{cm}}^3$, due to the significant discrepancies between the models and the experimental information found in most dense minerals [22,23,24,56,91,92,93,94], and higher than the density range of the fluid considered for the study. The optimal particle density was found to be at the maximum levels of this variable, which coincides with the particle density for copper sulphide minerals found in the literature [91], while the optimum fluid density level was 0.9 $\mathrm{g}/{\mathrm{cm}}^3$.
• It is widely accepted that particles stabilize the froth, improving the flotation performance. However, predicting the effect of particle addition on froth stability is challenging, since it is highly dependent on the lifetime of the bubbles (typically 0.1 to 60 s), which is analogous to variations in the air velocity in the flotation cells. Using the maximum bubble pressure method, the results of the work reported by Hadler and Cilliers [95] show that the addition of particles results in lower surface tension, both in the dynamic (i.e., short-term) lifetime of the bubble and toward equilibrium (i.e., the bubble at 60 s). This research concludes that increased particle loading at the air–water interface, either through higher surfactant concentrations or lower air velocities (longer bubble life), results in lower surface tension and greater froth stability [95]. In this work, the greatest recovery efficiency was achieved at medium levels (39.863 mN/m) of the range established for the surface tension parameter. Surface tension depends directly on the frothing agents used in the process [96], finding that the higher the surface density, the lower the concentration of the froth.
• With respect to the speed of the bubbles, the relationship between the bubble size and the elevation rate is directly proportional, that is, small bubbles rise more slowly than large bubbles, and small bubbles in the presence of many solutes rise more slowly than in water alone [97]. Controlling the bubble rate depends on both the bubble size and gas flow rate. Therefore, the simulation model developed indicates that the optimal levels of bubble velocity are reached at lower levels of the considered range of this variable.
• Another variable of interest in the optimization process of the operating dynamics of rougher cells is the kinematic viscosity of the fluid, defined as the ratio between the dynamic viscosity and the density of the fluid [98]. Viscosity is a measure of the fluid’s resistance to gradual deformation, either by shear stress or tensile stress. Therefore, higher the viscosity levels, the greater the probability of delay in the drainage of hydrophobic particles, and therefore, the lower the efficiency of mineral recovery [99]. At the current stage, to mitigate the negative impact of high viscosities on flotation efficiency, some polymeric dispersants and metal ions are effective in reducing this variable and improving flotation in some applications [100]. Therefore, the optimal kinematic viscosity is reached at low levels of the factor ($\mathsf{\nu }=0.008{ \mathrm{cm}}^2/\mathrm{s}$), considering the range defined in Table 2, which is in agreement with the results reported in the literature, while the Reynolds number of the bubble—the nondimensional ratio between the inertial and viscosity forces of the moving fluid—reaches its optimal level at medium levels of the factor ($R{e}_b\cong 42$); this can be explained by the impact of the variable on both the collision efficiency and the attachment efficiency, that is, increases in Reynolds number result in higher collision efficiency ($d{E}_c^{YL}/dR{e}_b>0$), but lower adhesion efficiency ($d{E}_a^{YL}/dR{e}_b<0$); conversely, decreases in Reynolds number lead to lower collision efficiency and higher attachment efficiency.
• Finally, the optimal level for bubble–particle contact angle is 39.739 (see the theoretical impact of the contact angle on the stability model, presented in Appendix A), while those for superficial gas velocity and gas hold-up in the froth zone are 1.8 and 0.95, respectively; that is, the optimum levels are reached at medium levels for these factors. Additionally, the volume of the flotation cell is not a parameter whose variation is efficient from a practical point of view (therefore, it was not sensitized in this work); however, the size and shape of the tank does have an impact on mineral recovery, as validated by studies considering the effects of the hydrodynamic phenomenon with respect to geometric and dynamic similarities [101,102]. On the other hand, the recovery versus time function increases monotonously, with an asymptotic behavior at higher levels of this factor, meaning that the optimal time for mineral recovery from the froth depends on an economic analysis. However, to calculate the optimum according to the defined domain, the time that maximizes the response is 60 min, that is, the upper bound of the established domain.

Complementing the simulation model and the corresponding optimization, the incorporation of machine learning techniques into conventional optimization could have the potential to identify operational variables with a considerable impact on mineral recovery in rougher cells or RCS circuits. This might provide a better understanding of the dynamics of the process through the application of algorithms such as neural networks [103,104], making it possible to obtain an abstraction of the process, or Bayesian networks [105], making it possible to identify causal dependency relationships between the operative variables.

5. Conclusions and Future Works

5.1. Conclusions

Mineral deposits tend to be heterogeneous, making it necessary for the processing parameters to evolve and to adjust in response to the feeding characteristics. In Chile, the production of copper from oxidized minerals is decreasing, giving way to the extraction (en masse) of copper from sulphide minerals through flotation processes. In this research, a mathematical model of a rougher cell was developed and a simulation was performed to study the impact of operational parameters on mineral recovery. The impacts of the variables gas flow rate, bubble size, particle size, bubble speed, turbulence dissipation rate, fluid kinematic viscosity, particle density, fluid density, surface tension, Reynolds number of the bubble, bubble–particle contact angle, superficial gas velocity, gas hold-up in froth zone and time in mineral recovery were examined; however, the framework could be extended to various types of minerals (using variables such as types of frothers or flocculants, for example), and could be applied in a variety of geological domains in mines.

The analytical models present a good fit to the sample data, and on the basis of the simulation analysis for low, medium, and high levels of the operational variables, it was evident that the variables with the greatest impact on mineral recovery were particle size, bubble velocity and particle density, while those parameters that had a moderate influence were the gas flow rate, bubble size, turbulence dissipation rate, fluid kinematic viscosity, fluid density and the Reynolds number of the bubble. Finally, surface tension, bubble–particle contact angle, superficial gas velocity and gas hold-up in the froth zone had a negligible impact on mineral recovery under the modeling, simulation and optimization scheme developed.

The recovery of the rougher cell is dependent on relatively high levels of gas flow rate, which, together with bubble size, determines the aeration capacity of the flotation equipment. This means that for low bubble sizes, there is less contact/adhesion surface, while the larger the bubble size, the greater the instability of the bubble; therefore, the optimal levels correspond to the values between low and medium considered in this study.

Particle size and surface tension, on the other hand, should be centered on the mean levels of the factor since recovery efficiency is negatively impacted when operating at the extremes. The turbulence dissipation rate, as in the case of the particle size, results in weak recovery at the extremes, because at low, the float velocity tends to remain constant, while at high levels, the drop in recovery is due to bubble–particle instability. The optimal particle density is at the maximum levels of this factor, while the optimal fluid density level is considerably lower. The optimal fluid kinematic viscosity is achieved at lower levels of this factor, since increases in viscosity resulting in lower efficiency of mineral recovery, while the Reynolds number of the bubbles should also be centered on medium values, due to the negative impact on both collision efficiency and attachment efficiency. Additionally, the optimal mineral recovery is achieved with high levels or medium levels of bubble–particle contact angle, superficial gas velocity and gas hold-up in the froth zone.

Finally, the recovery function has an asymptotic behavior with increasing values of time, and the optimum is achieved at the maximum value in this domain; however, it should be noted that in an industrial context, the optimal recovery time depends on an economic analysis that determines the marginality of this variable with respect to the response.

5.2. Future Works

To further advance the operational research of flotation processes, the following lines of research are being considered in order to expand the framework to incorporate rougher cells in series: the modeling and simulation of aggregated RCS circuits, and the inclusion of different modes of operation [106], depending on the mineralogical content of the feed. On the other hand, the development of analytical or machine learning models incorporating additional operational variables into the process, such as mineralogical characterization, types of frothing agents or flocculants, could be considered.