Simulation Algorithm for Water Elutriators: Model Calibration with Plant Data and Operational Simulations

Abstract

A dynamic simulation algorithm based on 1-D transient convection/diffusion transport per particle size class is proposed to simulate a hydraulic classifier operated to selectively remove quartz from an iron oxide concentrate produced by processing the ore from an iron ore mine in northeastern Canada. The calibrated model is used to simulate the operation of dense bed hydraulic classifiers of different sizes and/or under different operating conditions. The simulator predicts the behavior and characteristics of the pulp at different depths within the classifier as a function of time. The simulator is validated by confronting the simulation results to experimental data obtained from sampling industrial and laboratory classifiers. The simulator is then used to assess the role of the fluidization or teeter water and of bed density on the quality of the produced separation of quartz from the valuable iron oxide of the processed ore. The knowledge acquired in the noise-free environment of simulation provides clues on the way to manipulate the hydraulic classifier operating variables in a process control strategy for an industrial scale unit.

1. Introduction

Mineral concentration processes refer to the operations used to separate the valuable minerals from the non-valuable minerals (gangue) in mined ores or rocks. Mineral concentration equipment either use the physical (specific gravity, magnetic susceptibility, electrical conductivity, etc.) or chemical (surface composition) properties of the minerals to achieve their separation [1]. Reliance on the differences in the specific gravities between valuable and non-valuable minerals is certainly the oldest and most traditional way to carry out such separations. Gravity separation devices are inexpensive to operate and do not require the addition of chemical reagents to carry out the separation of valuable minerals from the gangue and are thus often considered as environmentally friendly (Burt, 1999 [2]). Das and Sakar (2018, [3]) present a review of gravimetric concentration processes. The present study focuses on a type of gravity separation equipment named hydraulic classifier (Hc) or water elutriator.

The different types of hydraulic classifiers are presented in the work of Tripathy et al. [4] and the one studied here is the floatex density separator. Hc units have been the subject of several studies and publications (Bazin et al. 2011–2012 [5,6], Das et al. 2007–2010 [7,8,9,10]; Sar et al. 2017 [11]; Sadeghi and Bazin 2020 [12], among others). Generally, the prevailing separation mechanisms in an Hc are hindered settling, elutriation and fluidization [4,13]. These separators consist of continuously-fed units that can process large amounts of solids and generate dense fluidized particle beds that can greatly improve the separation acuity (Galvin et al., 1999 [14]). Hc combines low capital installation and operating costs (Luttrell et al., [13]) and should therefore always be considered as an economically viable option for the separation of minerals.

The main objective of the study is to develop a phenomenological model that explicitly accounts for the common operating variables (flow rate of the teeter water, bed density, and so forth) of the classifier and the geometric characteristics of the apparatus. The proposed mathematical model can help the operators to improve the operation of existing processes by simulating the performances of different control strategies or by predicting the performances of separation equipment prior to its commissioning into a circuit consisting of several gravity separators. The separation of particles in an Hc results from a combination of particle size and specific gravity effects given the properties of the liquid (density, viscosity), which is often an aqueous base. Indeed, a coarse and low-density non-valuable mineral particle may behave similarly to a small and dense heavy valuable mineral particle, making it difficult to anticipate the classifier’s separation quality without pilot-scale testing of the ore under scrutiny. However, multiplication of pilot-scale testing may rapidly lead to an increase in the costs of a pre-feasibility study. Hence, the development and use of virtual mathematically-based models for industrial equipment are certainly viewed as the ‘voie royale’ to help reduce the risk associated with the commissioning of a circuit incorporating that equipment or to assist plant operators in designing a control strategy for the existing operating equipment.

The objective of this paper is to propose an Hc mathematical model that can serve the two indicated purposes, namely, design and optimization. To establish the material balance of the Hcs, the algorithm considers that the volume flow of material is the resultant of the forces generated from dispersion and convection (hindered settling or elutriation). This relationship is used, among others, for the simulation of hydraulic classifiers [15,16,17] and for the simulation of fluidized beds [18,19,20,21,22,23]. The use of the dispersion term allows the modification of the velocity profile of the solid particles, which results from a non-uniform flow profile of the liquid, a mixing effect caused by the turbulence of the flow and the rebalancing of the density gradient as a function of height by the Archimedes’ principle, while the convection term allows to model the effect of gravity forces and the movement of the surrounding medium. The method of solving the equation of motion, which is a partial derivative equation, is inspired by the research works of Kim [15], Lee [23] and Honaker [16]. The combination of their solution methods results in a computational algorithm with high stability when there is a large variation of the solid concentration according to the depth.

