Using DEM in Particulate Flow Simulations

There is no doubt that the dynamics of solid particles is the center of interest in the mineral processing industry, including crushing, grinding, classification, mineral separation, and leaching, just to name a few. But most processes involve fluid as a carrier and/or media, making the study of fluid dynamics and coupling between fluid dynamics and particle motion essential part for such a particulate flow. Various modeling and coupling approaches capable of considering particle behaviors, fluid dynamics and coupling effects, have been actively pursued, researched and developed in recent years. By name, particulate flows usually include one or more continuous fluid phase, and one or more type of particles, or say generally, discrete phases. The discrete particles/bubbles/ droplets are often dispersed in a continuous phase, so a discrete phase is also called a disperse phase in continuum multiphase modeling. There may be strong interactions between discrete phases and continuous phase, and strong interactions among discrete phases for dense particulate flows. The coupling physics pose a huge challenge to researchers, since coupling physics between fluid dynamics and particle motion requires coupling numerical modeling approaches. There is no one-fits-all solution for all applications especially after considering limitations, accuracy, computational costs of various numerical models. As computer technology for hardware and software advances so rapidly, it also push scientific and engineering simulations to high standards of requirement with respects to accuracy, fidelity, efficiency. There is increasing research activity of using Discrete Element Method (DEM) in particulate flow simulations. Discrete Element Method (DEM) (Cundall & Strack, 1979; Landry et al., 2003; Walton, 1992) is a Lagrangian model and is well accepted nowadays to model solid particle behavior. In principle, the DEM is based on the concept that individual particles, each of which is usually assumed to be semi-rigid, are considered to be separate and are connected at boundaries by appropriate contact laws. DEM naturally captures characteristics of each particle, therefore further dynamics like breakage and wear can be modeled locally at the small scale. Using DEM to track dynamics of particles, although the computing cost is high, eliminates the need of modeling fluid dynamics of particle phase, therefore improves fidelity of simulations. Interactions among discrete phases can be addressed more accurately inside DEM, thank to that microscopic physics has been clearly understood and described at most times. 2


