The Numerical Modelling of the Flow in Hydrocyclones †

Over the last decade, the progress in numerical modelling of hydrocyclones has advanced at a rapid pace. A summary of the recent advancements is undertaken in the paper to summarize the major directions and developments in hydrocyclone modelling. The important aspects that contribute to the accurate simulation of general internal swirling flows and other characteristic features occurring in hydrocyclones modelling are covered. The relevant issues are: the formulation of the governing equations (classical versus stress divergence form), the imposition of boundary conditions especially regarding the open part of the outlet boundary conditions and the problem of proper representation of the computational domain. The complex flow pattern that occurs in hydrocyclones invalidates several assumptions in the standard turbulence models. Challenges still remain in predicting accurately the effects of turbulence anisotropy. A performance review of selected equations closures is presented with insight on their potential and limitations. The flow predictions obtained using the Reynolds Stress Model (RSM) and Large Eddy Simulation (LES) turbulence models are examined. The predictions of the free surface flow and air core dimensions have yet to be fully resolved. The same concerns the problem of simulating large number of particles. The recent important developments related to the solid/gas and sold/liquid interaction problems are presented together with intriguing challenges one need to overcome in order to numerically handle multi-component flows. The Lagrangian and multicontinuum Eulerian techniques are considered. The attempt to systemize various approaches for multi-component flows according to their accurateness and relevance to the various flow regimes in hydrocyclones is undertaken. The issue of experimental validation is crucial to provide confidence in hydrocyclone flow-field models and the current state of experimental data is examined. Finally, suggestions for the future direction of hydrocyclone modelling are highlighted.


Introduction
Hydrocyclones are primarily used for the separation of particulate slurries as well as liquid-liquid systems.A wide variety of industries, including petrochemical, mineral processing, and pharmaceutical, have recognized the advantages of hydrocyclones in separation.The effectiveness of the hydrocyclone in diverse applications signifies the economic and operational value as a separational unit.The foremost advantage is the versatility of the hydrocyclone in a variety of applications ─ concentrating, dewatering, and separation.The economic cost of a hydrocyclone is significantly less compared with other common unit of operations.Due to absence of moving parts, a hydrocyclone has an extended operational life and low maintenance costs versus other separation technologies.These characteristic traits allow continuous operation without severely affecting the performance of the process.A significant drawback to hydrocyclones is the amount of waste that is produced from the separation.Another disadvantage is the large dependence of the separational ef ficiency on the mass flow rate, which is a difficult parameter to control.Although the affect can be reduced by certain measures, the efficiency of separation will vary with time in a process.Despite the disadvantages, hydrocyclones are still a widely used unit of separation.
A hydrocyclne operates by tangentially inserting a feed that causes a high centrifugal force to occur.The slurry is injected into the hydrocyclone via one or several tangential inlets.The confinement of the feed in the conical/cylindrical section of a hydrocyclone produces a high swirl flow, which is the primar y mechanism of separation for particles or droplets from liquid─solid or liquid-liquid mixtures.The combined axial and tangential components of velocity field characterize a flow streamline that rotates around the central axis.Bank and Gauvin (1977) revealed that the radial velocity has a minor contribution in the separation of the particles.The wall curvature directs the flow into a spiral path that travels downwards towards the apex, as shown in Fig. 1.While the wall friction reduces the tangential velocity, the radial contraction in the conical section maintains the magnitude of the velocity.The centrifugal force acting on the par ticles is directly propor tional to the mass of the particle.Larger particles tend to be driven towards the wall, while smaller particles tend to accumulate in the center.The heavy particles in the suspension travel towards the hydrocyclone wall and then down to the apex.In the conical contraction section, an axial reversal flow is formed and moves upwards carrying the smaller particles through the vortex finder.According to Rietema (1961), the total pressure drop across the hydrocyclone consists of the inlet velocity head, friction losses and the centrifugal head.The pressure drop determines the maximum volumetric throughput that can be processed.
Although hydrocyclones have been used in industry since the 1950s, the details of the flow behavior are still not well understood.Modelling and experimental investigation of the complexity and range of phenomena of fluid flow in a slurr y separator has been attempted for many years and is presented in numerous scientific contributions.The nature of fluid flow and in particular the swirling flow characteristic has been studied experimentally by Fisher and Flack (2002), theoretically using computational methods by Nowakowski et al. (2004), or simplified analytical approaches with the most extensive presented by Bloor and Ingham (1987).All these strands of the research confirm that a vortex constituting a swirl flow is composed of two separate vortices, a free and forced.A forced vortex resembles a solid body rotation.Near the wall, the swirling components become less dominant and the device geometry exerts more influence.The relationship between the tangential velocity and the radius is inversely proportional in the free vortex and proportional in the forced vortex.Although this phenomenon seems to be obvious, it is intuitively understood and is obser ved in prevalent studies, there are some uncertainties related to how the swirl flow could be controlled.Detailed investigation of the velocity patterns in various flow regimes may provide a fruitful avenue for understanding a number of open issues regarding internal swirling flow and the implications on separational efficiency in a hydrocyclone.Detailed knowledge of the flow structure is required if one is to consider such issues as energy saving, cost or product quality.Several benefits could arise from this knowledge, for instance: areas of high erosion may be identified and potentially minimized or accounted for in design; design modifications for improved separation or reverse engineering of a hydrocyclone's geometry could be obtained.
Along the hydrocyclone axis, an area of low pressure is created by high angular momentum.This may cause the formation of a rotating free liquid surface at the centre.If the hydrocyclone is open to the atmosphere, air is inhaled through the outlet orifices and forms a central air core.If the air core exists, the pressure at the air-liquid interface is equivalent to atmospheric pressure if we neglect both the surface tension and viscous forces.Air core formation in a hydrocyclone is not currently known.The prevalent thought, however, is that the air enters from both outlets and eventually reaches a critical point, which combines the two air vortices, Cullivan et al. (2004).
A majority of the current hydrocyclone modeling work has largely ignored the complex features of aircore development, by making several simplified assumptions.
The operational states of hydrocyclones have been defined by the type of discharge at the hydrocyclone apex.The three different operational states can be distinguished, which are spray, roping, and transitional discharge.Spray discharge is the occurrence of unobstructed flow through the apex and exits the apex in the form of an "umbrella" .The roping discharge states occurs when the apex has been impeded.The roping state reduces the separation efficiency by affecting the split ratio of the device.Neesse et al. (2004) show that the optimal hydrocyclone operation conditions occur in the transitional discharge state between roping and spray.Datta and Som (2000) present the effects of changing different operating parameters on the discharge angle.
A popular method of hydrocyclone design consists of var ying operational parameters, such as pressure drop, cut size, and volumetric throughput to determine the effects of the changes.For example Kelsall (1952) developed a design process that uses an empirical model to produce a classification curve at different flow rates.Using the classification curve, the hydrocyclone's geometrical parameters could be defined for different process streams.A classification curve is normally used to determine the d50 parameter, which is the particle size that has an equal chance of exiting either the overflow or underflow.
The emergence of high performance computing has allowed the numerical simulation of partial differential equations become a viable alternative to physical experiments and empirical methods.The problem associated with mathematical modelling of the flow patterns involves the solution of strongly coupled, non-linear partial differential equations of the conservation of mass and momentum and lies well beyond any analytical approach.Additionally in the context of modelling a flow in hydrocyclones the challenges are related to the complex flow behaviour arising from the three-dimensional flow entry; multi-phase interactions; and the mechanisms governing the formation of an air-core (when the hydrocyclone is open to the atmosphere).Consequently, computational studies have been in general, limited to low particleconcentration flows and to simplified geometries of the hydrocyclone entry region.
Even with cur rent limitations, computational methods of design have significant advantages over the empirical data, such as for example freedom to change the geometry quickly for verification of possible changes in separational efficiency.
The paper reviews the existing theoretical models and addresses how the computational algorithms could be extended to include complex interactions between continuous and dispersed phases.Specifically, the paper considers methodologies for velocity field prediction and modeling of the particle distribution.A strategy for future developments in these areas is outlined.The paper also briefly discusses essential measurement techniques for the validation of hydrocyclone flow field.

