Simulation of particles and sediment behaviour in centrifugal field by coupling CFD and DEM

Chemical Engineering Science 94 (2013) 7–19

Contents lists available at SciVerse ScienceDirect

Chemical Engineering Science journal homepage: www.elsevier.com/locate/ces

Simulation of particles and sediment behaviour in centrifugal field by coupling CFD and DEM Xiana Romanı´ Ferna´ndez n, Hermann Nirschl Institute of Mechanical Process Engineering and Mechanics, Karlsruhe Institute of Technology (KIT), Straße am Forum 8, 76131 Karlsruhe, Germany

H I G H L I G H T S c c c c

Simulation of the multiphase flow in a centrifuge using a coupling of CFD and DEM. Analysis of the solids volume concentration in the sediment formed. Analysis of the effect of the particles on the flow pattern. Sediment dependence on the rotational speed, particle size and friction coefficient.

a r t i c l e i n f o

abstract

Article history: Received 26 May 2012 Received in revised form 19 January 2013 Accepted 15 February 2013 Available online 26 February 2013

Centrifugal separation systems, such as solid bowl centrifuges, are used to effectively separate fine particles from industrial fluids. Knowledge of the streams and sediment behaviour inside such solid bowl centrifuges is required to determine the geometry and the process parameters that lead to an optimal performance. The multiphase flow of water, air and particles is calculated numerically in a cylindrical solid bowl centrifuge using a Computational Fluid Dynamics (CFD) software coupled with a Discrete Element Method (DEM) software. The CFD software calculates the flow of water and air and transfers the hydrodynamic forces to the DEM software in order to calculate the particle positions adding the hydrodynamic effect to the contact forces. Afterwards, the influence of the particles on the flow is considered. Thus, an unsteady and coupled simulation calculates flow and particles alternatively. Different sediment forms and solid volume concentrations are obtained when varying the geometry, operating conditions, particle size and contact model parameters such as the friction coefficient. & 2013 Elsevier Ltd. All rights reserved.

Keywords: Centrifugation Sediment Simulation Computational fluid dynamics Discrete element method

1. Introduction Centrifugal separation of particles in a suspension is one of the most common problems appearing in industry such as chemical processing, mining and mineral processing, solid-fuel industry, food processing, waste water treatment, pharmaceutical processes and biotechnology. Centrifuges are, together with sedimentation tanks and filters, one of the most commonly employed devices for mechanical solid–liquid separation. Sedimentation uses as its driving force a mass force based on the different densities of the phases. In sedimentation tanks this driving force is gravitational force while in centrifuges it is centrifugal force. Filters, instead, make use of differential pressure as a driving force. The choice of which technique should be applied for a

n

Corresponding author. Tel.: þ49 721 608 42404; fax: þ 49 721 608 42403. E-mail addresses: [email protected], [email protected], [email protected] (X. Romanı´ Ferna´ndez). 0009-2509/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ces.2013.02.039

specific separation task depends on the properties of the mixture, the purpose of the separation, whether the thickening of the solid phase or the clarification of the liquid, and the available time and space. The rapid development of the centrifugal technology nowadays allows high efficiencies separating very small particles, even in the nanometre range, from multiphase flows. Among other centrifugal devices as decanters and disk stack separators, which are more complex, solid bowl centrifuges can be used to reach the mandatory clarifying purity. Despite centrifugation being a well-established operation, the knowledge of the complex phenomena occurring inside of centrifuges is limited. The prediction of the separation efficiency of a centrifuge is often restricted to the cut size, the smallest particle size that can be separated. This prediction is usual based on a certain flow pattern assumed for the centrifugal device, which can be erroneous. Moreover, the flow pattern rarely accounts for the influence of already settled particles forming a sediment that can interact with the flow pattern changing it and thus changing the settling conditions for subsequent particles.

8

´ndez, H. Nirschl / Chemical Engineering Science 94 (2013) 7–19 X. Romanı´ Ferna

The investigation of the sediment behaviour is not only important during the centrifugation performance due to the influence in the separation efficiency but also for the discharge of the solids. In industrial centrifuges the discharge of the solids is carried out by scrolls, peelers or nozzles and its efficiency depends on the sediment behaviour. Due to the increasing computing resources and the variety of simulation models, Computational Fluid Dynamics (CFD) has become an important tool to study the flow within industrial machinery, and thus within centrifuges. The simulation of the multiphase flow inside a solid bowl centrifuge using CFD has already been object of research (Romani Fernandez and Nirschl, 2009,2010; Romani Fernandez et al., 2011). In these studies, after the continuous phase had converged, the particulate phase was simulated neglecting particle–particle interactions and taking into account only the hydrodynamic forces. Thus, it was possible to track the trajectory of the particles and then calculate the separation efficiency of the simulated centrifuge depending on the particle properties and the operating conditions. However, in order to simulate the sediment build-up and its behaviour, it is necessary to consider particle–particle interactions. For this purpose it is appropriate to apply the Discrete Element Method (DEM) (Cundall and Strack, 1979), originally developed for the simulation of granular materials and then extended to the whole particle processing technology. This method accounts for particle–particle and particle–wall interactions by means of various contact models. Other forces such as electrostatic or Brownian forces can be also implemented when necessary. In centrifugal field, the hydrodynamic forces (drag, lift) play a crucial role in the movement of the particles and must be included. Hence, a coupling of the DEM simulation with the CFD simulation is necessary. For the simulation of the fluids a commercial CFD software ANSYS FLUENTs (Version 12.0.16) was used. For each position in the CFD computing field, the velocity, pressure, density and viscosity of the flow must be transferred to DEM software EDEMs, produced by DEM Solutions Ltd. The particles settle on the bowl wall forming a sediment which influences the flow. Thus, the information regarding the particle positions must be passed on to the CFD software in order to recalculate the flow considering the pressure drop caused by the sediment. This complex and coupled system is calculated numerically. The advantage of this method is that the internal flow, the particle trajectories and the deposition of the particles can be represented accurately. It allows a detailed description of the flow and particle behaviour under various operating conditions, which is of high importance for the design, optimisation and laying-up of centrifuges. Nevertheless, the simulation results should be examined carefully from a scientific point of view and must be validated with experimental work. The simulation of the sediment in a centrifuge is of great importance to investigate how it influences the flow pattern and thereby the settling of the subsequent particles. The object of this study is to analyse to what extent is it possible to predict the build-up and the behaviour of a sediment inside a centrifuge using the coupling between CFD and DEM for the numerical simulation. It is also important to specify which parameters are necessary for the modelling as well as what is its influence on the results. Thus, a study regarding operating parameters, particle system properties and contact model parameters was performed. The calculations were made using the software ANSYS FLUENT 12.0.16 and EDEM 2.3, which include an additional DEM-CFD coupling module for FLUENT provided by DEM Solutions Ltd. FLUENT is responsible for the computation of the multiphase flow of liquid and air and EDEM calculates the particle positions including the hydrodynamic forces and contact forces. Some

