Hydrodynamic Modeling of Swirling Flow and Particle Classification in Large-Scale Hydrocyclones t

The fluid flow in a hydrocyclone is extremely challenging from the modeling point of view, because the flow field is a vortex flow with flow reversal. Furthermore, there exists a great need to predict the particle size classification efficiency in this device. The mining and chemical industries use hydrocyclones in large numbers wherever there is need to separate fines from coarse particles. Due to the complexity of the flow, early researchers took the empirical approach of relating key operating variables with classification efficiency. Although the empirical models were very successful in size classification predictions, they could not be used universally since such models are confined to the bounds of data from which they were developed and their applicability diminished even further when models developed at a given site were used in other situations. Consequently, models based on the physics of fluid flow began appearing in the literature, from those incorporating very simple ideas at first to very detailed models requiring numerical solutions in recent years. In this manuscript we present such a model which uses all the known principles of fluid flow as much as possible, leaving only the turbulence closure to empiricism, as this phenomenon is still unsolved in the larger context of fluid mechanics. The complete fluid-flow based model lends itself as a computer-aided design software for hydrocyclone, not merely confined to the traditional cylinder plus cone design but to explore newer geometries as well. In fact the traditional design suffers from a few inherent flaws. The numerical modeling approach is well suited to examining with ease the classification performance of conceptual designs. For this reason the model must be perfected to the point where it can predict traditional industrial hydrocyclones for different slurry properties and device designs.


Introduction
The fluid flow in a hydrocyclone is extremely challenging from the modeling point of view, because the flow field is a vortex flow with flow reversal. Furthermore, there exists a great need to predict the particle size classification efficiency in this device. The mining and chemical industries use hydrocyclones in large numbers wherever there is need to separate fines from coarse particles. Due to the complexity of the flow, early researchers took the empirical approach of relating key operating variables with classification efficiency. Although the empirical models were very successful in size classification predictions, they could not be used universally since such models are confined to the bounds of data from which they were developed and their applicability diminished even further when models developed at a given site were used in other situations. Consequently, models based on the physics of fluid flow began appearing in the literature, from those incorporating very simple ideas at first to very detailed models requiring numerical solutions in recent years. In this manuscript we present such a model which uses all the known principles of fluid flow as much as possible, leaving only the turbulence closure to empiricism, as this phenomenon is still unsolved in the larger context of fluid mechanics.
The complete fluid-flow based model lends itself as a computer-aided design software for hydrocyclone, not merely confined to the traditional cylinder plus cone design but to explore newer geometries as well. In fact the traditional design suffers from a few inherent flaws. The numerical modeling approach is well suited to examining with ease the classification performance of conceptual designs. For this reason the model must be perfected to the point where it can predict traditional industrial hydrocyclones for different slurry properties and device designs.

