A 3 D Simulation of a Moving Solid in Viscous Free-Surface Flows by Coupling SPH and DEM

This work presents a three-dimensional two-way coupled method to simulate moving solids in viscous free-surface flows.The fluid flows are solved by weakly compressible smoothed particle hydrodynamics (SPH) and the displacement and rotation of the solids are calculated using the multisphere discrete element method (DEM) allowing for the contact mechanics theories to be used in arbitrarily shaped solids. The fluid and the solid phases are coupled through Newton’s third law of motion. The proposed method does not require a computational mesh, nor does it rely on empirical models to couple the fluid and solid phases. To verify the numerical model, the floating and sinking processes of a rectangular block in a water tank are simulated, and the numerical results are compared with experimental results reported in published literatures. The results indicate that the method presented in this paper is accurate and is capable of modelling fluid-solid interactions with a free-surface.


Introduction
Fluid-solid interaction in free-surface flows has gained a lot of interests in recent years due to its involvement in a wide variety of industrial processes and engineering disciplines such as hydraulic, ocean, and environmental engineering.Accurate numerical simulations are helpful in obtaining detailed information of the interaction process and in improving the engineering design while avoiding tedious and timeconsuming physical experiments.There are generally two types of coupling model [1] applied to calculate the interaction between fluid and solids: one-way coupling and twoway coupling.Substantial efforts have been made in the past to develop numerical methods for two-way coupling of fluid and solid.For example, Singh et al. [2] used the distributed Lagrange multiplier method for particulate flows with collisions.Takizawa et al. [3] modelled fluid-solid interactions with free-surface flows using CIP method.Latham et al. [4] simulated discrete solids in motion within fluids by coupling of DEM and CFD based on the arbitrary Lagrangian-Eulerian (ALE) frame.Wu et al. [5] presented a two-way coupled Moving Solid Algorithm (MSA) to calculate the motion of the solids in free-surface flows.In contrast to the above-mentioned mesh-based techniques, meshless methods do not require an interface capturing scheme or the moving mesh technology, which is a clear advantage when the problem involves breaking waves and moving boundary.Of the meshless techniques now available, smoothed particle hydrodynamics (SPH) is proving to be popular and robust.It was originally developed for astrophysical simulations [6,7] and has been extended to simulate free-surface flows by Monaghan [8].Instead of using a mesh, the SPH method uses a set of interpolation nodes placed arbitrarily within the fluid.This gives several advantages in comparison to mesh-based methods when simulating nonlinear flow phenomena.More complete reviews on SPH can be found at [9,10].Attempts have been made in most recently to develop SPH based method for two-way coupling of solid and fluid in a freesurface flow.Hashemi et al. [11] modelled free floating bodies using a ghost particle method which is usually used for fairly regular boundaries.Liu et al. [12] simulated fluid-structure coupling problems, especially for moving structures, using incompressible smoothed particle hydrodynamics (ISPH) model.Ren et al. [13] developed a two-dimensional SPH-DEM model to simulate the wave-structure interaction by describing the movement of the solids based on multisphere DEM introduced by Favier et al. [14].Canelas et al. [15] coupled SPH and DEM to describe the motion of arbitrarily shaped solids in viscous fluids based on the concept of DCDEM introduced by Cummins and Cleary [16].Yang et al. [17] developed an improved SPH-EBG coupling method for modelling the interaction of fluid with free-surface and flexible structure with free and fixed ends.
This work presents a 3D coupled DEM-WCSPH method to simulate moving solids in viscous free-surface flows.The fluid flow was solved by a weakly compressible smoothed particle hydrodynamics (WCSPH) which has an advantage of avoiding solution of the pressure Poisson equation.The motion of the solids was determined by the discrete element method and does not need to be prescribed beforehand.The DEM is also a Lagrangian method which was developed by Cundall and Strack [18] to describe granular material and nowadays is widely used in particulate flows.The DEM model in our study has been implemented using a multisphere approach introduced by Favier et al. [14] to enable solids with different complex shapes to be modelled.The accuracy of the simulation results was verified by comparison with published laboratory experiments of falling and floating solids in water.

