Settling Suspensions Flow Modelling: A Review

In spite of the widespread application of settling suspensions, their inherent complexity has yet to be properly predicted by a unified numerical model or empirical correlation, and usually industries still possess customized charts or data for their particular suspension. This is, clearly, rather inefficient and can lead to oversized dimensioning, low energy efficiency and even operation limitations/difficulties. In this manuscript a review of empirical correlations, charts and numerical models that have been employed to predict the behaviour of settling suspensions is briefly described, providing information on the advantages and drawbacks of each method. Their evolution throughout the years: from Durand and Condolios correlations, to empirical models by Wasp, single phase simplifications with mixture properties by Shook and Roco, and to other Euler-Euler or Euler-Lagrangian numerical models, will be presented. Some considerations on recent particle migration and turbulence modification publications will be added. In addition, information about some current CFD application of LatticeBoltzmann and Discrete Element Method (DEM) will be given. Lastly, data from CFD modelling employed by the authors that is able to predict turbulence attenuation in settling flows with medium sized particles for different concentrations is reported.


Introduction
A solid-liquid settling suspension is a heterogeneous mixture of solid particles in a liquid, which is a subclass of multiphase flows.Different sized particles, ranging from micro to millimeters, and having diverse densities, can be considered when speaking of solid-liquid settling suspensions.When suspensions contain medium or coarser particles with density higher than the liquid they tend to settle and accumulate at the bottom of the vessel or pipe.These are called settling suspensions.One of the first recorded settling suspensions flow investigations was done in 1906 by Nora Blatch, where pressure drop as a function of flow, density, and solid concentration was accounted for, in a 25 mm (1 in) diameter horizontal pipe (Abulnaga, 2002).The mining industry was the first, amongst many other industries, to have dealt with settling suspensions flows in the mid-nineteenth century.Other examples include not only classical industrial sectors such as paints, oil, cement, coal, drugs and foodstuffs, but also emerging ones as those dealing with "intelligent" materials, biological systems, and also in applications related with environmental remediation processes.Also, in many industrial processes the concentrated solid-liquid mixtures, called pastes, are either subjected to molding as in the case of casting metals, or extrusion, as in the case of ceramics, polymers, or foods, such as pasta.These suspensions are of great practical interest having become ubiquitous in everyday life, either as a natural or formulated product (Abulnaga, 2002;Balachandar & Eaton, 2010).
In spite of the widespread application of settling suspensions, their inherent complexity has yet to be properly predicted by a unified numerical model or empirical correlation, and usually industries still possess custom charts or data for their particular suspension.This is, clearly, rather inefficient and can lead to oversized dimensioning, low energy efficiency and even operation limitations/difficulties.When working with settling suspensions a number of variables have to be accounted for such as flow patterns, transition velocities, the flow behavior in pipes of different geometries, and also particle concentration, shape, size, and size distribution.For concentrated set-tling suspensions, incorporating the modelling of phenomena such as particle-particle interaction, particle-wall interaction, shear-induced migration, turbulence attenuation and augmentation, and lift forces is paramount for a proper suspension behavior characterization.
As pointed by several authors, selecting the material to include in a review paper can become a daunting task, mainly because existing materials in the literature provide a thorough and much appreciated job on the matter.So, adopting a point of view similar to Shook in his review paper from 1976 (Shook, 1976), rather than providing an exhaustive literature review, the choice here was to try to add on existing materials and simultaneously present a different interpretation.Thus, with this manuscript a review of empirical correlations, charts and numerical models that have been employed to predict the behaviour of settling suspensions is done, providing information on the advantages and drawbacks of each method.Their evolution throughout the years: from Durand and Condolios correlations, to empirical models by Wasp, single phase simplifications with mixture properties by Shook and Roco, and to other Euler-Euler or Euler-Lagrangian numerical models, a historical review will be summarized.Some considerations on recent particle migration and turbulence modification publications will be added.In addition, information about some current CFD application of Lattice-Boltzmann and Discrete Element Method (DEM) will be presented.Lastly, results from CFD modelling employed by the authors that are able to predict turbulence attenuation in settling flows with medium sized particles for different concentrations is reported.