Background Literature
Fortunately, the physics of fluid flow in hydrocyclones has been the subject of investigation in the last decade as reviewed by Chakraborti and Miller (1992), and so the key ideas in modeling have been forthcoming steadily. Leaving aside the early investigations of measuring fluid-flow with mechanical aids and dyes, and also early works which had some semblance to the Navier-Stokes momentum balance equations, we come to the work of Bloor and Ingham (1987), who gave an analytical solution to the Navier-Stokes equations, albeit for overly simplifying assumptions. In the region near the central axis, the vortex conservation was applied, with inviscid and rotational flow assumptions, which yields axial and radial components. In the region along the wall the boundary-layer approach was used to derive velocities. Despite Bloor and Ingham's efforts in predicting velocity data measured by Kelsall (1952), a data set often referred to in numerous publications, the lack of an adequate turbulence description led others to revise the analytical solutions as necessary. The first successful work in predicting the fluid flow in hydrocyclones is that of Pericleous and Rhodes (1986), who used the PHOENICS computer code for the solution of the partial differential equations. Using the simple Prandtl mixing length model and the axisymmetry assumptions, the authors reported the velocity predictions in a 200-mm hydrocyclone. Later, Hsieh and Rajamani (1991) numerically solved the turbulent momentum equations to obtain the velocities and compared them with the Laser Doppler V elocimetry measurements in a 75-mm hydrocyclone. The success of this comprehensive work may be attributed to the Prandtl mixing-length model assumed in this work, to be described later in this manuscript. This work showed that, by a simple balance of forces on a particle present in the flow field, the trajectory of the particle can be traced inside the hydrocyclone, from which the entire size-classification efficiency is computed. In a sequel, Monredon, Hsieh and Rajamani (1992) showed that the same model is evidently applicable even if the vortex finder and spigot dimensions were changed drastically besides the operating conditions. All of these modeling works have been confmed to hydrocyclones processing slurries in the 5 to 10 o/o solids range, and so Rajamani and Ludovic (1992) took a further step in the direction of higher percent solids, reporting that volume concentration of a particle in a computational cell is proportional to the residence time of that particle within that cell. This idea extended the numerical model to an iterative scheme, with which the spatial concentration distributions were predicted. Notable among the numerical modeling works is that of Davidson (1994), who modeled the multiphase flow by treating the fluids and solids as interpenetrating continua, from the constitutive relations for solidliquid mixtures developed by Roco (1990). The treatment of solid particles as an interpenetrating continuum of liquid is a pioneering approach that has been proved successful in slurry flow in pipes. Davidson's use of these concepts in hydrocyclone modeling awaits further experimental confirmation. Nevertheless, it brings about the importance of Bagnold and turbulent diffusive forces exerted by the dense suspension of particles near the walls.
The modeling of turbulence is much needed because it disperses the classified particles. In the larger context of modeling thus far, turbulence has been dealt with a linear algebraic stress model, and as far as particle phase is concerned a diffusiontype model is used. Dyakowski and Williams (1992) have carried out a detailed analysis of turbulence in small diameter hydrocyclones. In the fluid mechanics literature it is very common to describe the turbulent quantities by the k-£ model, but Dyakowski and Williams (1992) pointed out that the equality of normal Reynolds stresses presumed in the k-£ model is questionable in small diameter hydrocyclones due to the extreme curvature of streamlines as well as the presence of high swirl. The authors proposed a revised k-£ model, incorporating both partial differential and algebraic equations for Reynolds stresses.
Indeed such cumbersome numerical solution of mean flow and turbulence is largely unnecessary for modeling size-classification in hydrocyclones. Barrientos and Concha (1992) used the analytical solutions of Bloor and Ingham (1987) with a set of data from an operating hydrocyclone to determine a few unknown constants. While the model has been confirmed in 150-mm diameter hydrocyclones, this work too must await more experimental con-96 firmation in larger diameter hydrocyclones. Hwang et a!. (1993) too used Bloor and Ingham's inviscid flow analytical solutions to verify velocities measured in a 60-mm hydrocyclone. There is some merit in these analytical solutions, but their applicability seems to be device specific, in other words limited to the data set from which the constants were estimated.
There is second major issue in modeling vortex flow in a hydrocyclone. It is the size and the stability of the air core. The air core occupies space along the central axis, and hence it inherently regulates flow out of the spigot and vortex finder. In the extreme, at very high inlet flow rates the air core can completely occupy the spigot opening, causing all the flow to reverse toward the vortex finder, resulting in substantial degradation in size-classification efficiency. Until recently, this issue has been set aside because solving the Navier-Stokes equations itself posed a challenge. Hsieh and Rajamani (1991) and Davidson (1994) specify the diameter of the air core in their numerical calculations. Pericleous and Rhodes (1986) intentionally keep numerous air bubbles along with the particles in the numerical solution. The air bubbles naturally migrate to the central axis forming the air core, and hence its dimension falls out of the numerical calculation. However, the stability of such a numerical procedure is questionable, according to Davidson (1994). Barrientos et a!. (1993) computed the air core diameter by balancing the viscous forces and surface tension force at the air/water interface. This is a plausible explanation but needs more experimental validation. An alternative is to solve the complete Navier-Stokes equations for air core diameter as done by Steffens et a!. (1993). These authors computed the air core sizes by simplifying the hydrocyclone geometry to a flat cylinder with a radial inlet and one central outlet. Until data on hydrocyclones operating with high concentration of solids appear, Steffens solution can be regarded as valid for water flows only. Even though the stability and the nature of the air core surface is of minimal importance in modeling, yet further understanding depends on detailed studies. In this direction Williams et a!. (1994) have used Electrical Impedance Tomography (EIT) to map the oscillating air core as well as its dimensions. In summary, the physics of air core is at its infancy but is expected to develop very rapidly. A rigorous method of computing the air core diameter while computing velocities is shown in this manuscript.
This paper describes the development of the fluid-flow model and provides model verification with experimental data collected in a 250-mm hydrocyclone. The experimental data includes both velocity measurements with a Laser Doppler Velocimetry and particle size classification measurements in a sump-pump recirculation system.