Methodology
2.1.SPH for Fluid Phase.For the fluid phase, the following continuity and momentum equations are employed and given by where  is the time, v is the velocity,  is the pressure, g is gravity acceleration,  is the kinetic viscosity, and  is the density.
The key idea of SPH is that the concerning fluid domain is discretized by use of the moving particles which carry all the properties, such as mass, density, velocity, and pressure.And any quantity  of a particle can be approximated by summing over the relevant quantities of its neighboring particles within its support domain [9,10]: where , , and r are, respectively, the mass, density, and position of a particle, ℎ is the smoothing length,   = |r  − r  | being the distance between particles  and , and  is the smoothing kernel.In general, the higher the order of the kernels, the greater the accuracy of the SPH scheme.The Wendland kernel has some advantages over the Gaussian and the Cubic Spline [19].So, it is used in this study and written by [20] where   is 7/4ℎ 2 in 2D and 7/8ℎ 3 in 3D,  =   /ℎ.Following Morris et al. [21] and Lo and Shao [22], (1) can be written in SPH discretization as where k  = k  − k  , and ∇  denotes the gradient taken with respect to the coordinates of particle .
As mentioned before, the fluid phase is treated as weakly compressible where the density in the fluid is allowed to vary slightly.The closure problem is solved by using an equation of state to relate changes in pressure to changes in density.Following Monaghan [8] and Batchelor [23], the relationship between pressure and density for a particle  is given by where constant  =  2 0  0 /,  0 is the reference density,  0 is the sound speed at the reference density, and  = 7 for a fluid like water.The speed of sound  0 is generally chosen as 10 times the maximum velocity in the fluid to ensure the fluctuations in density are less than 1%.
When solving any boundary value problem, the proper implementation of boundary conditions is crucial.Due to its Lagrangian description and to kernel truncation, it is challenging to implement solid boundaries in SPH.In this work, the so-called dynamic boundary conditions [24] are used for both the moving and fixed boundary.It consists of creating boundary particles that satisfy the same equations of continuity (4) and state (6) as the fluid particles, but their positions are not updated using the momentum equation ( 5) and they remain fixed in position (fixed boundaries) or move according to some externally imposed function (moving objects).
Finally, particle positions updated every time step using the XSPH variant [10]: where   = (  +   )/2 and  is a constant, whose values range between zero and unity, and  = 0.5 is often used.This method is a correction for the velocity of particle .This correction lets particles to be more organized and, for high fluid velocities, helps to avoid particle penetration.

DEM for Solid
Phase.DEM coupled with SPH has been used to numerically model the moving solids in viscous freesurface flows.The motion of each solid is tracked based on Newton's laws of motion as follows: where solid  possesses mass   , velocity U  , and angular velocity Ω  .g is gravity acceleration.Force F   is due to fluidsolid interaction and force F   represents any solid contact that might occur.Integrating (8) into time advances the linear motion of the solid, whereas (9) accounts for the rotational motion.
In order to calculate the fluid-solid interaction force F   , the boundary of a moving solid is represented by groups of SPH particles and these moving boundary particles have an interparticle spacing that is equal to half of the fluid particle spacing to prevent the fluid particles from penetrating the moving solid boundary.The fluid-solid interaction is ensured through the interactive force balance condition satisfying Newton's third law of motion.Following Ren et al. [13], the fluid-solid interaction force F   can be determined as where the inner summation means the total force on moving boundary particle  belonging to solid  due to the neighborhood fluid particle .Following Latham et al. [4], a multisphere approach originally proposed by Favier et al. [14] was used for modelling the complex-shaped multibody dynamics in which the surface of each solid is represented by a cluster of small spheres whose diameter is equivalent to the spacing of the moving boundary particles.The solids are allowed to interact via contact forces when the surface spheres of different solids overlap.No relative movement between spheres of the same body is allowed.Based on the multisphere approach, the resultant contact forces and torque acting upon solid  are evaluated as where   is the radius of surface sphere  belonging to  and n  represents the normal unit vector pointing from the center of sphere  belonging to  to the center of sphere  belonging to . f , and f , are the normal and tangential forces, respectively, on surface sphere  of solid  due to surface sphere  of solid .Based on the model by Cundall and Strack [18], the normal contact forces f , are calculated by where   is the normal spring stiffness, the overlap   = (  +  )−|x  −x  |, where   and   are the surface sphere radii, and x  and x  are the surface sphere centers.v , is the normal relative velocity, and   is the normal damping coefficient given by where   is the restitution coefficient.

Numerical Examples and Discussion
In this section, numerical examples are presented to validate both solid-solid and fluid-solid interactions.For the cases involving water, the initial velocity of the fluid particles is considered zero.The fluid particles are assigned an initial density based on hydrostatic pressure.Therefore, the density of particle  (located at depth ) is calculated taking in account of the water column height as where (, ) is the initial water depth at position (, ) and  is the vertical distance from the bottom.The pressure is calculated using equation of state (6) following the density value.In all the simulations presented herein, the modelled fluid is water with reference density  0 = 1000 kg/m 3 and kinetic viscosity  = 1 × 10 −6 m 2 /s.For all simulations, the smoothing length is defined as ℎ = 1.2(Δ 2 + Δ 2 + Δ 2 ) 1/2 , where Δ = Δ = Δ = Δ represent the initial particle spacing.

