Modelling of the Separation Process in a Ferrohydrostatic Separator Using Discrete Element Method

This paper presents a model of the separation process in a ferrohydrostatic separator (FHS) which has been designed and developed at DBGS, De Beers, South Africa. The model was developed using special discrete element method software package called Particle Flow Code in 3D (PFC). Special attention has been paid to the selection of the simulation parameters in order to achieve the required feed rates. The simulation was carried out using spherical particles and density tracers of different sizes and densities ranging between 0.004 and 0.008 m and 2700 and 3800 kg/m3, respectively. The tracers were used to set the apparent density of the ferrofluid (the cut-point) and to provide a measurement of the efficiency of the separation. The model is replacing the ferrofluid by imposing a drag force on the particle. The results of the simulation were presented in the form of the distribution of the density tracers into the sink fraction. These results are realistic and show the advantages of DEM to understand the complex flow behavior of granular materials.


INTRODUCTION
A ferrohydrostatic separator (FHS) features a stationary magnetic fluid, housed in a separation chamber placed between pole-pieces of a magnet.Magnetic fluid-based separations were pioneered in the early 1960s [1].More recently, numerous ferrohydrostatic separators have been designed for separation of nonferrous metals [2,3], scrap metal [4], coal preparation [5], for recovery of gold, platinumgroup metals and diamonds [6], biologically active compounds [7].
A single nonmagnetic particle suspended in a static ferrofluid is acted upon by force of gravity and by two buoyancy forces.The first is the classical Archimedes gravity-related force, and the other is the magnetically induced buoyancy force due to the magnetic "weight" of the magnetic fluid as introduced by Rosensweig [8].These forces acting on a nonmagnetic particle will force the particle to move along a certain path and the trajectory of the particle can be computed as it was described in [9].
A particle of radius b (small enough to assume that the magnetic field gradient is uniform within the particle), volume V p , and density ρ p immersed in a magnetic fluid of density ρ f , dynamic viscosity η and saturation magnetization M f in the presence of an external magnetic field gradient ∇H will experience the following forces (see Figure 1): (i) force of gravity (ii) buoyancy force (iii) the magnetically derived buoyancy force where μ 0 = 4π × 10 −7 Tm/A is the magnetic permeability of vacuum, (iv) drag force (written for the regime where Stokes Law is applicable and with the magnetic fluid at rest) Equations ( 1)-( 4) are written in SI units.Besides, this system of units is used throughout.The bouyancy forces can be added together and the socalled apparent density of the ferrofluid can be defined as:  α being the angle between the direction of the magnetic field gradient and the horizontal direction.In, for example, (5) the magnetization and the magnetic field intensity are expressed in A/m. Figure 2 depicts trajectories of nonmagnetic particles (of diameter 0.003 m) of different densities (given in kg/m 3 ) as they were computed in [9].The values of the parameters used in calculation are the apparent density of the ferrofluid ρ app = 3500 kg/m 3 , the ferrofluid density ρ f = 990 kg/m 3 , the ferrofluid viscosity η = 1.65 × 10 −3 kg m −1 s −1 , the angle α = 2 • and drop height h = 0.003 m.It can be observed that a particle having the density greater than the apparent density will follow a sink trajectory while a particle of density lower than the apparent density will float.Besides, a particle having the same density as the apparent density of the fluid is sinking.
It is obvious that this model based on a single particle is not suitable for an assembly of particles interacting with one another.The interactions can deflect the particle from its correct path and thus can lead to the misplacement of the particles.A single particle model can be used, eventually, for the calibration process, where single particles are introduced into the fluid.
A software package which can model the interaction between particles can be successfully used to model the separation process.Such a software is Particle Flow Code in 2D and 3D (PFC 2D and PFC 3D ) developed by Itasca Consulting Group Incorporation, USA.
The 2D model of the De Beers FHS unit has been presented in a previous paper [19] and has been developed for the vertical plane situated in the middle of the chamber as it is shown in Figure 3.In that paper, the results of the 2D model have been compared with the experimental findings, but in an average way and by using no density tracers of certain densities.In that way the theoretical results were close to the average experimental ones.A full validation of the model must take into account every single density tracer and thus, to monitor the behaviour of the tracers of each density and each size.
On the other hand, the 2D model has used a simple approach in terms of the feeding system.Particles were generated right under the surface of the ferrofluid eliminating thus the randomness of a feeding system consisting of a vibrating feeder tray.
In the 2D model the feed rate has been calculated using the same number of calculation steps between two consecutive particles generated under the surface of the ferrofluid.This time has not been constant knowing that the time step in PFC is not constant.
The most important fact that affects the separation process is the agglomeration of the particles above the baffle plate.In the 2D model this phenomenon cannot be very well modelled using the representation of the particles as discs.
The experimental work carried on the FHS unit has shown that the particles had the tendency to migrate towards the middle of the chamber.This migration can increase the agglomeration of the particles within the chamber and thus, affects the efficiency of the separation.This fact can be explained based on the magnetic field pattern generated by the electromagnet within the gap.This phenomenon could not be included in the 2D model of the separator, too.
This paper presents the 3D model of the separation process in the FHS machine using the PFC 3D software package.The 3D model has been developed in order to overcome the weak points of the 2D model presented above.The simulation results, presented as distribution curves, are realistic and show clearly the advantages of using DEM to model the behaviour of granular materials.
Once the model is developed and validated against experimental results, then it can be used to redesign the separator in order to improve the separation efficiency of the unit.