Transport Equations
The flow field is essentially symmetric except in the region of the inlet. Symmetric flow field is established about halfway in the cylindrical part. Hence axisymmetry (Pericleous and Rhodes, 1986;Davidson, 1994;Dyakowski and Williams, 1992) is a useful assumption that reduces the fluid-flow problem to two dimensions. In the computational model the fluid is assumed to enter through a full 360° ring inlet instead of the single tangential inlet tube. The governing equations are written in vorticity stream-function form to avoid the pressure balance during computations. For two-dimensional incompressible turbulent flow with constant properties, the dimensionless transport equations are as follows: Re ar2 +--;: ar -r2+ az2 (1) Stream function : Here Re is the Reynolds number defined as RcU 0 IP.

Turbulent Transport
The k-8 model is often found in recent literature, as it models the generation and dissipation of turbulent kinetic energy by means of partial differential equations akin to that of Navier-Stokes equations. However, this model in its standard form is not applicable in hydrocyclones, since the high swirling flows impose anisotropic turbulence. Therefore, the revised k-8 model of Dyakowski and Williams (1992) with an algebraic Reynolds stress model as proposed by Boysan et al. (1982) seems like the most appropriate in this direction. However, the use of the k-8 model entails solving three more partial differential equations which require at least five empirical constants. To avoid this extra numerical effort, one tries the simpler Prandtl mixing length model first. It relates the turbulent transport terms to the gradient of the mean flow quantities. For two-dimensional boundary-layer flows, particularly those developing over rigid boundaries, the mixinglength offers both accuracy and numerical simplicity. The mixing length model used in the present model uses the following expression for turbulent viscosity: uses a single unknown parameter, the mixing length A, which is the characteristic length scale of turbulence and whose dependence over the flow field must be prescribed only with empirical information. In this work we consider that the mixing length varies both radially and axially and that separate expressions should be prescribed for swirling and non)swirling components as follows (Hsieh and Rajamani, 1991): quantities in a 250-mm hydrocyclone, and its adequacy and validity is beyond the scope of this study.

Numerical Solution
The physical dimensions of the actual hydrocyclone are explicitly specified in the numerical solution, and the entire domain is divided into a rectangular grid system. A 180 x 98 grid system covering a half view of the hydrocyclone, as seen in Figure 2, was sufficient for the solution. Such a grid system results in a grid spacing of 1.25 mm in the radial direction and 4 mm in the axial direction.
The fluid dynamic equations are a set of coupled parabolic and elliptic differential equations. Equations 1 and 3 are parabolic and pose an initial value problem, while equation 2 is elliptic and poses a boundary value problem. The parabolic equations are solved with the Hopscotch method, and the elliptic equation is solved by the Successive Over-Relaxation method (Hsieh, 1988).

