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.
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 bodies A and B in space, the Minkowski sum of them can be defined aswhere x and y are coordinate points belong to A and B, respectively. If A and B 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.
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:where is the normal stiffness of two dilated polyhedral elements in contact; is the normal overlap vector; is the normal relative velocity vector; 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:where is the vector of tangential elastic deformation; is the sliding friction coefficient between elements. is the tangential stiffness, and . Here the ratio is defined as 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 Figure 2. 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 of different contact types can be determined, and thus the contact force can be calculated by the Hertzian contact model [28].
(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:where W 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 . It must satisfy the normalization condition and compact condition. In SPH, the particle approximation of a field function at particle i can be expressed by the adjacent particles in its support domain.where N is the number of adjacent particles in the support domain; j denotes an adjacent particle of particle i; is the mass of adjacent particle j; is the density of adjacent particle j; . The divergence of the field function is only related to the smoothing function and can be written as follows:The smoothing length h 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:where . 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:where ρ0 is the initial density of the fluid; is generally equal to 7 [19]; is the speed of sound, which is usually equal to 10 times of the maximum particle velocity, i.e., .
Ignoring the viscous stress of fluid, the momentum conservation equation can be written as follows:where and are the pressures of particle i and particle j, respectively. In order to eliminate the nonphysical oscillation in SPH simulations, the artificial viscous term by Monaghan is adopted:where and are standard constants and generally range from 0 to 1.0; in addition, there arewhere and . 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 aswhereHence, the final form of momentum equation yielded by particle approximation of SPH can be obtained, written asThe explicit leap-frog (LF) method is used to integrate the above discretized equations [29]. The time step of integration is determined as follows:In 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:where is the boundary repulsive force; is the unit normal vector of the wall boundary; is the pressure correction term; is the minimum distance between the SPH particle and boundary; is the repulsive force function, written aswhere is the normalized distance, and ; is related to the smoothing length and sound speed, and .
(a) Traditional particle boundary
(b) Geometric boundary
In order to balance the pressure in different depths of fluid, the pressure correction term are employed here, written aswhere is the initial height of fluid; is the height of SPH particle. Here note that the water level is zero plane, i.e., , 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
Figure 4 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 5 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 is employed, where is the gravity acceleration and a is the initial length of water. The simulation by the proposed SPH method is in good agreement with other approaches.
(a) t = 0.1s
(b) t = 0.3s
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 Figure 4. 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.
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.
(a) Single cube
(b) Three cubes
Figure 9 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 along x 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].
(a) t = 0 s
(b) t = 0.9 s
(c) t = 1.6 s
(d) t = 2.2 s
Figure 11 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.
(a) t = 0 s
(b) t = 1.4 s
(c) t = 1.4 s
(d) t = 2.3 s
(a) t= 0.98 s
(b) t = 1.28 s
In the simulation of dam break flow on three cubes, the displacements of bottom particles in x direction and top cubes in x direction and z 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.
(a) Displacement of bottom cube in x direction
(b) Displacement of top cube in x direction
(c) Displacement of top cube in z 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.
The simulation is shown in Figure 14. 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 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.
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.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This study is financially supported by the National Key Research and Development Program of China (Grants numbers 2018YFA0605902, 2016YFC1401505, and 2016YFC1402705) and the National Natural Science Foundation of China (Grants numbers 11572067 and 11772085).