GENERAL FEATURES OF THE PFC 3D PACKAGE
PFC 3D models, in 3D, the dynamic behaviour of assemblies of arbitrarily-sized spherical particles.A particle generator allows the statistical generation of specified distributions of particles to be done automatically.Particles radii may distribute uniformly or according to a Gaussian distribution.Properties are associated with individual particles (like initial velocity) or contacts.Any type of force can act on the particles, for example, drag force, gravity force, friction force, and so forth.Any number of arbitrarily oriented planes may be specified as walls, each with its own set of contact properties.These walls represent usually boundary conditions which force the particles to move only inside a certain volume.The most important feature of the software is the contact model.PFC 3D allows assigning a contact model to each particle or wall, which describes the physical behaviour occurring at a contact.

THE DEVELOPMENT OF THE MODEL
The simulation process has been developed in several steps, each step being described in the following sections.

Geometry of the chamber
The design of the chamber for F&L machine is shown in Figure 4.The model will be used to simulate the behaviour of the particles inside the chamber.Based on the chamber design and the position of the ferrofluid pool inside the chamber, the walls used in the modelling are the following: bottom of the chamber, the splitter plate, the baffle plate, and the side walls.Besides these walls, the ferrofluid surface will be another wall used in simulation.A visual representation of the walls is reproduced in Figure 3.
The entire geometry of the unit was built in 3D using POV-Ray (www.povray.org).POV-Ray creates threedimensional, photo-realistic images using a rendering technique called ray-tracing.This package is far more flexible and produces more quality images that built in tool within PFC 3D .A picture of the unit obtained with POV-Ray is shown in Figure 5(a).