The Formulation of the Governing Equations and the Velocity Field Prediction
Nowakowski et al. (2004) provided a compilation of many of the major contributions to the study of hydrocyclone fluid dynamics between the years: 1982-2003.The application of full three-dimensional modelling to the hydrocyclone has been attempted only during the last decade with early attempts being made by Concha et al. (1998).Authors such as He et al. (1999) demonstrated that full three-dimensional modelling is essential in order to accurately model the hydrocyclone flow field, as a result of its inherent axial asymmetry.The present authors have endeavoured to compute hydrocyclone performance using their own finite element code and the commercial computational fluid dynamics package Fluent.Some aspects of the numerical algorithm development have previously been described in detail in Nowakowski and Dyakowski (2003).The present review is focused on fully three-dimensional approach and introduces a summary of numerical studies from the year 2000 onwards.
The objective is to compute the flow problem in a hydrocyclone modelled by the Navier-Stokes equations in three-dimensional domain Ω with boundary ∂Ω = ΓDw+ΓDe+ ΓN.The ΓDw stands for solid walls and ΓDe for inlet boundary, while ΓN indicates the outflow parts of boundary.The underlying equations written in primitive variables and stress-divergence form are: (1) (2) where constitutive relations for incompressible fluids are given by (3) The standard notation for velocity u and pressure p is used and ρ is constant density; µ dynamic viscosity; τ is the total stress and I is the identity tensor.The initial velocity field is assumed to satisfy incompressibility constraint (2).The advantage of stress-divergence form is that it can readily treat problems with fluid constitutive equations in a more general way than the classical form where the stress tensor is decomposed.Therefore, the extensions of the application of numerical techniques to Non-Newtonian fluid can be considered.Secondly the stress-divergence form leads to the physical Neumann boundary conditions that represent force f (4) and Dirichlet boundary condition for velocity field u at the inlet and solid walls (5) where n is unit outer vector normal to ∂Ω, f is ap- plied vector valued force on the boundary, while û is specified velocity profile at the inlet.
The velocity ub must satisfy the global condition.
(6) that follows from integrating the continuity equation (2) over the domain Ω and using the divergence theorem.A serious difficulty encountered is the determination of the pressure field and the fulfilment of the incompressibility constraint.While the continuity equation represents a constraint for the velocity field, the pressure variable, which appears in the momentum equation through the term ∇p, provides a degree of freedom necessary to satisfy such a constraint.It is important to note that no boundary condition is prescribed for the pressure.On inflow and outflow boundaries, the pressure can only be supplemented as a part of the Neumann type of boundary condition.
The specific form of the governing equations depends on the choice of dependent variables.The primitive variable formulation (the equations written in terms of u and p) has been preferred by the authors and other hydrocyclone investigators.The primitive variable formulation expressed in the "stress-divergence" form (Gresho, 1991) is commonly used in the finite element numerical framework.
A numerical solution of the equations of motion in a hydrocyclone has also been obtained using a streamfunction-vorticity formulation by Ovalle et al. (2005).
The numerical treatment of even the simplest turbulent flows generates enormous difficulties.Primarily, this is due to the complicated hydrocyclone geometry generating swirling flow and the fact that the non-linearity in the Navier-Stokes equations, gives rise to a very broad range of spatial and temporal turbulent scales.
The Reynolds Averaged Navier-Stokes equations (RANS) modelling approach enables the approximation of the turbulence effects on the mean flow.The Reynolds equations have the same form as the instantaneous equation 1, with the exception of the secondmoment or correlation tensor, R = < u  ⊗ u  >. (7) The equations were solved using a finite element code introduced by Nowakowski and Dyakowski (2003) for laminar and turbulent regime.In order to clarify the three-dimensional structure of the flow and show its swirling features the particle trajectories were plotted in Fig. 2. The simulation results demonstrated that the short-circuit flow moves from the feed inlet and along the outer wall of the vortex finder tube until joining the rest of fluid in the overflow.The results also revealed the development of the double helical flows, the inner-upward and outerdownward helical flows (Fig. 2), and the change of the processing vortex core as it moves about the hydrocyclone axis as well as the characteristic pattern of the developing flow at the under flow and overflow outlets.The numerical simulations were able to reproduce the swirling character of the flow in both vertical and horizontal directions and confirmed the existence of asymmetr y as obser ved in labora- tor y experiments.At the outlets, the mathematical description of the problem (stress-divergence form) enabled imposition of natural boundary conditions in terms of forces, provided that the necessary data can be determined experimentally.In the simulations performed, zero stress boundary conditions were prescribed at the hydrocyclone outlets.This boundary choice is neutral in character, does not restrict the profile of velocity field, and is dictated by the need to truncate computational domain.Such an imposition of boundary conditions avoids any a priori assumptions concerning the velocity field at the outlets and most importantly, does not define the mass split-ratio.
The postulated descriptions of outflow boundar y conditions enable simulation of a characteristic spray profile velocity field at the spigot.Doby et al. (2005) investigated the interaction between the swirling flow and velocity profile at the outlet.The studies were carried out for fluids with different properties and it was shown that with the proper imposition of boundary conditions the flow develops naturally.
Owing to the complex hydrocyclone geometry an unstructured mesh generation technique has been used by a majority of researches studying the numerical flow in hydrocyclones.The key advantage of the unstructured grid methodology is the ability to handle not only the hydrocyclones of the same class but also with different geometry configuration.