Empirical Correlations
One of the first recorded empirical correlations for pressure drop estimation considering fully suspended heterogeneous flows of solid-liquid settling suspensions in horizontal pipes was developed by Durand and Condolios in 1952.This correlation was constructed based on a collection of pressure drop data associated with the flow of sand-water and gravel-water mixtures with particles of sizes ranging from 0.2 to 25 mm.and pipe diameters from 3.8 to 58 cm. with solids concentrations up to 60 % by volume (Aziz & Mohamed, 2013).These studies culminated with the establishment a relation between the pressure drops of water and slurry, given by Eqn. 1, where i and i w are the pressure drop of slurry and of water respectively, k is a constant, C D is the drag coefficient for the free falling particle at its terminal velocity, g is the gravity, D i is the pipe internal diameter, C o is the volumetric concentration of solids and V m is the average flow velocity.
Another important result was the classification of flow regimes, based on particle size and for particles having a specific gravity of 2.65 (Abulnaga, 2002;Aziz & Mohamed, 2013): 1) Particles of a size less than 40 μm are transported as a homogenous suspension; 2) Particles of a size between 40 μm and 0.15 mm are transported as suspension that is maintained by turbulence; 3) Particles of a size among 0.15 and 1.5 mm are transported by a suspension and saltation; 4) Particles of a size greater than 1.5 mm are transported by saltation.
Although being quite useful for narrow sized highly turbulent flows, it fails to account the effect of particle concentration, size and shape.
In 1967 Zandi and Govatos using an extensive number of data points, improved Durand's correlation to different solids and mixtures (Abulnaga, 2002) and defined an index number, Ne, that defined the limit between saltation and heterogeneous flows.While Durand and Condolios based their studies on the drag coefficient, Newit based its work on the terminal velocity as a means to determine the pressure drop (Eqn.2).In 1955 a comprehensive paper was published where a thorough study on solid-liquid flows resulted in several flow regime specific correlations.These correlations, which were not more than a set of criteria, allowed to define the flow regime and theirspecific set of equations where i and i w are the pressure drop of slurry and of water respectively, K 2 is a constant, ρ s is the density of the solids, ρ L is the density of the liquid, C o is the volumetric concentration of solids, V t is the terminal velocity of the particle and V m is average flow velocity.Wasp's correlation from 1977 is based on the assumption total pressure loss is a sum of contributions from both a homogeneous distribution of particles (vehicle part of the flow) and from the excess pressure drop resulting from a heterogeneously distribution of particles (bed formation part of the flow).This correlation can be applied to solid-liquid settling suspensions with varied sized particles, typically in industrial slurries, by splitting the particles' sizes into fractions.The procedure, which was thoroughly described in the literature (Crowe, 2005), for the pressure drop estimation with this correlation is an iterative one where the correlation's critical factors determine: i) the particle size, split between the homogeneous and heterogeneous parts of the flow; ii) the equivalent ho-mogeneous vehicle properties (i.e., density and viscosity) as a function of particle size and concentration.This correlation assumes that the homogeneous vehicle is a Newtonian fluid: in 1980 Hanks extended this approach to account for non-Newtonian properties of the vehicle.
Turian and Yuan's correlation (Crowe, 2005;Peker & Helvaci, 2011), in 1977, extended the pressure loss correlation scheme, by taking into account the fact that various flow regimes are observed depending upon the flow conditions.Their correlation (Eqn.3) utilizes regimespecific coefficients, K, m 1 , m 2 , m 3 , and m 4 to estimate the pressure drop.

 
f and f L are the friction factors for slurry and water, respectively, at the same mean velocity, V m , ϕ is the volumetric fraction of solids and ρ s is the specific gravity of the solids .
To decide on which coefficients to use for each regime the authors proposed a regime delineation scheme based on a regime transition discretization number (Eqn.4).

 
1 In Eqn. 4 K t , n 1 , n 2 and n 3 are coefficients for determining the regime number.
An exhaustive review on the empirical correlations for solid-liquid settling flows is beyond the scope of this manuscript and additional details have been presented in the literature (Crowe, 2005;Lahiri & Ghanta, 2008;Miedema, 2013;Peker & Helvaci, 2011).