Experimental Work
There are two kinds of experiments in this modeling study: Particle size classification experiments and liquid phase velocity measurements. In the size classification experiments, a KREBS D10B (250-mm diameter) hydrocyclone was installed in a sump-pump recirculation system. Table I and II shows the dimensions of the hydrocyclone and the operating conditions along with two sets of vortex finders and spigots used in this study. Experiments were done at various flow rates and solids concentrations Boththe overflow and underflow streams were sampled, and the size distributions as well as percent solids and respective flow rates were measured The size-classification efficiencies were determined from the test results by calculating the fraction of the feed reporting to the underflow for each size interval Duplicate tests exhibit a scatter of ± 5o/o error in the experimental results.
The velocity measurements were done with a TSI Laser Doppler V elocimetry (LDV) system. A schematic drawing of the LDV)pump)sump hydrocyclone rig is shown in Figure 1. A 100-mW argon-ion laser (ILT model5500A) generated a coherent light source. The laser beam passed through a beam splitter (TSI model 9115), which splits the laser beam into two parallel, equal intensity beams. Beam spacing was fixed by the prism design. One of the split beams 98 passes through a Bragg cell which shifts the frequency by 40 MHz. The flow direction is detected by downmixing the required frequency to 40 MHz with the help of the frequency shifter (TSI model 9186A). The frequency shifter also helps in the measurement of low-velocity flow, removes the pedestal frequency and reduces the fringe bias. The intersecting beam measurement volume is a 1.997mm x 0.1997mm ellipsoid. Back-scattering alignment was chosen because of the presence of the air core at the center of the hydrocyclone. The output signals from the photomultiplier are sent to the signal processor (TSI IF A 550) which filters the signal, digitizes and transmits the Doppler frequency to the computer. The computer software displays the velocity spectrum measured at that point. The entire LDV system is mounted on a three-dimensional traverse system which allowed the probe volume to be positioned with ± 0.5 mm resolution.

Experimental Results
A few qualitative observations can be made with LDV measurements. The tangential velocity profiles in four different planes along the axis are shown in Figure 2. These profiles always exhibit a forced vortex near the air core and a free vortex thereafter. The slopes of the forced vortex region are nearly the same throughout the axis, suggesting the inviscid nature of the flow in the central region. In all the four planes shown, the tangential velocity reaches a maximum of about 8 m/s, which corresponds to a centrifugal force of 210 g. Hence the centrifugal force overwhelms the gravitational force, even more in smaller diameter hydrocyclones. The gradient of the forced vortex in the conical section is nearly the same as that in the cylindrical section. Since the cone narrows downwards the increase in tangential velocity is compensated by frictional dissipation of momentum, hence the constant gradient. The measured axial velocities for Series I and II are shown in Figure 3. Not shown here is the asymmetry observed near the inlet as well as the apex region, but everywhere else symmetry prevails. In fact multiple flow reversals were observed in the horizontal planes over the vortex finder in a 75-mm hydrocyclone (Hsieh and Rajamani, 1991). The vortex flow spirals downward along the cone walls, and since the cone narrows, part of the flow in the central column spirals upwards. As two flows interact, multiple flow reversals occur. The axial velocities decrease rapidly toward the air core, presumably to zero at the air/water interface. But the LDV is unable to measure the velocities in the vicinity of the air core surface since the air core vibrates continuously. 100 not be investigated any further. One would expect the fluctuating components to diminish in magnitude as the streams gets weighed up by solid particles. Finally, the radial velocity could not be measured at all reliably with the LDV, since its magnitude is of the order of 10 em/sec. Hence this component was computed by applying the continuity condition.