attempts at simulating multiphase flows in filters using this ¨ coupling methodology have already been carried out (Muller, 2009; Schilling et al., 2009). Other authors have already simulated complex multiphase flows in other equipment, such as fluidised beds (Deen et al., 2007), filters (Dong et al., 2003; Ni et al., 2006) and hydrocyclones (Chu et al., 2009; Chu and Yu, 2008). Nevertheless, due to the high velocity gradients and complex interactions between the phases, the simulation of multiphase flow and in particular the sediment behaviour in centrifugal field is still a challenge. The following chapters explain, the general methodology used for particulate multiphase flows and specifically the equations applied here including the coupling between both simulation software. Subsequently, the centrifugal geometry used and the results of the obtained sediments are presented.

2. Model description There are several multiphase models available in order to simulate a particulate phase. A summary of the different simulation methods for multiphase flows in CFD can be found elsewhere (van Wachem and Almstedt, 2003). The multiphase model applied in this study for the particulate phase is a Euler– Lagrange method. The Euler–Lagrange methods solve the continuum conservation equation just for the characterisation of the continuous phase. The dispersed phase is described as mass points and their trajectories are determined by integrating the force balance for each individual particle. This way, instead of a partial differential equation, just an ordinary differential equation must be solved for the dispersed phase. The Euler–Lagrange methods are therefore numerically simpler than the Euler–Euler methods, but to represent a dispersed system statistically, a sufficiently large number of particle trajectories must be calculated. The weakness of the Euler–Lagrange method is that it is only valid for diluted suspensions. A maximum value for the volume concentration of the disperse phase lies between 5% ¨ et al., 2007) and 10%–12% (ANSYS, 2009), over this value (Schutz the interaction between particles should be considered. The nature of the problem and the information that should be obtained from the simulation usually decides which of the two methods should be chosen. The Lagrangian approach for the discrete phase applies the second Newton’s motion law on the particles. Different forces acting on the particles can be considered: drag forces, gravitational forces, centrifugal forces, electro-magnetic forces, Brownian forces, etc. In order to account for the mutual influence of the continuous and discrete phase during the simulation, a four-way coupling was applied. The four-way coupling involves the interaction within the disperse phase: impacts between particles, particle–wall interactions, van der Waals forces, electrostatic forces, etc. Thus, the restriction for the volume fraction is not valid here anymore. The Discrete Element Method (DEM) (Cundall and Strack, 1979) is a Lagrangian method originally developed for the simulation of bulk solids and then extended to the whole particle technology. This method solves the Newton’s equations of motion for translation and rotation of the particles. In these equations the impact forces of particle–particle and particle–wall interactions are included. Other forces such as electrostatic or Brownian forces can also be implemented if necessary. For describing the collisions of particles a soft sphere model was used. In the soft-sphere approach the particle interactions are modelled through a potential force. When two particles collide, they actually deform. This deformation is described by an overlap displacement of the particles. The repulsive force after a collision is then calculated as a function of overlap. This approach uses a

´ndez, H. Nirschl / Chemical Engineering Science 94 (2013) 7–19 X. Romanı´ Ferna

spring-damping system to calculate the behaviour of the particles during the contact time. The soft-sphere model allows contacts of a particle with more neighbouring particles during a time step. The net contact force acting on a particle is added in the case of multiple overlaps. This model for contact forces was first described by Cundall and Strack (1979), who proposed a parallel linear spring-damping model for the normal direction and a parallel linear spring-damping in a series with a slider for the tangential direction. A comparison of the performance of different soft-sphere models can be found in the literature (Di Renzo and Di Maio, 2005; Stevens and Hrenya, 2005). While the DEM model was originally used in environments absent of fluid, increasing efforts have been made to extend this model to the mixture of fluid and particles. For that, a coupling of the DEM calculations with the fluid flow simulation is necessary to be able to include the hydrodynamic forces in the particle simulation. An important aspect in the calculation of particles in a fluid is the energy dissipation into the surrounding fluid. In viscous fluids there are, in contrast to air, two additional effects that lead to energy losses. On the one hand, the drag force leads to lower particle velocities. On the other hand, when a particle is about to hit a surface or another particle a slowdown occurs. This is due to the drainage of the fluid, which is located in the gap between the two particles or the particle and the surface. Zhang et al. (1999) introduces a correction factor for the resistance force in the calculation of particle paths to reduce the impact velocity. Indirectly, this effect is also taken into account by Legendre et al. (2006), who considers the decrease of a particle velocity shortly before the collision happens, as part of the bouncing mechanism. In his work, as well as in the studies of Gondret et al. (2002), the coefficient of restitution e is presented as a function of Stokes number St, which describes the ratio of the dynamic response time of a particle to flow changes and the characteristic flow time. Both authors determined empirically a critical Stokes number of about 10. Below this value a particle do not bounce when it hits a wall and, instead, it remains in contact with it.

9

external volumetric force. All terms are discretised and calculated for each volume cell with the exception of tt, the tensor of turbulences, which is modelled via a k-e model: the k-e Renormalisation Group (RNG) model (Yakhot and Orszag, 1986). This model is more reliable for a wider class of flows than the standard k-e model. 2.2. Mathematical formulation of the particles The particle–particle interactions during settling and in the sediment are considered by using DEM, which calculates the translation, Eq. (4), and rotation, Eq. (5), of a particle based on the forces and momenta acting on them: ! dvp 1 X ¼ Fi , ð4Þ mp dt ! do p 1X ¼ Mi , Ip dt

ð5Þ

where vp and op are the linear and angular velocities of the particles, mp and Ip are the mass and the momentum of inertia of the particle respectively. Mi represents the momenta and Fi the different forces applied on each particle. The Hertz–Mindlin contact model (Mindlin, 1949), a non-linear soft-sphere model based on the Hertzian contact theory (Hertz, 1881), implemented in EDEM was used to model the particle– particle and particle–wall contacts (DEM Solutions, 2010b). Electrostatic and van der Waals forces were not considered here because no submicron particles were simulated with this method. Fig. 1 represents the forces in tangential and normal direction appearing when two particles collide. The normal force Fn is a function of normal overlap dn, equivalent Young’s modulus E* and equivalent radius R*: 4 n pffiffiffiffiffin 3=2 E ð6Þ R dn , Fn ¼ 3 2

ð1n2i Þ ð1nj Þ 1 þ , n ¼ Ei Ej E

ð7Þ

1 1 1 ¼ þ , Ri Rj Rn

ð8Þ

2.1. Mathematical formulation of the fluid The Volume of Fluid (VOF) method (Hirt and Nichols, 1981) is designed to track the interface of two phases that are not interpenetrating, and hence adequate to simulate the gas–liquid multiphase flow. The presence of solid particles was ignored for fluid simulation purposes, which is an acceptable assumption for low solid concentration. A continuity conservation equation, Eq. (1), is solved for each phase, and the volume fraction of each phase aq in any cell has to obey Eq. (2), @ ! ðaq rq Þ þ rUðaq rq v q Þ ¼ 0, @t n X

aq ¼ 1:

with Ei, ni, Ri and Ej, nj, Rj being the Young’s modulus, Poisson’s ratio and radius respectively of each sphere in contact.

ð1Þ

ð2Þ

q¼1

As the flow inside the bowl is non-laminar, the Reynolds Averaged Navier–Stokes equation, Eq. (3), was chosen to solve the velocity field. With the VOF method a single momentum equation is solved throughout the domain, and the resulting velocity field is shared among the phases. The dependency on the volume fraction is implemented by using volume averaged values for the density r and the viscosity m, ! @ ! ! ! ðr v Þ þ rðÞ v r v ¼ rp þ rt þ rtt þ F : @t

ð3Þ

In Eqs. (1) and (3), v represents the velocity, p the pressure, t the shear stress tensor, tt the tensor of turbulences and F an

Fig. 1. Contact forces between two spherical particles in a soft sphere model.

´ndez, H. Nirschl / Chemical Engineering Science 94 (2013) 7–19 X. Romanı´ Ferna

10

The damping force in the normal direction is given by Eq. (9): rffiffiffi 5 pffiffiffiffiffiffiffiffiffiffiffiffiffinffi rel F dn ¼ 2 b Sn m vn , ð9Þ 6 where m* is the equivalent mass, which can be calculated analogue to the equivalent radius, b a coefficient depending on the restitution coefficient e given by Eq. (10), Sn the normal stiffness given by Eq. (11) and vrel n the normal component of the relative velocity between the particles. lne

b ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2

ln e þ p2 pffiffiffiffiffiffiffiffiffiffi Sn ¼ 2 En Rn dn

ð10Þ Fig. 2. Coupling of CFD and DEM.

ð11Þ

The tangential force depends on the tangential overlap dt and the equivalent shear modulus G*: pffiffiffiffiffiffiffiffiffiffi ð12Þ F t ¼ 8 Gn Rn dn dt : The first terms of the tangential force can be grouped into the tangential stiffness St, analogue to the normal stiffness: pffiffiffiffiffiffiffiffiffiffi ð13Þ St ¼ 8 Gn Rn dn : For the damping force in the tangential direction the following equation was used: rffiffiffi 5 pffiffiffiffiffiffiffiffiffiffiffiffiffin rel b St m vt : ð14Þ F dt ¼ 2 6 being the tangential component of the relative velocity with vrel t between the particles. The tangential force is limited by the Coulomb friction ms  Fn where ms represents the static friction coefficient. If the net tangential force is smaller than the friction force, the particle would not move in the tangential direction. The rolling friction is accounted for by applying a torque to the contacting surfaces as a function of the normal force Fn and the coefficient of rolling friction mr:

ti ¼ mr F n Ri oi :

ð15Þ

A linear cohesion contact model (DEM Solutions, 2010b) was also included to account for the cohesiveness occurring between the particles and between the particles and the wall when external forces are applied on the sediment. The cohesion force is calculated with Eq. (16) as a function of contact area Acontact and cohesion energy density kcohesion and is added to the Hertz– Mindlin normal force. This force represents in a macroscopic way van der Waals forces, electrostatic forces and other interparticle forces occurring at a microscopic level through the constant kcohesion. F ¼ kcohesion Acontact

ð16Þ

This methodology alone considers the contact forces but not the hydrodynamic forces. Through coupling of CFD and DEM both, hydrodynamic and contact forces can be determined.

3. Simulation methodology 3.1. CFD-DEM coupling A DEM simulation has to be coupled with a CFD simulation to consider the hydrodynamic forces, which play a crucial role on the particle trajectories. Fluid flow and particle movement have to be computed in an alternative manner. A schematic representation of how the coupling between the CFD software FLUENT and the DEM software EDEM works can be seen in Fig. 2. First, the flow is determined with CFD. Afterwards, the velocity, pressure, density and viscosity of the flow at each position in the computational domain, are transferred to DEM in order to calculate the

hydrodynamic forces and include them in the balance. After the particle positions have been calculated, the fluid flow is computed again considering the influence of the particles as an additional sink term in the momentum equation. The coupling methodology used between FLUENT and EDEM is a Lagrangian formulation. This means that all particles are reduced to mass points without volume when the flow is calculated. The volume occupied by the particles is not considered in the continuity and momentum conservation equations. A sink term with the addition of the drag forces of the particles included in each cell accounts for the momentum exchange. This method is restricted to very low particle concentrations, which is a reasonable assumption in the whole centrifuge except for the region where the sediment is formed. The Euler coupling available in the coupling module would consider the volume of the particles in the calculation of the continuous phase. However, it is not compatible with the VOF model used to simulate the multiphase flow of liquid and air. In order to consider the effect of the particles volume on the fluid flow, the Ergun equation, Eq. (17), was implemented. This equation calculates the pressure drop Dp of a fluid flowing through a packed bed of length L as a function of the porosity e, the superficial velocity v and the averaged particle diameter Dp:

Dp L

¼ 150

ð1eÞ2 m v 3

e

D2p

þ 1,75

ð1eÞ r v2 : e 3 Dp

ð17Þ

The pressure drop caused by the sediment (discrete phase) is taken into account in the fluid flow (continuous phase) by a sink term in the momentum equation, Eq. (3). Thus, even if a Lagrangian coupling is set in the software, the volume of the particles forming the sediment is considered when calculating the fluid. The sink term of the momentum equation in the three dimensional directions is calculated as:   m 1 Si ¼  vi þ C 2 r9v9vi ð18Þ 2 pc where pc represents the permeability and C2 the resistance coefficient due to turbulent flow. These parameters can be calculated by the Ergun equation, Eq. (17), as a function of sediment porosity in each cell and the average particle diameter of the particle distribution forming the sediment: pc ¼

D2p e3 150 ð1eÞ2

ð19Þ

C2 ¼

3:5 ð1eÞ Dp e 3

ð20Þ

In order to corroborate this method to consider the influence of the sediment on the flow, some simulations and validating experiments of a sediment build-up in a pipe passed by a flow were performed successfully (Noerpel, 2009).

´ndez, H. Nirschl / Chemical Engineering Science 94 (2013) 7–19 X. Romanı´ Ferna

In the case of a centrifuge the sediment built at the wall rotates with the same angular velocity as the bowl. The first layer of particles is subjected to the non-slip condition at the wall. The next particles are attached to the particles of the first layer after the contacts had been calculated. Although a small relative velocity in tangential direction may occur between the fluid and the sediment, the main direction of the sediment impeding the flow is the axial one. Radial flow is only present in the feed accelerator, near the impinge point of the feed on the rotating pool and at the outlet weir. In the remaining domain the fluid flow does not have any radial component. Hence, only the sink term for the z-direction was implemented to calculate the permeability and resistance coefficients in this direction. To reduce calculation time, repercussion of the particles on the flow was only considered near the wall where the particles accumulate (radius40.04 m and porosityo0.9).

time steps in order to correctly capture the particle contact behaviour. The choice of the time steps for the DEM is of crucial importance in order to achieve stability in the calculation and to avoid very long simulations, but still represent the real particle behaviour. It is therefore necessary to analyse the forces acting on the particles to determine the different time scales. The Hertzian contact time represents the total time of contact in the Hertz theory of elastic collisions between two elements (Hertz, 1881). It can be calculated as: !1=5 mn2 t H ¼ 2,87 n n2 : ð22Þ R E vrel For a good numerical performance, there must be a minimum of 6 time steps for each contact, which leads to a maximum time step of Dt ¼tH/6. The particle sedimentation time is the time required for a particle to settle a distance equal to one particle diameter:

3.2. Boundary conditions and solution method t sed ¼ In the CFD simulation, the liquid enters the centrifuge through the inlet at a given velocity, flowing then into an accelerator. The value of the inlet velocity defined in this case corresponds to a throughput of 5.8 l/min. The accelerator, which consists of two rotating plates, and the rotor are defined as rotating no-slip walls with a certain rotational speed. The fluid leaves the centrifuge via boreholes (outlet in Fig. 3). Here atmospheric pressure is imposed as boundary condition. The drag force of the fluid, liquid or air, on the particles is calculated as a free stream drag force (Lord Rayleigh, 1913): ! ! ! ! ! F freestream ¼ 0:5 C D Ap ð v f  v p Þ9 v f  v p 9,

ð21Þ

where CD represents the drag coefficient for a spherical particle ¨ (Loffler and Raasch, 1992; Schlichting, 1960), Ap the particle projection area, and vf and vp the fluid and particle velocity respectively. During the coupled simulation, the CFD calculation must be unsteady. The DEM time steps are typically smaller than the CFD

11

18m : Dp o2 rðrp rl Þ

ð23Þ

The time for a particle to adapt to the surrounding fluid flow is the particle relaxation time. It is calculated as: 2

t relax ¼

1 Dp rp : 18 m

ð24Þ

Here, the Hertzian contact time is the smallest time scale and a value of DtDEM ¼10  6 s is chosen to capture all physical phenomena occurring in the process. A DEM:CFD time step ratio between 1:10 and 1:100 is recommended for a stable simulation (DEM Solutions, 2010a). Consequently, a time step DtCFD ¼ 10  4 s is set in FLUENT. A parameter setting is necessary to represent the real particle behaviour by a contact model. In this study the values of the physical parameters of the rough glass particles and the steel wall listed in Table 1 were taken from the EDEM User’s Guide (DEM Solutions, 2010a). The friction coefficient is an empirically determinable property which can be found in the literature for various material combinations. Some authors have already determined this friction coefficient for sand and glass particles for DEM simulations (Asaf et al., 2007; Li et al., 2005). In this work various values of the static friction coefficient were tested to analyse the influence of this parameter in the sediment behaviour. One of the values tested was obtained from shear experiments using a sediment made of quartz particles in a Jenike shear cell under normal pressure. This methodology has been already used by other ¨ authors for the calculation of the friction coefficient (Hartl and Ooi, 2008). Two different particle sizes are tested: 125 and 200 mm. For either case a logarithmic normal distribution with a standard deviation of 2.5% is defined. In order to reproduce the experimental concentration of approximately 0.4% v/v, 80,000–100,000 particles per second are given at the inlet pipe of the centrifuge. Table 1 Parameters used in the contact model for the DEM simulation.

Fig. 3. Geometry of the cylindrical centrifuge used for the simulation of the sediment.

Particle radius

R

125/200 mm

Particle density Shear modulus of the particles Shear modulus of the wall Poisson’s ratio of the particles Poisson’s ratio of the wall Restitution coefficient Static friction coefficient Rolling friction coefficient

rp Gp Gw up uw e ms mr

1016/2550 kg/m3 2.2  108 Pa 7.0  1010 Pa 0.25 0.30 0.010 0.1/0.787/1 0.100

´ndez, H. Nirschl / Chemical Engineering Science 94 (2013) 7–19 X. Romanı´ Ferna

12

eight boreholes distributed along the circumference at the top of the bowl. Fig. 5 shows the results of the phase distribution of air and water in the centrifuge. Depending on the operating conditions, the radius of the interface may change in a small radial range corresponding to the size of the weir. For high rotational speeds the interface tends to be a perpendicular and the amount of fluid in the centrifuge thus decreases. In Fig. 5 the radius of the interface changes by 35% along the height of the centrifuge from the inlet to the upper part of the centrifuge because the C-value is 174 and the influence of gravity is notable. The main flow occurs in the direction of the rotational speed of the bowl. Both the liquid pool and the air core rotate. The tangential velocity of the water jet coming from the inlet accelerator (z* ¼0.275) is lower than the tangential velocity of the bowl because the liquid has to travel through the air core and loses tangential velocity. This causes the entire rotating liquid pool to decelerate compared to the angular velocity of the bowl. The values of the tangential velocity of the liquid averaged over the circumference as a function of the radius are plotted for different heights in Fig. 6.

Fig. 4. Grid on a horizontal cross section of the centrifuge.

3.3. Geometry and mesh The geometry of the centrifuge where the sediment was simulated is shown in Fig. 3. Due to the simple cylindrical geometry, a structured grid with approximately 100,000 cells was created. This grid resolution was found to be a good compromise between accuracy and computational demand. In this case, the grid was refined at the bowl wall (Fig. 4) to better determine the porosity of the sediment and its influence on the flow, although the particle volume was much smaller than the cell volume (0.3–1.3%). The computations were performed for a system having the dimensions shown in Fig. 3. In order to facilitate the transfer of the numerical results to systems of other dimensions, the following non-dimensional parameters were used: axial position (z*¼z/H), radial position (r*¼r/R), radius of the feed accelerator (r*feed ¼rfeed_accelerator/rweir) and centrifugal acceleration ratio or C-value, which represents the ratio between centrifugal acceleration (o2  r) and acceleration due to gravity g: C¼

o2 Ur g

:

Fig. 5. Contours of volume fraction of air obtained by the fluid simulation in the centrifugal geometry with r* feed¼ 0.75 at 2550 rpm and 5.8 l/min flow rate.

ð25Þ

4. Results and discussion 4.1. Fluid flow Prior to discussing the particle simulation results, it is necessary to analyse the multiphase fluid flow of water and air inside the centrifuge. Water is fed in axially through the inlet at a certain velocity and passed to the accelerator consisting of two rotating discs. There, the water changes its direction and gains tangential velocity. The rotational movement is transferred just at the plate surfaces, defined as no-slip boundaries, via momentum diffusion. Then, the water exits the accelerator in the form of continuous liquid streams that travel through the air core and enter the rotating liquid pool. Due to the Coriolis force of the rotating system, the water streams describe a spiral path until they reach the rotating liquid pool. The water leaves the centrifuge through

Fig. 6. Liquid tangential velocity versus radial position for different heights at 2550 rpm and 5.8 l/min.

´ndez, H. Nirschl / Chemical Engineering Science 94 (2013) 7–19 X. Romanı´ Ferna

The average tangential velocity of the liquid in this geometry is approximately 29.6% below the velocity of a solid body rotation. The values at the wall have to reach the solid body rotation because of the no-slip boundary condition. Hence, these values are closer to the solid body rotation than the others. The steep decrease of the tangential velocity near the wall was found to be smoother for a 6 times finer grid. However, the computational effort for the multiphase simulation including the particles was unaffordable with the finer grid. The small error bars, which represent the standard deviation from the circumferential average, indicate that the values of the tangential velocity do not change in circumferential direction. The deviation from rigid body rotation decreases for lower flow rates because less liquid has to be accelerated. The simulation of this geometry shows how important it is to accelerate the feed properly to the angular velocity of the liquid pool before entering it. Only then can the angular velocity of the bowl be reached in the rotating liquid, thus achieving the highest centrifugal acceleration possible. Although the main flow occurs tangentially, a secondary flow in axial direction arises due to the throughput. As expected from the boundary layer theory developed by Bass for the axial flow in solid bowl centrifuges (Bass, 1959,1962), a layer of fast moving fluid with high axial velocity at the gas–liquid interface was observed in the results (Fig. 7). Near the impact height of the feed into the rotating pool, two whirls with opposed directions are formed. Above the height of the inlet, a homogeneous axial boundary layer is formed at the interface. The error bars represent the deviation between the different circumferential positions, which is about 10% near the interface and near the outlet weir, the positions with

13

the highest variation. The layer thickness and the axial velocity values change slightly with the axial position. A thickness of 9.5 mm was calculated averaging the thickness over the centrifuge height. This value is similar to the one proposed by the theoretical ¨ model of Gosele (1968) of 11.3 mm, but much thicker than the values of one 1.7 mm calculated following the theory of Reuter (1967). Along the boundary layer, a recirculation layer is formed due to the lower tangential velocities, which cause a small pressure gradient in radial direction reducing the stratification of the flow and allowing the fluid to move in axial direction. As the angular velocity of the rotating pool is lower than that of the bowl, layers are formed at the bowl wall (Greenspan and Howard, 1963). 4.2. Particles and sediment behaviour The particles at the inlet pipe have the same axial velocity as the fluid. They then enter the feed accelerator and gain radial and tangential velocity. Subsequently, they pass the air core until they enter the liquid pool. There, they are subjected to the centrifugal force of the rotating liquid pool until they reach the wall or already settled particles. The results of the EDEM simulation of 200 mm glass particles coupled with the multiphase flow of air and liquid calculated in FLUENT at 2550 rpm and 5.8 l/min can be seen in Figs. 8 and 9. 4.2.1. Influence of the centrifugal force The sediment was simulated at two different bowl rotational speeds: 1275 and 2550 rpm. Consequently, two different velocity and pressure gradients in the bowl and two different sediment shapes were obtained. The effective rotational speeds of the liquid and the normal pressures acting on the bowl wall are listed in

Fig. 9. Particles (200 mm) coloured by velocity (see scale Fig. 8) in a vertical slice through the feed accelerator and the sediment at the rotor wall obtained with a friction coefficient of 0.787 at 2550 rpm and 5.8 l/min.

Table 2 Rotational speed and normal pressure achieved in the coupled simulations.

Fig. 7. Liquid axial velocity versus radial position for different heights at 2550 rpm and 5.8 l/min.

nbowl (rpm)

nliquid (rpm)

C-value

pnormal (kPa)

1275 2550

1080 1870

58 174

12.9 42.9

Fig. 8. Particles (200 mm) and path lines simulated with EDEM coupled with a two phase simulation of a centrifuge in FLUENT obtained with a friction coefficient of 0.787 at 2550 rpm and 5.8 l/min.

14

´ndez, H. Nirschl / Chemical Engineering Science 94 (2013) 7–19 X. Romanı´ Ferna

Fig. 10. Sediment formation behaviour of 200 mm particles coloured by velocity (scale similar to Fig. 8 but with values divided by 2) with a friction coefficient of 0.787 at 1275 rpm and 5.8 l/min.

Fig. 11. Sediment formation behaviour of 200 mm particles coloured by velocity (see scale Fig. 8) with a friction coefficient of 0.787 at 2550 rpm and 5.8 l/min.

Fig. 12. Solids volume fraction versus radial position for different numbers of particles obtained with 200 mm particles at 1275 rpm and 5.8 l/min with a static friction coefficient of 0.787.

Table 2 for both cases. Due to the slip of the rotating pool, the effective rotational speed of the liquid is clearly lower than the rotational speed of the bowl. In both cases, the evolution of the sediment can be observed in Figs. 10 and 11 as a function of the number of particles calculated. Both figures represent a vertical cross section of the centrifuge with the particles in the region near the feed accelerator where the sediment is built. Due to the low tangential velocity and hence low particle settling velocity, the particle trajectory is influenced by gravity when the particles leave the feed accelerator in the simulation at 1275 rpm. Then, thanks to the drag force of the flow in the rotating liquid pool, the particles travel in radial direction until they settle on the wall or on already settled particles. Meanwhile, they are affected slightly by the axial flow. The low normal pressure acting on the settled particles yields a sediment having the form of a heap. In contrast to this, the sediment obtained at the higher rotational speed of 2550 rpm, and hence, at the higher normal pressure acting on the bowl wall is flatter and the particles are not affected by the axial flow present in the liquid pool.

The solids volume fraction of the two different sediments is presented in Figs. 12 and 13. The values are averaged over the circumferential position and the error bars indicate the standard deviation. In both cases the volume fraction is close to zero outside the sediment and the concentration increases where the sediment is located. At the rotational speed of 1275 rpm, the sediment is thicker and grows in radial direction as the solids volume fraction increases with time at smaller radius. This is due to the lower normal pressure acting on the sediment. In the case of 2550 rpm, where angular velocity is higher in the rotating pool, the sediment is flat from the beginning of its formation. After a certain number of particles have settled (approximately 32,000), the solids volume fraction remains constant in radial direction with time. The sediment grows in axial direction. The highest solids volume fractions reached, about 0.54 in both cases, are lower than the dense random packing of mono-dispersed spheres of 0.63, but agree with the value of about 0.56 achieved for very loose random packing produced by settling spheres (Dullien, 1992).

´ndez, H. Nirschl / Chemical Engineering Science 94 (2013) 7–19 X. Romanı´ Ferna

15

Fig. 13. Solids volume fraction versus radial position for different numbers of particles obtained with 200 mm particles at 2550 rpm and 5.8 l/min with a static friction coefficient of 0.787.

Fig. 14. Sediment formation behaviour of 125 mm particles coloured by velocity (see scale Fig. 8) at 2550 rpm and 5.8 l/min.

Fig. 15. Solids volume fraction of 125 mm particles versus radial position for different heights of the sediment at 2550 rpm and 5.8 l/min.

A simulation made under changed flow rate of 2.5 l/min led to similar results. The influence of the flow rate is reflected by the effective angular velocity achieved of the rotating liquid pool. With higher flow rates, the deviation from the rigid body rotation is larger due to the larger amount of liquid to be accelerated, thus leading to a lower centrifugal force acting on the particles.

4.2.2. Influence of the particle size The particle size influences the time step used for the calculation to capture all the physical phenomena occurring in the

process. Smaller particle sizes require smaller time steps and, hence, longer simulation times. The coupled CFD and DEM simulations were performed for two different average particle diameters of 200 and 125 mm with a standard deviation of 2.5% in both cases. The particle trajectories and depositions are very similar for 200 mm and 125 mm particles under the same operating conditions. The difference lies in the size of the sediment. It is smaller when simulating the 125 mm particles, as it can be seen in Fig. 14. Fig. 15 represents the solids volume fraction for various heights of the sediment after the calculation of 100,000 particles. There, it can be seen that a similar maximal solids volume fraction

16

´ndez, H. Nirschl / Chemical Engineering Science 94 (2013) 7–19 X. Romanı´ Ferna

Fig. 16. Sediment build-up of 200 mm particles coloured by velocity (see scale Fig. 8) with a friction coefficient of 0.1 at 2550 rpm and 5.8 l/min.

Fig. 17. Sediment build-up of 200 mm particles coloured by velocity (see scale Fig. 8) with a friction coefficient of 1at 2550 rpm and 5.8 l/min.

Fig. 18. Solids volume fraction of 200 mm particles versus radial position at different times with a friction coefficient of 0.1, at 2550 rpm and 5.8 l/min.

Fig. 19. Solids volume fraction of 200 mm particles versus radial position at different times with a friction coefficient of 1, at 2550 rpm and 5.8 l/min.

of 0.55 is reached in the middle of the sediment. The solids volume concentration and, hence, the porosity do not strongly depend on the particle size in the range of 100–1000 mm, as was explained in literature (Yang et al., 2000).

2007; Li et al., 2005). The values of 0.1, 0.787 and 1 were chosen for the parameter study. The value of 0.787 was measured for a quartz suspension under normal pressure using a Jenike shear cell ¨ as it was done by other authors (Hartl and Ooi, 2008). The sediments obtained were very similar after 80,000 particles had been calculated. However, the evolution of the sediment as a function of time differs depending on the friction coefficient. Figs. 16 and 17 show the particles coming from the inlet accelerator and settling at the bowl wall for two different friction coefficients. Figs. 18 and 19 represent the solids volume fraction in the middle of the sediment as a function of the radius. In the case of ms ¼0.1, the sediment is thicker and the solids volume concentration is low at the beginning of the formation. After approximately 50,000 particles have settled and the pressure over the sediment has increased due to the mass of these particles, the sediment flattens and spreads over the bowl wall. This can be explained by the lower value of the friction coefficient, which

