DEM modelling of 3D polyhedra with applications to gabion rockfall barriers

Rock particle shape plays a crucial role in shear resistance and energy consumption during transient loads. The dynamics of such granular materials are complex and cannot be properly described using closed-form solutions when the problem involves more than a few particles. Thus, for the sake of computational efficiency, it is common practice to implement simplified numerical models that involve a limited number of particle interactions. In this study, a novel approach is used to capture realistic particle shapes while maintaining a relatively high simulation efficiency. The geometry of the particles is determined by a Delaunay triangulation which operates on a set of vertices and returns the corresponding network of facets and grid connections. Inertial and material properties are assigned to the rock prototype which are representative of realistic gravel particles. The algorithm is validated by performing a series of numerical simulations for various particle configurations, demonstrating that mass and momentum are conserved. A potential application of this work is related to rockfall barriers and their response to rigid boulder impacts. This innovative model, based on the discrete element method, is shown to be capable of simulating rock particles with realistic shapes and complex physical interactions.


Introduction
The modelling of granular materials involves characterizing interactions between a significant number of particles. During collisions, mechanical contact effects between particles play a key role in the observable physics and are a strong function of particle geometry. A common tool for numerically simulating particle interactions is the discrete element method (DEM) [1]. One such instance of this approach for the simulation of particle collisions is presented in [2], where the DEM is used to capture damage behaviours observed during rockfalls with clump-simulated boulders.
An and Tannant [3] developed a contact model for the fully dynamic simulation of rock impacts and rebound processes. The model allows for the controlled removal of kinetic energy that occurs during an impactrebound event in a DEM simulation. The parameters needed for the contact model are the normal contact stiffness, the transition force, and an exponent for the power law which determines the contact force upon unloading. The transition force is the normal contact force at which the material transitions from a linearelastic to an elastoplastic stress-strain relationship. The model has been successful in simulating rock impact processes that result in kinetic energy loss. The new contact model was tested in PFC 2D for a simple DEM simulation of a rock falling vertically onto a horizontal surface. The rock was modelled as a sphere with radius *Corresponding author: yma@ltu.edu 0.05 m and density 2600 kg/m 3 . The Spherical Discrete Element Code (SDEC) was used in [4] to investigate boulder impact events on reinforced soil embankments with the goal of developing more specific regulations for the design of rockfall protection structures. From the numerical experiments, it was found that the maximum impact force from the boulder is related to the kinetic energy by a power law with an exponent less than one.
In a number of these models, homogenous spherical bodies are used to simulate rock particles. However, granular materials are more likely to possess a polyhedral geometry with nonhomogeneous characteristics. Herein we present the Polyhedral Reinforced Interior Shell Model (PRISM), which facilitates the generation of convex angular bodies from a Delaunay triangulation acting on a point cloud of vertices. This approach differs from other polyhedral particle models because it allows the user to assign internal stiffness properties that are independent of the surface elasticity. Furthermore, the interior stiffness can vary in a controlled, yet anisotropic manner by implementing a relatively small amount of cylindrical particles that adjoin the vertices to the centroid of the polyhedron. In other models, such as LMGC90, polyhedral shape classes exist, but meshes are often necessary in assigning various material properties to irregular structures. Yade [5] also has a polyhedron shape class, which is approximately ten times faster than the presented PRISM model but is not compatible with several other classes within the Yade domain, such as wire meshes made from grid connections. Furthermore, the polyhedral particles are modelled as ideally rigid bodies. The PRISM model is advantageous due to its flexibility in allocating stiffness properties proportionally to the mass of the particle components, thus allowing for larger stable time-steps. This feature may be beneficial in numerical experiments, such as the triaxial tests performed by Lee et al. [6] that require large numbers of polyhedral particles. Overall, this proposed model offers a simplified solution for generating bodies which possess physical inertial and material properties equivalent to their solid counterparts. The process for constructing these particles is outlined below, and a validation of the model is presented by comparing simulations with analytical results. To demonstrate the suitability of the PRISM method for engineering applications, a simulation is presented of a gabion barrier subjected to a boulder impact.

