# Coupled DEM-SPH Method for Interaction between Dilated Polyhedral Particles and Fluid

Authors: Shunying Ji; Xiaodong Chen; Lu Liu
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

The discrete element method (DEM) and smoothed particle hydrodynamics (SPH) can be adopted to simulate the granular materials and fluid media respectively. The DEM-SPH coupling algorithm can be developed for the dynamic interaction between the two media. When the particle material is simulated by polyhedral element, a fluid-solid coupling interface would lead to the complex geometry between the granular particle and the fluid. The boundary particle method is traditionally used for the fluid-solid interface but with low computational efficiency. In this paper, the dilated polyhedral element is constructed based on Minkowski sum theory, while the contact force between the elements is calculated by Hertzian contact model. Accordingly the dilated polyhedra based DEM is established. The weakly compressible SPH is adopted to simulate the fluid medium, while the interaction on the geometrically complex fluid-solid interface is evaluated with the repulsive force model which can be determined by the contact detection between SPH particles and solid particles in geometry. This method avoids the storage and calculation of a large number of boundary particles, which can be potentially applied for the complex fluid-solid boundary. In order to improve the computational efficiency, a GPU-based parallel algorithm is employed to achieve high performance computation of SPH. The acceleration of the parallel algorithm is evaluated by the cases of dam break. The numerical simulation of the impact of dam break on cubes is implemented. The simulation results are verified with the corresponding experimental and simulation results. Therefore, the rationality and accuracy of the DEM-SPH coupling method for numerical simulation of the interaction between granular materials and fluid media are illustrated. This method is then adopted for the impact of falling rocks on underwater pipeline. The force of water and rocks on the pipeline is analyzed. This method can be further applied for real engineering problems.


---

## Body

## 1. Introduction

The coupling of granular materials and fluid media is widely existed in geological hazards, geotechnical engineering, ocean engineering, and process engineering [1–3]. The discrete element method (DEM) is usually adopted to simulate solid particles, while the computational fluid dynamics (CFD) is adopted to simulate the fluid. Most of the coupling effects of the two methods are based on the interface coupling, i.e., coupling force on the interface balances.The dilated polyhedral element is generated by the Minkowski sum of a sphere and an arbitrary polyhedron, which can be used effectively for the numerical simulation of irregular granular materials [4, 5]. Because the dilated polyhedral element has the geometric characteristics of polyhedron, the dilated polyhedron can be employed to construct the geometrically continuous mesh grid. Then the bonding and breaking model based on the dilated polyhedron can be developed to simulate the failure process of the continuum materials [6]. Meanwhile, the lattice Boltzmann method (LBM) can be used to establish the DEM-LBM method for the fluid-solid coupling to analyze the hydraulic fracturing phenomena of brittle materials [7, 8]. Since the dilated polyhedral element transforms sharp corners and edges into smooth spherical and cylindrical surfaces, it avoids the disadvantage of pure polyhedron in contact detection that it is difficult to detect contact directly by geometric features. So the dilated polyhedron improves the efficiency and stability of contact detection. In addition, the Hertzian contact force model can be employed to calculate the contact force between smooth surfaces. Therefore, the dilated polyhedral element has a good application in the numerical simulation of the mechanical behavior of irregular granular materials and the failure process of brittle materials.Traditionally the mesh-grid CFD methods have the advantages of solid theoretical basis and mature application in engineering. A lot of open source and commercial software can be used for the DEM-CFD coupling simulations [9]. However, mesh-grid CFD method is difficult to deal with free surface problems, and the mesh distortion can easily occur because of the large deformation at the fluid-solid interface. Besides, the computational efficiency of tracking and identifying free surface and coupling complex shaped solid is also cumbersome, which restricts the large-scale coupling simulation with DEM. In recent years, the LBM based on mesoscale statistical theory develops rapidly for fluid simulations. Because of the microparticle characteristics of LBM, it is convenient to describe the interaction between multiphase flows, and the DEM-LBM coupling method is developed [10]. Because the number of lattices required by LBM for real engineering simulations is quite large, the computational efficiency of LBM in real engineering applications is low. In addition, LBM usually utilizes gas and liquid simultaneously to deal with the free surface, which also reduces the computational efficiency. Therefore, the DEM-LBM coupling is mostly used for the basic research such as multiphysical field and multiphase flow, while it is relatively rare for large-scale simulations in engineering.The smoothed particle hydrodynamics (SPH) based on Lagrangian particle method was first proposed by Lucy, Gingold and Managhan in 1977 [11, 12]. It originated from the research of astrophysics and then employed in the fields of electromagnetics, heat conduction, nuclear physics, and computational fluid dynamics [13, 14]. This method can effectively avoid mesh distortion, adaptively simulate free surface, and naturally track moving interfaces. It has advantages in simulations of large deformation and free surface. In SPH, the dynamic relaxation approach based on the weakly compressible SPH (WCSPH) can solve the basic continuity, momentum and energy equations respectively. The method is widely used in the fluid-solid coupling analysis for it is simple and easy to be implemented [15]. In addition, the Incompressible SPH (ISPH) method based on Poisson pressure equation (PPE) can calculate pressure more accurately, which is more stable than WCSPH. But it costs a lot of computational time to solve the matrix equation of pressure, and thus the computational efficiency of ISPH is much lower [16]. An explicit scheme of ISPH, also named as EISPH, is developed to solve the PPE using SPH summations [17], which avoid the calculation of matrix equation. Another merit of EISPH is the larger time step can also be used for the time integration.In the fluid-solid coupling simulation of granular materials and fluid media, SPH and DEM have highly consistent data structure and operation logic in the implementation of computer program, which can be employed to simulate the fluid-solid coupling problem effectively and suitable for parallel computational environment [18–20]. Robinson et al. (2014) developed the DEM-SPH coupling method for particles of different sizes by the local average method [21]. The simulation results are verified with those of other numerical methods. Sun et al. (2015) used theδ-SPH method to simulate the fluid-structure interaction between the cylinder and the eccentric block [22], while the solid boundary is constructed through the boundary particles. The simulation results are in agreement with the experimental data. In addition, the integral value of the SPH particle smoothing function near the particle boundary is fitted into formulas. Based on the formulas, the interaction between SPH particle and complex boundary can be calculated without boundary particle in simulations [23]. At present, the boundary particle method is mainly used to deal with the fluid-solid boundary with simple and constant shape, while the local average method is mainly for solid particles with spherical shape. So it is difficult to construct complex geometric models. The boundary approach based on the boundary integral method adopts pure geometry boundary instead of traditional particle boundary. This boundary treatment can be used to establish DEM-SPH coupling method with the local average method [24]. It should be noted that the above methods utilize spherical elements in DEM, which cannot describe the real particle morphology. Therefore, the irregular shaped discrete element models for real granular materials, such as rock and broken ice, should be developed for the fluid-solid coupling simulation and avoiding the empirical selection of mesocalculation parameters as far as possible [25].In this paper, the dilated polyhedra based discrete element method according to the Minkowski sum theory is established. The weakly compressible SPH is adopted for the fluid simulation. The interaction between DEM and SPH is calculated by the repulsive force model. The DEM-SPH coupling algorithm is implemented in the GPU-based parallel environment. This method is validated with results of several corresponding examples to prove the reliability and accuracy.


