A contact detection algorithm for triangle boundary in GPU-based DEM and its application in a large-scale landslide
Introduction
Numerical methods are widely used in engineering and science for the simulations of physical and mechanical processes. The discrete element method (DEM), firstly proposed by Cundall and Strack (Cundall and Strack 1979), is one of the representative numerical methods that can simulate the meso-mechanical behaviors, dynamics, and failure processes of geomaterials and other granular materials. DEM has proven effective to understand the discontinuous behaviors of granular materials, such as particle breakage (Eliáš, 2014, Liu et al., 2020b), rock failure (Liu et al. 2020a), direct shear test (Yang et al. 2017), grinding mill (Govender et al. 2015b), landslides (Zhang et al. 2020). In recent years, to simulate increasingly complex physical processes, coupling algorithms between DEM and other numerical methods, such as SPH (Tan and Chen, 2017, Wu et al., 2016, Xu et al., 2019), LBM (Ding and Xu, 2018, Li et al., 2018), FEM (Guo and Zhao, 2016, Shen et al., 2019, Zheng et al., 2018) have been developed in which the DEM is used to simulate the discontinuous part. High performance computing architecture and efficient parallel algorithms are needed to perform a large-scale and robust analysis for the physical processes of granular systems.
In DEM, the materials are divided into a set of particles (or blocks) that interact with each other. Using an explicit scheme in the time domain, the particles’ velocities and accelerations are calculated based on Newton’s second law. Collision detection, contact forces computation, and numerical integration are the main processes within each time step. In these processes, the collision detection is often the largest time-consuming process and may account for over 70% of the total computing time (J.R. and R. 1995). To avoid the time complexity of order O(N2) over N particles in the computed domain, and to make the contact detection suitable for complex particle geometry, bounding volumes for particles are used that split the contact detection into broad phase and narrow phase.
There is a wide range of algorithms for broad phase detection to find the collisions between bounding volumes. Munjiza (Munjiza et al. 2006) proposed a CPU based spatial partitioning algorithm that uses uniform grids and applied simple improvements to increase the efficiency. The sweep and pure method was proposed by Jonathan(Cohen et al. 1995), which is an efficient collision detection method for bounding volumes with varied sizes and detects the inversion of the bounding geometries using an insertion sort algorithm. Hierarchical methods (Hubbard, 1996, Klosowski et al., 1998) separate geometries in a tree structure and traverses them in a bottom-up or top-down manner to perform collision detection. In most circumstances, the uniform grid method has higher efficiency simulating evenly sized geometries and is simpler to implement in parallel programming while other methods can better adapt to complex computing domains, which is specifically compared by Retief (Lubbe et al. 2020).
In recent years, with the advances in Graphics Processing Unit (GPU) computing, implementations for parallel algorithms like the uniform grid method (Anderson et al., 2008, Govender et al., 2015b) and the LBVH method (Garanzha et al., 2011, Karras and Aila, 2013, Lubbe et al., 2020) for collision detection have gradually been developed and proven efficient in large scale simulations. In the uniform grid method, the cell size should not be smaller than the bounding volume of the largest particle to ensure that possible collisions are not neglected. A cell size that is much larger than the size of a particle being detected results in a significant increase in computational cost. Therefore, for a particle system with various particle sizes (especially when the particle size differs greatly), the computational efficiency is usually low in the uniform grid method. Unfortunately, in many natural systems, particles with varied sizes cannot be avoided, therefore improved grid methods must be developed to solve this problem. A discrete function representation (DFR) method was developed by Williams (Williams and O’Connor 1999) to represent a large circle or sphere by small particles. Moreover, clumping a certain number of spheres is a popular method used in DEM to describe arbitrarily shaped particles (Xu et al. 2016), and is an approach that is also effective in decreasing the computational cost in the grid method. Furthermore, multi-level grid methods (Huang et al. 2020) and hierarchical grid methods (Weinhart et al. 2020) are developed for particles with varied sizes and exhibits high efficiency. However, complex boundaries are not included in the two methods above.
When using DEM to simulate the physical processes of granular materials in industry or nature, the simulation of the interactions between particles and complex boundaries is common. For example, DEM is widely used for the simulation of landslide dynamics, and the terrain is usually described by a series of triangles (Wu et al., 2017, Wu et al., 2018). Due to the complexity of the terrain, the quality of triangle meshes is often not refined enough to represent the terrain accurately. In this case, a large number of small triangles will be contained in some cells because the cell size in the grid method is required to be larger than the biggest triangle, which will increase the computational cost during the broad phase traversal. Furthermore, dynamic boundaries such as in the milling process of granular materials in a mixer, a dynamic grid algorithm may be needed to increase the computing efficiency.
In this study, to improve the detection efficiency of the contact between particles and complex boundaries, a spatial discretization for triangles in the uniform grid method and the corresponding detection algorithm are provided, and the algorithm based on the contact area between the triangle and sphere in the narrow phase detection is used. Then, the algorithms for the translatable and rotatable grid methods are provided and implemented on DEM based on GPU. To verify the algorithms, two benchmarks are used, one is a sphere moving on a slope and the other is the milling process of granular material in a mixer. In the end, the failure process of a large-scale landslide is efficiently simulated based on the improved algorithms.
Section snippets
Discrete element method (DEM)
In this study, the cohesive fracture model (CFM) (Liu et al. 2020a) is used for the DEM contact law. In each step, when two particles A and B are in contact, with the penetration depth of and shear increment , the resultant force for particle A in response to particle B can be decomposed to the normal force and the shear force which are governed by normal stiffness () and the shear stiffness (),
where, and between
a sphere sliding along a slope
To validate the contact algorithm between the spherical particle and the triangular boundary, a sphere sliding on a slope composed by a series of triangles (Fig. 6) is firstly used in this section. The inclined angle of the slope is 31° and three slope models S1, S2, S3 (Fig. 6) are used in the simulation. During the simulation, the gravity of 9.8 m/s2 in vertical direction (in the negative direction of z axis) is applied. The mechanical parameters of the sphere and the slope after it has been
Application for the simulation of a large-scale landslide
Generally, the terrain of a slope in nature is very complex, and it will have a great influence on the dynamic process of the landslide. So, it will be important for a numerical simulation of the landslide dynamics to precisely represent the geometry of the terrain. The triangular meshes are usually used for the representation of the geometry of a terrain. In this section, as an application for the large sale and complex example, the algorithms provided in this study will be used to the
Conclusion
The interaction between particles and complex boundaries are common in DEM simulation, and the corresponding contact detection will greatly influence the computation efficiency. To solve this, a GPU accelerated improved contact detection algorithm based on the uniform grid method for complex triangle boundary is proposed. In the improved grid method, the large-sized triangles are separated into relative grid cells so that it simulates geometries of various sizes more efficiently. The proposed
Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Acknowledgements
This work was supported by the projects of “Natural Science Foundation of China, China (41941019, 52079067)”, “Research Fund Program of the State Key Laboratory of Hydroscience and Engineering (2020-KY-04)”
References: (36)
- et al.
General purpose molecular dynamics simulations fully implemented on graphics processing units
Journal of Computational Physics
(2008) - et al.
The 2000 Yigong landslide (Tibetan Plateau), rockslide-dammed lake and outburst flood: Review, remote sensing analysis, and process modelling
Geomorphology
(2015) - et al.
Study on the multiphase fluid-solid interaction in granular materials based on an LBM-DEM coupled method
Powder Technology, Elsevier B.V.
(2018) Simulation of railway ballast using crushable polyhedral particles
Powder Technology
(2014)- et al.
Discrete element simulation of mill charge in 3D using the BLAZE-DEM GPU framework
Minerals Engineering
(2015) - et al.
Parallel hierarchical multiscale modelling of hydro-mechanical problems for saturated granular soils
Computer Methods in Applied Mechanics and Engineering, Elsevier Ltd
(2016) - et al.
An improved contact detection algorithm for bonded particles based on multi-level grid and bounding box in DEM simulation
Powder Technology, Elsevier B.V.
(2020) - et al.
Coupled DEM-LBM simulation of saturated flow velocity characteristics in column leaching
Minerals Engineering
(2018) - et al.
A cohesive fracture model for discrete element method based on polyhedral blocks
Powder Technology, Elsevier B.V.
(2020) - et al.
Study on the particle breakage of ballast based on a GPU accelerated discrete element method
Geoscience Frontiers, Elsevier Ltd
(2020)
Analysis of parallel spatial partitioning algorithms for GPU based DEM
Computers and Geotechnics, Elsevier
A super-large landslide in Tibet in 2000: Background, occurrence, disaster, and origin
Geomorphology
Application of open source FEM and DEM simulations for dynamic belt deflection modelling
Powder Technology
A hybrid DEM-SPH model for deformable landslide and its generated surge waves
Advances in Water Resources, Elsevier Ltd
Fast, flexible particle simulations — An introduction to MercuryDPM
Computer Physics Communications, North-Holland
Morpho-sedimentary and stratigraphic characteristics of the 2000 Yigong River landslide dam outburst flood deposits, eastern Tibetan Plateau
Geomorphology
Assessing the impacts of a large slope failure using 3DEC: The Chiu-fen-erh-shan residual slope
Computers and Geotechnics
Post-failure simulations of a large slope failure using 3DEC: The Hsien-du-shan slope
Engineering Geology
Cited by (14)
Numerical study of the effects of loading method on mixing of two kinds of pebbles in HTGR: A GPU-DEM simulation
2023, Progress in Nuclear EnergyModelling large-scale landslide using a GPU-accelerated 3D MPM with an efficient terrain contact algorithm
2023, Computers and GeotechnicsDEM contact model for spherical and polyhedral particles based on energy conservation
2023, Computers and GeotechnicsCitation Excerpt :The contact relationship between intersection circle and the target triangle can be obtained easily. All the possible intersection types between the target triangle and the intersection circle in the intersection plane can be classified into 11 cases based on the number of edges that outside of the intersection circle (Ne) and the number of vertex that outside of the intersection circle (Nv) (Zhou et al. 2021b). The contact point can be calculated based on the area weighted method.
A Coupled FEM-MPM GPU-based algorithm and applications in geomechanics
2022, Computers and Geotechnics