Turbulence Modelling
All models, discussed below, are based on the time-averaged equations and the turbulence effects are incorporated through the Reynolds stress tensor.
To obtain a quantitative solution several assumptions concerning the physics of the turbulent flow have to be made.
In RANS, the right choice of turbulence closure that addresses the unknown correlation terms of the velocity fluctuations (Reynolds stresses) is vital for correct separation prediction.The highly rotational fluid flow in the hydrocyclone makes other wise successful RANS turbulence models fail to produce practically reasonable results.It is well known, for example, that for swirl flows the k-ε model (Harlow and Nakayama, 1968) predicts a too rapid decay of the swirl because it assumes isotropic viscosity, when in fact Reynolds stresses in hydrocyclone are highly anisotropic (Petty and Parks, 2001).The k-ε turbulence model is a two-equation model in which transport equations are solved for the turbulent kinetic energy, k, and its dissipation rate, ε.There are many forms of the k-ε model, which have been used for several decades, and is currently the most widely utilized model for industrial applications.The k-ε model provides accurate solutions for fully turbulent flows.Although Dai et al., (1999) modified the k-ε model to calculate the anisotropic character of turbulence flow within a hydrocyclone by developing the model constants Cε1 , Cε2 and Cµ..These constants in the standard k-ε model were determined from experimental measuremrnts in a wind tunnel with almost uniform velocity fields.It is apparent that these constants, determined in this way, are not applicable to hydrocyclone flow field.
Many research groups used the computational fluid dynamics package Fluent, with its finite volume platform, to investigate the capability of different turbulence models for hydrocyclone simulation.Cullivan et al. (2003) highlighted that modelling of hydrocyclone flow-field requires comprehensive modelling of turbulence to account for significant anisotropy and turbulence generation mechanisms.The Fluent full differential stress model (DSM) with highorder treatment of the pressure-strain term has been shown to be a good compromise between accuracy and computational cost.The DSM does not represent an entirely accurate description of the turbulence, in particular sub-grid scale modeling and equilibrium turbulence is assumed.The equilibrium assumption presumes that the rate of transfer of turbulence energy down through the energy containing length- scales is constant, which is not expected to hold for the rapidly developing and short residence time hydrocyclone flow-field.Even so, the DSM turbulence model offers practicality in terms of computational resource.
Earlier turbulence models were purely dissipative and did not account for Reynolds-stress relaxation effects.In present practice of modelling turbulence in hydrocyclones Reynolds stresses are calculated from a set of partial differential equations.Each Reynolds stress component can be viewed as a 'property' of fluid.Its progress is determined by its rate of production, destruction, diffusive spread and redistribution between components.The redistribution mechanism is governed by the "pressure-strain" tensor which is zero for an isotropic turbulence.The RSM model consists of a set of six partial differential equations for each component and one partial differential equation for the dissipation rate (Appendix A).Reynolds Stress Model (RSM) has been applied as turbulence closure in the study presented by Kraipech-Evans et al. (2008).The simulations using RSM model and supporting validation studies in the context of hydrocyclones were also presented in recent papers of Bhaskar et al. (2007) and Wang and Yu (2006).
The RSM model was found to predict well anisotropic turbulence.The RSM provides information about all the stress components and contains exact terms for swirling effects in the stress transport equations.
The Algebraic Stress Model (ASM) model was developed to account for anisotropy without solving for Reynolds stress transport equations.The ASM is based on the assumption that the convective and diffusive transport terms are removed or modelled, thus reducing the equations of Reynolds stress transport to algebraic equations.Chu and Chen (1999) introduced the algebraic stress turbulence model in simulation of a hydrocyclone flow, which is computationally less expensive than the Reynolds stress model, but still possesses a few of the RSM characteristics.
Alternatives such as large-eddy simulation (LES) and proper-orthogonal decomposition are still subject to active research and these require considerable additional computational cost.Slack et al. (2000) have demonstrated the DSM to perform remarkably well in comparison to LES without such prohibitively high computational costs.
The direct numerical simulation (DNS) and large eddy simulation (LES) are becoming more important and viable options for the solution of transitional and turbulent flows governed by the system of equations ( 1) and ( 7).The number of grid points required for a DNS is determined by the ratio of the respective length scales of the energy carrying dissipative eddies.This ratio is a very rapidly growing function of the Reynolds number and consequently the application of direct simulation is at present limited to flows at moderate Reynolds numbers and rather simple geometries.In the LES simulation approach, the severe Reynolds number restrictions are bypassed by directly simulating the large scales only and modelling the effects of missng small scales by using a subgrid model.Unlike the Reynolds averaged Navier-Stokes equations approach (RANS), large-eddy simulations retain the full three-dimensionality and time dependence of the fluctuating turbulent field.As a result of recent progress the LES technique has been used for cyclone (Slack et al., 2000) or hydrocyclone (Delgadillo and Rajamani, 2005) simulation.It seems that the computationally more expensive LES provides the best solution for capture of time dependent vortex oscillations and non-equilibrium turbulence, which will potentially impact the separation efficiency.Challenging tasks associated with LES are the implementation of a sub-grid scale model accounting for the effects of particles and the correlation of different time scales related to collisions and aggregation of particles.
Although many of simpler turbulence models are often incapable of handling the strong vortical features that are present in hydrocyclones, the emphasis on simplicity, computational speed and robustness in an industrial context does not always allow the adoption of more advanced models.Using the averaged Navier-Stokes equations for modelling requires empirical or semiempirical information on the structure of the turbulence and the relationship to the mean flow.The information needed for the modelling is not easily obtained in hydrocyclones.