## 2. DEM-SPH Coupling Method for Dilated Polyhedral Elements

### 2.1. Dilated Polyhedron Element Based DEM

For any two geometric bodiesA andB in space, the Minkowski sum of them can be defined as(1)A⊕B=x+y∣x∈A,y∈Bwhere x and y are coordinate points belong toA andB, respectively. IfA andB are set as arbitrary shaped polyhedron and sphere respectively, the Minkowski sum of them is the dilated polyhedron, as shown in Figure 1. The corner and edge of pure polyhedron can be transformed into spherical and cylindrical surface by the Minkowski sum theory. Accordingly, the contact pattern can be classified as six types, including sphere-sphere, sphere-cylinder, sphere-face, cylinder-cylinder, cylinder-face, and face-face. The contact force model in different contact pattern could be established by the Hertzian model [26], which overcomes the limitation of the linear contact force model and the aimlessness of the choice of stiffness parameters. Meanwhile, the dilated polyhedron avoids the singular phenomena caused by the uncertain contact normal when calculating the contact force at the corner or edge of polyhedron.Figure 1
Dilated polyhedral element based on the Minkowski sum theory.Because the geometry of dilated polyhedron is complex, it is complicated to calculate its mechanical parameters precisely, such as mass, centroid, moment of inertia, etc. Therefore, the minimum envelope polyhedron approximation of the dilated polyhedron is used instead of the dilated polyhedron itself to calculate the mechanical parameters.In the contact force model between dilated polyhedral elements, the normal contact force considers both elastic force and viscous force. It can be written as follows:(2)Fn=knδnκ-Cnδnκ-1δ˙nwhere kn is the normal stiffness of two dilated polyhedral elements in contact; δn is the normal overlap vector; δ˙n is the normal relative velocity vector; Cn is the effective normal damping coefficient, which can be calculated by referring to the contact force model of spherical elements [27];κ is usually 1.5 in the Hertzian model.The tangential force considers the Mindline model and the Mohr-Coulomb friction law. It can be expressed as follows:(3)Fs∗=ksδnκ-1δs(4)Fs=min⁡Fs∗,μFnδswhere δs is the vector of tangential elastic deformation; μ is the sliding friction coefficient between elements. ks is the tangential stiffness, and ks=rsnkn. Here the ratio rsn is defined as 1/2(1+ν) where ν is the Poisson's ratio, according to the relationship between the elastic modulus and shear modulus of the isotropic material [28].According to the surface geometric characteristics of dilated polyhedron, the contact between dilated polyhedrons is divided into six categories: sphere-sphere, sphere-face, sphere-cylinder, cylinder-cylinder, cylinder-face, and face-face contact, as shown in Figure2. Therefore, the contact points can be calculated by geometric detection of different contact types, from which the curvature radius of the two elements at the contact points can be obtained. Accordingly, the contact stiffness kn of different contact types can be determined, and thus the contact force can be calculated by the Hertzian contact model [28].Figure 2
Contact types of dilated polyhedral elements.
(a)
Sphere-sphere (b)
Sphere-face (c)
Sphere-cylinder (d)
Cylinder-cylinder (e)
Cylinder-face (f)
Face-face