Introduction
There is no doubt that the dynamics of solid particles is the center of interest in the mineral processing industry, including crushing, grinding, classification, mineral separation, and leaching, just to name a few.But most processes involve fluid as a carrier and/or media, making the study of fluid dynamics and coupling between fluid dynamics and particle motion essential part for such a particulate flow.Various modeling and coupling approaches capable of considering particle behaviors, fluid dynamics and coupling effects, have been actively pursued, researched and developed in recent years.By name, particulate flows usually include one or more continuous fluid phase, and one or more type of particles, or say generally, discrete phases.The discrete particles/bubbles/ droplets are often dispersed in a continuous phase, so a discrete phase is also called a disperse phase in continuum multiphase modeling.There may be strong interactions between discrete phases and continuous phase, and strong interactions among discrete phases for dense particulate flows.The coupling physics pose a huge challenge to researchers, since coupling physics between fluid dynamics and particle motion requires coupling numerical modeling approaches.There is no one-fits-all solution for all applications especially after considering limitations, accuracy, computational costs of various numerical models.As computer technology for hardware and software advances so rapidly, it also push scientific and engineering simulations to high standards of requirement with respects to accuracy, fidelity, efficiency.There is increasing research activity of using Discrete Element Method (DEM) in particulate flow simulations.Discrete Element Method (DEM) (Cundall & Strack, 1979;Landry et al., 2003;Walton, 1992) is a Lagrangian model and is well accepted nowadays to model solid particle behavior.In principle, the DEM is based on the concept that individual particles, each of which is usually assumed to be semi-rigid, are considered to be separate and are connected at boundaries by appropriate contact laws.DEM naturally captures characteristics of each particle, therefore further dynamics like breakage and wear can be modeled locally at the small scale.Using DEM to track dynamics of particles, although the computing cost is high, eliminates the need of modeling fluid dynamics of particle phase, therefore improves fidelity of simulations.Interactions among discrete phases can be addressed more accurately inside DEM, thank to that microscopic physics has been clearly understood and described at most times.
For fluid dynamics modeling, the terminology Computational Fluid Dynamics (CFD) is well-known dedicated to it.A conventional CFD is usually based on continuum mechanics principles and control volume methods.After decades of intensive research, conventional CFD is a well-developed technology with a series of mature well-defined numerical and physical models for single phase, turbulence and multiphase flows.It is basically a grid-based Eulerian model and is usually computationally efficient, especially for single phase flows.In practice, using CFD is not as easy as expected, since most mineral processing applications involve complex geometry and free surface flows.Generating appropriate volume mesh for complex geometry is a challenge even with help of commercial programs.Free surface flow also adds cost and difficulty to simulations.In the commonly used VOF method, the entire possible physical domain, even the space that is to be occupied by fluid occasionally, has to be meshed, and interface capturing and reconstruction scheme have to be implemented (Gao et al., 2003).Among CFD approaches, the one named Eulerian-Eulerian model, or say, multi-fluid model, has been extensively studied, implemented in MFIX (Syamlal, 1998;Syamlal et al., 1993), CFX, and FLUENT, and applied to simulations of fluidized beds (Gera et al., 2004;Sun & Battaglia, 2006a).The formulation of this model is essentially based on the continuum fluid dynamics.It considers both the fluid phase and solid particles to be interpenetrating continua whose dynamics are governed by the Navier-Stokes equations (Goldhirsch, 2003;Huilin et al., 2003;Savage, 1998).The particle mixture can be divided into several disperse phases with different properties.Closure of the model requires formulation of constitutive equations for each phases and inter-phase momentum transfer models, where often the most difficulties are encountered and approximations are made (Jenkins & Savage, 1983;Srivastava & Sundaresan, 2003).From the difficulty of building continuum models for granular flows (Gao et al., 2006;Savage, 1998), people realized that many of the physics-based governing equations work well at small scale, but non-linear physics makes derivations of those equations at a larger scale based on simplifications and assumptions no longer valid.For particulate flows with a wide property distribution, modeling errors are easily accumulated and computation costs are largely amplified.The major drawback of the Eulerian-Eulerian approach is that it cannot capture essential characteristics of individual solid particles regarding size and shape, and thus cannot effectively identify influence of these characteristics on process performances.Contact of individual particles with structure is often the major source of wear and erosion.Sizeand shape-change processes, such as breakage and chemical reaction of individual particles, usually are core features of the mineral processing industry.However, the extensive Eulerian-Eulerian research laid solid foundation for coupling DEM with CFD.Lots of ideas and equations can be adopted in DEM-CFD full coupling.DEM-CFD assumes the high theoretical fidelity, since each phase is kept to have its natural properties.We treat fluid as continuous continuum, and particles as discrete entities.The basic concepts of interpenetrating phases for multiphase flows still hold, although we only need to model and compute fluid phase.Instead of modeling several disperse particle phases in kinetic theory (Savage, 1998), we derive particle motions directly from DEM, therefore improves modeling accuracy.Smooth Particle Hydrodynamics (SPH) has been used to simulate fluid dynamics for years (Monaghan, 1988;1994;Morris et al., 1997).In SPH, a fluid field is represented by particles, each of which is associated with a mass, density, velocity, viscosity, pressure and position.Particles are moved by averaging (smoothing) their interaction with spatial 30 Hydrodynamics -Optimizing Methods and Tools www.intechopen.comneighbors based on the theory of integral interpolants using kernel functions which can be differentiated without use of the grid.SPH, as a Lagrangian particle-based method, has its particular characteristics.It has some special advantages over conventional grid-based CFD.The most significant one is the meshfree feature.SPH does not require a pre-defined mesh to provide connections between particles when solving the governing equations.The SPH particles themselves are adaptive to geometry and free surface confinement.From a numerical implementation point of view, DEM-SPH, a Lagrangian-Lagrangian model, is the best incorporation for particulate flows because it can totally eliminate the need for a volume mesh.The meshfree feature is very attractive to the mineral processing industry, where geometry is complex and free surface flow is typical in many applications.However, SPH is not without problems.While Eulerian-Eulerian approaches, modeling naturally discrete particles as continua, have difficulty to give correct constitutive equations, similarly, SPH, modeling naturally continuous media as particles, compromises accuracy in some aspects.It resolves the dissipative term poorly in comparison with grid-based methods.SPH has a limited ability to deal with steep density gradient or other large property changes.Boundary conditions do not fit naturally in the particle approach, so they are difficult to implement in SPH.It is hard to capture fluid dynamics where complex boundary conditions are of critical importance.In summary, there is no single one-fits-all solution.Every model has strengths on some aspects and weaknesses on others, especially considering accuracy and cost factors.People have to be able to, based on preliminary understanding of the physical characteristics of a system of interest, pick up the right models, combine them together, and develop/use appropriate models for the specific system to capture major phenomena to discover and investigate the controlling mechanisms behind them.In this work we present three numerical coupling approaches to capture the physics of interest: one-way coupling with CFD, DEM-CFD coupling, DEM-SPH coupling.A one-way coupling is basically to run fluid dynamic solver separately from DEM, then import fluid flow solution to DEM, where the fluid effect on particles is considered.The one-way approach is practically important for industry applications, because at complex situations full coupling modelings are hard to converge, if not impossible.It has the advantages of using commercial package.This advantage may become very attractive in industry applications where the flow condition is very complex and density of particles is not very high.A one-way coupling can be extended to so-called 1.5-way coupling if multiphase fluid solver is used instead of single phase solver.We applied one-way coupling to a slurry pump, where solid particles and fluid are well mixed so that it is appropriate to treat slurry as a kind of single phase mixture, and the FLUENT is used to solve the flow field.But the CFD solver cannot give direct answers to our concerns: wear effect of particles on pump structure and particle breaking probability, therefore CFD results are imported into our DEM code to simulate the detailed behavior of individual particles.For the strong coupling physics, the full coupling of DEM-CFD or DEM-SPH is necessary.In the DEM-CFD coupling, we employ a lot of widely accepted Eulerian-Eulerian multi-fluid models that have been intensively studied in the continuum multiphase fluid dynamics.Convergence of the coupling models is usually a huge challenge.The numerical methods are discussed in each section of model descriptions.The segregation of different sizes of particles in a fluidization bed is controlled by both particles motion and fluid dynamics.Due to the simple geometry of a bed, DEM-CFD is the best candidate for this application.