Air-core Modelling
The air core is an interesting and not well-understood phenomenon that occurs in hydrocyclones and drastically influences the efficiency of separation.Predictions of the air-core diameter have been provided through: application of Bernoulli's equation applied with a minimisation procedure (Davidson, 1988); through the use of an effective air-core interface viscosity (Dyakowski and Williams, 1995); as well as through application of Young-Laplace's relation (Concha et al., 1998).The later two approaches offer the potential to account for asymmetric air-core geometry.Dyakowski and Williams (1993) showed that the air core is not stationary, but has an oscillatory nature.The oscillations could be caused either by fluctuations in the feed or induced inertial waves (Ovalle and Concha, 2005).Delgadillo and Rajamani (2005) investigated the dynamics of the air core using a LES turbulence model.Doby et al. (2008) numerically examined different flows at varying viscosities and revealed that low viscosity feeds tend to develop a larger low-pressure area in the center region of the hydrocyclone than feeds at high viscosity (Fig. 3).As shown in Fig. 3, air from the atmosphere pushes the fluid due to the low pressure formation in the center area and upper section in the vicinity of the vortex finder.The total area of the low pressure reduces with increasing viscosity.Based on the contours of the apex and upper section, the air core formation appears to begin in both the upper and at the apex.
The air core starts from each of the outlets and combines together in the central part of the cyclone.The numerical and experimentally obser ved trends for hydrocyclones operating with an air core and with an inserted metal-rod were comparatively discussed in Kraipech et al. (2008) (Fig. 4).The air-core characteristics obtained in this work are similar to the experimental observation as can be seen in Fig. 4. The air core was found to be unstable and its size, shape and position are unfixed because of the instability of the gas-liquid interface.
Calculation of the position of the air-core can be done using the interface-tracking algorithm, which is based on the volume of fluid method (VOF) of Hirt and Nichols (1981).The equation of motion (1) and (2) are solved for the mixture to obtain the velocity field which is shared between the two phases.The volume fraction averaged density ρ is given by ( 8) The fraction of fluid in each cell is defined as αk, which varies between 0 and 1.An additional transport equation for the volume fraction has to be solved in order to track the air-core.It has the following form: (9) Delgadillo and Rajamani (2005) used the interface tracking algorithm based on the volume of fluid method (VOF) for calculating the position of the air-core.It was found however that the RSM model predicts an air-core, which is irregular and does not agree with experiment.Alternatively, if VOF or other multi-continuum approaches are not used in numerical simulations, where the air core is considered, the Fig. 3 The effect of viscosity on the formation of the air core (Doby et al., 2008).
Fig. 4 The obtained pressure (Pascal) of fluid projected on a vertical plane of the hydrocyclone: a) operating with an air core, b) experimental observation (Kraipech et al., 2008).
interfacial boundary conditions are required.Natural boundary conditions, as postulated by Nowakowski and Dyakowski (2003), relating forces could be advantageous for such simulations.
The other approach that could be proposed to tackle air-core is the capture of the liquid-gas interface by means of the level set method of Osher and Sethian (1988).This is relatively simple and versatile approach for computing and analyzing the motion of an interface in three dimensions.Caiden et al. (2001) further developed this approach for the two-phase flow regime addressing the boundary conditions and coupling at the compressible/incompressible interface.
The accurate modelling of the air core still remains a challenging task and it is expected that various new approaches will be adapted to hydrocylone simulation as a result of the recent progress in the development of numerical methods for tackling free surface flows.The important issue is to address an air-core creation phenomenon as a strongly coupled problem to the slurry concentration and swirling flow pattern, which governs the flow split between the outlets.