Semi-Empirical Model
Acknowledging the limitations of purely empirical methods, researchers devoted their attention to other methods that incorporated both theoretical and semiempirical knowledge as found in the work carried by Bagnold (Bagnold, 1966;Shook & Daniel, 1965).These works had diverse outcomes: one of the most relevant was an equation for energy loss based on the dispersive stress defined by Bagnold, in an attempt to describe solid-liquid settling suspension flow.Some studies (Shook et al., 1968) added on the work done by Bagnold mechanisms describing particle suspension by dispersive stress, incorporating the influence of turbulence suspension of particles using the eddy diffusivity concept together with Richardson-Zaki equation for settling velocity (Shook et al., 1968), which allowed to derive an equation for concentration distribution in steady state.
This correlation is a modification of the individual particle settling velocity, u o , based on an empirical parameter, n, dependent on the flow regime, represented by the terminal Reynolds number (Felice & Kehlenbeck, 2000;Crowe, 2005;Peker & Helvaci, 2011), and also on the ratio between particle and vessel diameter, dp/D.The empirical parameter, n, also known as the Richardson-Zaki exponent has been the subject of several publications.Traditionally this parameter was determined using a set of equations, each defined for a different range of terminal Reynolds number, however: however, these equations can be cumbersome to use, since there are regions where overlapping occurs, and a continuous function was presented as an alternative by Rowe (Rowe, 1987) for the determination of the Richardson-Zaki exponent.Batchelor (Batchelor, 1982;Batchelor & Wen, 1982) extended the work of Richardson and Zaki extending the application of the previous correlation to dilute suspensions (Eqn.6) The hindered settling function is given by h(ϕ) = 1 -nϕ and Eqn.6 provides the settling velocity of randomly dispersed spheres in suspensions in Stokesian regime (Batchelor & Wen, 1982;Peker & Helvaci, 2011).This expression is valid for dilute suspensions where the volumetric fraction of solids is low enough for flocculation to occur.In the studies conducted by this author, the empirical parameter, n, was suggested, for suspensions with negligible interparticular forces, to be 5.5 when Pe number is large and 6.5 for suspensions with very small Pe (Batchelor & Wen, 1982).Additional studies, that the author encourages as further reading, were published by Di Felice (Di Felice, 1999) where a detailed revision on the Richardson-Zaki and other equations for the settling velocity and their dependence with the Reynolds Number is presented.Buscall (Buscall & Goodwin, 1982) conducted several experiments for dilute and concentrated dispersions of polystyrene latex and compared the results with theoretical results from the Richardson-Zaki (Richardson & Zaki, 1954) and Batchelor (Batchelor & Wen, 1982) equations.This study further investigated the empirical parameter, n.Peker & Helvaci (Peker & Helvaci, 2011) also provided a thorough assembly of the values for the Richardson-Zaki exponent in their book.The most impressive aspect of Richardson-Zaki equation is how it effortlessly depicts the complex phenomena in particle-fluid interaction forces and how that is simply compressed into only two parameters: nevertheless, that over a half cen-tury has passed since its inception, it still is used presently by researchers (Kaushal, Seshadri & Singh, 2002;Baldock et al., 2004;Bargieł & Tory, 2013).
It is well known that the most notorious shortcoming of the Richardson-Zaki equation, based on experiments with highly concentrated particles, is failing to take into account the fact that the settling velocity should tend toward zero at the maximum concentration, ϕ max , (Bürger & Wendland, 2001;Kusuda et al., 2007).To circumvent this shortcoming, in some works (Kusuda et al., 2007) the Richardson-Zaki equation has been corrected (Eqn.7).

 
Additionally this equation has been extended to the study of more complex solid-liquid systems behaviors, such as the settling of cohesive particles (Kusuda et al., 2007;Bürger & Wendland, 2001).
For cohesive particles sediments, when flocculation occurs and assuming the size and density of the floc are constant, since some of the effects of concentration on the settling velocity are due to the size of the particle, the volume concentration of flocs, φ, as well as the maximum concentration of flocs, φ max , should be used.Eqn. 8, which can be reduced to the Richardson-Zaki relationship when ϕ < 0.1, is given for the hindered settling velocity of cohesive sediments (Kusuda et al., 2007).
Another model (Eqn.9) for the hindered settling velocity of large concentrations of cohesive sediments is by Winterwerp (Kusuda et al., 2007).

  
Additional models for the settling velocity of aggregated solid particles in a fluid are given in the literature by Bürger and Concha (Garrido, Bürger & Concha, 2000;Bürger, Concha & Tiller, 2000) who employ a one-dimensional phenomenological model (Eqn.10 and 11) for the prediction of the batch sedimentation behavior of solid-liquid suspensions where flocculation of fine particles occurs.This model is dependent on two constitutive material-specific functions, the Kynch batch flux density function, f bk (ϕ), and the effective solids stress, σ e (ϕ), both of which depend only on the local solids concentration and are calculated based on Kynch's (Kynch, 1952) and Richardson-Zaki's (Richardson & Zaki, 1954) works.In Eqn. 10 and 11 Δρ is the solid-fluid mass density difference and p e is the excess pore pressure; other authors employ fractal based equations, for the settling velocity of aggregates (Eqn.12), based on empirical drag correlations (Smith & Friedrichs, 2011;Johnson, Li & Logan, 1996), where θ is a particle shape factor (1 for spherical particles), ρ s the particle density, ρ l the fluid density d p the particle diameter, D f the f loc diameter, η f the fractal dimension, µ the fluid dynamic viscosity and Re p the floc Reynolds number.
Further developments on the use of the Richardson-Zaki equation in depicting the behavior of agglomerates composed of fine micrometric sized particles in a fluidized bed have been published by Valverde & Castellanos (Valverde & Castellanos, 2006).Assuming that the fluidized bed is composed of simple agglomerates, they modified the Richardson-Zaki equation for the settling velocity: Where * o u is the settling velocity of an individual agglom- erate and ϕ * is the volumetric fraction of agglomerates in the fluidized bed, which are defined in Eqn. 14 and Eqn. 15, respectively The number of primary particles is represented by N, u o is the Stokes settling velocity of a single particle and k is defined as R G being the radius of gyration, and d p is the particle size.
Using Eqns.14 to 16, Eqn. 13 can be rewritten as Finally, for the value of n, since fluidized beds of fine particles are usually operated in the small Reynolds number regime, the best choice, in our opinion, is to fix n as 5.6 in agreement with the theoretical derivation in the dilute limit (Valverde & Castellanos, 2006).
Presently, the number of publications regarding solidliquid pipe flow of cohesive particles is scarce by comparison to gas-solid pipe flow of cohesive particles and mostly for laminar flow (Vaezi G, Sanders & Masliyah, 2011;Grof et al., 2009).These works involve fundamental studies of floc and cluster structure interaction.For turbulence studies Toorman (Toorman et al., 2002) published a very interesting paper on turbulence modulation due to the presence of suspended cohesive sediment, and more recently Chai (Chai, Yang & Wang, 2014) presented a model which takes flocculation, sedimentation and turbulent diffusion into account, to analyze the vertical transport of cohesive fine sediment.
Numerical studies are also present in the literature for solid-liquid sedimentation of cohesive particles.Publications can be found which are dedicated to the simulation of flocculation processes for differential settling of cohesive sediments, which were simulated via the Lattice Boltzmann method (LBM) in which the hydrodynamics and attractive van der Waals forces were considered (Zhang & Zhang, 2011).
A more detailed review on cohesive particles is beyond the scope of this manuscript, however, for further readings on the subject the authors recommend the following works found in the literature (Concha & Bürger, 2002;Kusuda et al., 2007;Mantovanelli & Ridd, 2006;Bürger & Wendland, 2001;Garrido, Bürger & Concha, 2000).Most studies in the literature are related to the description of tidal and marine sediments.
Other authors obtained the pressure drop using the dispersion coefficient as a function of local distribution of solids, also describing the settling phenomena making use of the Richardson-Zaki equation (Rasteiro, Figueiredo & Franco, 1993;Rasteiro, Rebola & Scarlet, 1988).Although the coupling of these concepts provided somewhat accurate good comparisons with experimental data, they were only possible to implement through significant simplifications.
An engineer or researcher examining the literature will be overwhelmed by the sheer amount of publications on empirical correlations based on dimension analysis for the critical velocity and pressure loss in settling and nonsettling suspension pipe flows and listing all of them is beyond the scope of this manuscript.Each one of these correlations assumes an enhancement on the quality of results compared to existing publications.Traditionally, empirical correlations have been used to effectively design pipelines; nevertheless, these successful predictions are limited to specific ranges of variables and lack universality since outside the specified range these correlations produce disappointing results.Moreover, for empirical correlations to be effective predictive tools, their coefficients need to be fine-tuned using experimental data from tests in the pipeline system.This is in itself, a logical fallacy, since accurate predictions from empirical correlations to properly design a pipeline need data from that same pipeline.
Some thorough reviews and books on earlier iterations on empirical and semi-empirical correlations for both pressure drop and critical deposition velocity can be found in the literature (Abulnaga, 2002;Crowe, 2005;Peker & Helvaci, 2011;Shook, 1976).