More recently, Payenzo et al. (2011, [5]) analyzed the operation of an Hc for iron ores using the partition or efficiency curves of the minerals. However, their investigation was limited to simply observing the actual operation of the classifier. The authors concluded on the limited applicability of their model, as it was not possible to use their approach to simulate the response of the classifier to actions implemented on the two main operating variables of the classifier, namely, the density of the bed and the flow rate of the fluidization water. Most of the mathematical models proposed so far for the hydraulic classifiers do not explicitly account for all operating variables of the classifier (feed water and underflow), a weakness that the model proposed in this study attempts to overcome. Unlike many previous studies dealing with hydraulic classifiers that analyze the behavior of spheres in water, this study is concerned with the modeling of an industrial unit processing an actual ore. Despite the fact that gravity separation devices have been operating for eons, there is still no well-established mathematical representation for an industrial Hc. To improve the modeling of the behavior of industrial Hc, it was necessary to introduce new concepts, such as the convection generated by the underflow, supervision and control of simulation variables (concentration, underflow rate and pulp density) and the addition of methods for estimating the parameters of the model (bed porosity and axial dispersion coefficient). The developed mathematical model does not require physical assembly to be simulated and can help operators to improve designs using hydraulic classifiers, allowing the simulation of different control strategies and the design of equipment for particular applications. In addition, the simulation tool allows the prediction of the states of unmeasured and/or unmeasurable variables, such as local density, overflow composition and interface location, thus allowing it to be used as a virtual sensor.

2. Materials and Methods

2.1. Sampling the Industrial Hydraulic Classifier

Most of the data used for the development of the model for the Hc come from the sampling of the industrial Hc units. A sampling campaign on a hydraulic classifier usually lasts less than one hour. The hydraulic classifier was first allowed to operate under steady state during one hour before starting to collect samples from the feed, underflow and overflow streams. Each stream was cut three times during the one-hour sampling campaign and the cuts were combined to give the composite samples used for the analysis.

The samples were analyzed to obtain the slurry solids’ concentration, the elemental (Fe, SiO₂, Al₂O₃, CaO, Na₂O, MgO, etc.) composition and the particle size distribution of the solids. The material within each size fraction was assayed for the same elements as the head sample. The data were finally processed using a material balance program that calculates the solids split between the Hc underflow and overflow streams. The materiel balance algorithm was also built to calculate the mineral (hematite, magnetite, quartz, mica, garnet, ilmenite, etc.) concentrations from the elemental assays for the three sampled streams of the Hc [24,25]. The estimated mineral concentrations and the solid split were used to calculate the previously defined performance indices.

Samples were also collected during steady-state operation at various depths inside the Hc in order to evaluate the type of mixing along the vertical axis. These samples were analyzed to obtain the slurry solids’ concentration and the elemental composition. The results of this sampling, shown in Figure 1, clearly demonstrate [24] that the concentration of solids and the density of the pulp vary with depth and that the mixing within a Hc is clearly not that of a perfect mixer justifying the discretization approach used here for the model development.

The operating principle of the mine, concentrator and pellet plant is detailed in Appendix A.1. Appendix A.2 explains the operation of the Hcs in the pellet plant, while Appendix A.3 discusses the performance indices used to describe the separation efficiency.

Analytical Methods

The solid samples were assayed using a Pananalytical Epsilon 1 Energy Dispersive X-Ray Fluorescence (ED-XRF) (Malvern Panalytical, Westborough, MA, USA). XRF samples were fused with lithium Borax to prepare glasses for the chemical analysis. The mineral composition of the samples obtained from the tests conducted on the laboratory set-up was back calculated from the density of the sample, as measured using a Gas pycnometer (HumiPyc model 2 from Instruquest, Boca Raton, FL, USA). The size distribution of a sample was measured by using Tyler screens ranging from 1200 to 38 µm.

2.2. Laboratory Set-Up

In order to verify the versatility of the simulation algorithm in various situations (geometry, pulp characteristics, flow regime, etc.), experiments were also carried out using a laboratory unit [24].

Table 1. Operating conditions of experiments on industrial and laboratory set-ups.

Characteristics	Units	Type of Hc
		Laboratory	Industrials
			Primary	Secondary	Tertiary
Height	m	1.5	4	4	4
Width	m	0.06	2.13	2.13	3.15
Area	m2	0.003	4.5	4.5	7.8
Shape	-	Circular	Square	Square	Circular
Bed density	t/m3	1.1–1.4	1.7–2.2	1.7–1.8	1.7–1.8
Teeter water	l/min	1–2	~2000	~1000	~1000
Feed (solid)	t/h	0.05–0.1	150–300	80–100	90–110
Feed (water)	l/min	~4	3000–4000	~2000	~4000

gives the characteristics and operating conditions for the laboratory set-up and presents an average value of the operating parameters of the Hcs of the pellet plant in Port-Cartier. The first experimental step in the laboratory set-up was to mix the solid and water in the feed tank. Once this step was completed, the hydraulic classifier was filled with water using the cone water and teeter water injectors. Afterwards, the feed, cone water and teeter water flow rates were set to the selected value for the experiment. Once the system was started, it was operated for 30 min to allow a steady state to be reached, before taking the first samples. The obtained samples were then weighed, dried and sieved. The last step consists in determining the density of the solid using a gas pycnometer (HumiPyc). This last step allows to estimate the proportion of heavy and light particles of the sample.