Yade
The proposed rock particle model is implemented in Yade. The utility of the open-source DEM software is rooted in its engines which execute algorithms depending on the geometry, material properties, and constitutive laws of the particles. A sweep-and-prune method is used for collision detection. The rigid-body kinematics of the model are decomposed into their translational and rotational counterparts. The former is addressed by using a "leapfrog" scheme, and the latter requires an extension of this method, outlined in [7], which numerically solves Euler's equation: for rigid bodies in a non-inertial reference frame. represents the net torque, is the diagonalized inertia tensor, and is the angular velocity. This approach is warranted on account of the irregular convex polyhedral geometry of the rock particles and is implemented by selecting the "aspherical particle" option in Yade, which assigns the proper inertial values to a particle.

Constructing a polyhedron
The physical and geometric structure of the rock particles are composed of grid nodes, grid connections, and pfacets (i.e., deformable particle facets in Yade). The grid nodes are placed at the vertices of the polyhedron and are made from a cohesive frictional material which are assigned a Young's modulus, normal cohesion, shear cohesion, angle of friction, bending yield stress, and density. The grid connections form bonds between grid nodes and are represented by a cylindrical geometry with physical properties derived from the frictional material class (i.e., frictMat) in Yade. The pfacets, also made from the frictMat material, form a mosaic of tessellations about the triangulated faces. The construction process is illustrated in Figure 1. The topologies of these three particle classes (grid node, grid connection, and pfacet) are assigned by a Delaunay triangulation. This algorithm acts on a set of points and returns a collection of partitioned surfaces (triangles) and volumes (tetrahedrons) along with their respective vertices located by = ( , , ) , = ( , , ) , = ( , , ) , and = ( , , ) . Grid connections are then formed between neighbouring nodes based on the Delaunay triangulation. It is assumed that the polyhedron is homogeneous, so the centroid is found by performing a volumetrically weighted average over the tetrahedrons: in which the vector = ( , , ) locates the centre of mass. and = ( , , ) are the volume and centroid of the i th tetrahedron, respectively. The vector is defined in terms of the vertices: and the volume of each tetrahedron is calculated from, A grid node is then placed at the location . The inertia tensor of the solid polyhedron must be calculated to determine the principal axes of the rigid body. It will be useful to introduce the relative distances between the centres of mass of the tetrahedrons and the polyhedron ( ′ , ′ , ′ ) = ( − , − , − ) . The inertia tensor is given by: in which the entries of the matrix are: ( ′ 2 + ′ ′ + ′ 2 + ′ ′ + ′ ′ + ′ 2 + ′ ′ + ′ ′ + ′ ′ + ′ 2 + ′ 2 + ′ ′ + ′ 2 + ′ ′ + ′ ′ + ′ 2 + ′ ′ + ′ ′ + ′ ′ + ′ 2 )/10 (7) Similarly, for the products of inertia we have: , = (2 ′ ′ + ′ ′ + ′ ′ + ′ ′ + ′ ′ + 2 ′ ′ + ′ ′ + ′ ′ + ′ ′ + ′ ′ + 2 ′ ′ + ′ ′ + ′ ′ + ′ ′ + ′ ′ + 2 ′ ′ )/20 (8) The mass of the tetrahedron is denoted by . The terms associated with the parallel axis theorem are defined as , 2 = ′ 2 + ′ 2 and , = ′ ′ such that where is a principal moment of inertia and the eigenvector = ( , , ) is the corresponding principal axis. The aspherical grid node, located at the centroid, is rotated into the coordinate system formed by the principal axes, given by the following matrix: Note that Yade modifies the orientation of an object by the implementation of quaternions. Equations for converting rotation matrices to quaternions are readily available in existing literature (e.g. [9]).

Model validation
The