Mechanistic Models
The study of solid-liquid settling suspension flows where a non-homogeneous distribution of particles exists, has provided us with one certainty: any model or correlation accuracy in predicting flow characteristics is intrinsically related to its capability of incorporating the flow regime mechanisms (Crowe, 2005).

Two-Layer Model
In 1970, Wilson (Miedema, Riet & Matoušek, 1995;Wilson, 1970) developed a mechanistic model in which the flow is divided in two layers.In the first layer of the Two-Layer model, the upper part of the flow, the suspended particles linger while in the bottom layer, the second layer, the particles have settled.This model development started with experimental results obtained for narrow particle size distribution of solid-liquid settling suspensions; however it is not suited for cases where there is a low contact of the particles in the bottom layer.In such cases a homogenous model is preferred.One of the issues with the Two-Layer model is that the existence of two layers and an interface inside the pipe is purely conceptual and used only for the sake of numerical representation purposes.
The fundamental bases of this model are as follows (Crowe, 2005): -The flow is divided in two hypothetical layers, an upper layer of particles less than 74 μm and a lower layer containing all particle sizes in the slurry; -Each layer has its own uniform velocity and volumetric solids concentration and there is no slip between the solids and the liquid within the layers; -Since the suspension in the upper layer behaves essentially as a liquid, as far as the wall shear stress is concerned, then the wall shear stress in the upper layer is kinematic, i.e., velocity-dependent; In the lower layer the particles experience a Coulombic friction force; Several iterations of the Two-Layer model have been proposed in the literature that are either simplifications or modifications of the original model in which a stationary or moving bed, is in the bottom layer, and in the upper layer, a heterogeneous suspension with a particle concentration gradient is present (Crowe, 2006).Additional developments were later introduced by other authors (Crowe, 2005;Gillies, Shook & Wilson, 1991;Shook et al., 1982;Shook et al., 1981).A more detailed description of the Two-Layer model, which is beyond the scope of this review, can be found in the literature (Crowe, 2005;Peker & Helvaci, 2011).

Three-Layer Model
The Three-Layer model was introduced in 1995 (Doron & Barnea, 1995) and was developed by joining experimental information with the Two-Layer model.Since the Two-Layer model fails to accurately predict the suspension behavior for low flowrates where a stationary bed is present, this model states that in suspension pipe flow three different flow regimes occur at the same time.In the top layer a heterogeneous flow, in the middle layer a moving bed and in the bottom layer a stationary bed.The additional complexity of the Three-Layer model equips it with the capability of predicting flow patterns transitions, however, due to the aforementioned complexity, supplementary expressions and constitutive relations are required for closure of the equation set (Crowe, 2005;Doron & Barnea, 1995;Ramadan, Skalle & Saasen, 2005).
Mechanistic methodologies are a considerable improvement to empirical correlations in the depiction of settling and non-settling suspension flows.Still, layered models present difficulties in predicting the flow of particles between layers (Roco & Shook, 1984) and, additionally, these models' predictions are still obtained with the help of parameters that require accurate experimental data.Again, this defeats the purpose of the "predictive" aspect intended for a model.In addition, some of the assumptions required for the successful application of mechanistic models may not hold in this case since the assumption that velocity in each layer is uniform in both Two and Three-layer models (Crowe, 2005).

Deterministic Models
Traditional approaches for predicting the behavior of multiphase flows were based on empirical correlations and mechanistic approaches, as seen in previous sections, which resulted from extensive experimental data compiled by equipment designers.These methods had the drawback of being case specific, i.e., they failed to produce accurate predictions if any of the conditions, such as particle data, inlet conditions or geometry, were altered.With the advent of computational modelling techniques and ever evolving computer hardware, the traditional approaches have been refined or replaced, providing scientists, engineers and equipment designers with an enhanced predictive capability and lack of restrictions to adjust process conditions to better suit their demands (Massoudi, 2010).
The number of CFD codes and software, either proprietary or open source, has grown considerably through recent years.Although, single phase CFD codes are well established in the literature, for multiphase flows they are still an open problem (Balachandar & Eaton, 2010;Borhani, 2010;Sommerfeld, Wachem & Oliemans, 2008), in spite of extensive research.When categorizing CFD codes for multiphase flow, more precisely for solid-liquid settling suspension flows, the following approaches have been considered:

Single-Phase Numerical Models
This approach is only suitable for solid-liquid settling suspension flows were the solids concentration is quite low and there is one-way coupling, i.e., where the presence of the particles has little or no impact on the overall properties of the liquid phase (Crowe, Troutt & Chung, 1996).Earlier works derived a turbulent model that used the properties of the mixture for the calculations of settling suspension flow behavior (Roco & Shook, 1984).Two equations turbulence models (Wilcox, 2006) became increasingly popular in two-phase applications, and some recent works have employed this approach to highly concentrated solid-liquid settling suspension flows with turbulence modulation (Bartosik, 2010;Bartosik, 2011).
In other approaches with two-equation single phase turbulence models for solid-liquid settling suspension the authors introduced additional parameters into the equations to incorporate the particle influence on the carrier phase (Hsu, 2003;Jha & Bombardelli, 2009;Rizk & Elghobashi, 1989).
Single-phase numerical models offer a computational inexpensive tool for predicting pressure drop, velocity and turbulence profiles, while providing some insight into particle-boundary layer information, although in a limited fashion, namely now when more complex numerical models are available in software packages, either commercial or open-source.Recently these models have been employed in the study of turbulence modification for concentrated solid-liquid settling suspensions flow (Bartosik, 2011;Bartosik, 2010).

Euler-Euler Numerical Models
In the "Dense Phase approach", Eulerian or even "twofluid" approach, the two components are interacting with each other in a way that the behaviour of each phase in-fluences the other and are considered to be at the same location at the same time.The volumetric fraction is of paramount importance as this variable will dictate the amount of each phase at a given time and place.The Eulerian models provide an averaged depiction of a multiphase system, and in the literature a wide range of averaging processes can be found, namely time, volume or ensemble based averaging (Ishii & Hibiki, 2011).With this averaging methodology two advantages arise, one being that with the averaging process all the forces are inherently present in the model, the other advantage is that the computational cost is not dependent on the number of particles, making the Eulerian approach more suited for large systems with a great number of particles.A drawback of the Eulerian modelling approach, also a consequence of the averaging process, is the loss of detail, which creates the need for closure equations for the turbulence and interaction forces.This approach is widely employed in fluidization, gas-solid flows, pneumatic and hydraulic conveying, and suspension flows (Balachandar & Eaton, 2010;Sommerfeld, Wachem & Oliemans, 2008).The Eulerian approach has become increasingly popular for concentrated or dense suspension flows, either using a single fluid approximation (Mixture Model, Volume of Fluid, cavitation models, etc.) or two-fluid approximation (Euler-Euler Model or Euler-Granular Models) (Brennen, 2005).
The Mixture Model (Manninen, Taivassalo & Kallio, 1996) was used to perform a series of numerical studies on pipe flows of both zircon-water and silica-water mixtures to a maximum of 20 % solid volumetric fraction: all of them showed good agreement with the experimental data (Ling et al., 2003); this same model was also used together with the Standard/High Reynolds k-ε Turbulence Model for highly concentrated solid-liquid flow in pipes (Kaushal et al., 2012) with results far from satisfactory due to an over-prediction of the pressure drop in the pipe section, which increased with solids concentration.Another approach, using the Mixture Model and a Low Reynolds Turbulence closure, was employed to describe highly concentrated flows of solid-liquid suspensions (Silva et al., 2014;Silva et al., 2013).A new photocatalytic reactor (XiaoWei & LieJin, 2010) with solar concentrator for hydrogen production was simulated using an Algebraic Slip Mixture model (ASM) with a catalyst volumetric fraction up to 15 % .
An increasing number of publications where the Two-Fluid approach incorporating the Kinetic Theory of Granular Flow is employed to characterize particle-particle interaction can also be found in the literature to study highly concentrated solid-liquid settling suspensions pipe flow (Kaushal et al., 2012;Lahiri & Ghanta, 2008;Lahiri & Ghanta, 2010), slush nitrogen (Jiang & Zhang, 2012), ice slurry (Wang et al., 2013), with good results in repro-ducing experimental data.The Kinetic Theory of Granular Flow is an adaptation form the Kinetic Theory of Gases, and in this way the particle-particle interactions are quantified in the flow.However, although it provides some good results for concentrated solid-liquid flows for different size and density of the particles, it is a very complex numerical model, with a high computation requirement and often with boundary conditions issues that require some simplifying assumptions.
The predominant applications of the existing twofluid models exhibit problems hindering their use for more complex flows of engineering interest.Amongst the main issues one can point out numerical instabilities, very timeconsuming, difficulty in dealing with the complex geometries since the calculation time becomes prohibitively expensive and none of the existing models has shown to be able to determine the minimum in the pressure gradient versus slurry velocity, which characterizes the transition to bed flows (Messa, Malin & Malavasi, 2014).