Modelling of the magnetic field in the electromagnet gap and calculation of the apparent density of the fluid
In order to calculate the apparent density of the ferrofluid, a model of the electromagnet was developed.The electromagnet was modelled using Faraday, software for electromagnetic designs produced by Integrated Engineering Software, Canada.The FHS electromagnet design is presented in Figure 5(b), as it was modelled in Faraday.
Unlike the 2D model of the separation, where the separation has been modelled within a plane situated in the middle of the gap, the 3D model must take into account the third component of the magnetic field, oriented from one pole to the other.It has been observed experimentally that this component is responsible for the movement of the particles toward the middle of the electromagnet gap.The basis for this behaviour is that the apparent density of the ferrofluid is greater close to the pole tips than at the middle of the gap.
In order to quantify the contribution of the third component of the magnetic field to the separation process, it is necessary to map the magnetic field within a plane that contains the cross-section area of the gap.This plane is presented in Figure 6.
The mapping of the magnetic field was done by dividing the area between the pole pieces in small squares Physical Separation in Science and Engineering (0.005 × 0.005 m).Within each square the magnetic field intensity is considered constant.In other words, within each small square the apparent density of the ferrofluid is constant.The resulting grid is also shown in Figure 6.
The Faraday model of the electromagnet can provide the values of the magnetic field intensity and its gradient along the vertical grid lines shown in Figure 6.The apparent density of the fluid ρ app can be then calculated using the following equation: where ρ ff is ferrofluid density, μ 0 represents the magnetic permeability of vacuum, M ff the magnetization of the ferrofluid, ∇H the gradient of the magnetic field intensity and g the gravitational acceleration.
The apparent density can be used in the DEM model to calculate the magnetic buoyancy force acting on the particles.
The profile of the apparent density along the vertical grid lines is shown in the plot presented in Figure 7.The first plot represents the variation of the apparent density along the vertical line situated in the middle of the gap.The other plots represent the variation of the apparent density along vertical lines situated at different distances from the middle of the gap.
The modelling and the calculation of the apparent density prove that the apparent density is constant along the ver- tical lines only in a limited region.It can be also noticed on the plot that the apparent density increases sharply for the last four lines.This is due to the fact that the vertical lines are crossing the pole pieces of the electromagnet.At the same time, the Faraday model of the electromagnet shows that the apparent density increases in regions situated close to the pole pieces of the electromagnet.This is consistent with the observation that the particles are pushed toward the middle of the gap where the apparent density is lower than close to the pole pieces.
The points on each plot are fitted using an equation of the following form (see Figure 8): The values of the parameters y 0 , a, b, and c are different for each vertical line situated at different distances from the middle of the gap.Equation (7) combined with the values of y 0 , a, b and c were used to calculate the magnetic buoyancy force within the cross-section area of the gap.The apparent density inside the electromagnet gap also changes as a function of the current density in the coils.Using different values of the current density in the coils it is possible to change the cut-point density of the ferrofluid.By changing the current density in the coils the apparent density will still be given by (7), although the values of the parameter y 0 will be different for different current densities.This procedure will be used to set the cut-point density (apparent density) in the model.

Forces used in the model
The following forces were used to model the separation process: (i) gravity force (ii) buoyancy force (iii) magnetic buoyancy force (iv) drag force In the above equations ρ p represent the density of the particle, V p is the volume of the particle, η is the dynamic viscosity of the fluid, b is the radius of the particle and v p is the velocity of the particle.

Modelling of contact forces
In order to calculate contact forces between contacting entities, the model used the PFC embedded Hertz-Mindlin contact model for particle-particle contact [20,21] and the spring-dashpot contact model for particle-wall contact [22].The Hertz-Mindlin contact is fully described by the shear modulus G and Poisson's ratio ν of the material the particles are made off.The contact between particles and walls is defined by the normal and tangential stiffness of the string k n and k s and the normal and tangential damping parameters of the dashpot c n and c s .Both contacts models include a slip model and is defined by the friction coefficient at the contact μ.Another parameter used was the local damping or particle damping α [23].This parameter is used as a measure of the energy lost through frictional sliding.Its default value was 0.7 and it has been proved to affect significantly the vibration movement of the particles on the feeder tray.The values of these parameters as they were used in the model are shown in Table 1.
It must be noticed here that the values of the damp parameters were chosen from a range of a values as a function of the cut-point density.These values were modified until the model replicated the experimental cut-points.This will be explained in detail in the next section of the paper.

METHOD AND RESULTS
The model has been used to simulate the separation process in a FHS unit at two different cut-points and five different feed rates.The simulation was carried out in three steps.
(1) The cut-point setting, where two different sets of particles of certain densities were dropped in ferrofluid to check the apparent density of the ferrofluid.(2) The feed-rate setting, where the hopper was filled with particles, the feeder tray was vibrating and the discharge time was computed.(3) The separation process itself.
The above steps are described in detail in the following sections.