Solid Dropping onto a Flat
Floor.In this simulation, a sphere is dropped under the influence of gravity onto a flat  floor.The purpose of this simulation is to verify that the momentum exchange between a solid and a flat boundary agrees with theory based on the principle of impulse and also to see if the motion of the solid subject to a gravitational force is accurate.The problem contains one sphere whose bottom is    flow.The simulation involves a cube falling into calm water in a tank.The trajectory of the cube is compared to laboratory experiments proposed by Wu et al. [5].In their experiments, a cube was released by a clamp at the water surface and the trajectory of the cube was recorded by a high speed camera.In order to reduce the water splash caused by the cube falling into water, the bottom of the cube was immersed in the water before the cube was released.The clamp was controlled by an electromagnetic device so the cube could be released instantly (see Figure 3).
The problem domain contains a water tank of 150 mm × 140 mm × 140 mm, and the water depth in the tank was  131 mm.The cube with a side length of 20 mm was assumed be a rigid as in the experiment.The density of the cube was 2120 kg/m 3 , the normal spring stiffness was 1 × 10 5 N/m, the coefficient of friction  = 0, and the restitution coefficient  = 0.5.A time step of 1.0 × 10 −4 s and an initial fluid particle spacing Δ = 2 mm are used in the simulation.Figure 4 compares the vertical displacements of the cube to the experimental results in Wu et al. [5].The simulation results are in good agreement compared with the measured trajectory in the vertical direction.Figure 5 presents some snapshots at different instants of falling cube evolution with the velocity magnitude (V = √V 2  + V 2  + V 2  ) on the central cutting section.The color bar is common to all snapshots and velocities are in meters per second.

Floating Solid in Water.
In this simulation, a wooden block rising from the bottom of a water tank was considered.Figure 6 shows the initial configuration of the problem.The numerical results are also compared with the experimental data provided by Wu et al. [5].The size of the water tank was 150 mm × 140 mm × 140 mm, and the water depth in the tank was 52 mm.A rectangular wooden block (width = 48 mm, length = 49 mm, height = 24 mm, and density 800.52 kg/m 3 ) was released from the bottom of the water tank.In our simulation, the normal spring stiffness was 1 × 10 5 N/m, the coefficient of friction  = 0, and the restitution coefficient  = 0.5.A time step of 1.0 × 10 −4 s and an initial fluid particle spacing Δ = 2 mm are used in the simulation.Figure 7 shows the comparison of measured and predicted vertical displacements of the block.It shows a good agreement between the simulated results and experimental data by Wu et al. [5]. Figure 8 shows the snapshots of the configuration of the floating solid at different time and the calculated freesurface profiles are very similar to that of the laboratory results.

Conclusions
A three-dimensional numerical model has been developed to simulate moving solids interacting with a free-surface flow by coupling the smoothed particle hydrodynamics and the discrete element method.The motions of the solids do not need to be prescribed beforehand and the mesh-free method is used for the computational fluid dynamics.The DEM model has been implemented using a multisphere approach to enable solids with arbitrarily complex shapes to be modelled.The comparison of results from the coupled approach with known analytical solutions and experimental data indicates that the proposed methods could be an efficient way of simulating complex moving solids in viscous freesurface flows.Our numerical examples focus mainly on fluidsolid interaction; however, multibody systems with solidsolid contacts are also within the capabilities of the model.Further investigations might be needed to consider the solid deformable and turbulent effects.

Figure 4 :
Figure 4: Comparison of measured and predicted vertical displacements.

Figure 5 :
Figure 5: The falling cube in the water tank with plotting the velocity magnitude.

Figure 6 :
Figure 6: Initial configuration of the floating solid problem.

Figure 7 :
Figure 7: Comparison of measured and predicted vertical displacement of the floating solid.

Figure 8 :
Figure 8: Snapshots of numerical (left) and experimental (right) results of the floating block.
The effective mass  * = (  +   )/    , where   and   are the mass of surface sphere.The tangential component of the contact force f , is calculated by a Coulomb friction law using a coefficient of friction   .It can be expressed as    t  −   k ,      f ,      <        f ,      −       f ,      t       f ,      >        f ,      , (14) where   ,   , and   are the tangential spring stiffness, tangential overlap, and tangential damping coefficient, respectively.The tangential relative velocity k , = k  − k , and the tangential unit vector t  = k , /|k , |.The tangential damping coefficient given by