Euler-Lagrange Numerical Models
The Lagrangian approach, also known as "Dilute Phase approach", is employed when the amount of the dispersed phase is small and does not disturb the motion of the continuous phase.This approach is predominant in case studies of sprays, atomization and flows with bubbles, where droplets and particles which are treated as the dispersed phase.Amongst the Lagrangian approach three major modelling techniques can be outlined: "Point-Particle Direct Numerical Simulation (DNS)", "Point-Particle Large Eddy Simulations (LES)" and "Point-Particle Reynolds Averaged Navier-Stokes (RANS)".The DNS modelling approach requires that the particles must be smaller than the Kolmogrov scales, i.e., the time scales of the particle have to be smaller than the time scales of the smaller scales of the fluid.This requirement limits the DNS application to very low Reynolds numbers or to very small particles.To overcome this limitation LES modelling can be used.Both DNS and LES application are limited to dilute systems where collisions and hydrodynamic interactions can be neglected and a one-way coupling between the dispersed and carrier phases is assumed (Balachandar & Eaton, 2010;Fox, 2012;Hiltunen et al., 2009;Mashayek & Pandya, 2003).
In recent studies (Adams, Fairweather & Yao, 2011;Soldati & Marchioli, 2012) one-way coupled Eulerian-Lagrangian models where employed in the study of dilute solid-liquid flows of suspension and re-suspension of particles; a solid-liquid settling suspension flow in horizontal pipes was investigated (Capecelatro & Desjardins, 2013) for operating conditions above and below the critical deposition velocity.A high-fidelity large eddy simulation framework is combined with a Lagrangian particle tracking solver to account for polydispersed settling particles in a fully developed turbulent flow.Two cases were simulated, the first having a Reynolds number of 85 000 and the second considers a lower Reynolds number of 42 660.Since most studies of the Lagrangian properties of turbulence trace point-like particles (Toschi & Bodenschatz, 2009) they cannot still be generally applied to all types of particles, due to the high computational demand, which hinders their application for processes with a large number of particles.They will, nevertheless, most likely become standard tools in the future.
Some excellent reviews on Lagrangian-Eulerian Methods on multiphase flows were given by Subramanian (Subramaniam, 2013) and Zhou (Zhou, 2010).

Lattice-Boltzmann Numerical Models
In the last two decades the lattice Boltzmann method (LBM) has been developed into an established CFD approach for solving fluid flow problems.Important developments have been done in LBM's capability for several flow problems, containing multiphase flows, turbulence, and microfluidics.Amongst the numerous areas, the solidliquid systems have received special emphasis considering the unique advantage of LBM in its computational efficiency and parallel scalability.Traditionally, conventional CFD numerical schemes are based on discretisation's of macroscopic continuum equations, like finitedifference, finite-element or finite-volume methods, have been used to solve the velocity and pressure field from Navier-Stokes equations: on the other hand LBM is based on microscopic models and mesoscopic kinetic equations in which the fluid is described by a group of discrete particles that propagate along a regular lattice and collide with each other.This scheme is particularly successful in fluid flow applications involving interfacial dynamics and complex boundaries (Aidun & Clausen, 2010;Chen et al., 2010;Yu & Fan, 2010).
The LBM can serve as an alternative flow solver for different types of incompressible flows.The incompressible Navier-Stokes (NS) equations can be obtained in the nearly incompressible limit of the LBM.Pressure, in the LBM, is calculated using an equation of state.In contrast, in the direct numerical simulation of the incompressible NS equations, the pressure satisfies a Poisson equation with velocity strains acting as sources and solving this equation often produces numerical difficulties requiring special treatment, such as iteration or relaxation (Aidun & Clausen, 2010;Chen et al., 2010;Yu & Fan, 2010).For the modelling of solid-liquid systems, the LBM, due to its simple implementation, becomes particularly appropriate for simulations involving large numbers of particles.Furthermore, it can be coupled, if it's regarded only as a solver for the fluid flow, with various methods for particles such as Discrete Element Method (DEM) or Lagrangian tracking (Aidun & Clausen, 2010;Balachandar & Eaton, 2010;Borhani, 2010).
Among recent publications on LBM application in solid-liquid settling suspensions flows, it is important to refer: Shardt & Derksen (Shardt & Derksen, 2012) simulations of up to 45 % solids volume fraction of rigid non-spherical particles with low density ratios at moderate Reynolds numbers (< 1) using the LBM coupled with DNS studies; two-dimensional (2D) and threedimensional (3D) CFD studies of solid-liquid settling suspensions flows, by Kromkamp et al. (Kromkamp et al., 2006) where Couette flows of single, two and multi-particle systems were conducted; Gao et al. (Gao, Li & Wang, 2013) particle-resolved simulation method for turbulent flow laden with finite size particles where the method was based on the multiple-relaxation-time lattice Boltzmann equation.In this case, a maximum of 51200 particles in 3D have been considered in their simulations and the authors note that particle-laden turbulent flow is a multiscale problem that requires state of the art computers to include all relevant scales into the simulations with realistic physical parameters.
The Lattice-Boltzmann approach, due to its relatively simple implementation for parallel computing and hybrid combinations of the Eulerian lattice with a Lagrangian grid system (Aidun & Clausen, 2010;Subramaniam, 2013), shows great promise relatively to traditional approaches, still it is only possible to employ it in the simulations of suspension flows with dilute concentrations.