31
Using DEM in Particulate Flow Simulations www.intechopen.com In the DEM-SPH coupling, a multiphase DEM-SPH model is proposed and described in detail.The coupling between the solid particles and the fluid phase are considered via volume fraction, pressure and drag force.In mills, particles are so dense that the behaviors of solid particles dominate the overall physics, and separation of solid particles from fluid is apparent.But fluid dynamics cannot be neglected because fluid damping/drag may change particle-particle contacts and particle-boundary contacts.Furthermore, the rotating nature, presence of a free surface flow, and liner geometry make CFD simulations of mills prohibitively expensive.Numerical scheme is simple at this stage, time stepping scheme could be either explicit Euler or lower order Runge-Kutta for more accuracy.The DEM-SPH is applied to a preliminary study of a mill.

DEM model with flow effect
The DEM simulation is based on a 3D soft particle model (Cundall & Strack, 1979;Silbert et al., 2001;Walton & Braun, 1986) where small deformations and multiple contacts on a particle are allowed, and friction and rotation are also taken into account.Each particle has six degrees of freedom of motion.For simplification of the description of the DEM model, spherical particles are assumed, although in our real application and implementation, we use spherical, tetrahedral and irregular convex shape particles together.The complex shapes definitely add difficulty to the computational geometric and solving for particle rotation, but the basic theory behind it is the same as for all spherical particles.Although our program is fully parallelized for a distributive memory machine, this work is not intended to cover the parallelization scheme and readers who are interested are referred to Plimpton (1995), which is the framework we follow.The movement of a particle with mass m, moment of inertia I can be described by Newton's law and the kinematic relation: where subscripts i and j are for identifying particles, u is the velocity of the mass center, ω the angular velocity, r the position, R the radius of the particle, F c,ij the contact force of particle j acting on i, F sf,i the fluid-solid interaction which is assumed to act at the mass center of particles, and n ij the normal direction of the contact pointing to particle i from j.
The governing equations are simply ordinary differential equations in time.Each particle is evolved by integrating the governing equations and applying the initial condition.The major task of modeling and simulation thus becomes actually formulating and calculating the force terms.
The implementation of contact forces is essentially a reduced version of that employed by Walton & Braun (1986), developed earlier by Cundall & Strack (1979).Contact force } is first calculated from the deformation through the spring-dashpot model,

32
Hydrodynamics -Optimizing Methods and Tools www.intechopen.comassuming soft particles: where m eff =( m i m j )/(m i + m j ), subscript n denotes the normal direction and t the tangential direction, k is the spring coefficient, λ the damping coefficient, and δ the elastic displacement for a contact, respectively.The damping effect in the tangential direction can be neglected.The tangential displacement δ t between particles is obtained by integrating surface relative velocities over time during deformation of the contact.Actually just this history dependent feature makes the computation, especially parallel computation, more expensive.The magnitude of δ t is truncated as necessary to satisfy a local Coulomb yield criterion |F t |≤µ|F n |.
In the above DEM model, the normal compression between two particles is easily written as: here, The normal direction and tangential direction are defined as: The value of the spring constant should be large enough to avoid particle interpenetration, but not too large to require an unreasonably small time step ∆t, since an accurate simulation typically requires ∆t ∼ t c /50, here t c is the characteristic contact time during a collision process between particles.The amount of energy lost in collisions is characterized by the inelasticity through the coefficient of restitution e.For the linear spring-dashpot model, the following relations can be taken as guidance to find the damping coefficient λ n e n = exp(−λ n t c /2) , ( 10) After contact force is calculated, the equations of motion, which are ordinary differential equations, can be numerically integrated to get the particle trajectories.The boundary surfaces are represented by triangles.Any meshing tool generating surface triangles can be used here.
In comparison with particle-particle contact, the only difference is the geometry resolution for the particle-triangle contact.The overlap δ n is equal to a particle radius minus the distance to a triangle.The advantage of DEM is that it can capture behaviors of individual particles, and collect detailed contact information, such as velocity, contact force, shear and impact energy spectrum, so that it can model wear more accurately.Most severe wear happens in particulate flows when solid particles collide on a surface.We model the particle wear effect on a boundary triangle as: where subscript i denotes the particle groups, ∆h is the thickness to be worn off the surface, A is the triangle area, E sh,i is the cumulative shear energy over the time period of interest, and C wr,i is the wear coefficient that is a function of the triangle and particle material properties and the particle size.Similarly, we model the particle breakage as where Br i is the breakage percentage (probability) of the group i of particles when passing through a process, E imp,i is the cumulative impact energy over the retention time of the process, and C br,i is the breakage coefficient that is a function of the material properties and the particle size.The constants C wr,i and C br,i are to be determined from the calibration with experiment or operational data in this work.See the reference for details and recent developments about the wear model (Hollow & Herbst, 2006;Qiu et al., 2001) and particle breakage model (Herbst & Potapov, 2004;Potapov et al., 2007).

