Grain Sedimentation with SPH-DEM and its Validation

. Our mesoscale simulation method [M. Robinson, S. Luding, and M. Ramaioli, submitted (2013)] for multiphase ﬂuid-particle ﬂows couples Smoothed Particle Hydrodynamics (SPH) and the Discrete Element Method (DEM) and enjoys the ﬂexibility of meshless methods, such as being capable to handling free surface ﬂows or ﬂow around complex and/or moving geometries. We use this method to simulate three different sedimentation test cases and compare the results to existing analytical solutions. The grain velocity in Single Particle Sedimentation compares well (< 2% error) with the analytical solution as long as the ﬂuid resolution is coarser than two times the particle diameter. The multiple particle sedimentation problem and Rayleigh Taylor Instability (RTI) also perform well against the theory, but it was found that the method is susceptible to ﬂuid velocity ﬂuctuations in the presence of high porosity gradients. These ﬂuctuations can be damped by the addition of a dissipation term, which has no effect on the terminal velocity but can lead to slower growth rates for the RTI.


INTRODUCTION
Fluid-particle systems are ubiquitous in nature and industry, occurring in sediment transport and erosion, the rheology of avalanches, slurry flows and soils, industrial fluidized beds and the dispersion and mixing of particles in the food, chemical and painting industries.
The length-scale of interest determines the method of simulation for fluid-particle systems. For very small scale processes it is feasible to fully resolve the interstitial fluid between the particles, but for many applications this is infeasible and it becomes necessary to use unresolved, or mesoscale fluid simulations. Such a mesoscale type method is the focus of this paper and the domain of applicability for the SPH-DEM method.
Fluid-particle simulations at the mesoscale are often given the term Discrete Particle Models (DPM). These models fully resolve the individual solid particles using a Lagrangian model for the solid phase. The fluid phase does not resolve the interstitial fluid, but instead models the locally averaged Navier-Stokes equations and is coupled to the solid particles using appropriate drag closures. Most of the prior work on DPMs have been done using grid-based methods for the fluid phase, and a good review of the area can be seen in Zhu et al. [1] We present a DPM method based on the coupling of Smoothed Particle Hydrodynamics (SPH) for the fluid phase and DEM for the solid particles. This results in a purely particle-based solution method and therefore enjoys the flexibility that is inherent in these methods. This is the primary advantage of this method over existing grid-based DPMs. In particular, the model described in this paper is well suited for applications involving a free surface, including (but not limited to) debris flows, avalanches, landslides, sediment transport or erosion in rivers and beaches, slurry transport in industrial processes (e.g. SAG mills) and liquid-powder dispersion and mixing in the food processing industry.

SPH FLUID PHASE
Here we briefly describe the governing SPH equations for the fluid phase, based on the locally averaged Navier-Stokes equations (LANSEs) derived by Anderson and Jackson [2]. For more details see Robinson et al. [3], Robinson et al. [4].
We define a smooth porosity field by smoothing out the DEM particle's volumes according to the SPH interpolation kernel W a j (h) = W (r a − r j , h) where V j is the volume of DEM particle j. For readability, sums over SPH particles use the subscript b, while sums over surrounding DEM particles use the subscript j.
To calculate the continuity and momentum equations in the LANSEs, we first define a superficial fluid density ρ equal to the intrinsic fluid density scaled by the local porosity ρ = ερ f .
Substituting the superficial fluid density into the averaged continuity and momentum equations reduces them to the normal Navier-Stokes equations. Therefore, our approach is to use the weakly compressible SPH equations with variable h (resolution/smoothing length) terms [5,6,7] and adding fluid-particle drag terms (as specified below).
The rate of change of superficial density becomes , The SPH acceleration equation is given by where f a is the coupling force on the SPH particle a due to the DEM particles. The viscous term Π ab models the divergence of the viscous stress tensor and is calculated here using the term proposed by Monaghan [8] where α = 10μ/(ρhc s ), u sig = c s + u n /|r ab | and ρ ab = 0.5(ρ a + ρ b ).
The fluid pressure in Eq. (3) is calculated using the weakly compressible equation of state where the reference density ρ 0 is scaled by the local porosity to ensure that the pressure is slowly varying with porosity as: The smoothing length h a varies according to the superficial density (and hence with the porosity) and is calculated by h a = 1.5(m a /ρ a ) 1/3 .