Discrete Element Method (DEM) Numerical Models
One of the main challenges in simulating settling suspensions flows derives from their intricate behaviour brought about by the complex interactions between individual particles and their interactions with surrounding liquid and wall.Understanding the underlying mechanisms has been the aim of particle scale research which in recent years has grown worldwide, a result from the intense development of both discrete particle simulation techniques and computational capabilities.An important discrete model is the discrete element method (DEM) originally developed in 1979 by Cundall and Strack (Cundall & Strack, 1979).This method uses the Newton's equation of motion to contemplate a finite number of discrete particles interacting through contact and noncontact forces moving translationally and rotationally.Both trajectories and transient forces acting on individual particles are extremely difficult to obtain by experimental techniques, is the type of information provided by DEM simulations (Zhu, Zhou & Yang, 2008).
In recent publications DEM has been combined with CFD techniques to describe solid-liquid settling suspensions: modelling solid-liquid suspension flows in the density-driven segregation of a binary particulate suspension incorporating 10 000 particles in a closed container, using a hybrid combination of the discrete element method (DEM) with computational fluid dynamics (CFD) (Qiu & Wu, 2014); using a new Lagrangian-Lagrangian algorithm, also referred to as the DEM-SPH method, for solid-liquid flows in both a dam break problem and a quasi-steady solid-liquid flow in a cylindrical tank (Sun, Sakai & Yamada, 2013); simulating dense medium cyclone separation (DMC), combining DEM with CFD in coal preparation, which is a process with a simple design but where the flow pattern within is complex, due to the size and density distributions of the feed as well as the turbulent vortex formed (Chu et al., 2009).In this case study, DEM is used to model the motion of discrete particles by applying Newton's laws of motion and CFD is used to model the motion of the slurry medium by numerically solving the local-averaged Navier-Stokes equations using the Volume of Fluid (VOF) and Mixture multiphase flow models; transported in a fluid, for predicting the location of the puncture point location of the particles in an elbow, DEM was used to describe the kinematics and trajectory of the discrete particles as well as the particleparticle interaction while the hydrodynamic modelling of the fluid phase was based on the volume-averaged Navier-Stokes equations, and a fluid density-based buoyancy model was adopted to calculate the solid-fluid interaction force (Zhang et al., 2012).
Regarding numerical studies involving DEM for solidliquid settling suspensions in horizontal pipe there seems to be a lack or complete absence of publications on the subject, which can be attributed to the limitation on the number of particles that it is possible to simulate even with this method that offers an alternative to DNS due to its parallel computing capability (Zhang et al., 2012;Zhu, Zhou & Yang, 2008).A review on DEM application in particulate systems was given by Zhu et al. (Zhu, Zhou & Yang, 2008).

Turbulence Modulation
As pointed out in section 2.4.1, with low solids volumetric fraction the usual assumption is that the turbulence of the fluid phase is equal or very similar to the singlephase flow.Yet, as the solids volumetric fraction increases additional phenomena appear where turbulence augmentation, dissipation and distortion become significant.For solid-liquid settling suspensions flows the phenomena of turbulence attenuation is a rather interesting one for design engineers, since this would allow solids conveying of concentrated suspensions at energy expenditures similar to those of single-phase f lows (Balachandar & Eaton, 2010;Crowe, Troutt & Chung, 1996).
Crowe and Elghobashi have done extensive work on turbulence modification, but mainly for gas-liquid and gassolid flows (Crowe, 2000;Kenning & Crowe, 1997;Elghobashi & Truesdell, 1993;Truesdell & Elghobashi, 1994).In his studies, the following conclusion was attained "small particles will attenuate the turbulence while large particles will generate turbulence" (Crowe, Troutt & Chung, 1996).And while this seems to hold true for gassolid and gas-liquid suspensions flows, recent studies (Kaushal et al., 2012;Lahiri & Ghanta, 2010;Matoušek, 2005) seem to contradict this statement for solid-liquid settling suspensions flows.In fact, quite the opposite seems to be the case for solid-liquid settling suspensions but only for highly concentrated solids volumetric fractions.
In a recent publication Tanaka (Tanaka & Eaton, 2008) presents a dimensionless parameter, the particle moment number, Pa, that was derived using dimensional analysis of the particle-laden Navier-Stokes equations.This analysis was based on a set of 80 experimental measurements where the turbulent kinetic energy was modified by particles.Data for the turbulent kinetic energy augmentation in air and water was included as well as data for the turbulent kinetic energy attenuation in air, but there is a void of data for the turbulent kinetic energy attenuation for water.This is a very thorough study that in spite of the absence of information on turbulence attenuation when the medium is water represents a significant step towards predicting turbulence modification in particle laden flows.
Searching the literature for current numerical studies trying to characterize on turbulence attenuation for solidliquid settling suspensions, some manuscripts are found where: drag correlations are modified in an attempt to reproduce experimental data where turbulence modulation occurs (Hadinoto, 2010;Hadinoto & Chew, 2010); singlephase Low Reynolds turbulence models are modified to incorporate turbulence modulation (Bartosik, 2010;Bartosik, 2011); Mixture Models with a Low Reynolds Turbulence closure for turbulence modulation in highly concentrated solid-liquid flows (Silva et al., 2013); Euler-Euler model is used for pipe flow of concentrated slurries (Kaushal et al., 2012;Lahiri & Ghanta, 2010).