3. Mathematical Model

3.1. 1-D Unsteady-State Material Balance Equations

To develop a procedure for simulating Hcs, the overall material balance as well as the species-specific material balances (nature, density and particle size class) must be established. The first step in solving the material balance is to make the following assumptions:

Assumption 1.

The pulp is evenly distributed crosswise so that lateral gradients of any property are inexistent.

Therefore, the solid will be uniformly distributed horizontally regardless of the vertical position, making the material balance unidirectional in $z$. Equation (1) and Figure 2 illustrate that the local variation of material volume flux ($\partial J$) over a differential axial slice ($\partial z$) reduces to the accumulation of materials in terms of concentration ($\partial C$) over the time interval ($\partial t$):

$$-\partial J/\partial z=\partial C/\partial t.$$

(1)

The particles’ movement along the vertical direction in the Hc is accounted for in a series of species-labeled transient one-dimensional volume balances, which relate the above-defined accumulation term to the linear convection transport and axial dispersion transport for each species, as shown in Equation (2). This relationship indicates that the variation of the local concentration (m³_solid/m³_element) of a species $k$, at the particle-size class $j$, in the interval $\partial z$ varies over time $\partial t$ due to the additive contributions of axial dispersion and convection terms:

$$\frac{\partial C_{j,k}(z,t)}{\partial t}=\underset{Dispersion}{\underbrace{D\frac{{\partial }^2{C}_{j,k}(z,t)}{\partial Z^2}}}-\stackrel{Convection}{\overbrace{\frac{\partial (V_{j,k}(z,t,\phi )C_{j,k}(z,t))}{\partial Z}}}$$

(2)

where $D$, $C$, $V$, $\phi$, $z$ and $t$ are, respectively, the axial dispersion coefficient, concentration, velocity, total concentration, the height and time. The next step prior to solving these equations consists in discretizing the variables of Equation (2). To do this, the Hc is discretized into $N$ elements of equal height ($\Delta z$), as shown in Figure 2. By convention, element 1 is located just above the underflow valve, element $N_f$ corresponds to the element where the feed is introduced, and element $N$ is just below the unit overflow. It is also necessary to discretize the time variable by selecting an appropriate interval ($\Delta t$). The solution of Equation (2) is further simplified by make the following assumptions:

Assumption 2.

The axial dispersion coefficient is constant for each particle-size class ($j$), species ($k$) and the characteristics of the pulp in depth ($z$).

The initial and boundary conditions of Equation (2) are:

$$C_{j,k}\left(z,0\right)=0 \forall 0<z<H$$

(3)

$$\sum C_{j,k}\left(z,t\right)\le {\phi }_{max} \forall 0<z<H \& t\ge 0$$

(4)

The height, H, of the classifier extends from the underflow to the overflow (Figure 2). The initial condition (Equation (3)) assumes that the system is free of solids at the beginning of the calculations. Boundary condition 4 states that the concentration cannot exceed the maximum concentration that can be obtained experimentally or estimated (fluidized bed) with the method presented in Section 3.3.

The partial temporal derivative term on the LHS in Equation (2) is approximated with a scheme of (Euler) forward differences. For the RHS terms, a centered difference for the dispersive-term second-order spatial derivatives (Leapfrog) [15] and a combination of forward and backward differences for the convective first-order spatial derivatives (Upwind/Downwind) are used [16]. Thus, the following finite-difference simplification of Equation (2) is obtained:

$$\begin{gathered} \frac{{\left(C_{j,k}\right)}_n^{t_{i+1}}-{\left(C_{j,k}\right)}_n^{t_i}}{\Delta t}=D\frac{{\left(C_{j,k}\right)}_{n+1}^{t_i}-2{\left(C_{j,k}\right)}_n^{t_i}+{\left(C_{j,k}\right)}_{n-1}^{t_i}}{{\left(\Delta z\right)}^2} \\ +\frac{{\left(V_{j,k}^{\downarrow }{C}_{j,k}\right)}_{n+1}^{t_i}+{\left(V_{j,k}^{\uparrow }{C}_{j,k}\right)}_{n-1}^{t_i}-{\left(V_{j,k}^{\downarrow }{C}_{j,k}\right)}_n^{t_i}-{\left(V_{j,k}^{\uparrow }{C}_{j,k}\right)}_n^{t_i}}{\Delta z} \end{gathered}$$

(5)

where $V_{j,k}^{\uparrow \&\downarrow }$ is the ascending or descending velocity of the particles in absolute value ($\left|V_{j,k}\right|$). Simplifying and rearranging leads to:

$$\begin{gathered} {\left(C_{j,k}\right)}_n^{t_{i+1}}={\left(C_{j,k}\right)}_n^{t_i}+\alpha \left({\left(C_{j,k}\right)}_{n+1}^{t_i}-2{\left(C_{j,k}\right)}_n^{t_i}+{\left(C_{j,k}\right)}_{n-1}^{t_i}\right) \\ +\beta \left({\left(V_{j,k}^{\downarrow }{C}_{j,k}\right)}_{n+1}^{t_i}+{\left(V_{j,k}^{\uparrow }{C}_{j,k}\right)}_{n-1}^{t_i}-{\left(V_{j,k}^{\downarrow }{C}_{j,k}+V_{j,k}^{\uparrow }{C}_{j,k}\right)}_n^{t_i}\right) \end{gathered}$$

(6)

with:

$$\alpha =\frac{D\mathsf{\Delta }t}{\Delta z^2} \& \beta =\frac{\mathsf{\Delta }t}{\Delta z}$$

where the value of the variables at time $i$ is known and where $i+1$ must be found from a forecast of the value at the next time step. The previous equation is valid only if the elements have identical volumes. When the elutriator has elements with different volumes (Figure 2), Equation (6) takes the more general form:

$$\begin{gathered} {\left(C_{j,k}\right)}_n^{t_{i+1}}={\left(C_{j,k}\right)}_n^{t_i}+\alpha \left(R_h{\left(C_{j,k}\right)}_{n+1}^{t_i}-2{\left(C_{j,k}\right)}_n^{t_i}+R_l{\left(C_{j,k}\right)}_{n-1}^{t_i}\right) \\ +\beta \left(R_h{\left(V_{j,k}^{\downarrow }{C}_{j,k}\right)}_{n+1}^{t_i}+R_l{\left(V_{j,k}^{\uparrow }{C}_{j,k}\right)}_{n-1}^{t_i}\right)-\beta \left({\left(V_{j,k}^{\downarrow }{C}_{j,k}+V_{j,k}^{\uparrow }{C}_{j,k}\right)}_n^{t_i}\right) \end{gathered}$$

(7)

with:

$$R_h=\frac{v_{n+1}}{v_n} \& R_l=\frac{v_{n-1}}{v_n}$$

where v_n represents the volume of element n. The next step in solving this set of algebraic equations is to make the following assumption:

Assumption 3.

Anything that leaves the hydraulic classifier cannot return.

Therefore, the material balances for the boundary elements are:

$${\left(C_{j,k}\right)}_1^{t_{i+1}}={\left(C_{j,k}\right)}_1^{t_i}+\alpha \left({\left(C_{j,k}\right)}_2^{t_i}{R}_h-2{\left(C_{j,k}\right)}_1^{t_i}\right)+\beta \left({\left(V_{j,k}^{\downarrow }{C}_{j,k}\right)}_2^{t_i}{R}_h-{\left(V_{j,k}^{\uparrow }{C}_{j,k}+V_{j,k}^{\downarrow }{C}_{j,k}\right)}_1^{t_i}\right),$$

(8)

$${\left(C_{j,k}\right)}_N^{t_{i+1}}={\left(C_{j,k}\right)}_N^{t_i}-\alpha \left(2{\left(C_{j,k}\right)}_N^{t_i}-{\left(C_{j,k}\right)}_{N-1}^{t_i}{R}_l\right)+\beta \left({\left(V_{j,k}^{\uparrow }{C}_{j,k}\right)}_{N-1}^{t_i}{R}_l-{\left(V_{j,k}^{\uparrow }{C}_{j,k}+V_{j,k}^{\downarrow }{C}_{j,k}\right)}_N^{t_i}\right),$$

(9)

$${\left(C_{j,k}\right)}_u^{t_{i+1}}=\alpha R_h{\left(C_{j,k}\right)}_1^{t_i}+\beta R_h{\left(V_{j,k}^{\downarrow }{C}_{j,k}\right)}_1^{t_i},$$

(10)

$${\left(C_{j,k}\right)}_o^{t_{i+1}}=\alpha {\left(C_{j,k}\right)}_N^{t_i}{R}_b+\beta R_b{\left(V_{j,k}^{\uparrow }{C}_{j,k}\right)}_N^{t_i},$$

(11)

In these equations the subscripts u and o indicate the underflow and overflow streams of the Hc, respectively. The solid balance is completed with addition of the feed rate. The method consists of introducing the volume of feed solids to be added at the beginning of each time step according to:

$${\left(C_{j,k}\right)}_{N_f}^{t_i}={\left(C_{j,k}\right)}_{N_f}^{t_i}+\frac{Q_{j,k,a}\mathsf{\Delta }t}{v_{N_f}}$$

(12)

where Q is the volumetric flow rate and the index α represents the feed stream. As shown in the work of Roy [24], particle settling might occur from the pulp before being introduced into the classifier so that the solid must be introduced in zone 2 and the water in zone 3 (Figure 3). In the industrial simulation case, the solid and the fluid are introduced in the same differential element.