DEM SOLID PHASE
Given a DEM particle i with position r i , the equation of motion is where m i is the mass of particle i, c i j is the contact force between particles i and j (acting from j to i) and f i is the fluid-particle coupling force on particle i. For the simulations presented below, we have used the linear spring dashpot contact model where δ is the overlap between the two particles and n i j is the unit normal vector pointing from j to i.
The force on each solid particle by the fluid is [2] where V i is the volume of particle i. The first two terms model the effect of the resolved fluid forces (buoyancy and shear-stress) on the particle. The fluid pressure gradient and the divergence of the stress tensor can be obtained from the SPH momentum equation given in Eq.
(3), and are evaluated at each solid particle, using a Shepard corrected [9] SPH interpolation.
The force f d models the drag effects of the unresolved (i.e. smoothed) variations in the fluid variables and is calculated from the local porosity ε i and the superficial velocity u s = ε i (u f −u i ). These two values are calculated at each DEM particle position, again using a Shepard corrected SPH interpolation. For the results in this paper we use both the simple Stokes drag force and a more general drag law proposed by Di Felice [10] The coupling force on SPH particle a is determined by a weighted average of the fluid-particle coupling force on the surrounding DEM particles.
where f j is the coupling force calculated for each DEM particle using Eq. (8) and is a correction factor to guarantee equal and opposite forces between the two phases.

SINGLE PARTICLE SEDIMENTATION
The first test case models a single particle sedimenting (SPS) in a 3D fluid column under gravity. The water column has a height of h = 0.006m and the bottom is a no-slip boundary. The boundaries in the x and y directions are periodic with a width of w = 0.004 m and gravity acts in the negative z direction. The single DEM particle is initialised at z = 0.8h. It has a diameter d = 10 −4 m and a density ρ p = 2500 kg/m 3 .
For the initial conditions of the simulation, the position of the DEM particle is fixed and the SPH fluid is allowed to reach hydrostatic equilibrium. The particle is then released at t = 0 s. In Fig. 1 the evolution of a DEM particle's vertical speed in water is shown for one-way and two-way coupling and for a reference fluid with parameters corresponding to water. We have performed similar simulations (data not shown) with fluids corresponding to air and a water-glycerol mixture [3]. The SPH-DEM results reproduce the analytical velocity curve within 0.3-1% error besides short-lived higher deviations at the initial onset of motion (approx 5%). One of the key assumptions of the SPH-DEM method (and any fluid-particle method that uses an unresolved fluid phase) is that the fluid resolution is sufficiently greater than the DEM particle diameter d. This ensures a smooth porosity field calculated via Eq. (1). By varying the fluid resolution h, we have found that accurate results are achieved as long as h ≥ 2d.