Numerical Approach to Multi-component Flows
In modelling the hydrocyclone performance, the influence of the particles on the flow is significant, particularly in the dense slurr y flow, when the exchange of momentum from the particle─fluid, particle-particle and particle wall interactions affect the velocity of the fluid.This may cause inefficiency in separation performance.The frequently reported but anomalous 'fish-hook' effect, which results in an excess of fines reporting to the underflow, has not been categorically explained.The presence of the vortex-ring behind particles seems to be worthy of consideration in terms of hydrocyclone performance.
It may be responsible an additional mechanism for finer particles reporting to the underflow in the wake behind the larger particles.Such a mechanism may explain the shape of the selectivity curve, and the fact that the bypass value is higher than the water recovery to the underflow.This selectivity curve does not have a sigmoidal shape, but exhibits a dip in regions of finer particle size.This dip is known as the fish─ hook effect.Kraipech et al. (2005) studied three impor tant par ticle interactions, which are the liquid-solid (drag) interaction, particle-particle lubrication, and particle-particle collisions interactions.The study of the importance of these particle interactions was conducted by applying the time scale analysis.This method is based on the force balance equation and can determine the predominant interaction of the flow within a hydrocyclone.The result shows that the particle-particle interactions due to the lubrication and collision mechanisms only play an important role in the area near the hydrocyclone wall, and near the air core.The majority of the hydrocyclone's area is dominated by the particle-fluid interaction.Therefore such a estimation can only be used to obtain a rough information of the predominant particle interaction and the concentration map.The idea behind this study is to reduce the complexity of modelling the particle motion in a hydrocyclone, since the particleparticle interaction can be ignored in most areas of a hydrocyclone.
Solid-liquid interaction problems are still difficult to solve with existing computational methods.The most popular numerical techniques for simulation of fluid-solid two-phase flows can be classified into several different categories.

A multi-continuum model
The most natural approach is to use a continuum theory that treats the solid and the fluid as interpenetrating continua, each governed by conser vation laws, either postulated or derived by averaging.The work of Drew and Passman (1999) focuses on the fundamental aspects of such systems of equations.This Eulerian-Eulerian approach results in field equations for the flow proper ties for all phases in the system, but leads to unknown terms representing the interactions between the phases.These terms must be modelled to close the description of the system (constitutive relationships).The nature of the detailed interactions between the solids and the fluid is not the consequence of mixture theories alone and hence the derivation of such constitutive equations is an active research area (see Marchioro et al., 2001).Experimental data, recently published in the literature, showed that during the sedimentation process smaller particles move at almost identical velocities to the larger particles, from which it can be concluded that the smaller particles appear to be influenced by larger particles wake.Such observations are crucial for construction of the set of constitutive equations required for the multi-continuum model.Thus, many ad hoc approaches, based on phenomenological con-siderations and physical intuition, have been applied but with only limited success.The Eulerian-Eulerian approach has been used in practical engineering multiphase flow simulations.The two-dimensional multi-continuum approach of Pericleous and Rhodes (1986;1987) was well advanced for its time and represents to some degree the direction that modern day multiphase hydrocyclone modelling could take.However, modelling a distribution of types and sizes of particles complicates the continuum formulation because separate continuity and momentum equations have to be solved for each particle size and type.Nowakowski et al (2000) introduces the concepts and principles of multi-continuum framework for studying of the hydrodynamics of the flow in a hydrocyclone.
In this model the carrying liquid is described as one continuum, and each particle fraction, with its characteristic size is described as a separate continuum.Particle-particle and particle-liquid interactions derived from a lubrication theory and a collision theory are described.