Cut-point setting
In order to set the cut-point, particles of 0.006 m size were used.The setting was done by generating forty density tracers of the same density (±3% deviation from standard value due to the manufacturing process) onto the vibrating feeder tray.The behaviour of the tracers was monitored and the percentage of the tracers reported into the sink fraction was calculated and printed on the screen.The two cut-points set were 3250 kg/m 3 and 3400 kg/m 3 , respectively.In the case of the first one, the tracers of density 3200 kg/m 3 must float while those of density 3300 kg/m 3 must sink.Three images taken during the simulation are presented in Figure 9.The first one shows the majority of the tracers of density 3200 kg/m 3 floating, the second one shows most of the 3300 kg/m 3 Physical Separation in Science and Engineering sinking and the third one the 3400 kg/m 3 split between sink and float fractions.In this case, the values of the viscous damp parameters were 1.4 for the normal one and 0.8 for the shear one and the particle damp inside the ferrofluid 0.7.The apparent density of the ferrofluid was increased by increasing the y 0 value in (7).
The above figures show that the desire cut-point can be achieved by a proper selection of the DEM parameters.This can be achieved by a trial and error procedure.The parameters such as: particle damp, y 0 (for apparent density in ( 7)), normal and shear viscous damp are the main parameters which can be modified to set the cut-point.

Feed-rate setting
Once the cut-point has been set, the next step was to establish how different parameters affect the setting of the feed rates.In order for the model to be validated, in the future, against experimental results and then to be used as an optimization tool it will be required to set a certain feed rate which is close to the feed rates used on the real machine.The main parameters which can affect the feed rate are frequency of the vibration of the feeder tray, the amplitudes of the vibration (horizontal and vertical), and the particle damp.
The study of the influence of each of these parameters on the feed rate was conducted using 10 000 gravel of size between 0.004 and 0.008 m and density between 2700 and 3800 kg/m 3 .The particles were generated randomly within the hopper.The total mass of the particles generated was 4.5 kg.The particles were left to settle for few seconds, then the gate of the hopper was opened and the feeder tray started to vibrate.The movement of the feeder tray was modelled as combined oscillations along the vertical and horizontal directions.
The feed rate was computed as the ratio of the total mass of the particles and the time required to empty the hopper It was noticed from these simulations that the particle damp was the most sensitive parameter in terms of feed-rate setting.Therefore, was decided to keep constant the other three parameters (frequency and amplitudes) and find a relationship between feed rate and particle damp.The particle damp is used to dissipate energy by effectively damping the equation of motion.Nine simulations were run at different values of the particle damp and the values of the feed rates are presented in Figure 11 as a function of the particle damp.It can be noticed that a greater value of the particle damping implies a lower feed rate as less energy is transferred from the vibrating feeder tray to the particles.The constant parame- ters were set at the following values: (i) frequency of the vibration of the feeder tray f = 5 Hz, (ii) the amplitude of the vibration along the horizontal direction A x = 0.0005 m, (iii) the amplitude of the vibration along the vertical direction A z = 0.0007 m.
The points were fitted with a curve of the equation y = y 0 + ax + bx 2 and the values of the parameters obtained were The curve obtained can be used to set the desired feed rate.Table 2 shows five required values of the feed rate and the corresponding values of the particle damp computed based on the fitted curve.It must be noticed that this equation is valid only for f = 5 Hz, A x = 0.0005 m and A z = 0.0007 m.
After the value of the particle damp was obtained using the fitted curve, a simulation was run at that specified particle damp and the theoretical feed rate computed.These values are also presented in Table 2.
The last feed rate, namely, 522 kg/h, was not achieved under above conditions even for a very low particle damp 0.01.In order to reach this feed rate, the particle damp was set at 0.01 and the other parameters related to the vibration of the feeder were set as follows: (i) frequency of the vibration of the feeder tray f = 15 Hz, (ii) amplitude of the vibration along the horizontal direction A x = 0.0008 m, (iii) amplitude of the vibration along the vertical direction A z = 0.0012 m.
Using the above values, the theoretical feed rate was 536 kg/h.It must be specified here that these values of the particle damp were used only for the particles situated on the feeder tray or inside the hopper.As soon as the particles leave the feeder tray the damp was changed to zero and then once the particles penetrate the ferrofluid surface the damp value was changed again in order to take into account the viscosity of the ferrofluid.