3.2. Particle and Pulp Interstitial Velocities

The interstitial velocity can be estimated by considering it to be approximately equal to the hindered-settling velocity of the particles ($U_{j,k}$), to which the interstitial velocity of the pulp ($V_p$) is added:

$$V_{j,k}=U_{j,k}+V_p.$$

(13)

Several equations have been developed to estimate the hindered-settling rate of the particles. The method used in the algorithm to calculate the settling velocity of a single particle in water is that of Hartman et al. [26]. The conversion from free to hindered-settling velocity is performed with the equations of Masliyah [27] and the Richardson–Zaki index adapted as suggested by Rowe [28]. The three methods are explained in Appendix B.

Before calculating the motion generated by the displacement of the pulp, another assumption must be made. Thus:

Assumption 4.

To calculate the velocities of the water and the solid from the flow, we consider that the water flows mostly in the space not occupied by the solid, which is inversely the case for the solid too.

Therefore, the interstitial velocity of the pulp ($V_p$) is calculated by adding the interstitial velocity of water ($V_{w,n}$) to the interstitial velocity of the solid ($V_{s,n}$) as follows:

$$V_p=V_{w,n}+V_{s,n}$$

(14)

with:

$$V_{p,n}=\frac{Q_{w,n}}{a_n{\epsilon }_n}+\frac{Q_{s,n}}{a_n\left(1-{\epsilon }_n\right)}$$

(15)

and

$${\epsilon }_n=\left(1-{\phi }_n\right)=\left(1-{{\displaystyle \sum }}_{j=1}^J{{\displaystyle \sum }}_{k=1}^K{C}_{j,k,n}\right).$$

(16)

In these equations, the variables $a$, $\epsilon$ and $\phi$ represent, respectively, the area (normal to flow), the void fraction and the total solid volumetric concentration. The subscripts$w$ and $s$ designate water and solid, respectively. When the solid concentration reaches the maximum compaction concentration (see Section 3.3.2), the algorithm estimates the net particle velocity from the pulp volumetric flow rate as follows:

$$V_{j,k,n}=V_{p,n}=Q_p/a_n.$$

(17)

The next step in the resolution is to obtain an estimate of the overall motion of the pulp surrounding the particles ($V_{p,n}$). For this purpose, the hydraulic classifier is separated into three distinct zones according to the direction of the pulp flow, as shown in Figure 3.

The flow rates for each zone are obtained as follows:

$$Q_{zone 1}=Q_u-Q_c$$

(18)

$$Q_{zone 2}=Q_u-Q_c-Q_f=Q_a-Q_o$$

(19)

and

$$Q_{zone 3}=Q_o$$

(20)

To find the contribution of the solid and liquid phase, it is possible to estimate the volume concentration of underflow and overflow solids during the time step as being equal to the volume concentration of solids in the adjacent element at the beginning of the time step. This yields:

$$Q_{o,w}=Q_o\left(1-{{\displaystyle \sum }}_{j=1}^J{{\displaystyle \sum }}_{k=1}^K{C}_{j,k,N}\right)$$

(21)

and

$$Q_{u,w}=Q_u\left(1-{{\displaystyle \sum }}_{j=1}^J{{\displaystyle \sum }}_{k=1}^K{C}_{j,k,1 }\right)$$

(22)

where $Q_{o,w}$ and $Q_{u,w}$ are estimates of the water flow rate in the overflow and in the underflow.

3.3. Methods to Estimate Model Parameters

The use of a phenomenological simulation model of hydraulic classifiers presented in this study implies prior estimation of the axial dispersion coefficient ($D)$ and the porosity in the fluidized bed.

3.3.1. Dispersion Coefficient

Figure 4 shows the variation of the hematite partition curve as a function of $D$ (axial dispersion coefficient) during the simulations of an experiment for the industrial (Figure 4a) and laboratory (Figure 4b) scales. These figures suggest that a value of $D$ between 15–25 cm²/s yields the best fit for the industrial data, while a value of 4 cm²/s is more appropriate for the laboratory set-up. This latter value is very close to the 4–6 cm²/s reported by Kim [15]. The values obtained from the laboratory set-up also agree with the values found in the literature [18,19,20,21,22,29,30,31]. The higher $D$ describing the industrial set-up is expectable because a larger back-mixing prevailing in the larger sizes as well as the flow rate of treated solid per unit area is larger in the industrial unit (e.g., ~10 m/h) compared to that of the laboratory one (~4 m/h). From the simulation results, a method for estimating $D$ has been developed based on the value of the flow rate of treated solid per unit area. This method is presented in the work of Roy [24] and allows the simulation of the Hc when no experimental set-up is available.