### 2.2. Weakly Compressible SPH

SPH approximates the field function by the smoothing function in a certain range of space domain. The related physical quantities are represented by the particle approximation. Thus, the partial differential equations can be transformed into time-dependent ordinary differential equations. The variation of the field variables of each particle with time can be obtained by the time integration. For a continuous function in any fieldΩ, it can be approximated by the smoothing function as follows:(5)fx≈∫Ωfx′Wx-x′,hdx′whereW is the smoothing function;h is a smoothing length defining the effective range of the smoothing function. The smoothing function is an even function related to x. It must satisfy the normalization condition and compact condition. In SPH, the particle approximation of a field function at particlei can be expressed by the adjacent particles in its support domain.(6)fxi=∑j=1NmjρjfxjWijwhereN is the number of adjacent particles in the support domain;j denotes an adjacent particle of particlei; mj is the mass of adjacent particlej; ρj is the density of adjacent particlej; Wij=W(xi-xj,h). The divergence of the field function fxi is only related to the smoothing function and can be written as follows:(7)∇∙fxi=-∑j=1Nmjρjfxj∇iWij∇iWij=xijrij∂Wij∂rij=xi-xjrij∂Wij∂rijThe smoothing lengthh has an important influence on the computational efficiency and accuracy. Its purpose is to keep enough adjacent particles in the support domain of SPH particles to ensure the approximate validity of continuous variables. According to the above basic equations, the governing equations of fluid can be expressed by the particle approximation, including the continuity (mass conservation) equation and the momentum conservation equation. The continuity equation for the particle density can be written as follows:(8)dρidt=∑j=1Nmjvij∙∇iWijwhere vij=vi-vj. Accordingly the density of each particle can be determined by the time integration. The pressure of each particle can be calculated by the equation of state (EOS), which can be written as follows:(9)p=c02ρ0γρρ0γ-1whereρ0 is the initial density of the fluid; γ is generally equal to 7 [19]; c0 is the speed of sound, which is usually equal to 10 times of the maximum particle velocity, i.e., c0=10vmax.Ignoring the viscous stress of fluid, the momentum conservation equation can be written as follows:(10)dvidt=-∑j=1Nmjpiρi2+pjρj2∇iWijwhere pi and pj are the pressures of particlei and particlej, respectively. In order to eliminate the nonphysical oscillation in SPH simulations, the artificial viscous term by Monaghan is adopted:(11)Πij=-αIIc-ijϕij+βIIϕij2ρ-ij,vij∙xij<00,vij∙xij≥0where αII and βII are standard constants and generally range from 0 to 1.0; in addition, there are(12)ϕij=hijvij∙xijxij2+φ2(13)c-ij=12ci+cj(14)ρ-ij=12ρi+ρj(15)hij=12hi+hjwhere ci=cj=c0 and φ=0.1hij. Moreover, in order to avoid the tension instability caused by the negative pressure of particles near the free surface and local particle aggregation, the tension correction term is employed to improve the simulation [23, 24], written as(16)Rijfij4=RijWijW∆x,h4where(17)Rij=0.01piρi2+piρj2,pi>0andpj>00.2piρi2+pjρj2,elseHence, the final form of momentum equation yielded by particle approximation of SPH can be obtained, written as(18)dvidt=-∑j=1Nmjpiρi2+pjρj2+Πij+Rijfij4∇iWijThe explicit leap-frog (LF) method is used to integrate the above discretized equations [29]. The time step of integration is determined as follows:(19)∆t=0.25hc0In fact, the time step is closely related to the state change of fluid. So the time step considering the viscous dissipation property and particle acceleration can also be used [30, 31].


### 2.3. DEM-SPH Coupling Algorithms for Dilated Polyhedral Elements