4.2.3. Influence of the friction coefficient Simulations with three different static friction coefficient values were performed to study its influence on sediment formation and behaviour. The friction coefficient characterises the tangential force necessary to move a particle under a certain normal force. This coefficient is a macroscopic representation of the interparticle forces that lead to sticking or gliding in a certain medium. The friction coefficient can be determined empirically and it can be found in literature for various material combinations. Some authors have already determined this friction coefficient for sand and glass particles for DEM simulations (Asaf et al.,

´ndez, H. Nirschl / Chemical Engineering Science 94 (2013) 7–19 X. Romanı´ Ferna

means that a higher normal pressure is needed to make the particles move. In the case of ms ¼1, the sediment is flat from the very beginning and the solids concentration profile remains constant. This kind of sediment evolution was also observed with ms ¼0.787. 4.3. Influence of the sediment on the flow The location and quantity of particles are recorded as a variable in the CFD software, i.e. the volume fraction of the particles, which is used to calculate the pressure drop of the flow through the sediment in axial direction. The effect of the sediment on the flow pattern is shown in Figs. 20 and 21 representing

17

the tangential and axial velocities after calculation of 60,000 particles at 2550 rpm and 5.8 l/min. The tangential velocity values, averaged over the circumference, are lower than those of a rigid body rotation as seen in all the centrifuges without baffles and with this kind of inlet accelerator. In this case, the deviation is 27.2%, which leads to an effective rotational speed of the liquid of 1870 rpm. A simulation of the same geometry under the same conditions but without particles led to a slightly higher deviation of 29.6%. The effect of the sediment is noticeable at heights of z* ¼0.125 and z* ¼0.250 in the region near the wall. The sediment of approximately 1 cm in thickness influences the layers of liquid rotating above; its angular velocity is increased. The axial velocity values shown in Fig. 21 correspond to a certain angular position, 2701, and could not be averaged over the circumference. This is due to the jets coming from the inlet accelerator which impact at different angles (at 2701 and z* ¼0.250) and to the outlet flow, predominating near the outlet weir (at 2701 and z*¼ 0.975). Both caused large error bars during averaging. The form of the axial flow pattern changes only qualitatively compared to the simulation without particles. At the wall and approximately 1 cm inside in radial direction, where the sediment is located, hardly any axial fluid flow occurs because of the resistance of the sediment. Upsides, the axial velocities obtained are larger in the recirculation layer. Thus, the sediment influences the recirculation layer at the bowl wall impeding it through the sediment but forcing it to increase externally.

5. Conclusions

Fig. 20. Tangential velocity values averaged over the circumferential position as a function of the radius obtained with the FLUENT-EDEM coupling using a friction coefficient of 0.787 at 2550 rpm and 5.8 l/min.

Fig. 21. Axial velocity values at the angular position of 2701 as a function of the radius obtained with the FLUENT-EDEM coupling using a friction coefficient of 0.787 at 2550 rpm and 5.8 l/min.

For the first time, a numerical simulation of the multiphase flow of gas, liquid and particles was performed in a cylindrical solid bowl centrifuge using a combination of CFD and DEM. The flow of gas and liquid was calculated with the multiphase model Volume of Fluid model (VOF) in order to obtain a sharp gas–liquid interface in the FLUENT software. This software was coupled with the EDEM software used to simulate particles of different sizes taking the hydrodynamic and contact forces into account. The volume occupied by the particles is not considered in the calculation of the fluid, but the effect of the particles is taken into account by a sink term in the momentum equation. The results reveal a main flow of water in tangential direction inside the centrifuge bowl. A secondary axial flow appears as a thin boundary layer at the gas–liquid interface rapidly moving towards the weir with a recirculation layer alongside. Particles travel through the air core and enter the liquid pool. Then, they move on spiral paths until they reach the wall or already settled particles and form a sediment. The particle trajectories and the sediment build-up agree qualitatively with the expected results. This proves that the coupling method can be applied to rotating geometries, such as centrifuges. When maintaining the contact model parameters and particle size constant, different sediment behaviour can be observed for different angular velocities in the rotating liquid pool. In the case of low tangential velocity, the particles form a sediment heap due to the low normal pressure exerted on it. In contrast to this, the sediment built at higher angular velocity and normal pressure acting on the bowl wall is much flatter and extends along the wall rather than in radial direction. The particle diameter was reduced from 200 mm to 125 mm, and no significant changes in the sediment behaviour were noted. However, some effort is still required to reduce the particle size even further without significantly increasing computation time. The value of the friction coefficient influences the behaviour of the sediment. Low friction coefficients lead to a sediment heap at

´ndez, H. Nirschl / Chemical Engineering Science 94 (2013) 7–19 X. Romanı´ Ferna

18

the beginning of the formation. After a certain amount of particles has settled; however, the normal pressure acting on the basis of the sediment increases and it spreads along the bowl wall. High values of the friction coefficient lead to a flat sediment with high solids concentrations from the very beginning of formation. Nevertheless, further studies in this field, including the experimental determination of the friction coefficient for different materials and particle size distributions, are needed. Moreover, experimental validation of the simulations using the same particle system is strongly recommended.

rel t w

relative tangential direction wall

Acknowledgements The authors owe special thanks to the AIF as part of the German Ministry of Economy BMWI for financial support.

References Nomenclature Acontact Ap C C2 CD Dp e E F g G H I kcohesion L m n p pc r r* R S t v z z*

2

area, m projection area, m2 centrifugal acceleration ratio or C-value, dimensionless resistant coefficient, dimensionless drag coefficient, dimensionless particle diameter, m restitution coefficient, dimensionless Young’s modulus, Pa force, N gravity acceleration, m s  2 shear modulus, Pa centrifuge height, m momentum of inertia, kg m2 cohesion energy density, J m  3 length, m mass, kg rotational speed, min  1 pressure, Pa permeability, m2 radial position, m non-dimensional radial position, dimensionless centrifuge radius, m stiffness, Pa m time, s velocity, m s  1 height, m non-dimensional height, dimensionless

Greek letters

a b Dp

e m mr ms u

r t tt ti $

volume fraction, dimensionless contact model coefficient, dimensionless pressure difference, Pa porosity, dimensionless dynamic viscosity, kg m  1 s  1 rolling friction coefficient, dimensionless static friction coefficient, dimensionless Poisson’s ratio, dimensionless density, kg m  3 shear stress tensor, Pa turbulent stress tensor, Pa torque, N m angular velocity, s  1

Index i, j l n p q

discrete elements liquid normal direction particle fluid phase

ANSYS, 2009. Modeling Discrete Phase, User’s Guide Ansys Fluent 12.0, Chapter 23, p. 1207. Asaf, Z., Rubinstein, D., Shmulevich, I., 2007. Determination of discrete element model parameters required for soil tillage. Soil Tillage Res. 92, 227–242. ¨ Bass, E., 1959. Stromungen im Fliehkraftfeld I. Periodica Polytechics chem. Ingenieurwes 3, 321–340. ¨ ¨ ¨ Bass, E., 1962. Stromungsund Absetzvorgange in Rohrenzentrifugen 3, 249–264. Chu, K.W., Wang, B., Yu, A.B., Vince, A., 2009. CFD-DEM modelling of multiphase flow in dense medium cyclones. Powder Technol. 193, 235–247. Chu, K.W., Yu, A.B., 2008. Numerical simulation of complex particle-fluid flows. Powder Technol. 179, 104–114. Cundall, P.A., Strack, O.D.L., 1979. Discrete numerical model for granular assemblies. Geotechnique 29, 47–65. Deen, N.G., Van Sint Annaland, M., van der Hoef, M.A., Kuipers, J.A.M., 2007. Reviewof discrete particle modeling of fluidized beds. Chemi. Eng. Sci. 62, 28–44. DEM Solutions, L., 2010a. EDEM-CFD Coupling for FLUENT User Guide 2.3. DEM Solutions, L., 2010b. EDEM 2.3 User Guide-Appendix A: Contact model theory. Di Renzo, A., Di Maio, F.P., 2005. An improved integral non-linear model for the contact of particles in distinct element simulations. Chem. Eng. Sci. 60, 1303–1312. Dong, J.K., Yang, R.Y., Zou, R.P., Yu, A.B., Roach, G., Jamieson, E., 2003. Simulation of the cake formation and growth in sedimentation and filtration, Third International Conference on CFD in the Minerals and Process Industries CSIRO, Melbourne, Australia. Dullien, F.A.L., 1992. Porous Media. Fluid Transport and Pore Structure, 2nd ed. Academic Press Inc., San Diego. Gondret, P., Lance, M., Petit, L., 2002. Bouncing motion of spherical particles in fluids. Phys. Fluids 14, 643–652. ¨ ¨ ¨ Gosele, W., 1968. Schichtstromung in Rohrenzentrifugen. Chemie IngenieurTechnik 40, 657–659. Greenspan, H.P., Howard, L.N., 1963. On a time-dependent motion of a rotating fluid. J. Fluid Mech. 17, 385–404. ¨ Hartl, J., Ooi, J.Y., 2008. Experiments and simulations of direct shear tests: porosity, contact friction and bulk friction. Granular Mat. 10, 263–271. ¨ ¨ ¨ ¨ die reine Hertz, H., 1881. Uber die Beruhrung fester elastischer Korper. Journal fur und angewandte Mathematik 171, 156–171. Hirt, C.W., Nichols, B.D., 1981. Volume of Fluid (VOF) method for the dynamics of free boundaries. J. Comput. Phys. 39, 201–225. Legendre, D., Zenit, R., Daniel, C., Guiraud, P., 2006. A note on the modelling of the bouncing of spherical drops or solid spheres on a wall in viscous fluid. Chem. Eng. Sci. 61, 3543–3549. Li, Y.J., Xu, Y., Thornton, C., 2005. A comparison of discrete element simulations and experiments for ‘sandpiles’ composed of spherical particles. Powder Technol. 160, 219–228. ¨ Loffler, F., Raasch, F., 1992. Grundlagen der Mechanischen Verfahrenstechnik. Vieweg, Braunschweig. Lord Rayleigh, 1913. The resistance of spheres in the movement of air. Comptes rendus hebdomadaires des seances de L’Acade´mie des Sciences 159, 109–110. Mindlin, R.D., 1949. Compliance of elastic bodies in contact. J. Appl. Mech.-Trans. Asme 16, 259–268. ¨ Muller, T., 2009. Modellbildung und Simulation der Partikelabscheidung an Fasern ¨ Mechanische Verfahrenstechnik. Universitat ¨ mittels CFD und DEM, Institut fur Stuttgart. Ni, L.A., Yu, A.B., Lu, G.Q., Howes, T., 2006. Simulation of the cake formation and growth in cake filtration. Miner. Eng. 19, 1084–1097. Noerpel, S., 2009. Numerische Modellierung des Sedimentaufbaus in Computa¨ Karlsruhe. tional Fluid Dynamics. Diploma Thesis, Universitat ¨ ¨ Reuter, H., 1967. Stromungen und Sedimentation in der Uberlaufzentrifuge. Chemie-Ing. Tech. 39, 311–318. Romani Fernandez, X., Nirschl, H., 2009. Multiphase CFD simulation of a solid bowl centrifuge. Chem. Eng. Technol. 32, 719–725. Romani Fernandez, X., Nirschl, H., 2010. A numerical study of the impact of radial baffles in solid bowl centrifuges using computational fluid dynamics. Phys. Sep. Sci. Eng., 10. Romani Fernandez, X., Spelter, L.E., Nirschl, H., 2011. Computational Fluid Dynamics and Discrete Element Method in centrifuges, Fluid Dynamics/ Book 2. InTech Open Acces Publisher. ¨ Schilling, M., Schutz, S., Piesche, M., 2009. Numerical simulation of the transport and deposition behaviour of particles on filter fibres using Euler–Lagrange and

´ndez, H. Nirschl / Chemical Engineering Science 94 (2013) 7–19 X. Romanı´ Ferna

coupling of CFD and DEM. In: Proceedings of the 6th International Symposium on Multiphase Flow, Heat Mass Transfer and Energy Conversion, Xi’an, China. Schlichting, H., 1960. Boundary Layer Theory, 4th ed. McGraw-Hill, New York. ¨ Schutz, S., Piesche, M., Gorbach, G., Schilling, M., Seyfert, C., Kopf, P., Deuschle, T., Sautter, N., Popp, E., Warth, T., 2007. CFD in der mechanischen Trenntechnik. Chemie Ingenieur Technik 79, 1777–1796. Stevens, A.B., Hrenya, C.M., 2005. Comparison of soft-sphere models to measurements of collision properties during normal impacts. Powder Technol. 154, 99–109.

19

van Wachem, B.G.M., Almstedt, A.E., 2003. Methods for multiphase computational fluid dynamics. Chem. Eng. J. 96, 81–98. Yakhot, V., Orszag, S.A., 1986. Renormalisation-group analysis of turbulence. Phys. Rev. Lett. 57, 1722–1724. Yang, R.Y., Zou, R.P., Yu, A.B., 2000. Computer simulation of the packing of fine particles. Am. Phy. Soc. 62, 3900–3908. Zhang, J.P., Fan, L.S., Zhu, C., Pfeffer, R., Qi, D.W., 1999. Dynamic behaviour of collision of elastic spheres in viscous fluids. Powder Technol. 106, 98–109.