CFD coupling with DEM
The basic conceptual theory for CFD coupling with DEM comes from the Eulerian-Eulerian multi-fluid model, where the fluid phase and the solid particle mixture are described as interpenetrating continua.The particle mixture can be divided into a discrete number of phases, each of which can have different physical properties.Generally, n sets of governing equations have to be solved for a multiphase flow with n phases, and an exponentially increasing number of constitutive equations are required for closure of the model.The approach of CFD coupling with DEM is proposed to overcome the closure difficulty.The Eulerian control-volume multiphase CFD governing equation is used to describe the fluid dynamics, while the DEM is used to model the solid mixture dynamic behaviors.
The governing equations for incompressible flow, continuity and momentum equations, for the fluid phase, are: where ρ is the fluid density, u f the fluid velocity, θ f the fluid volume fraction, and P the pressure, and F d the drag force, and T the viscous stress tensor.
The drag force F d between particles and fluid phase is generally defined as where β is the drag force coefficient, u s the velocity vector of solid particles, and (u s − u f ) is the slip velocity between the two phases.The volume fraction θ s and velocity fields u s of the solid phase are obtained through averaging particle data from DEM simulation, thus θ f can be obtained from Without losing generality, in the DEM-CFD coupling, the drag correlation by Syamlal et al.
(1993) is adopted: where dp is the average diameter of particles, and V r is the ratio of the falling velocity to the terminal velocity of a single particle.The following form for V r is used From the above governing equations and drag force calculation, we can see volume fractions appear allover everywhere.This is the core modeling concept for interpenetrating disperse flow, whereas the exact particle shape has not been followed, rather the volume/mass conservation must be maintained.This multiphase modeling approach is a perfect compromise for large scale particulate flows with a huge number of small (relative to flow characteristic length) particles.
The volume fraction and velocity of the solid phase are needed for each cell.It sounds simple at the first glance.It may be simple for regular mesh, not for irregular mesh and irregular particles.As we mentioned, particle shape is not of importance at modeling particle fluid interaction, nor should be mesh geometry.We compromise the geometry details for the efficiency of numerical calculation.Many researchers use kernel function to calculate the volume fraction to avoid handling complex mesh and particle geometry.The volume fraction of each cell is obtained from the particle spatial distribution and the volume of each particle through the averaging processes: where the subscript c denotes a cell and i denotes a particle.N p is the total number of particles in the system.K is a kernel function, which should be bell-shape function with local support, such as Gaussian function K(ξ)=exp[−(ξ/w) 2 ], where ξ =( |x i − x c |)/w and w is the bandwidth.In this study, the following kernel function is used: This function is very close to Gaussian function (with scaled bandwidth w = 0.45), and it is more efficient to compute than Gaussian function, because it only requires to loop over the particles in the neighbor cells.Recently, Xiao & Sun (2011) dedicated a big portion of their paper to the particle volume fraction calculation.Please refer to that for details.Similarly, the solid phase velocity is obtained via the following coarse-graining procedure:

35
Using DEM in Particulate Flow Simulations www.intechopen.com On the DEM side, two major solid-fluid interaction terms acts on a particle.In Equation 1force F sf,i equals the sum of buoyancy force and drag force acting on particle i: The drag of Equation 16is the aggregate drag force acting on all the particles in a cell.For the motion equation 1 of a solid particle, the drag force on each particle is needed.There are many ways of redistributing the aggregate force back to individual particles.This can be obtained from the cell whose center is nearest to the particle, or from the neighboring cells whose contribution portion is determined by distances between particles and cell centers.
Each particle shares the drag force proportional to its surface area, or proportional to its volume, or as functions of other properties.The uncertainty also occurs with buoyancy force.
There is no sure answer to the question of where pressure values are used in calculation of buoyancy force on a particle.Here in the application of DEM-CFD full coupling, we take the simple scheme by distributing the nearest cell force according to surface areas of particles.
For the one-way coupling applications, particles are small and well mixed with fluid phase, so that one-way coupling is justified, and Wen-Yu drag relation (Li & Kuipers, 2002;Rong & Horio, 1999) is used.Thus, the drag force calculation can be as simple as: where V i is the volume, and D i the hydraulic diameter of particle i.For the drag coefficient C d , please refer to DEM-SPH coupling section.
For DEM-CFD coupling, the fluid equations are solved using a solver provided by a library in the OpenFOAM (Open Field Operation and Manipulation) toolbox (Rusche, 2002).This solver is used in the current study with modifications to accommodate the fact that only the fluid phase is solved and the disperse phase is tracked in the Lagrangian framework.Finite volume method is used to discretize the equations on an unstructured mesh.For the time integration, Euler implicit scheme is used, which has only first order accuracy but is unconditionally stable.The convection and diffusion terms are discretized with a blender of central differencing (second-order accurate) and upwind differencing (first order accurate).The advantage of blended differencing is that high accuracy is achieved while the boundness of the solution is ensured.A sophisticate stepping control and interpolation over time was brought up by (Xiao & Sun, 2011) to enhance accuracy and convergence of the DEM-CFD coupling solver.The velocity-pressure coupling is handled with the modified PISO algorithm (Rhie & Chow, 1983).In this algorithm, momentum equation is first solved to get a predicted velocity field, and then the pressure equation (obtained by combining momentum equation into continuity equation) is solved for a corrected velocity field.This process is repeated until the velocity field satisfies continuity equation.The PISO algorithm prevents from the decoupling of pressure-velocity and the oscillation in the solution, eliminating the necessity of a stagger grid.Therefore, a collocated grid is used in the models, where all the variables are stored in the cell centers, thus it is a significant simplification over a stagger grid.

SPH coupling with DEM
SPH equations are obtained by interpolating fluid dynamics governing equations over disordered mass points in the influence range of an interpolation kernel function (Monaghan, 1988).The kernels are analytical functions which can be differentiated without using a mesh.Although control volume CFD can be tuned to get accurate solutions of physical problems, it requires tremendous work, including generating the mesh, in order to couple with DEM to account for the multiphase flows.On the other hand, SPH is a method which gives reasonable accuracy and couples well with the particle method DEM without requiring a mesh.Monaghan & Kocharyan (1995) originally built SPH multiphase models for interpenetrating multi-fluids.In the later work, Monaghan (1997) improved the multiphase SPH solver by using an implicit drag technique.In this work, we will modify the interpenetrating multiphase SPH model to couple SPH for fluid dynamics with DEM for solid particles.SPH is only used to model fluid phase, and particles are represented and evolved by DEM The governing equations, continuity and momentum equations, for the fluid phase are: where β is the drag force coefficient, u s the velocity vector of solid particles and T the viscous stress tensor.The fluid density ρ in the multiphase model is related to the fluid volume fraction θ and actual fluid density ρ by The equation of state has to be defined in order to fully describe the dynamics of the fluid.The actual equation of state of incompressible flow is very stiff, requiring extremely small time steps.In SPH, the fluid pressure is an explicit function of local fluid density.Therefore, it is necessary to use a quasi-compressible equation of state and an artificial speed of sound as a reference value.Monaghan (2000) used the equation of state similar to that defined for water: where ρ 0 is the reference density, c 0 the speed of sound at ρ 0 , and γ a constant with physical meaning of the ratio of specific heat for ideal gases.Monaghan (2000) took γ = 7 for incompressible flows.The choice γ = 7 results in large changes in pressure from small changes or perturbations in density.We confirmed that Equation 31 works well for single phase flow SPH.However, in this work considering multiphase free surface flow, the density changes could also be exaggerated by changes or errors in calculation of volume fraction (see below).We also follow suggestions of Morris et al. (1997) regarding choosing appropriate values of γ and c 0 .The speed of sound is:

37
Using DEM in Particulate Flow Simulations www.intechopen.comSPH is based on the theory of integral interpolants.If the kernel functions are some types of delta functions, then a field variable A(r) can be approximated by the weighted averaging over a limited range of neighboring particles as: where m a denotes mass of SPH particle a at the position r a , and similar notations for ρ a , A a .W(r − r a , h) is the kernel which is a function of smoothing length h and distance between positions r and r a .
For clarification of the description, the subscripts a and b are used for the SPH fluid particles, and i and j for the DEM solid particles.Integrating the governing equations and simplifying the integrals as above, we can get the overall SPH governing equations: where V j is the volume of DEM particle j, θ j the solid phase volume fraction at the position of DEM particle j, θ a the fluid phase volume fraction at the position of SPH particle a, and τ ab the viscous term.The solid-fluid inter-phase interaction force is represented by the second and third group of terms on the right hand side of Equation 36.The second term is the pressure gradient on a solid particle, and third term is the drag force.The drag force between two particles acts along their center line, working like dashpot damping in the DEM models.Thus the inter-phase interaction in Equation 1 can be written as: In the above equations, we have used the notation for a vector, and for a kernel and a kernel derivative, where r ab is the distance between particle a and b.Similar notation is used for other terms.
The kernel function is the commonly used cubic spline function where q = r/h and r is the distance between particles.The Gidaspow drag correlation, which combines the Wen-Yu relation and the Ergun equation, is commonly used in CFD multiphase modeling (Gera et al., 1998;2004;Li & Kuipers, 2002;Rong & Horio, 1999) and it is used here: where D j is solid particle effective hydraulic diameter.The ρ a , θ a and θ j terms are intentionally grouped together as factors, because these terms will be canceled when β aj is substituted back to Equations 36 and 37.
From the comparison of governing equations of DEM-CFD to those of DEM-SPH, one apparent difference is that DEM-CFD calculates aggregate coupling force first, and distributes back to particles, while the DEM-SPH calculates coupling force for individual particles and collects the total force on SPH particle.Due to the high cost the drag force calculation, the cost of DEM-SPH may be much higher than DEM-CFD.Further modification to the DEM-SPH by adopting DEM-CFD approach is the future work of authors.
No matter which way of coupling, we have to obtain the collective volume fractions θ j and θ a .Using kernel function for volume fraction calculation is the most natural way in the SPH framework.Please refer to Monaghan (1997); Monaghan & Kocharyan (1995) for more theory support.The most simple form of fluid fraction in the SPH type of kernel functions can be defined as follows: where W * is a kernel function that can be different from W, and a different smoothing length should be used.The smoothing length should be at least twice the maximum particle size.
After that, θ j can be obtained by smoothing over SPH particles and using the same smoothing length as in SPH equations.One can use the same way as in the above DEM-CFD, Equation 23, for volume fraction calculation, however, using a imaginary sphere as V c , because there is no real mesh in SPH.
The imaginary sphere radius should be at least twice the maximum particle size.In this work,

39
Using DEM in Particulate Flow Simulations www.intechopen.comwe use a combination of both approaches.We first obtain the intermediate volume fraction in a way like Equation 23, then smooth it over SPH particles in the normal way as We have to go back to formulate the viscous term in order to complete the set of modeling equations.This work employs a SPH viscous diffusion model used by Morris et al. (1997): where µ is the dynamic viscosity.This expression uses only the first order kernel derivative.It conserves translation momentum accurately, while angular momentum is only approximately conserved.For applications with low fluid velocities like in this work, this formulation is appropriate.
In this work, several practical approaches suggested in a series of SPH publications by Monaghan (1989;1994;2000) are employed in order to build a stable, robust solver.First, the so-called XSPH velocity correction is employed: The XSPH variant is used to move the particles and also in the continuity equation for consistency.The adjustment is important for free surface flows and useful for high speed flows.It basically keeps the particles more ordered and moves particles in a velocity similar to the average velocity in their neighborhood to prevent fluid penetration.Second, artificial pressure is introduced into the momentum equation to avoid SPH particle clustering.See Monaghan (2000) for details.The overall time advancing scheme for a DEM and SPH coupling system is: 1) calculate the solid and fluid coupling terms based on the old field values, 2) calculate the fluid-fluid interactions and integrating SPH particles.3) calculate the solid-solid contact forces and move solid particles.Two different time stepping schemes for SPH particle integration have been applied and tested in this work.One scheme is analogous to the explicit method of control volume CFD for incompressible flow, where the source term of pressure and force terms are calculated based on the old velocity field at t n , then the pressure equation is solved, and finally the new pressure field is used to update the velocity to t n+1 .Similarly in the SPH solver, at the first sweep, the density changes and force terms are calculated based on fields at t n ; at the second sweep, the pressure and density are updated; at the third sweep, the velocity fields are updated and the particles are moved.The method is enhanced by incorporation of the leap-frog approach: velocities are updated at intervals midway, t n+1/2 , between time steps t n and t n+1 .The other scheme is the simple predictor-corrector (Monaghan & Kocharyan, 1995), or, the second order Runge-Kutta method.First, values (velocity, density, position) at t n+1/2 are predicted from t n and t n−1/2 .Then, force and other changes under the predicted conditions are calculated.After that, field values at t n+1/2 are corrected using the new changes.Finally, 40 Hydrodynamics -Optimizing Methods and Tools www.intechopen.comvalues at t n+1 can be obtained using the trapezoidal rule.This method can achieve second order accuracy.Boundary conditions for SPH particles are similar to that for solid DEM particles: spring-dashpot contact models are applied to the particles assuming spherical particles with diameter equal to the initial separation distance.The contact coefficients are taken from DEM parameters for the same size of DEM particle.