The fluid-solid boundary is regarded as the boundary condition of SPH in this paper. Meanwhile, the SPH particles are subject to the repulsive force from the boundary, while the repulsive force will react on the DEM elements at the same time. Thus the coupling between the two media can be achieved. The traditional SPH repulsive force boundary model is composed of particles, which is not conducive to the complex shaped boundary [15]. In this article, the particle boundary is transformed into the geometric surface boundary as shown in Figure 3. According to the traditional SPH repulsive force boundary model, the repulsive force on SPH particles is calculated by the distance between particles and the boundary and the normal direction of the boundary, which can be expressed as follows:(20)fB=nBεzRdwhere fB is the boundary repulsive force; nB is the unit normal vector of the wall boundary; εz is the pressure correction term; d is the minimum distance between the SPH particle and boundary; R(d) is the repulsive force function, written as(21)Rd=A1q1-qwhere q is the normalized distance, and q=d/2h; A is related to the smoothing length and sound speed, and A=0.01c02/h.Figure 3
The repulsive force model at the boundary between solid and fluid.
(a)
Traditional particle boundary (b)
Geometric boundaryIn order to balance the pressure in different depths of fluid, the pressure correction termεz are employed here, written as(22)εz=0.02,z≥0zh0+0.02,-h0≤z<01,zh0>1where h0 is the initial height of fluid; z is the height of SPH particle. Here note that the water level is zero plane, i.e., z=0, and thus it is negative below the zero plane.In fact, the interaction force between dilated polyhedral particle and SPH particle only considers the normal direction ignoring the tangential direction. It is more like free-slip condition. The no-slip condition can be conducted by considering the tangential force between solid and fluid particles in the future.Generally, the polyhedral surface is treated as SPH boundary in DEM-SPH coupling. The force between SPH particles and DEM elements satisfies Newton's third law. The detection between SPH particles and DEM elements is solved in geometry to determine the minimum distance between them which is utilized to calculate the repulsive force. By this simplified boundary treatment approach, the fluid-solid coupling of DEM-SPH algorithm for dilated polyhedral elements and fluid can be established.


## 2.1. Dilated Polyhedron Element Based DEM

For any two geometric bodiesA andB in space, the Minkowski sum of them can be defined as(1)A⊕B=x+y∣x∈A,y∈Bwhere x and y are coordinate points belong toA andB, respectively. IfA andB are set as arbitrary shaped polyhedron and sphere respectively, the Minkowski sum of them is the dilated polyhedron, as shown in Figure 1. The corner and edge of pure polyhedron can be transformed into spherical and cylindrical surface by the Minkowski sum theory. Accordingly, the contact pattern can be classified as six types, including sphere-sphere, sphere-cylinder, sphere-face, cylinder-cylinder, cylinder-face, and face-face. The contact force model in different contact pattern could be established by the Hertzian model [26], which overcomes the limitation of the linear contact force model and the aimlessness of the choice of stiffness parameters. Meanwhile, the dilated polyhedron avoids the singular phenomena caused by the uncertain contact normal when calculating the contact force at the corner or edge of polyhedron.Figure 1
Dilated polyhedral element based on the Minkowski sum theory.Because the geometry of dilated polyhedron is complex, it is complicated to calculate its mechanical parameters precisely, such as mass, centroid, moment of inertia, etc. Therefore, the minimum envelope polyhedron approximation of the dilated polyhedron is used instead of the dilated polyhedron itself to calculate the mechanical parameters.In the contact force model between dilated polyhedral elements, the normal contact force considers both elastic force and viscous force. It can be written as follows:(2)Fn=knδnκ-Cnδnκ-1δ˙nwhere kn is the normal stiffness of two dilated polyhedral elements in contact; δn is the normal overlap vector; δ˙n is the normal relative velocity vector; Cn is the effective normal damping coefficient, which can be calculated by referring to the contact force model of spherical elements [27];κ is usually 1.5 in the Hertzian model.The tangential force considers the Mindline model and the Mohr-Coulomb friction law. It can be expressed as follows:(3)Fs∗=ksδnκ-1δs(4)Fs=min⁡Fs∗,μFnδswhere δs is the vector of tangential elastic deformation; μ is the sliding friction coefficient between elements. ks is the tangential stiffness, and ks=rsnkn. Here the ratio rsn is defined as 1/2(1+ν) where ν is the Poisson's ratio, according to the relationship between the elastic modulus and shear modulus of the isotropic material [28].According to the surface geometric characteristics of dilated polyhedron, the contact between dilated polyhedrons is divided into six categories: sphere-sphere, sphere-face, sphere-cylinder, cylinder-cylinder, cylinder-face, and face-face contact, as shown in Figure2. Therefore, the contact points can be calculated by geometric detection of different contact types, from which the curvature radius of the two elements at the contact points can be obtained. Accordingly, the contact stiffness kn of different contact types can be determined, and thus the contact force can be calculated by the Hertzian contact model [28].Figure 2
Contact types of dilated polyhedral elements.
(a)
Sphere-sphere (b)
Sphere-face (c)
Sphere-cylinder (d)
Cylinder-cylinder (e)
Cylinder-face (f)
Face-face


## 2.2. Weakly Compressible SPH