Predicted Velocity Profiles
The purpose of predicting velocity profiles is to verify the assumed form of the Prandtl mixing-length model and check the relevancy of the boundary conditions. The ability to predict the flow field leads to accurately predicting solid particle trajectories. The LDV data and predicted tangential velocities are shown in Figure 5. It is seen there is good quantitative agreement both in the forced-and free-vortex regions. The gradients and the peak velocities are also accurately predicted at each horizontal level.  Figure 6 shows the LDV data and predictions of axial velocities. It is evident that the model is able to track the upward flow in the center and downward flow along the wall. There is a slight discrepancy in the conical part: There is a small lateral shift in the loci of zero axial velocity, and the predicted maximum axial velocity is shifted toward the air core more than the LDV data would indicate. The flow patterns near the air core, where the axial velocities rapidly tend to zero, are also well predicted.
The mixing length model. This model was arrived at after examining many others in the work done on a 75-mm hydrocyclone (Hsieh, 1988). Indeed the same formulation is used here in the 250-mm hydrocyclone, which confirms that the basic structure of turbulence generation and dissipation remains same in both smaller and larger hydrocyclones. It is reasonable to presume this is true also in much larger hydrocyclones of 625-mm diameter. Predicting velocities is the key to the hydrocyclone modeling, and it hinges squarely on the turbulent closure as is true in many other fluid dynamic applications.

Air Core Diameter
The diameter of the air core measured in Series I and II were about 70 mm just below the vortex finder and 60 mm at the apex. The entire envelope of the air core was specified in the numerical calculations to solve for the velocities. In the design and evaluation of industrial hydrocyclones operating with the slurries, the size of the air core is unknown and must be computed. Instead of taking highly simplified methods proposed by Barrientos et al. (1993) or Steffens et al. (1993), we use the Navier-Stokes equation combined with the turbulence closure model itself to determine air core diameter. The Navier-Stokes equation itself is a balance for viscous and body forces; hence the solution must satisfy zero pressure condition at the air core when the inlet pressure is specified. The only quantity that is missing in such an approach is the pressure jump across the air/liquid interface. In addition, the hydrocyclone entrains the central column of air through the apex and discharging through the vortex finder. Therefore there is a small but nonzero pressure gradient in the air column which is neglected here.
Once the velocity set has been solved for at all nodal points, the gradient of the pressure in the two directions can be easily determined by substituting the computed velocities in the Navier-Stokes equations. To make matters even easier the SIMPLER algorithm (Patankar, 1980) is used to calculate the pressure. The consistency of the pressure calculations can be verified by measuring pressure along the wall in addition to the inlet pressure. Figure 7 shows the measured and predicted pressures. The agreement is consistent with the excellent velocity predictions reported in Figures 5 and 6. Curiously the pressure drops to zero at the air c?re boundary. Hence this pressure prediction can be used in the reverse manner to determine the envelope of the air core. First the calculations are begun with as small an air core as the grid spacing would permit. Then fixing the measured inlet pressure, the entire pressure profile is mapped over the whole domain. If the pressure at the presumed air core surface is greater than zero, the core must be widened in the next iteration and vice versa. The results of such calculations would be reported in the future.

