Mathematical Problems in Engineering

Volume 2017 (2017), Article ID 3174904, 7 pages

https://doi.org/10.1155/2017/3174904

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

^{1}College of Water Resources & Civil Engineering, China Agricultural University, Beijing 100083, China^{2}State Key Laboratory of Simulation and Regulation of Water Cycle in River Basin, China Institute of Water Resources and Hydropower Research, Beijing 100038, China

Correspondence should be addressed to Liu-Chao Qiu and Yu Han

Received 30 September 2016; Revised 14 November 2016; Accepted 19 December 2016; Published 10 January 2017

Academic Editor: Payman Jalali

Copyright © 2017 Liu-Chao Qiu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

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.

#### 1. 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 time-consuming physical experiments. There are generally two types of coupling model [1] applied to calculate the interaction between fluid and solids: one-way coupling and two-way 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 free-surface 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.

#### 2. Methodology

##### 2.1. SPH for Fluid Phase

For the fluid phase, the following continuity and momentum equations are employed and given bywhere is the time, is the velocity, is the pressure, 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 are, respectively, the mass, density, and position of a particle, is the smoothing length, 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 in 2D and in 3D, .

Following Morris et al. [21] and Lo and Shao [22], (1) can be written in SPH discretization aswhere , 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 bywhere constant , is the reference density, is the sound speed at the reference density, and for a fluid like water. The speed of sound 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 and is a constant, whose values range between zero and unity, and 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.

##### 2.2. DEM for Solid Phase

DEM coupled with SPH has been used to numerically model the moving solids in viscous free-surface flows. The motion of each solid is tracked based on Newton’s laws of motion as follows:where solid possesses mass , velocity , and angular velocity . is gravity acceleration. Force is due to fluid-solid interaction and force 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 , 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 can be determined aswhere 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 aswhere is the radius of surface sphere belonging to and represents the normal unit vector pointing from the center of sphere belonging to to the center of sphere belonging to . and 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 are calculated bywhere is the normal spring stiffness, the overlap , where and are the surface sphere radii, and and are the surface sphere centers. is the normal relative velocity, and is the normal damping coefficient given bywhere is the restitution coefficient. The effective mass , where and are the mass of surface sphere.

The tangential component of the contact force is calculated by a Coulomb friction law using a coefficient of friction . It can be expressed aswhere , , and are the tangential spring stiffness, tangential overlap, and tangential damping coefficient, respectively. The tangential relative velocity and the tangential unit vector . The tangential damping coefficient given bywhere and . The value of causes the normal and tangential springs to react on similar time scales [24].

#### 3. 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 kg/m^{3} and kinetic viscosity m^{2}/s. For all simulations, the smoothing length is defined as , where represent the initial particle spacing.

##### 3.1. 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 0.5 m above a flat floor. The solid sphere has a density of 2.6 × 10^{3} kg/m^{3}, a diameter of 0.1 m, and a normal spring stiffness of 1 × 10^{5} N/m. The collision is assumed to be frictionless and purely elastic. So, coefficient of friction and the restitution coefficient . A time step of 1.0 × 10^{−4} s and an initial fluid particle spacing m are used in the simulation. The sphere is dropped and falls under the force of gravity until it contacts the level surface. The translational motion of the sphere can be described in three stages: free fall, contact, and rebound. Figure 1 shows snapshots of the falling sphere at different time. Figure 2 compares the vertical velocity and the height of the solid as a function of time. The results agree very well with the analytical solution given in the work of Chen et al. [25].