SPH approximates the field function by the smoothing function in a certain range of space domain. The related physical quantities are represented by the particle approximation. Thus, the partial differential equations can be transformed into time-dependent ordinary differential equations. The variation of the field variables of each particle with time can be obtained by the time integration. For a continuous function in any fieldΩ, it can be approximated by the smoothing function as follows:(5)fx≈∫Ωfx′Wx-x′,hdx′whereW is the smoothing function;h is a smoothing length defining the effective range of the smoothing function. The smoothing function is an even function related to x. It must satisfy the normalization condition and compact condition. In SPH, the particle approximation of a field function at particlei can be expressed by the adjacent particles in its support domain.(6)fxi=∑j=1NmjρjfxjWijwhereN is the number of adjacent particles in the support domain;j denotes an adjacent particle of particlei; mj is the mass of adjacent particlej; ρj is the density of adjacent particlej; Wij=W(xi-xj,h). The divergence of the field function fxi is only related to the smoothing function and can be written as follows:(7)∇∙fxi=-∑j=1Nmjρjfxj∇iWij∇iWij=xijrij∂Wij∂rij=xi-xjrij∂Wij∂rijThe smoothing lengthh has an important influence on the computational efficiency and accuracy. Its purpose is to keep enough adjacent particles in the support domain of SPH particles to ensure the approximate validity of continuous variables. According to the above basic equations, the governing equations of fluid can be expressed by the particle approximation, including the continuity (mass conservation) equation and the momentum conservation equation. The continuity equation for the particle density can be written as follows:(8)dρidt=∑j=1Nmjvij∙∇iWijwhere vij=vi-vj. Accordingly the density of each particle can be determined by the time integration. The pressure of each particle can be calculated by the equation of state (EOS), which can be written as follows:(9)p=c02ρ0γρρ0γ-1whereρ0 is the initial density of the fluid; γ is generally equal to 7 [19]; c0 is the speed of sound, which is usually equal to 10 times of the maximum particle velocity, i.e., c0=10vmax.Ignoring the viscous stress of fluid, the momentum conservation equation can be written as follows:(10)dvidt=-∑j=1Nmjpiρi2+pjρj2∇iWijwhere pi and pj are the pressures of particlei and particlej, respectively. In order to eliminate the nonphysical oscillation in SPH simulations, the artificial viscous term by Monaghan is adopted:(11)Πij=-αIIc-ijϕij+βIIϕij2ρ-ij,vij∙xij<00,vij∙xij≥0where αII and βII are standard constants and generally range from 0 to 1.0; in addition, there are(12)ϕij=hijvij∙xijxij2+φ2(13)c-ij=12ci+cj(14)ρ-ij=12ρi+ρj(15)hij=12hi+hjwhere ci=cj=c0 and φ=0.1hij. Moreover, in order to avoid the tension instability caused by the negative pressure of particles near the free surface and local particle aggregation, the tension correction term is employed to improve the simulation [23, 24], written as(16)Rijfij4=RijWijW∆x,h4where(17)Rij=0.01piρi2+piρj2,pi>0andpj>00.2piρi2+pjρj2,elseHence, the final form of momentum equation yielded by particle approximation of SPH can be obtained, written as(18)dvidt=-∑j=1Nmjpiρi2+pjρj2+Πij+Rijfij4∇iWijThe explicit leap-frog (LF) method is used to integrate the above discretized equations [29]. The time step of integration is determined as follows:(19)∆t=0.25hc0In fact, the time step is closely related to the state change of fluid. So the time step considering the viscous dissipation property and particle acceleration can also be used [30, 31].


## 2.3. DEM-SPH Coupling Algorithms for Dilated Polyhedral Elements

The fluid-solid boundary is regarded as the boundary condition of SPH in this paper. Meanwhile, the SPH particles are subject to the repulsive force from the boundary, while the repulsive force will react on the DEM elements at the same time. Thus the coupling between the two media can be achieved. The traditional SPH repulsive force boundary model is composed of particles, which is not conducive to the complex shaped boundary [15]. In this article, the particle boundary is transformed into the geometric surface boundary as shown in Figure 3. According to the traditional SPH repulsive force boundary model, the repulsive force on SPH particles is calculated by the distance between particles and the boundary and the normal direction of the boundary, which can be expressed as follows:(20)fB=nBεzRdwhere fB is the boundary repulsive force; nB is the unit normal vector of the wall boundary; εz is the pressure correction term; d is the minimum distance between the SPH particle and boundary; R(d) is the repulsive force function, written as(21)Rd=A1q1-qwhere q is the normalized distance, and q=d/2h; A is related to the smoothing length and sound speed, and A=0.01c02/h.Figure 3
The repulsive force model at the boundary between solid and fluid.
(a)
Traditional particle boundary (b)
Geometric boundaryIn order to balance the pressure in different depths of fluid, the pressure correction termεz are employed here, written as(22)εz=0.02,z≥0zh0+0.02,-h0≤z<01,zh0>1where h0 is the initial height of fluid; z is the height of SPH particle. Here note that the water level is zero plane, i.e., z=0, and thus it is negative below the zero plane.In fact, the interaction force between dilated polyhedral particle and SPH particle only considers the normal direction ignoring the tangential direction. It is more like free-slip condition. The no-slip condition can be conducted by considering the tangential force between solid and fluid particles in the future.Generally, the polyhedral surface is treated as SPH boundary in DEM-SPH coupling. The force between SPH particles and DEM elements satisfies Newton's third law. The detection between SPH particles and DEM elements is solved in geometry to determine the minimum distance between them which is utilized to calculate the repulsive force. By this simplified boundary treatment approach, the fluid-solid coupling of DEM-SPH algorithm for dilated polyhedral elements and fluid can be established.