SEDIMENTATION OF A CONSTANT POROSITY BLOCK (CPB)
The second test case follows the sedimentation of a rigid porous block with constant porosity ε, here called a Constant Porosity Block (CPB). It has a cuboid shape with a width and depth equal to the water column and a height of h/2. The CPB is modelled by a regular grid of DEM particle that cannot move relative to each other, with a separation determined by ε. As the CPB falls, the fluid is displaced by its volume and flows upward through the DEM particles, affecting the terminal velocity. All the simulations use the Di Felice drag law, which is necessary to incorporate the effects of neighbouring particles (lower porosity) on the drag force.
Varying the block porosity ε allows us to evaluate the accuracy of the SPH-DEM model at different porosities. Figure 2 shows the average terminal velocity of the CPB over a range of porosities from ε = 0.6 to 1.0. Results using both water as the interstitial fluid are shown on the same plot by scaling the y-axis by the expected terminal velocity of a single DEM particle. The average terminal velocity is taken after the CPB has reached a steady terminal velocity.
The results for two different fluid resolutions are shown, h/d = 2 and 6. The lower resolution results show a systematically increased terminal velocity due to reduced drag at the edges of the block. This is caused by an excessive smoothing of the porosity discontinuity by the large width of the smoothing kernel, a feature not restricted to SPH-DEM but common to any fluid-particle method that uses an unresolved fluid phase. However, reducing the fluid resolution to h/d = 2 gives results better than 5% deviation (for the smallest ε) over the range of porosities tested.
At lower porosities the vertical velocity of the block suffers from increasing fluctuation around the mean. This is a consequence of fluctuations in the fluid velocity near the edges of the block. We have found that large porosity gradients cause an instability in the SPH particles. For this test case it was sufficient to add a small amount of artificial viscosity (we increased α by α art = 0.1 in Eq. (4)) to the simulations to ensure accurate results.

RAYLEIGH-TAYLOR INSTABILITY (RTI)
The third test case is identical to the CPB, but now the particles are allowed to move freely. This setup is similar in nature to the classical Rayleigh-Taylor (RT) instability, where a dense fluid is accelerated (normally via gravity) into a less dense fluid. The combination of particles and fluid can be modelled as a two-fluid system with the upper "fluid" (suspension) having an effective density ρ d , and an effective viscosity μ d , both higher than the properties of the lower fluid without particles. From this an expected growth rate can be calculated for the instability and compared with the simulated growth rate.
In Fig. 3, the growth of the RT instability versus time for ε = 0.8, fluid resolution h/d = 2 is shown using water as the surrounding fluid. The symbols show the vertical position of the lowest DEM particle, which provides an approximate measure of the instability amplitude. The two-fluid model [11] is included here as a benchmark, but it should be noted that this model contains some significant approximations and is not necessarily more accurate than the SPH-DEM results. While a constant porosity of 0.8 is used for the two-fluid RTI model, the porosity of the DEM particles ranges from 0.8 ≤ ε ≤ 0.93 during the growth of the instability. To account for this variation in porosity, we instead use the analytical model to obtain upper (using ε = 0.93) and lower (using ε = 0.8) bounds to the expected instability growth.
The SPH-DEM results are shown for the cases where the artificial viscosity is either applied (α art = 0.1) or not used (α art = 0.0). In both cases there is a clear exponential growth of the RT instability, but while the α art = 0.0 case shows an accurate growth rate of the instability, the addition of artificial viscosity slows down the growth rate erroneously. Therefore, while for the majority of applications (such as the CPB test case) the addition of a small amount of artificial viscosity has no significant effect on the results and is successful in eliminating the problematic velocity fluctuations seen near high porosity gradients, care must be taken when the results are sensitive to the fluid viscosity.

CONCLUSION
We have presented an SPH implementation of the locally averaged Navier Stokes equations and coupled this with a DEM model in order to provide a simulation tool for oneor two-way coupled fluid-particle systems. One notable property of the resulting method is that it avoids the use of a mesh and is completely particle-based. It is therefore suitable for those applications where a mesh presents additional problems, for example, free surface flow or flow around complex and/or moving geometries.
The SPH-DEM formulation was applied to three 3D sedimentation test cases and compares very favorably with known analytical solutions. We are currently applying this method to the dispersion of solids in fluid or fluid-gas environments. Other relevant directions for future developments are: removal of SPH velocity fluctuations near high porosity gradients; effects of different drag laws and the inclusion of the added mass and lift forces; the effects of different DEM particle contact forces and the inclusion of friction and lubrication forces; and the inclusion of surface tension effects.