A contact detection algorithm for triangle boundary in GPU-based DEM and its application in a large-scale landslide

https://doi.org/10.1016/j.compgeo.2021.104371Get rights and content

Abstract

Triangle boundaries are crucial geometric objects in the discrete element method (DEM). In this study, an improvement for the uniform grid method is proposed and implemented on the GPU. The improved algorithm includes the spatial discretization of triangles in which triangles are classified by their size and motion. The large and relative static triangles are separated into grid cells such that it relaxes the restriction of cell size in the uniform grid method. To verify the accuracy and efficiency of the algorithm, the improved grid method is compared with the normal uniform grid method and the linear bounding volume hierarchy (LBVH) method using a simulation of free-falling grains in a mixer. The results show that the improved grid method does significantly increase the simulation efficiency compared with the normal uniform grid method. Then, taking a large-scale landslide as an example, the improved grid method and the LBVH method are used to simulate the failure process of the landslide. The results indicate that the improved grid algorithm has higher efficiency for simulations of large-scale landslides than the LBVH method but at the cost of more memory usage.

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 unand shear increment Δus, the resultant force FAC for particle A in response to particle B can be decomposed to the normal force FCn and the shear force FCs which are governed by normal stiffness (kn) and the shear stiffness (ks),FAC=FCn+FCsFCn=knunnFCs={FCs}update+ksΔus

where, kn and ks 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)

Cited by (14)

  • DEM contact model for spherical and polyhedral particles based on energy conservation

    2023, Computers and Geotechnics
    Citation 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.

View all citing articles on Scopus
View full text