## 3. Validation Examples for DEM-SPH Coupling Method

The dam break flow is simulated by sequential and GPU-based parallel programs, respectively, and the acceleration ratio of GPU parallel computing is analyzed. In order to verify the reliability of DEM-SPH coupling method, the dam break flow impacting against single cube and three cubes is simulated, while the displacement of cubes is validated with the corresponding tests and simulations.


### 3.1. Acceleration of GPU-Based Parallel Algorithm

Figure4 shows the size of initial water body and water tank. The initial water body size is 0.5 m in length, 1 m in height, and 0.5 m in width. The size of water tank is 2 m in length, 1.5 m in height, and 0.5 m in width.Figure 4
Sketch of dam break.Figure5 shows the SPH simulation of dam break. The wave front is generated after the water impact on the boundary, and then flow back to the other side. The position of water front versus time is obtained to verify with other approaches. As shown in Figure 6, the result is compared with VOF simulation result and experimental result [32, 33]. Here a dimensionless time 2g/a  is employed, where g is the gravity acceleration anda is the initial length of water. The simulation by the proposed SPH method is in good agreement with other approaches.Figure 5
SPH simulation of dam break.
(a)
t = 0.1s (b)
t = 0.3sFigure 6
Position of flow front versus time and the comparison with other approaches [32, 33].To study the computational efficiency of GPU, the program is compiled in CPU serial and GPU-based parallel frameworks to simulate the dam break process of fluid with different number of SPH particles. The physical model is the same with that in Figure4. The hardware environment of CPU serial program is mainly Intel Core i7-4790 (3.6 GHz); the hardware environment of GPU parallel program is Intel Xeon 4114 (2.2 GHz) plus Nvidia Tesla P100. Figure 7 shows the acceleration ratio of CPU serial and GPU parallel computing for 10 k steps with different number of SPH particles. The acceleration ratio of GPU increases fast in this hardware environment if the particle number is less than 20 k. When the number of particles exceeds 40k, the acceleration efficiency gradually stabilizes and approaches 10. Therefore, the GPU-based parallel algorithm can significantly improve the computational efficiency of SPH simulation.Figure 7
Acceleration ratio between GPU parallel and CPU serial algorithm.


### 3.2. DEM-SPH Coupling Simulation of Dam Break Flow Against Cubes

In order to validate the SPH-DEM coupling method, the fluid-solid interaction process of dam break flow impacting against cubes is simulated. Referring to the work of Canelas et al. (2016) [34], the cubes is placed at 1.7 m to the initial water, as shown in Figure 8. The main simulation parameters are listed in Table 1.Table 1
Simulation parameters of dam break flow in DEM-SPH method.
Parameter Value Water tank: length × width × height 8 m × 0.7 m× 1 m Cube size 0.15 m × 0.15 m × 0.15 m Dilated radius of cube 0.01 m Density of cube 800 kg/m3 Distance between cube to water 1.7 m Friction coefficient of water tank 0.35 Friction coefficient of cube 0.45 Poisson’s ratio of cube 0.3 Water size 3.5 m×0.7 m×0.4 m Smoothing length 0.02 m Number of SPH particles 119 000Figure 8
Setup of DEM-SPH numerical simulations of cubes in break flow.
(a)
Single cube (b)
Three cubesFigure9 is the simulation of the dam break flow against single cube, from which we can see that the cube slides under the impact of flow without rolling. Figure 10 is the time history of cube displacement alongx direction, which is compared with the results of corresponding experiment and numerical simulation [34]. The results show that the cube moves at a constant acceleration under the action of water impingement after the fluid impacts on the cube at about 0.8 s. The motion of the cube is basically consistent with the SPH-DCDEM result and experimental result by Canelas et al. (2016). It should be noticed that the experimental data are the average result of several tests [34].Figure 9
Dam break flow impacts on single cube with the proposed DEM-SPH method.
(a)
t = 0 s (b)
t = 0.9 s (c)
t = 1.6 s (d)
t = 2.2 sFigure 10
The displacement inx direction of single cube simulated with the DEM-SPH method and its comparison with the SPH-DCDEM and experimental results [34].Figure11 is the dam break flow impacting on three cubes simulated by the proposed DEM-SPH method, from which we can see that the cubes collapse under the impact of flow, and then the three cubes slide under the flow. Figure 12 is the photo comparison with the numerical and experimental results by Canelas et al. (2016) at two different times [34]. It can be found that the fluid-solid coupling phenomenon has little difference with the numerical and experimental results of Canelas et al. (2016). In the discrete element simulation of Canelas et al. (2016), the cube is constructed by spherical bonding element, so that the contact model of the cube element is different from that in this paper, which results in differences in simulation results.Figure 11
Dam break flow impacts on three cubes simulated with the DEM-SPH method.
(a)
t = 0 s (b)
t = 1.4 s (c)
t = 1.4 s (d)
t = 2.3 sFigure 12
The photo of dam break flow on three cubes with the proposed DEM-SPH method and its comparison with the DCDEM-SPH and experimental results: DEM-SPH simulation, experimental results, and DCDEM-SPH simulation from left to right [34].
(a)
t= 0.98 s (b)
t = 1.28 sIn the simulation of dam break flow on three cubes, the displacements of bottom particles inx direction and top cubes inx direction andz direction are compared with the corresponding simulation and experimental results, respectively, as shown in Figure 13. Similar to the previous example, the experimental data are the results of statistical average of several tests; i.e., each group of experimental results would fluctuate near the average value. In the simulation of three cubes, the motion of the cube is affected by many factors, such as the impact of water flow, the interaction with other cubes, and the friction with the bottom of the flume. The cube motion may be subject to more interference, which will produce randomness in the experiments. Therefore, the simulation result, including that by us and Canelas et al. (2016), and experiment result have little difference and the trends are the same. In general, the results in this paper are in good agreement with the corresponding simulation and experimental results.Figure 13
Displacements of the three cubes inx andz directions with the proposed DEM-SPH method and its comparison with the DCDEM-SPH and experimental results [34].
(a)
Displacement of bottom cube inx direction (b)
Displacement of top cube inx direction (c)
Displacement of top cube inz direction