Concentrated Solid-liquid Settling suspensions Flow in a Horizontal Pipe: Mixture Model Studies
The authors of this review have also been working on numerical studies of solid-liquid settling suspensions flows in horizontal pipes, with highly concentrated solids volumetric fractions up to 40 % (v/v), using data from the literature for validation (Lahiri & Ghanta, 2010).CFD numerical studies where a Mixture Model, incorporating a Low Reynolds turbulence closure (Costa, Oliveira & Blay, 1999;Hrenya & Bolio, 1995) and using a Schiller-Naumann drag correlation (Pang & Wei, 2011), were used to describe experimental data from the literature on solid-liquid suspension flow of medium sized particles where wall turbulence attenuation phenomena is observed (Fig. 1).The focus of our studies was mainly on higher flow velocities and concentrations, since in an industrial environment working with flow velocities where a stationary or moving bed is observed is not desirable, and working in a regime where turbulence attenuation occurs is preferred, since there is an energetic advantage.A recent publication on energy consumption of solid-liquid settling suspensions flows demonstrates the interest on this subject (Jafari, Tanguy & Chaouki, 2012).
Previous numerical studies employing the Mixture Model with a High/Standard k-ε Turbulence Model including a traditional wall function, which is an empirical approximation employed in single-phase models, proved to be inadequate for the numerical study of highly concentrated solid-f lows as also seen in the literature (Ekambara et al., 2009;Troshko & Hassan, 2001).Thus, incorporating a Jones-Launder Low Reynolds Turbulence closure in the Mixture Model, for the settling particles studied, can circumvent this issue, since the Low Reynolds Closures solve the model through the integration of the turbulence equations until the wall.The numerical and experimental pressure drops for different concentrations and flow velocities are shown in Fig. 1, as well as the results using Durand and Condolios correlation.
There is a good agreement with the experimental pressure drop data for flow velocities from 3 up to 5 m/s for all solid concentrations.This improved agreement is more notorious for the highest flow velocities, which is the area of interest, since for 5 m/s there is a similarity between this pressure drop and the monophasic pressure drops (the region where turbulence attenuation is expected as explained in the literature (Lahiri & Ghanta, 2010;Tanaka & Eaton, 2008).A possible explanation for this phenomenon is related with lift-forces resulting from the viscousturbulent interface at the bottom layer of particles closest to the pipe wall (Kaushal & Tomita, 2007;Matoušek, 2005).
Overall, the tendencies detected in the numerical and experimental vertical solid volumetric fraction profiles (Fig. 2) are in concordance.The general tendency of the solid concentration lines is followed by the numerical ones, with the exception of the lower flow velocities, where there is a slight deviation, especially in the bottom region of the pipe cross-section.This can be again attributed to the moving bed regime that the Mixture Model does not represent accurately (Ling et al., 2003).With the inclusion of the Jones-Launder Low Reynolds Turbulence closure in the Mixture Model, the deviations in the solids volumetric fraction profiles for the highest velocities decrease significantly and overall the deviations are smaller when compared with the High Reynolds Turbulence closure, especially for the highest average concentrations.For the lowest solid volumetric concentration the High Reynolds Turbulence closure behaves well for all the velocities, with a better fit for velocities up to 3 m/s.For these lower velocities the Mixture Model with the inclusion of the Jones-Launder Low Reynolds Turbulence closure behaves poorly by comparison with the High Reynolds Turbulence closure, which can be attributed to the fact that the suspension viscosity is higher at the pipe bottom and since the Low Reynolds Turbulence closure integrates the turbulence equations all the way to the wall, pressure drop will rise accordingly, due to an over- shoot of the suspension and turbulent viscosity with this closure, which, for this concentration, is not compensated by the presence of strong lift forces resulting from particle/particle interactions in that region, as happens for the higher concentrations.

Future Directions
In the near future, with the continuous development of computational architectures, more complex DNS numerical studies will be possible allowing the widest application of current Eulerian and Lagrangian numerical models and providing further information on turbulence modula- tion and particle migration phenomena.Additionally, interphase correlations (drag, lift, etc.) can be further improved by DNS data.Euler-Lagrange numerical methods will become more widespread with enhanced computational structures and algorithms (LBM and DPM) for more concentrated suspensions and complex flows.
Thus, with the improvement of DNS and Lagrangian methods, providing microstructural data for settling suspensions flows, as well as further development of computational capability, Euler models coupled with that information will become more reliable and faster, allowing design engineers to obtain fast and accurate representations of the flow, thereby reducing the cost of empirical data acquisition and optimizing the energy consumption in pipe flow, thereafter minimizing the need for empirical correlations and charts.
Finally, it is paramount to have quality experimental data for settling suspensions, regarding turbulence data for both dispersed and carrier phases, which requires further development of experimental techniques.

Fig. 1
Fig. 1 Pressure drop comparison between experimental, Durand & Condolios correlation and simulated results, using the Mixture Model with a Low Reynolds k-ε Turbulence model, for settling particles with diameter of 0.44 mm, density of 2470 kg.m -3 and volumetric fractions from 0.1-0.4.

Fig. 2
Fig. 2 Comparison of normalized (using the efflux volumetric fraction, C VF ) experimental and numerical (Mixture Model with the Low Reynolds k-ε Turbulence model) vertical solid concentration profiles for settling particles, (a) 0.1 volume fraction, (b) 0.3 volume fraction, (c) 0.4 volume fraction.

CN
VF efflux volumetric concentration D i pipe internal diameter (m) D f the floc diameter (m) d p particle diameter or size (m) on the flow regime n 1 coefficients for determining the regime number in Eqn. 4 n 2 coefficients for determining the regime number in Eqn. 4 n 3 coefficients for determining the regime number in Eqn. 4 n 4 coefficients for determining the regime number in Eqn. 4 an individual agglomerate (m.s -1 )