Figure 5 shows the simulated axial profiles of the solid concentration for different values of $D$ against the axial profiles’ concentration of the plant measured data. The feed element is located at 50 cm from the top of the Hc and the teeter injector at 200 cm from the top. Figure 5 confirms that a low $D$ produces a variation in solid concentration as a function of depth that is not consistent with the observations made around the Hc feed zone. It also can be observed in Figure 5 that an increase in the value of $D$ decreases the slope of the solid concentration around the feed but increases the slope of the solid concentration in the fluidized bed zone. The difference between simulation and observation could be due to the inclusion of the dispersion term $(D\left({\partial }^{2}C/\partial z^2\right)$ in the computational algorithm to represent the pulp density re-equilibrium in the feed area due to Archimedes’ principle and that $D$ is kept constant over the height of the Hc (Assumption 2).

3.3.2. Methods to Estimate the Porosity of the Fluidized Bed

The Hc can be divided into 6 zones according to the characteristics of the pulp flow [16]. This division is shown in Figure 6. The water added to the feed increases the upward velocity of the fluid in zones A, B and C. This velocity is different from that in zones D, E and F. The teeter water causes the same type of variation (fluid velocity) between zones D and E and the collection zone F. As a result, the velocity of the water in the classifier varies according to the height in the device. This is also valid for the solid concentration (Figure 1), which also varies along the height of the apparatus, implying a variation in the density and viscosity of the pulp.

The minimum porosity or maximum volume concentration of compaction of the ore (zone F, Figure 6), before the underflow, is estimated from the results of sampling campaigns showing that the solid concentration in the underflow is always in the vicinity of 47% v/v. If this value is adjusted to account for the fact that the cone water exits with the slurry in the underflow, then the volumetric solid concentration in zone F is 67%.

The bed porosity in zones D and E is estimated according to a minimum fluidization velocity analysis for multi-species and multi-sized solid-containing pulps similar to that of Asif and Ibrahim in 2002 [32]:

$${\left(\overline{V_{mf}}\right)}^r={\displaystyle \sum }_{k=1}^K{X}_k{\left(V_{mf,k}\right)}^r$$

(23)

with $V_{mf,k}$ and $X_k$ which represent, respectively, the minimum fluidization velocity and the water-free solid volume fraction of species $k$. The exponent r can vary from $-1$ to $1,$ thus passing from harmonic to arithmetic mean. Asif and Ibrahim [32,33] obtained better results with $r=-1$, Rincon et al. [34] with $r=-0.5$, while in the present analysis the best estimation of bed porosity was obtained with $r=1$.

To determine the value of the porosity of the fluidized bed, Equation (24) was used. The simulator calculates the arithmetic averaging particle velocity ($V_{mean,n}^{\phi }$) as a function of total concentration ($\phi$) over a range of values from 30% to 70%. From these results, the simulator considers the maximal concentration in bed of particles to be the one where the average velocity is closest to zero.

$$\overline{V_n^{\phi }}\cong {\displaystyle \sum }_{j=1}^J{\displaystyle \sum }_{k=1}^K{X}_{j,k,n}{V}_{j,k,n}^{\phi }\approx 0$$

(24)

where $V_{j,k,n}^{\phi }$ is the hindered settling velocity of particles of dimensions $j$ of species $k$ when the total concentration of the element is equal to $\phi$.

4. Results and Discussion

4.1. Evaluation of the Calculation Algorithm

The simulation model will be validated with performance indices (Equation (A2), Appendix A.2) on a set of experimental data, including a total of 23 industrial experiments (1 to 23) and 8 laboratory experiments (L-1 to L-8). The experimental conditions are presented in

Table 2. Operating conditions of experiments on industrial, and laboratory set-ups.

Test	Type	Solid	Water	Teeter	Water	Pulp
		Feed Rate	Feed Rate	Water	Cone	Density
		[t/h]	[L/min]	[L/min]	[L/min]	[t/m3]
1	Indutrial	188	2998	1948	250	1.92
2	Ind.	188	3004	2191	250	1.85
3	Ind.	187	2997	1712	250	1.91
4	Ind.	188	3002	1943	250	1.89
5	Ind.	188	3000	1950	250	1.91
6	Ind.	188	3761	1950	250	1.88
7	Ind.	188	2334	1950	250	1.94
8	Ind.	188	2891	1950	250	1.96
9	Ind.	188	3003	1949	250	1.9
10	Ind.	173	2998	1949	250	1.76
11	Ind.	203	3001	1950	250	2.01
12	Ind.	188	3002	1950	250	1.97
13	Ind.	188	3003	1949	250	1.89
14	Ind.	188	3002	1950	250	1.77
15	Ind.	188	2998	1950	250	2
16	Ind.	188	2998	1950	250	1.96
17	Ind.	171	2961	1851	203	2
18	Ind.	171	2980	2199	201	2
19	Ind.	144	2510	2200	188	2
20	Ind.	144	2509	1701	186	2
21	Ind.	265	4672	2024	250	2
22	Ind.	90	2340	1000	237	1.77
23	Ind.	102	4000	960	250	1.77
L1	Laboratory	0.07	2.7	1.5	0.5	1.27
L2	Lab.	0.07	3	1.5	0.5	1.33
L3	Lab.	0.07	3.6	1.5	0.5	1.35
L4	Lab.	0.07	4	1.5	0.5	1.41
L5	Lab.	0.07	3.5	1	0.5	1.43
L6	Lab.	0.07	3.5	1.25	0.5	1.4
L7	Lab.	0.07	3.6	1.5	0.5	1.35
L8	Lab.	0.07	3.5	2	0.5	1.32

(and Figure A5 in Appendix C). Figure 7 shows the global underflow recovery of hematite (Fe₂O₃) and quartz (SiO₂) for all experiments. The simulation model provides a good estimate of the experimental values. It also demonstrates that the simulation algorithm can be adapted to different flow regimes (laminar/turbulent) and Hc geometry (size and shape). Furthermore, the simulation algorithm remains stable despite the presence of a discontinuity in the solid concentration of the pulp as a function of depth. This stability comes from the addition of an analysis loop that separates the particle velocity caused by convection into two variables ($V_{j,k}^{\downarrow }$, $V_{j,k}^{\uparrow }$) depending on their direction (Equation (5)).

4.2. Influence of Pulp Density

Hcs as physical separation devices generate a dense particle bed that offer the exploitation of the gravitational force and the density differences between solid particles and their surrounding environment. This makes quartz particles (SiO₂; specific gravity = 2.7) more dependent on the density variations of the pulp in zones D, E and F (Figure 6) than the hematite particles (Fe₂O₃; specific gravity = 5.2). However, the pulp density in these zones is not known as only the density of the slurry in the upper part of the Hc is measured, which includes zones A to D. Hence, lack of representation of the local density tends to inflate the error in simulating the behavior of light particles as their behavior in the pulp largely depends on the density of these zones.

Figure 8 shows the pulp density as a function of depth for three industrial-scale tests. Taken together, the experiments show that heavy particles behave significantly more predictably than light particles. This difficulty in adequately representing the displacement of light particles comes partly from the fact that the density of the pulp in the fluidized bed is close to the density of the quartz particles making all the particles of this species more dependent on the motion of medium then the dispersion force. Furthermore, as shown in Figure 5, the dispersion coefficient will be stronger in the feed zone than in the fluidized bed. Nevertheless, for simplicity, the dispersion coefficient has been assumed to be constant over the height of the apparatus (Assumption 2, Section 3.1). The dependence on the dispersion force of quartz makes its behavior slightly more difficult to predict.

4.3. Control Variables

The two main operating variables used to control the operation of an Hc against the natural disturbances (iron oxide composition and size distribution) are:

• The upper bed density as measured by a pressure sensor and manipulated via the opening of the Hc underflow valve;
• The flow rate of teeter water flow rate that is used to correct an underflow composition (e.g., SiO₂ content) that is off-spec.

Although these control variables are widely accepted as those that need to be manipulated to maintain a target operating performance for the Hc, it is somewhat surprising to find the void on the role of these variables in the literature dealing with the modeling of Hc. The proposed model provides a noise-free environment to study the impact of these variables on the performance of the Hc.

Figure 9 and Figure 10 show the values of the hematite and quartz weight recovery fraction at the Hcs underflow obtained experimentally and with the simulation model as a function of the different simulated operating variables, and as a function of the particle size. The figures show that to increase the underflow recovery, it is possible to decrease the flow of teeter water (experiments 1–4 and L1–L4); decrease the flow of feed water (experiments 5–8). When the unit is operated with a constant underflow valve opening value, it is possible to decrease the flow of the solids’ flow rate of the feed (experiments 9–12) or increase the opening of the underflow valve (experiments 13–16). While when controlling the bed density set point, it is necessary to increase the solids’ flow rate of the feed (experiments 17–20) or decrease the pulp density set point. The amount of light particle found in the underflow can be decreased by increasing the density of the pulp in the fluidized bed to inhibit the passage of light particles (experiments 5–8 and 13–16) or with decreasing the water flow rate at the top of the apparatus to reduce the loss of heavy particles in the overflow.

In summary, Figure 9 and Figure 10 show that the effect of the fluidization water is well captured by the simulator as well the effect of the bed density or valve opening, and this capability is used in [24] for studying control strategies that can reduce the variability reported in Figure 3 and as a way to define strategies to minimize iron losses while maintaining a low silica content in the Hc underflow.

4.4. Simulation of Primary, Secondary and Tertiary Hc

One of the objectives of the model development was to establish a mathematical model for the Hc that can be used to simulate an Hc without requiring calibration if the characteristics of the feed material change. The ability of the model to represent the operation of the primary Hc is demonstrated in Section 4.1, but whether this model could be used or not to predict the performances of secondary and tertiary Hcs (see Figure A2, Appendix A.2) remains to be checked. Finally,

Table 3. Observed and simulated performance indices for the primary, secondary and tertiary Hcs of the ArcelorMittal pellet plant.

Hc	Test	Solid	Water	Teeter	Water	Pulp	Weight Recovery (%)		% SiO2 in Hc Underflow
		Feed	Feed Rate	Water	Cone	Density	Hematite	Quartz
	#	[t/h]	[L/min]	[L/min]	[L/min]	[t/m3]	Obs.	Sim.	Obs.	Sim.	Obs.	Sim.
Primary	21	265	4672	2024	250	2.00	83.7	83.0	35.2	35.0	2.1	2.1
Secondary	22	90	2340	1000	237	1.77	25.6	26.9	1.9	3.6	2.6	4.5
tertiary	23	102	4000	960	250	1.77	40.0	43.9	2.3	1.4	3.0	1.7

gives the observed and simulated global performance indices for the primary, secondary and tertiary Hcs (see Figure A2, Appendix A.2). These results were obtained without altering calibration of the original model. The only input variables that were changed from one simulation to the other are the solid feed rate, the Hc feed mineral contents and the particle size distribution, and the dimensions (diameter and height) of the Hc. The Hc operating conditions (teeter water flow rate and bed density) applied for the simulation are those used in the plant. Despite the differences in the applications of the three Hcs, the model predictions are close to the Hc actual performances confirming that the fundamental approach used to develop the model ensures robustness and flexibility to simulate the behavior of various materials when processed in an Hc.

4.5. Dynamic Simulations

The model is built to allow the dynamic (time variation) simulation of the Hc response to step changes in the main Hc input variables. Details and other dynamic simulations can be found in [24]. An example of a dynamic simulation in response to step changes in the fluidization water flow rate is shown in Figure 11. This dynamic simulation example (Figure 11a) follows the same operating conditions as the industrial experiments 1 to 4 (Figure 11b). The iron oxide content in the Hc overflow is used as the output variable. The curve superimposed on the sampled values in Figure 11b was produced with an empirical simulation model proposed in the work of Payenzo [25], while that in Figure 11b was produced with the phenomenological simulation model presented in the present work. The simulated response is consistent with the observed response shown in Figure 11b where it can be seen that the simulator is able to predict the observed trend on the step response. This good prediction ability was possible without any adjustment to the calibrated model, confirming that the simulator could be used to virtually assess control strategies. Part of the observed error between the simulation and the experimental values comes from the fact that the experiments were carried out with a constant valve opening and that the simulation was carried out by modifying the underflow rate to follow the bed density obtained in these experiments. In addition, an error is contributed by the presence of particles that are not fully released and therefore have intermediate densities between heavy and light particles.

5. Conclusions

A one-dimensional dynamic model based on phenomenological convection/diffusion transport equations has been developed to simulate the operation of industrial Hc and of a laboratory for the separation of hematite from quartz. The model is able to predict both the steady state and dynamic responses of an industrial Hc to changes in teeter water addition and bed density.

The use of the assumptions presented in Section 3 allows to simplify the solution of the system of equations (Equations (2), (13), (21) and (22)) and allows to obtain the versatile simulation model. However, the study shows that assumption 2 is not always representative of the real behavior of the pulp (Section 3.3.1, Figure 5) and would require further investigation of the dispersion behavior as a function of depth under ranges of operating conditions similar to those of the experiments presented in Table 2.

In this research, it was possible to observe that the simulation of quartz is more subject to prediction errors than hematite. In addition to the error generated by the very low content of quartz in certain size classes and its non-uniform size distribution discussed in the work of Roy [24], it is possible to observe in Section 4.2 that the density of quartz is very close to the density of the fluidized bed, which increases the chances of error in the prediction of behavior and increases the influence of the dispersion at the expense of convection. Moreover, part of the error generated by the simulation model comes from the fact that the simulation algorithm does not take into account the particles that have not undergone a complete liberation, thus an intermediate density.

Nevertheless, the simulator has helped, for instance, to verify the following patterns and trends.

• To increase the underflow recovery, it is possible to:
  • Decrease the flow of teeter water or feed water;
  • Increase the opening of the underflow valve or decrease the pulp density set point;
  • Decrease the solids flow rate of the feed with a constant underflow valve opening;
  • Increase the solids flow rate of the feed while maintaining a constant slurry feed pulp density.
• To decrease the yield of light particle in the underflow, it is possible to increase the density of the pulp in the fluidized bed to help hinder their buildup in the underflow.
• To decrease the loss of heavy particles in the overflow, it is possible by decreasing the water flow rate at the top of the apparatus without disturbing the fluidized bed under the feed by reducing the amount of water introduced by the feed (when it is allowed).

The simulation algorithm proves to be an asset to help in the choice of operating conditions by allowing the testing of various adjustments without disturbing the industrial process. The model is flexible as it allows the representation of variations in the partition curves of hematite and quartz for industrial and laboratory applications.