## 3.1. Acceleration of GPU-Based Parallel Algorithm

Figure4 shows the size of initial water body and water tank. The initial water body size is 0.5 m in length, 1 m in height, and 0.5 m in width. The size of water tank is 2 m in length, 1.5 m in height, and 0.5 m in width.Figure 4
Sketch of dam break.Figure5 shows the SPH simulation of dam break. The wave front is generated after the water impact on the boundary, and then flow back to the other side. The position of water front versus time is obtained to verify with other approaches. As shown in Figure 6, the result is compared with VOF simulation result and experimental result [32, 33]. Here a dimensionless time 2g/a  is employed, where g is the gravity acceleration anda is the initial length of water. The simulation by the proposed SPH method is in good agreement with other approaches.Figure 5
SPH simulation of dam break.
(a)
t = 0.1s (b)
t = 0.3sFigure 6
Position of flow front versus time and the comparison with other approaches [32, 33].To study the computational efficiency of GPU, the program is compiled in CPU serial and GPU-based parallel frameworks to simulate the dam break process of fluid with different number of SPH particles. The physical model is the same with that in Figure4. The hardware environment of CPU serial program is mainly Intel Core i7-4790 (3.6 GHz); the hardware environment of GPU parallel program is Intel Xeon 4114 (2.2 GHz) plus Nvidia Tesla P100. Figure 7 shows the acceleration ratio of CPU serial and GPU parallel computing for 10 k steps with different number of SPH particles. The acceleration ratio of GPU increases fast in this hardware environment if the particle number is less than 20 k. When the number of particles exceeds 40k, the acceleration efficiency gradually stabilizes and approaches 10. Therefore, the GPU-based parallel algorithm can significantly improve the computational efficiency of SPH simulation.Figure 7
Acceleration ratio between GPU parallel and CPU serial algorithm.


## 3.2. DEM-SPH Coupling Simulation of Dam Break Flow Against Cubes