Simulation of the separation process
The previous two steps were necessary to set the apparent density of the ferrofluid and the desired feed rate.The previous used sample of 10 000 particles was mixed with density tracers.The tracers were of three different sizes (0.004, 0.006, and 0.008 m) and density range from 2800 to 3800 kg/m 3 in steps of 100 kg/m 3 .There were 50 tracers for each size and each density, the total number of tracers being 1650.Thus, the total mass of the particles and tracers was 4.81 kg.The main roll of the tracers was to provide an indication of the separation efficiency and whether the separation process is size dependent.
The efficiency of the separation can be quantified by plotting the distribution curve of the density tracers into the sink fraction.An example of such curve is presented in Figure 12.Two parameters can be computed from the distribution curves and their values can give a measurement of the efficiency of the separation process.These parameters are (i) actual cut-point CP 50 defined as the density at which the split is 50%-50%, (ii) Ep value defined as the ratio: Ep = (ρ 75 − ρ 25 )/2, where ρ 25 and ρ 75 are shown on Figure 12.
The meaning of these parameters are as follows.
(i) The value CP 50 gives a measurement of the shift in the actual cut-point compared to the set one.Any significant shift of the cut-point leads to a different masssplit and an increase number of misplaced particles.(ii) The E p value represents actually the inclination of the curve and gives an indication of the sharpness of the separation.The lower this value, the better the accuracy of the separation.
Figure 13 shows snap-shots of the separation process at two different cut-points and two different feed rates.

Simulation results and discussions
The simulations were run on a PC equipped with an AMD Athlon64 3.4 GHz processor, 1.5 GB RAM and one simulation took about one day.The results of the simulation were saved as percentages of the tracers reported into the sink fraction grouped by their sizes and densities.In order to examine the influence of the size of the particles on the efficiency of the separation, a distribution curve was plotted for each size of the tracers.At the same time, five different feed rates were used to study the influence of the feed rate on the separation process.Figure 14 presents the distribution curves obtained at CP = 3250 kg/m 3 and different feed rates for all the density tracers.Figure 15 shows the same results but at cut-point CP = 3400 kg/m 3 .
In order to compare the efficiency of separation among different size of the tracers and also to study the influence of the feed rate on the quality of the separation, the two parameters CP 50 and E p are computed from above theoretical curves.These values are shown in Tables 3 and 4.
It can be noticed from the curves and the tables presented above that the following hold.
(1) The separation process is size-dependent-at the same feed rate, the distribution curves are different for different sizes of the tracers.This can be seen by comparing, for instance, the curves at feed rate of 85 kg/h in Figure 14 or the corresponding values in Table 3  (2) The efficiency of separation deteriorates as the feed rate increases.Comparing the curves at 85 kg/h with the ones at 536 kg/h in Figure 15, it can be observed that at lower feed rate all the larger tracers of 3800 kg/m 3 are reported into sink fraction (which means no misplacements), while at greater feed rate the percentage of these tracers reported into sink fraction is below 100%.At the same time, the values of CP 50 and E p show an ascendant trend up to the point where these values cannot be computed due to a complete inefficient sep-aration.The increase in the value of CP 50 is an indication of increased number of the misplaced less dense tracers into the sink fraction.
(3) The separation is more effective at CP = 3250 kg/m 3 than at CP = 3400 kg/m 3 .Comparing the curves at the same feed rate, for instance 344 kg/h, but different cut-points it can be noticed that at CP = 3250 kg/m 3 the larger tracers of density 3800 kg/m 3 are all reported into the sink fraction, while at CP = 3400 kg/m 3 the same tracers are only 90% reported into the sink fraction.The tables also show that, in the case of the CP = 3400 kg/m 3 at 536 kg/h, the parameters CP 50 and E p could not be computed, while at the same feed rate but at CP = 3250 kg/m 3 these parameters were calculated.This proves that at CP = 3400 kg/m 3 and 536 kg/h the separation is nonexistent.

Physical Separation in Science and Engineering
All these observation have a very good explanation based on the design of the separation chamber and the physics of the separation process itself.The size-dependence can be explained by the fact that the forces acting on the particles are volumetric forces.A larger particle will experience a larger buoyancy magnetic force and therefore will find its path easier than the smaller ones.At the same time, a larger particle will knock the smaller ones and will force them to move into the wrong fraction.
The feed rate has a huge impact on the efficiency of separation mainly because of the agglomeration of the particles on the baffle plate.This agglomeration will lead to an increase number of misplaced particles as the feed rate increases.
The better separation at CP = 3250 kg/m 3 compared to CP = 3400 kg/m 3 can be explained by the mass split of the sample at these two cut-points.Thus, at CP = 3250 kg/m 3 the mass split is 50%-50% into sink and float fractions (considering that the density distribution of the particles is uniform), while at CP = 3400 kg/m 3 the mass split is 65% into float fraction and 35% into sink fraction.In the first case, the probability of the misplacement of a particle due to entrapment is lower compared to the second case when majority of the particles are floating.This upward stream of particles can easily entrapped heavier particles and forced them to report into float fraction instead of sink fraction.The distribution curves in Figure 15 show very well that even at very low feed rates the percentage of the 3800 kg/m 3 tracers into sink fraction is well below 100%.