One-way coupling for a pump
For slurry in pump operation, solid particles are well mixed with fluid, the particle volume fraction is usually not very high, velocity is high but streamlines do not interweave together, so that the particle-particle collision probability is low, thus the particle-boundary collision is the main interest here since particle-boundary interaction is the major source of wear and the reason for particle breakage under these conditions.DEM simulations were also performed to estimate the breakage percentage and liberation percentage of the particles during the operation of the pump at different conditions.The commercial software FLUENT is used to obtain the fluid dynamics solutions assuming single phase mixture flow.
The CFD simulation of a pump is not easy and straightforward like pipe flow.First of all, for high fidelity we use the manufacturer's geometry files in IGES format.These geometry files need a lot of cleaning, merging, and smoothing work in Gambit before appropriate CFD meshes can be generated.Moreover, due to the rotating part, the multiple reference frame technique has been used.Due to the unsymmetrical geometry, time dependent solutions have to be pursued, thus the sliding mesh method has been used.Fluid field solutions including pressure and velocity fields are imported into our DEM program, where DEM tracer particles are created and evolved as described above.During the simulations, particles are first set at the inlet segment of the pump, then particles flow along with fluid to enter the pump.Particles leaving the pump from the outlet are re-inserted at the inlet as in semi-periodic boundary condition, so that particle dynamics gradually reaches steady state.Once steady state is reached, shear and impact energy data are collected and recorded for at least three revolutions of the pump.
We performed the wear study on the Metso slurry pump HM300.Ultrasonic measurements of lining thickness as a function of location number is marked on Fig. 1.One hundred thousand (100,000) ore particles from size 0.5 mm to 3 mm are generated for DEM.Shear energy spectra are used in wear analysis as described by Equation 12.The prediction error was found to be less than 4%.The comparison of the wear prediction with measurements is shown in Fig. 2. We performed the ore breakage and valuable particle liberation study on Metso slurry pump MM300.We employed 5 different sizes (1.4 mm, 5.2 mm, 8.17 mm, 10.6 mm and 15 mm) of DEM particles to represent valuable particles, respectively.There are a total of 400,000 particles for ore, and 1000 particles for each of the 5 groups of valuable particles.Based on the impact energy spectra, the particle breakage percentage can be calculated according to Equation 13.The liberation percentage is proportional to be the breakage rate of the ores.The valuable particle breakage percentages and liberation percentages are estimated as shown in Table 1.The valuable particle liberation percentages are under 5% and the valuable particle breakage percentages are very small (< 0.1%) at the pump conditions under consideration.The trend of liberation percentages changing with flow rate and with solid concentration meet the qualitative expectation.As the volume flowrate increases, although the increased velocity enhances the collision of particles, but the retention time decreases, therefore, the liberation percentage decreases from the combination of two effects.This result may show us that the probability of collision, which increases with retention time, outweighs the collision strength increase with velocity at the current conditions.

DEM-CFD coupling for a fluidization bed
The DEM-CFD coupling has been applied to a study of particle segregation due to size.Here, we simulated and presented the segregation of two different sizes of particles in a model fluidization bed with inlet velocity of 1.5 m/s, which is larger than the minimum fluidization of the large particles reported in the experiments (Goldschmidt, 2001).To quantify the bed expansion, an average particle height is defined as here, m i is the mass of particle i.
To reduce computational cost while still keeping the essential physics of the process, a smaller bed as shown in Fig. 3, thus a smaller number of particles, is used, while the aspect ratio of the bed is kept the same as in the experiment (Goldschmidt, 2001).Two types of particles with the same density but different sizes were perfectly mixed initially.
The simulation was conducted for 10 seconds.The results in Fig. 4 show the close to zero segregation percentages, which agrees with the experimental results.The initial bed configuration as well as those at t = 6 s are presented in Fig. 5.The left panel shows that initially the particles are perfectly mixed and randomly packed.At t = 6 s, the two types of particles are still well mixed without visible segregation.The results are better than the multi-fluid modeling results (Sun & Battaglia, 2006b).