Particle Size Classification Experiments
The ultimate aim of the modeling work is to predict the size classification performance of the hydrocyclone. The efficiency of the hydrocyclone is measured by the "size classification curve," which is a plot of particle size versus the fraction of the amount of the feed (in each particle size interval) reporting to the underflow. Now, having solved the fluid dynamic equations for the water phase only, what is needed is a method of how particles migrate in the centrifugal field and exit via the apex or the vortex finder. Once particles are introduced in the fluid phase, the properties of the fluid phase change and there is some question about the validity of the computed velocity field. There are two ways to remedy this problem: The first one is to assume that the fluid is "dilute," especially when the concentration of solids in the feed is less than 15 weight percent, in which case the fluid-flow problem and particle classification problem can each be solved independently of the other. To a large extent this approach works because, even if the feed concentration is high, particles quickly migrate to the wall, and so the central body of fluid surrounding the vortex finder and below it becomes very dilute. Therefore particle separation occurs in this dilute regime. In fact industrial hydrocyclones processing 40-55% by weight of solids discharge a thick slurry of 65-70% by wt. in the underflow and a mere 20% by wt. in the overflow. Since the overflow is made up of fluid migrating from all along the air core, it can be concluded that the central column is "dilute" within the hydrocyclone. Anyhow, the second and more rigorous method is to start with the water phase only, solve for particle concentrations, and iterate between the Navier-Stokes solution and the particle concentration calculations. A logical way of calculating the concentration is to determine the residence time of particles in a control volume and proportion the feed concentration to the nodal point with respect to its residence time among all the particles in the inlet (Rajamani and Milin, 1992).
Under the assumption of dilute inlet slurry the path of the particles from the inlet to one of the outlets is easily determined by computing the particle slip velocities. The balance between the centrifugal force and the radial drag force yields the radial slip velocity: Similarly the balance between gravity and the axial drag force yields the axial slip velocity: 112 3 e1 CD Since the axisymmetry assumption was used, the 102 particle is considered to move along with the fluid in the azimuthal direction. With equations 8 and 9 the slip velocities at every grid point are calculated, and hence the trajectory of the particle phase is precisely calculated. Calculating the fraction of particles (in the feed) reporting to the underflow and repeating the above procedure for all sizes in the feed stream yields the classification curve. The experimental and predicted classification curves for a feed solid concentration of 10o/o by wt. are shown in Figures 8 and 9. There is generally a good agreement, especially in the vicinity of the size at which 50% reports to underflow. The overall flow split to the underflow is also predicted well. However, there is a slight discrepancy in the fine particle range, which is due to the computed drag force in the fine particle size range being higher than what might have been the actual drag force. It is believed that the use of Reynolds number versus drag coefficient commonly found in the literature may be in error because these correlations are determined under   unit gravitational field, while the centrifugal field inside the hydrocyclone is of the order of a hundred times gravity. Also, turbulent diffusion which is not accounted for in this modeling work may contribute to fine particle migration. In the coarse size range too there is a slight discrepancy: This is due to the axisymmetry assumption in which the model does not account for the coarse particles passing via the boundary layer on the outer wall of the vortex finder, especially on the surface 180° from the inlet flow direction. Nevertheless, the predictions are more than satisfactory for industrial purposes.

Concluding Remarks
The fluid-flow model fully takes into account the design of the hydrocyclone as well as the slurry properties. The design of the hydrocyclone, such as the cone angle, vortex-finder diameter etc., are specified, and the initial/boundary conditions of the numerical procedure are assigned, which in turn influences the computed vortex flow. The concentration of solid particles influences the the slurry viscosity which in turn affects computed flows. This influence is more so in the iterative particle concentration model. If a correlation of particle size and particle concentration versus slurry viscosity is available, then the numerical model reported here accounts for the effect of particle size distribution also. Turbulent diffusion of particles in slurries is still an emerging field, and as correlations become available, they could be incorporated into the model. The only empiricism involved is in the formulation of the turbulence model. A fully mechanistic explanation of turbulence, such as direct simulation approach, is not anticipated in the immediate future. Being so, the model which holds well in the 75-mm and 250-mm hydrocyclones may in fact correctly describe even larger hydrocyclones.

Acknowledgments
This research has been supported by the Department of the Interior's Mineral Institute program administered by the U.S. Bureau of Mines throughthe Generic Mineral Technology Center for Comminution under grant numbers G1115149 and G1125249.
Nomenclature CD drag coefficient dp particle diameter, em KONA No. 12 (1994) R radial distance from the axis of symmetry in cylindrical coordinates, em r radius of the hydrocyclone, em Reynolds number of the hydrocyclone defined as RcUolv radius of the hydrocyclone at the horizontal level in consideration, em dimensionless radial distance from the axis of symmetry in cylindrical coordinates dimensionless time mean inlet velocity, cm/s dimensionless radial velocity of the fluid particle radial slip velocity, cm/s v w tangential velocity of the fluid, cm/s dimensionless tangential velocity of the fluid axial velocity of the fluid, cm/s dimensionless axial velocity of the fluid particle axial slip velocity, cm/s dimensionless axial distance from the roof of the hydrocyclone in cylindrical coordinates Prandtl mixing length, em density of the fluid, g/cm3 density of the slurry, g/cm3 density of the particle, g/cm3 viscosity of water, g/cm-s viscosity of the slurry, g/cm-s