CONCLUSIONS
This paper has presented the 3D DEM model of the FHS unit developed with Particulate Flow Code.The previous model, developed in 2D, has been validated against the average experimental results and using only gravel in the model.The 3D model presented in this paper has overcome most of the 2D shortcomings by (i) including a proper feeding system consisting of a hopper and a vibrating feeder tray; (ii) simulating a mixture of gravel and density tracer of different sizes.In this way the simulation could investigate the influence of the size of the particles on the separation efficiency; (iii) including the third component of the magnetic force.
This force pushes the particles towards the middle of the chamber generating agglomeration of the particles with direct impact on the efficiency of separation; (iv) improving the feed-rate calculation.
The 3D model has shown that the larger tracers and implicitly the larger gravel are better separated based on their densities than the smaller ones.These particles are experiencing larger buoyancy magnetic force that the smaller ones and therefore are finding their correct path faster than the small ones.
The simulations run at different feed rates and the same cut-point helped us to study the influence of the feed rate on the efficiency of the separation.The separation is efficient only up to a certain feed rate.This threshold value can be calculated once a criterion is established.For instance, such a criterion could be the E p value smaller than 0.1, which is quite often used in densimetric separation.
The above observation can be explained by the fact that the particles start to agglomerate above the baffle plate and will not have enough space to follow the correct path.The model has shown this effect occurring at higher feed rates.
Another interesting observation was that the separation is more efficient at the cut-point equals to 3250 kg/m 3 that at the cut-point 3400 kg/m 3 .This can be explained based on the mass split of the material.Due to the fact that the density distribution among the gravel was a randomly uniform distribution between 2700 and 3800 kg/m 3 , at CP = 3250 kg/m 3 the mass split was around 50%, while at CP = 3400 kg/m 3 the split was uneven creating a predominate flux of particles moving upwards.This flux entrapped many heavily dense particles which were misplaced into the float fraction.This can be easily observed on the last plot in Figure 15 where only 40% of the 0.004 m tracers of density 3800 kg/m 3 were reported into the sink fraction.
The model does not take into account the shape of the particles and the coupling between the fluid and particles.In order to prove whether these factors affect the results of the simulation the model must be validated against experimental results.The validation exercise will be the subject of a future paper.

2 PhysicalFigure 1 :z
Figure 1: Section through the magnetic fluid bridge of depth δ, inclined with an angle α.

Figure 2 :
Figure 2: Trajectories of particles of different densities for h = 0.003 m.The dashed straight lines represent the trajectories of particle of ρ p = 3500 and 3550 kg/m 3 for h = 0.

Figure 3 :
Figure 3: The electromagnet and the plane in which the 2D simulation took place.

Figure 4 :
Figure 4: Chamber design and the walls for the model.

Figure 5 :Figure 6 :
Figure 5: (a) POV-Ray image of the FHS unit; (b) Faraday model of the FHS electromagnet.

Figure 7 :
Figure 7: Profiles of apparent density along the vertical grid lines within the cross-section area of the gap.

3 ) 3 Figure 8 :
Figure 8: Apparent density along a vertical line and the fitted curve.

Figure 11 :
Figure 11: Variation of the feed rate as a function of particle damp when f = 5 Hz, A x = 0.0005 m and A z = 0.0007 m.

Figure 12 :
Figure 12: Definition of the parameters computed from the distribution curve.

Figure 14 :Figure 15 :
Figure 14: Distribution curves of the density tracers at different feed rates and CP = 3250 kg/m 3 .

Table 1 :
Material properties and contact parameters.

Table 2 :
Particle damp values and the corresponding feed rates.