DEM-SPH coupling for a mill
We applied the DEM-SPH coupling multiphase solver described in section 2 and 3 to the modeling of a SAG mill.The mill operating conditions are: mill diameter 10 m; mill effective length 4 m; mill volume 386 m 3 ; total filling 31.7%;ball filling 15.4%; ore density 2700 kg/m 3 ; ball density 7850 kg/m 3 ; rotation speed 10.1 RPM; solid feed 1977 mt/h (mt = metric tons); solid mass percentage in feed 70%.We first ran the PBM model (Herbst & Pate, 2001) to estimate the size distribution under steady operation conditions.Due to the high cost of the particle methods, we use a truncated size distribution list that cuts the PBM size range to keep only the 5 coarsest sizes.We also assume that 25% solid fine particles are totally mixed with water to form a dilute slurry.The slurry has density 1380 kg/m 3 and dynamics viscosity 5.0E-3 kg/m s.Based on these practical treatments, we employ 6 groups of ore particles: 43 mm (in diameter) sphere, and 43 mm, 62 mm, 88 mm, 124 mm, 176 mm tetrahedron, and 4 groups of balls with diameter equal to 78 mm, 100 mm, 125 mm.At initialization, a total of 750,024 ore particles from the 6 ore groups according to the truncated size distribution are generated in the mill, and a total of 48,000 balls, 12,000 in each ball group, are set in the mill.Assuming 0.4 voidage of dense packing, we set up 868,158 SPH particles, with spacing distance 50 mm to make the fluid charge level about the same as that of the solid particles.Fig. 6 shows the initial setup of the particles in the mill.Due to the high density, all particles are initially packed regularly at the mill belly part.

43
Using DEM in Particulate Flow Simulations www.intechopen.comThe computation is expensive.The case has been running for 20 days using 16 processes but it only simulates 6 seconds of physical time, still not reaching steady state.We also realized the simulation is unstable due to explicit particle stepping, so time-step need to be controlled to be small.Here, we only can present some preliminary results.Fig. 7.A view with particle shapes at the mill center cross section.
Figure 7 shows the DEM and SPH particles at 5.68 seconds after the mill starts rotating from the initial setup.The light blue color represents SPH particles, red or brown denotes DEM ore particles, and black demotes steel balls.We can see the profiles of DEM and SPH particles spatial distribution that meet our expectation.In the figure DEM particles overshadow the SPH particles, so that it creates false impression that there are no SPH particles in most area.
The initial setup as shown in Fig. 6 may cause the segregation shown in Fig. 7. Please note at 5.68 seconds, the mill has finished less than one revolution.The potential power draw can be calculated as ∑ i,a g × r, where the sum is over all DEM and SPH particles.The change of this potential power with time is shown in Fig. 8, reflecting how close a mill is reaching the steady-state.The fluid power draw is less than 10% of the total power draw.It is understandable, at the first half revolution, potential power draw increases to its highest value because particles have to be lifted up by a mill.Unfortunately we do not have operational data to compare with.In another way, the power draw is calculated by ∑ F • u, where the sum is over all contact work acting on the boundary.At the time 5.68 seconds, power draw calculated in the second way is 9.18 MW.This value can be thought to be gross power draw, because it includes energy loss between inelastic contacts and viscous dissipation.The advantage of the DEM-SPH solver, a total particle method, can be shown in Fig. 9, where small solid and SPH particles pass through the grid.We can imagine the difficulty of meshing a mill with a grid for a conventional CFD solver, considering the multi-scale challenge here: the mill main body diameter is 10 m while grid size is 75 mm.But with SPH, it is flexible to control the solver by assigning SPH particle probability of passing through, or by applying different sets of triangles to SPH and DEM particles.

Conclusions
Three approaches to couple solid particle behavior with fluid dynamics have been described and three applications have been provided.For full coupling approaches DEM-CFD and DEM-SPH, they are physically equivalent, but may appear in different forms of equations.The governing equations have been carefully formulated.Numerical methods, difficulties and possible problems have been discussed in detail.The one-way coupling of CFD with DEM has been used in analysis of wear on lining structure and particle breaking probability during a pump operation.The DEM-CFD coupling has been applied to modeling fluidization bed.The multiphase DEM-SPH solver has been used in a wet grinding mill simulation.Each numerical approach has its strength and weakness with respect to modeling accuracy and computation cost.The final choice of best models should be made by application specialists on a case by case basis based on dominant features of physical phenomena and numerical models.
) 33 Using DEM in Particulate Flow Simulations www.intechopen.com

Fig. 5 .
Fig. 5. Two-dimensional snapshots of particle mixtures.Dark dots denote large particles and light ones denote small particles.Left: Initial configuration.Right: configuration after 6 seconds.

Fig. 6 .
Fig. 6.The initial setting of the particles

Fig. 8 .
Fig. 8.The calculated power draw variation with time

Fig. 9 .
Fig. 9.A view of the mill discharge end.

Table 1 .
Valuable particle liberation and breaking analysis in slurry pump.