Torque test
A cube, generated with the PRISM algorithm, is momentarily subjected to a couple force, resulting in a uniform angular motion about its centroid. Note that the relevant physics take place in the -plane, but the simulations involve three-dimensional cubes. Theoretically, the final angular speed of the object is given by = √2 / , where is the side length of the square, is the applied force, and is the time over which the torque acts on the body. See Figure 2 for a diagram of the setup. The body has a mass of = 1000 kg and a side length of = 1 m, corresponding to a moment of inertia of = 166.67 kg-m 2 . A force of = 10 6 N is applied to the cube for a time of = 10 −4 s, yielding a theoretical angular speed of = 0.84852 s -1 in the counterclockwise direction. The numerical result for the same input values corresponds to an angular speed of = 0.84848 s -1 , yielding a relative error of 0.005%. Note that the numerical model is deformable, so a Young's modulus of = 10 12 Pa is used with a grid node radius of = 0.08 m. The time step is Δ = 10 −9 s.

Collision test
Conservation of linear and angular momentum are verified simultaneously through a numerical simulation involving the collision of two cubes which are vertically offset from one another, as illustrated in Figure 3. The initial velocity of block "a" is purely in the horizontal direction, with a magnitude of = 10 m/s. Cube "b" is initially at rest and has its centre of mass offset by a distance of /4 in the -direction. Both bodies have a mass of = 1000 kg and side length = 1 m. In this experiment, conservation of linear momentum requires that = + . After the interaction, cubes "a" and "b" obtain final horizontal velocities of = 1.5426 m/s and = 8.4573 m/s, respectively, resulting in negligible error.

Fig. 3 Definition sketch of initial setup for collision test.
Cube "a" on left and cube "b" on right.
Conservation of angular momentum is considered about the stationary point , requiring the following condition to be satisfied: In equation (12) and denote the final angular velocities, and and denote the moment arms originating from point and extending to the centroids of cubes "a" and "b", respectively. The problem is 2D, so we only consider the component of angular momentum orthogonal to the -plane. The sum of the terms on the right-hand side yields a value of 0.00086 kg-m 2 /s. Note that | | = 2084.972 kg-m 2 /s, and the moment of inertia of cube "b" is = 166.67 kg-m 2 . Therefore, the error in angular momentum (0.00086 kgm 2 /s) would create a difference of 0.00052% in the angular velocity of square "b". From these simulations, it is evident that the model is capable of producing accurate results which compare well with theoretical predictions. Note that in this simulation, we use a time step of Δ = 10 −10 s, a Young's modulus of = 10 12 Pa, and a grid node radius of = 0.08 m. The authors performed similar numerical simulations as those described above, with larger time steps (Δ = 10 −5 s) that produced stable results but are not presented here for the sake of brevity.

Gabion barrier application
Results are presented for a numerical simulation in which a 75 kg boulder with an initial horizontal velocity of 10 m/s impacts a gabion barrier with dimensions 1.0 m × 0.5 m × 1.0 m. Gravity is set at 9.81 m/s 2 . The boulder and gravel content of the gabion are generated using the PRISM algorithm described previously. The rock particles have a density of 2500 kg/m 3 and a Young's modulus of 10 11 Pa. The grid connections belonging to the polyhedra have a radius of 0.02 m. The flexible wire mesh basket has a density of 7900 kg/m 3 and a Young's modulus of 2 × 10 11 Pa with a grid connection diameter of 0.01 m. The ground is made from pfacet elements with a density of 2500 kg/m 3 , a Young's modulus of 5×10 11 Pa, and a grid connection radius of 0.05 m. The angle of friction on all contacts is set at 30 o . A time step of Δ = 10 −7 s was used. The initial conditions are presented in Figure 4a. The gabion barrier is presented in Figure 4b and 4c, illustrating the state of the structure 0.05 s after the impact occurs. The maximum horizontal deformation was measured to be approximately 0.154 m, resulting in apparent plastic behaviour of the basket. A quantitative comparison was not performed, but the numerical results appear to be qualitatively consistent with postboulder-impact gabion photos inspected by the authors.

Conclusion
A novel method for generating convex polyhedral rock particles with Delaunay triangulated topologies has been proposed and verified in Yade through simple collision tests. This method is computationally advantageous, compared for example with "clump" methods, as it requires significantly fewer particles to accurately model the inertial characteristics of a given rock particle. The internal grid-connection structure allows for the reallocation of greater stiffness values to the massive centroid, thus permitting a larger time-step. The presented gabion-boulder-impact simulation suggests that this model is suitable to engineering applications involving rockfall protection structures.