Particle tracking technique
A second approach to multiphase flow simulations is the Eulerian-Lagrangian model for fluid-solid simulation (see e.g.Patankar et al., 2001).The approach provides a direct description of the particulate flow by tracking the motion of individual particles.Newton's second law, with empirical or semi-empirical forms for the hydrodynamic forces, governs particle motion.The particles do not perturb the flow field (unless computationally expensive coupling is undertaken) and the fluid satisfies the continuum equations that are solved on a fixed field in the usual Eulerian way.In hydrocyclones (see e.g.Ma et al., 2000) this one-way coupling between the phases has been employed mainly for dilute systems with a maximum solids volume fraction of 5%.For such cases the volume occupied by the particles in a computational cell may be neglected.However, the accumulation of solid particles in regions of high fluid strain rate and low vorticity can result in high values of the local particle concentration, indicating the presence of a significant (local) coupling of the two phases.When solids concentration exceeds 5% by volume, the presence of particles changes the viscosity stresses and results in the generation of extra inertial stresses.The latter, known as the Bagnold dispersive stresses, result from particle-particle collisions, which are important when particle concentration exceeds 10% by volume.
In case of shearing particles of mixed size, the larger particles drift towards the zone of least shear strain, e.g.towards the hydrocyclone axis, and the smaller particles towards that of greater shear strain, e.g. to the wall.Consideration of the Bagnold stresses might prove useful for explaining the 'fish-hook' effect (see e.g.Kraipech et al. 2002).
The equation of motion of a spherical particle in a fluid neglecting the interactions with other particles can be written as: where x p is the particle position, m p is the mass and ρ p is the density of the particle, u p the instanta- neous velocity of the particle and g the body acceleration.The term F p on the right hand side of the equa- tion represent the force caused by par ticle-fluid interaction.In the most general case it could constitute of the following contributors modelled using empirical formulas (see e.g.Kraipech et al., 2005): The term F D is the Stokes drag law that was de- rived with the assumption of uniform velocity and pressure field (13) The drag coefficient C D depends on variety parame- ters such as the particle size and shape, local Reynolds number, and local fluid density.It is assumed that particle density is sufficiently low that the collective effects on particle drag are unimportant, and the common assumption is made that the particles are spherical.Therefore these assumptions do not allow simulation of dragging mechanism of fine particles by large particles.Such an important mechanism may explain the shape of the selectivity curve in hydrocyclones, and the fact that the bypass value is higher than the water recovery to the underflow.The selectivity cur ve does not have a sigmoidal shape, but exhibits a dip in regions of finer particle size.The force F AM is due to the apparent mass, which is the force to accelerate the apparent mass of the particle relative to the fluid.
The force F BF is due to the temporal delay in the boundary layer around the particle as the relative velocity changes with time.This force represents viscous effects due to the acceleration of the relative velocity. (15) The F LS is the Saffman lift force produced by the pressure distribution developed on the particle due to the rotation induced by a fluid velocity gradient ( γ is the rate of fluid deformation).
( 16) The Magnus lift force F LM is due to the rotation of the particle.This force is caused by a pressure difference between both sides of the particle resulting from the velocity difference due to the rotation. ( The force F PG is due to the pressure gradient in the fluid surrounding the particle. ( The dense flow is characterized by high collision frequencies between particles, and hence their motion is dominantly influenced by particle-particle collisions.Interactions between the fluid and particles are of minor importance (Sommerfeld, 2001).

Other approaches
The third approach for simulating the motion of both the fluid and the solid particles is termed direct numerical simulation (DNS).To extract information implicit in the equations of motion for solidliquid flows, it is necessary to numerically solve the coupled system of differential equations consisting of the equations of fluid motion, and the equations of rigid-body motion (governing the particle motions), together with suitable initial and boundar y conditions.These equations are coupled through the noslip boundary condition on the particle surfaces and through the hydrodynamic forces and torques exerted by the fluid on particles (see e.g.Glowinski et al., 2001).Two classes of methods have been developed in DNS simulations: Arbitrary Lagrange-Euler (ALE) methods on body fitted moving unstructured grids and Fictitious Domain Methods, which allow the flow calculation on a fixed grid.The first approach (Hu et al., 2001) is enormously computationally expensive and at present it is not feasible to model 3D particulate flows with large number of particles of different sizes, shapes and densities.Computing acceleration based on Newton's second law with forces given by the local fluid forces on the particle would require a very fine, moving and deforming mesh for each particle.This is the obvious disadvantage of the approach, which on the other hand may be the only theoretical tool capable of studying the nonlinear and geometrically complicated phenomena of par ticle-par ticle and particle-wall interactions.Glowinski et al. (2001) developed the fictitious domain approach called the Distributed-Lagrange-Multiplier (DLM) method.The key idea is that fluid fills the whole computational domain inside as well as outside the particles boundaries, which eliminates the need for repeated remeshing.The Lagrange multiplier (physically a pseudo body force) is introduced to enforce the interior (fictitious) fluids to satisfy the constraint of rigid body motion, much like the pressure in incompressible fluid flow whose gradient is the force required to maintain the constraint of incompressibility.The method was used to investigate the sedimentation of 6400 disks settling down in a rectangular cavity and the fluidization of 1204 nylon spheres in a slit bed (Pan et al., 2002).
Other recent computational approaches to solidliquid flows, inspired by molecular dynamics, are cellular automata and the Lattice Boltzmann Method (LBM).Originating from discrete kinetic theor y, the latter approach has shown considerable promise and is being pursued by a number of investigators.Fully resolved simulations of colliding spheres were presented by Ten Cate et al. (2004) with a lubrication force used to account for subgrid hydrodynamic interaction between approaching particles.An interesting combination of the LBM and immersed boundary approaches is presented in Feng and Michealides (2005).There are still challenges to overcome in order to demonstrate the LBE schemes as a competitive methodology, compared to the direct solution of the conservation equations of continuum mechanics.These, apart from a limitation on admissible time and length scales, include a consistent modelling of energy transport and elimination of excessive numerical discretization errors (Nourgaliev et al., 2003).The interesting results produced by the LBM methods are not yet sufficiently practicable, for implementation for industrial scale hydrocyclones.
Experimental studies recently published on the influence of the fluid on the solid-solid collisions, investigated the mechanism of the rebound of a particle colliding with a fully immersed wall in a viscous fluid.This has been analysed using a high-speed video camera by Gondret et al. (2002) and Joseph and Hunt (2004).As a result the parameter such as a minimum particle Stokes number, below that no rebound occurs, normal and tangential coefficients of restitutions have been defined.This knowledge could prove to be useful for numerical simulations since it accounts for the combined effects of the interstitial fluid and the inelasticity of contact.Currently, several simulation methods do not allow for particle-wall collisions, because contact would break the lattice modelling of the fluid.
Although many theoretical and computational studies have helped to increase our physical understanding of solid-liquid flows, accurate and reliable predictions of such flows are still limited.The race to create state of the art numerical technique for solid-liquid flows has just acquired "momentum" and will have no shortage of intellectual stimulus due to the wide range of applications.The commercial CFD codes have reached a reasonable level of maturity and are accepted as an analysis tool for one-phase simulations.Despite the enormous advances in the development of numerical techniques a similar potential for solid-liquid flows, and multiphase flows generally, has yet to be fully exploited.

Experimental Validation
The confidence in hydrocyclone flow-field predictions can only be obtained by experimental validation.The validation study should concern the internal flow field.Doby et al. (2007) investigated the swirl flow mechanism in a cylindrical container.They constructed apparatus to enable the lids of a cylinder to be rotated.Due to shear forces, the rotation of the lid induces a flow field on the fluid.The test environment both numerical and experimental was created to conduct a flow pattern analysis.The study allowed validation of the numerical results and verification of experimental LDA technique before replacing the cylindrical container with specially constructed hydrocyclone.The container geometry was used as a benchmark to identify errors in measurements and data analysis.This simple swirl flow generator led to good predictions of a complicated flow.The study provided confidence in using the established tools for examining the flow field in hydrocyclones.The setup enabled incremental changes in the vortex structure to be studied and helped understanding when the fluid develops into unsteady 3D flow.Such an experimental setup could be used to study the performance of second-order closures in order to establish their potential and limitations in the context of swirling flows and possibly to contribute to numerical model improvement.
Laser doppler velocimetr y (LDV) technique has been widely used to measure the velocity patterns in a hydrocyclone.Chine and Concha (2000) used LDV techniques to measure the velocity profile in conical and cylindrical hydrocyclones.The majority of work has been done to measure the 2D flow.Fisher and Flack (2002) have provided an extensive amount of 3D LDA measurements, which could be used for numerical validation of velocity pattern.
Equally impor tant is validation of the hydrocyclone separation performance.The discussion of lowsolid flow and high-solid flow field measurements techniques is presented in Nowakowski et al. (2004).
The high-solids flow measurement precludes nonintrusive application of laser methods due to light extinction.There exist a wide range of applicable techniques, such as electrical impedance tomography, X-ray, Positron-emission, radioactive particletracking, ultrasound and nuclear magnetic resonance to mention a few.In case of electrical impedence tomography knowledge of the inter-relation between the conductivity and concentration fields permits reconstruction of the solids concentration.The challenge however is to form constitutive relations between the particle concentration and size distribution and the bulk conductivity.For example, during the hydrocyclone start-up phase or when operating with air-core instability, temporal conductivity changes are induced by air-core motion.

Conclusions
Advanced computational and experimental techniques are still needed to obtain a more in depth understanding of the complex physical phenomena affecting hydrocyclone performance.The knowledge of the different phenomena, such as particle-particle, particle-fluid, and particle-wall interactions would open the way to describe the particle effects on the suppression or generation of turbulence and on non-Newtonian slurry flows.
The appropriate use of turbulence models in numerical simulations is necessar y for a correct prediction of both the fluid flow field, and the separation performance in a hydrocyclone.This becomes particularly important due to the existence of highly rotational fluid flow of the anisotropic nature with the turbulence not fully developed.Such a situation often makes some of the otherwise very successful turbulence models fail to produce a practical result for the fluid flow in hydrocyclones.Nevertheless the standard computational methods used for hydrocyclone design have significant advantages over the empirical approaches, such as the possibility of changing the geometr y quickly to evaluate possible changes in the separation efficiency.Numerical methods should not be used without experimental verification and the understanding of their limitations.Future work should investigate an appropriate method to model turbulent flow within a hydrocyclone.Large Eddy Simulation (LES) is suggested technique that could be used as this approach solves for the largest scale motions of the flow while approximating or modelling only the small scale motions.It can be considered as a compromise between the Reynolds-Averaged Navier-Stokes equations (RANS) approach and the Direct Numerical Simulation (DNS) approach.LES approach, although still very costly computationally, is becoming more viable option for the solution of transitional and turbulent flow.
Apart from a fully 3D flow, the complex flow behaviour in a hydrocyclone arises from multiphase interactions and the mechanisms governing the formation of an air core.Further investigation on advance modelling, which allows simulation of an air core and multi-phase flow, will be required.For modelling multi-phase flow, an application of a multi-continuum model is particularly attractive because complex interactions between phases can be implemented.The complementar y time scale analysis can provide a quantitative discussion concerning the importance of various types of forces affecting the particle motion.
The accuracy of the prediction of the particle flow could be improved by taking into account the time dependency of fluid flow.
During the last decade, many papers related to the numerical simulation of a hydrocyclone flow were published.As a result, a framework for a 3D flow simulation of complex flow within a hydrocyclone has been laid out, which hopefully will allow the computational techniques to be improved as more research is performed and computing power increases.The primar y purpose of the future research should be focused on providing a numerical approach for simulation of the flow in hydrocyclones, which should be flexible enough to be used for multi-objective design purposes.To develop automatic design methods, the numerical simulations need to be combined with automatic search and optimisation procedures (Nowakowski et al., 2004;Slack, et al., 2004).

Appendix A
The Reynolds Stress Model (RSM) is a second-order closure model.The model represents the solution of the transport equations of the six Reynolds stress components together with the equations for the three mean velocities and the dissipation equation for the dissipation rate, ε, of the turbulent kinetic energy.The variables ρ, p and µ below represent density, pressure and molecular viscosity respectively.The velocity ui is decomposed into its mean and fluctuating components: The Reynolds stress term −ρ u i u j includes the tur- bulence closure, which is modelled in order to close momentum and continuity equations.Reynolds stresses are modelled by the following equation: (A1) The right hand side terms are: The P ij term represents the stress production, D T,ij is the turbulent diffusion term, D L,ij is the molecular diffusion term, Φ ij is the pressure-strain term, F ij is the rotation production term and ε ij is the dissipation term.Here k is the rotation vector, δ is the Kronecker delta and eijk = 1 if i, j, k are different and in cyclic order; eijk = -1.if i, j, k are different and in anti-cyclic order, and eijk = 0 if any two indices are the same.
The dissipation, pressure-strain, and turbulent diffusion terms cannot be computed exactly in terms of the other terms in the equations and therefore are modeled.The pressure-strain term is obtained from a linear-strain model: . (A7) The slow pressure-strain, Φ ij,1 is modelled as Here, variable k represents the turbulent kinetic energy: (A9) The rapid pressure-strain term, Φ ij,2 is modelled as . (A10) The adjustable constants C 1 , C 2 , have the standard values and are equal to 1.8, 0.6, respectively.
The close relation postulated for viscous dissipation term: .

(A11)
The adjustable constants Cµ, C1ε, C2ε σε are 0.09, 1.44, 1.92, 1.0 and 1.3 respectively.The mean strain rate is Eij expressed as .(A12) In the near-wall region of hydrocyclones the solution variables change with large gradients.The popular near wall models (wall functions) do not resolve accurately the near-wall region for strongly swirling flows.Due to the significant stream-wise pressure gradients, body forces and strong secondar y flows the assumption of local equilibrium does not hold.Therefore the Fluent two-layer zonal model was found to be more suitable.The two-layer near wall models divide the whole domain into two regions, a viscosity affected region and a fully turbulent region.The demarcation of the two region is determined by a wall distance-based, turbulent Reynolds number where y is the normal distance from the wall at the cell centres.In the near-wall region when Rey < 200, the one-equation model of Wolfstein (1969) is employed.For the values Rey > 200 the RSM model is used.

Fig. 1
Fig.1 Schematic diagram of a hydrocyclone explaining the principle of operation.

Fig. 2
Fig.2Simulation results obtained from the developed finite element software.Trajectories of mass less particles released from the hydrocyclone inlet(Nowakowski and Dyakowski, 2003).