In order to validate the SPH-DEM coupling method, the fluid-solid interaction process of dam break flow impacting against cubes is simulated. Referring to the work of Canelas et al. (2016) [34], the cubes is placed at 1.7 m to the initial water, as shown in Figure 8. The main simulation parameters are listed in Table 1.Table 1
Simulation parameters of dam break flow in DEM-SPH method.
Parameter Value Water tank: length × width × height 8 m × 0.7 m× 1 m Cube size 0.15 m × 0.15 m × 0.15 m Dilated radius of cube 0.01 m Density of cube 800 kg/m3 Distance between cube to water 1.7 m Friction coefficient of water tank 0.35 Friction coefficient of cube 0.45 Poisson’s ratio of cube 0.3 Water size 3.5 m×0.7 m×0.4 m Smoothing length 0.02 m Number of SPH particles 119 000Figure 8
Setup of DEM-SPH numerical simulations of cubes in break flow.
(a)
Single cube (b)
Three cubesFigure9 is the simulation of the dam break flow against single cube, from which we can see that the cube slides under the impact of flow without rolling. Figure 10 is the time history of cube displacement alongx direction, which is compared with the results of corresponding experiment and numerical simulation [34]. The results show that the cube moves at a constant acceleration under the action of water impingement after the fluid impacts on the cube at about 0.8 s. The motion of the cube is basically consistent with the SPH-DCDEM result and experimental result by Canelas et al. (2016). It should be noticed that the experimental data are the average result of several tests [34].Figure 9
Dam break flow impacts on single cube with the proposed DEM-SPH method.
(a)
t = 0 s (b)
t = 0.9 s (c)
t = 1.6 s (d)
t = 2.2 sFigure 10
The displacement inx direction of single cube simulated with the DEM-SPH method and its comparison with the SPH-DCDEM and experimental results [34].Figure11 is the dam break flow impacting on three cubes simulated by the proposed DEM-SPH method, from which we can see that the cubes collapse under the impact of flow, and then the three cubes slide under the flow. Figure 12 is the photo comparison with the numerical and experimental results by Canelas et al. (2016) at two different times [34]. It can be found that the fluid-solid coupling phenomenon has little difference with the numerical and experimental results of Canelas et al. (2016). In the discrete element simulation of Canelas et al. (2016), the cube is constructed by spherical bonding element, so that the contact model of the cube element is different from that in this paper, which results in differences in simulation results.Figure 11
Dam break flow impacts on three cubes simulated with the DEM-SPH method.
(a)
t = 0 s (b)
t = 1.4 s (c)
t = 1.4 s (d)
t = 2.3 sFigure 12
The photo of dam break flow on three cubes with the proposed DEM-SPH method and its comparison with the DCDEM-SPH and experimental results: DEM-SPH simulation, experimental results, and DCDEM-SPH simulation from left to right [34].
(a)
t= 0.98 s (b)
t = 1.28 sIn the simulation of dam break flow on three cubes, the displacements of bottom particles inx direction and top cubes inx direction andz direction are compared with the corresponding simulation and experimental results, respectively, as shown in Figure 13. Similar to the previous example, the experimental data are the results of statistical average of several tests; i.e., each group of experimental results would fluctuate near the average value. In the simulation of three cubes, the motion of the cube is affected by many factors, such as the impact of water flow, the interaction with other cubes, and the friction with the bottom of the flume. The cube motion may be subject to more interference, which will produce randomness in the experiments. Therefore, the simulation result, including that by us and Canelas et al. (2016), and experiment result have little difference and the trends are the same. In general, the results in this paper are in good agreement with the corresponding simulation and experimental results.Figure 13
Displacements of the three cubes inx andz directions with the proposed DEM-SPH method and its comparison with the DCDEM-SPH and experimental results [34].
(a)
Displacement of bottom cube inx direction (b)
Displacement of top cube inx direction (c)
Displacement of top cube inz direction


## 4. Numerical Simulation of Rock Falling on Underwater Pipeline

The rocks are usually used to fix the pipeline for oil transportation under sea [35]. The impact of falling rocks on underwater pipeline can be simulated by the proposed DEM-SPH coupling method to analyze the impact force on the pipeline. The main simulation parameters are listed in Table 2. The rocks are initially set above the water to fall into the water and impact the pipeline.Table 2
Simulation parameters of impact of falling rocks on underwater pipeline.
Parameter Value Water tank: length × width × height 11.2 m × 6 m × 2 m Diameter of pipeline 0.508 m Cube length 0.15 m Density of cube 3000 kg/m3 Elastic modulus of cube 4.76 GPa Poisson’s ratio of cube 0.3 Number of cubes 600 Density of sea water 1035 kg/m3 Height of sea water 1.1 m Smoothing length 0.05 m Number of SPH particles 560 000The simulation is shown in Figure14. The rocks fall into the water, and thus the free surface of water fluctuates. Because of the coupling interactions between rocks and water, the force on pipeline includes the components from water and rocks. Since the water is disturbed, the force from water is different from that in static water. Figure 15 illustrates the forces on the pipeline from water and rocks. The result indicates that the water force caused by disturbance is larger than the rock force. On the other hand, the force on the pipeline in static water is the buoyance force. It can be calculated by Fb=ρgV=23.025kN according to the size parameters in Table 2. The buoyance force is shown as solid straight line, which is beneath but close to the force of water. This potentially proves the simulation accuracy of the proposed method.Figure 14
Impact of falling rocks on underwater pipeline simulated by the proposed DEM-SPH method.Figure 15
The impact force on the pipeline from the rocks and water.


## 5. Summary

In this paper, the discrete element method (DEM) of complex particles is realized by developing the dilated polyhedral element and its contact model. The hydrodynamic simulation is implemented by the weakly compressible SHP method. A simplified boundary repulsive force model is employed to regard the coupling between dilated polyhedral element and SPH particles as a boundary problem of SPH. Hence, the DEM-SPH coupling method for the coupling simulation of granular materials and fluid is established. This method saves a lot of computing resources and improves the computational efficiency because it does not need to construct solid boundary particles.The proposed DEM-SPH coupling method is used to simulate the dam break flow and the impact process on cube. The results are in good agreement with the corresponding numerical simulation and experimental results. Meanwhile, the method is applied to simulate the impact of falling rocks on underwater pipeline, while the forces from water and rocks are analyzed. Generally, the proposed method is validated in computational accuracy and can be applied for analyzing real engineering problems.



---

*Source: 4987801-2019-08-21.xml*