In geotechnical engineering and rock mechanics, capturing the kinematic response of a particulate medium to loading is fundamental to modeling naturally deposited sediments, tunneling, and rock slope stability, among others. Fundamentally, in discontinuous media such as jointed rock masses, the displacement is governed by pre-existing discontinuities (bedding planes, faults, and joints) and continuum-based methods generally are not able to describe the highly localized nature of displacements along these structural features. Particulate-based numerical models, such as the Discrete Element Method (DEM) and Discontinuous Deformation Analysis (DDA), have been developed to model particulate systems, including fractured rock masses. These approaches are able to directly model the shape of individual particles within a particulate system, and can intrinsically accommodate large displacements associated with discontinuities at particle and block margins.

In this research, we first present a formulation of a new method to calculate the contact point between polyhedral particles that can be used to complement different particulate methods. Contact detection and contact geometry are among the most important steps in DEM and DDA simulations. Accurately representing the contact between two particles is crucial, and when modeling fractured rock using polyhedral particles, the accuracy of the contact point calculation is essential for obtaining realistic and reliable simulation results. We present a new algorithm for accurately calculating the contact point between two colliding polyhedral particles that uses the topology of the particles to assert their interaction with the plane of contact and not directly with each other. The new algorithm provides improved performance in terms of global stability of DEM models by mitigating numerically induced instability associated with errors and sporadic movement in the contact point calculation. We then present a comparative study of contact detection algorithms routinely used for simulating convex polyhedra in particulate methods. We critically assess, in terms of accuracy and computational efficiency, the advantages and limitations of each algorithm tested in the scope of this research. These algorithms are implemented within the same open-source software framework to allow for an objective assessment of their performance. In this case, DEM is used as the particulate solver; however, the characteristics of the selected contact detection algorithms are independent of this choice. The verification examples presented allow for informed decisions on which class of contact detection algorithm may be most appropriate for specific types of particulate simulations involving convex polyhedra. Lastly, we compare solutions obtained using Limit-Equilibrium (LE) and the Discrete Element Method (DEM), both in 2D and 3D, for analyzing rock slope stability and runout. Specifically, we illustrate the importance of kinematics in modeling of rock slopes and describing the progressive nature by which they displace and ultimately fail. While 3D LE methods provide a measure of the factor of safety against failure, the failure surface and mode are assumed, and the rock mass is typically represented by vertical slices (2D) or columns (3D) in the analysis. Thus, the kinematic response of the rock mass is artificially constrained since the natural structure within the rock is largely ignored, the failure mode is assumed a priori, and thus the factor of safety calculated using LE approaches may not necessarily describe the governing response. In contrast, in DEM a specific mode of failure is not assumed, since natural discontinuities, joints, shears, and fractures, as observed in the field can be used to create a more realistic representation of the rock mass. We also show that with an increasing number of rock blocks in the model (tighter spacing of the joints) without changing the intrinsic material properties, the rock mass is less stable. This has implications for rock slope stability evaluations, as rock that is more fractured will be less kinematically constrained and require more mechanical strength to remain stable. Additionally, during rock slide initiation and propagation, the rock within the sliding mass may fracture and disintegrate, such that it becomes less constrained as it deforms. The outcome is a progressive rock slope failure and accelerating displacements as the rock blocks within the sliding mass become more fractured. Thus, we show that in analyzing the response of fractured rock, it is important to not only consider the mechanical properties required for stability, but also the amount of displacement the sliding mass may undergo, as breakup of the rock due to displacement will have a destabilizing effect.