A surface mesh represented discrete element method (SMR-DEM) for particles of arbitrary shape

A surface mesh represented discrete element method (SMR-DEM) for granular systems with arbitrarily shaped particles is presented. The particle surfaces are approximated using contact nodes obtained from surface mesh. A hybrid contact method which combines the beneﬁts of the sphere-to-sphere and shpere-to-surface approaches is proposed for contact detection and force computation. The simple formulation and implementation render SMR-DEM suitable for three-dimensional simulations. Furthermore, GPU parallelization is employed to achieve higher e (cid:14) - ciency. Several numerical examples are presented to show the performance of SMR-DEM. It is found that on the particle level the method is accurate and convergent, while on the system level SMR-DEM can e (cid:11) ectively model particle assemblies of various regular and complex irregular shapes.


INTRODUCTION
Granular systems consisting of macroscopic particles widely present in nature, industries, and engineering.The discrete element method (DEM), originally proposed by Cundall and Strack [1], has been extensively applied to the simulations of granular assemblies.By simply considering inter-particle contact and Newton's equations of motion, DEM can provide detailed information of particle systems such as particle trajectories, flow dynamics, shear stress, contact forces, and energy dissipation.It has achieved notable success in a wide range of applications, e.g.granular flows [2][3][4], natural hazards [5,6], industrial processes [7][8][9][10], and geotechnical problems [11][12][13][14], among others.
In most DEM analysis, granular materials are modelled using spherical particles, which allow fast and straightforward contact detection and force computation.However, the majority of particles involved in natural and engineering processes are non-spherical and irregular, some examples are sand grains, rocks, ores, and pharmaceutical tablets.The shape of particles has a significant influence on macroscopic properties of granular assemblies, e.g.capacity of packing [15,16], resistance to shear and dilation under shear [17], ability to form arches [4,18], and particle segregation [19].Spherical representation of granular particles seems inappropriate for problems where particle shape plays a major role.Thus, accurate shape descriptions are required.
The main challenges in developing non-spherical DEM are the shape representation and contact detection.To date, many approaches have been developed to model non-spherical particles in DEM.An incomplete list includes: polygon/polyhedron [20][21][22], multi-spheres [23,24], superquadrics [10,25,26], probability-based method [27], discrete function representation [28,29], and level-set DEM [30,31].Detailed reviews on these methods can be found in Lu et al. [32] and Zhong et al [33], where accuracy and efficiency were extensively discussed.Many methods are being developed in parallel but no method is to be favoured in terms of accuracy and efficiency.For instance, the multi-sphere method is probably the most popular one owing to its capability of modelling particles with arbitrary shape.However, it was reported that the multi-sphere method has convergence issues, i.e. even with a large number of glued spheres, the accuracy of the method is not always guaranteed [17,34,35].
Apart from accuracy and efficiency, we believe that versatility and simplicity are also important criteria to assess methods for non-spherical particles.In fact, many methods are developed for certain particle shapes, i.e. [36] for ellipsoids, [37] for tablets, [38] for cubes, and [39,40] for cylinders.The polygon/polyhedron method [20,21] is not suitable for particles with smoothly curved surfaces.On the other hand, the multi-spheres method, although highly valued for its versatility, has difficulties in modelling particles with flat surfaces.The Gilbert-Johanson-Keethi algorithm-based DEM (GJK-DEM) can model convex shapes [17] and by gluing multiple particles it is able to model non-convex shapes [41].Its application to granular assemblies of arbitrary shape is however not yet reported.Generally, modelling particles of arbitrary shape is still a challenge in DEM analysis.Simplicity in implementation and data preparation should also be taken into consideration.Among the proposed methods, some need dedicated algorithms to represent non-spherical particles and detect contacts.Examples are the level-set function and the Gilbert-Johanson-Keethi algorithm.Others may require special treatment for data preparation, e.g. the multi-spheres method needs special algorithms [42][43][44][45] to fill irregular particles to achieve better shape representation with fewer glued spheres.
The modelling of arbitrarily shaped particles still remains a subject of intensive research.An ideal method should possess the following features: (1) it is efficient, or at least can be conveniently and efficiently parallelized; (2) its formulations and implementation are simple in threedimensions; (3) it should be stable and accurate in contact detection and force calculation, or at least have convergent behaviour; (4) it can model particles of any shape, e.g.regular, polyhedral, and arbitrarily irregular, using the same algorithm; (5) the data preparation should be simple.In this paper, we present an alternative surface mesh represented discrete element method (SMR-DEM) for three-dimensional analysis of granular systems.The presented SMR-DEM is highly versatile and can model any particle shapes, whether they are regular, polyhedral, convex, non-convex, or arbitrarily irregular.A hybrid contact method is proposed for contact detection and force computation.Except the surface mesh of grains, the SMR-DEM requires no additional algorithms for data preparation.Furthermore, high-performance computing of the proposed SMR-DEM framework is achieved using GPU acceleration.The features and performance of the SMR-DEM is demonstrated using several numerical examples involving regular and arbitrarily irregular particles.
This paper is organized as follows.In Section 2, we present the framework of SMR-DEM, including surface discretization, motion equations, and contact computation.Numerical implementation and GPU parallelization are presented in Section 3. Five numerical examples are performed to show the performance of SMR-DEM in Section 4. A discussion on numerical efficiency and limitations is provided in Section 5. Finally, some conclusions are given in Section 6.

Particle surface discretization
The time-driven soft approach is employed in SMR-DEM, where the overlap between particles, regarded as the particle deformation of the neighbouring contact surfaces, is used to calculate interaction forces.As shown in Figure 1, in SMR-DEM, particle surfaces are discretized using contact nodes, which can be conveniently determined in preprocessing as the vertices of the surface mesh.Note that the surface mesh is only used to determine the positions of contact nodes during preprocessing, it is not required in contact detection and force computation.There are mainly two benefits of describing particles using surface mesh representation.First, triangular mesh is the most common way to describe three-dimensional surfaces.Regardless of how the particle shape is obtained, either from X-ray CT [31,46], scanning electron microscope (SEM) [47], 3D laser ranging [48], the Structure-from-Motion technique [49], or computer-aided design (CAD), mostly the digitalized particles are in the form of triangular surface meshes.By directly utilizing the triangular meshes, additional algorithms for describing particle shapes, e.g. the filling algorithms in the multi-sphere method [42][43][44][45], super-quadric equations [28], or level-set functions [30,31], are no longer needed.This greatly simplifies the formulation, implementation, and data preparation.Second, with triangular surface mesh, arbitrarily complex shapes can be represented without difficulty.All regular and arbitrarily irregular particle shapes can be described using the same unified approach.Consequently, the particle representation is simple and SMR-DEM has a high degree of versatility.Additionally, various techniques and tools for triangular mesh generation, refinement, and optimization are readily available if improvement of initial input surface meshes is desired.
In SMR-DEM, all contact detection and force calculation are based on contact nodes, which carry the unit normals of the vicinity of particle surfaces.Let A i , i = 1, ..., N be the N triangles on the surface of particle P. For each triangle, the outward unit normal n i can be defined as , with a i , b i and c i being the vertices of triangle A i in a counter-clockwise order.Consequently, the unnormalized vector at contact node j is obtained by averaging the normal vectors of associated triangles where N j is the number of triangles associated with the contact node.The unit normal vector at contact node j is obtained as At each contact node, a nodal shape factor j is defined as which quantifies the maximum difference between the nodal unit normal and the unit normals of the associated triangles.A larger nodal shape factor indicates greater change in the local surface geometry.If the nodal shape factor is larger than a threshold λ, the contact node is labelled as a sharp corner.Its contact treatment is different from non-corner nodes and the details will be given in Section 2.3.
In SMR-DEM, the centroid and moment of inertia of a particle can be computed directly using the surface mesh.The mass M of a host particle can be expressed as where ρ(x) is the material density, which is assumed to be constant in one particle.Ω is the spatial domain bounded by the particle surface S .Using the divergence theorem, the mass is given by where N is the number of triangles on the surface of particle.As x • n i is constant across A i , and | ni | equals two times the area of A i , the mass of a host particle can be computed as where a i is the position of the first vertex of triangle A i .Similarly, the centroid X of a particle can be calculated from where V is the particle volume.Applying the divergence theorem once more, the above equation can be rewritten as where l is the index of spatial coordinates, X l , x l , n l denote the l-th component of the vectors X, x, and n, respectively.n il is the l-th component of the normal vector of the i-th triangle.In the equation above the Einstein summation convention is not applied.With the midpoint quadrature rule [50], the integration over triangles can be transferred to the vertices.The final expression of the coordinate of centroid is written as where a il , b il and c il are the l-th coordinate component of vertices a i , b i and c i , respectively.It should be noted that the expressions of mass and centroid are valid for arbitrarily shaped particles even if they are non-convex and not simply connected.
Once the centroid is determined, the host particle can be partitioned into N tetrahedra by connecting the centroid with the surface triangles.It is easy to calculate the moment of inertia for each tetrahedron with respect to the centroid of the host particle.The total moment of inertia of the host particle is thus the summation of contributions from all the constituent tetrahedra where Ω i is the integral domain of the i-th tetrahedron, and r = x − X is the relative position vector between a point x in the i-th tetrahedron and the center of mass X of the particle.Note that this method for computing moment of inertia is not universally applicable for irregular particles.For example, if the spatial domain of a particle is not simply connected, or if the centroid locates outside of the particle, the method above gives incorrect results.Under those circumstances, one should seek recourse from other approaches for computing the moment of inertia.

Equations of motion
In DEM, the motion of a granular system are solved by applying Newton's second law to individual particle.Because granular particles are assumed to be rigid, the movements of a particle can be decomposed into the translation of its centroid and the rotation around the centroid.The equations of translation are written as where M k , V k , and X k denote the mass, translational velocity and position of the centroid of the k-th particle, respectively.F k is the total force acting on the particle, which includes the body force and the inter-particle contact forces.
The angular velocity and rotated angle are determined as the following where J k , Ω k and Θ k are the inertial matrix, angular velocity, and rotated angle of particle k, respectively.T k is the total torque with respect to its centroid.The mass and inertial matrix are calculated using Eqs.( 5) and ( 9), respectively.Note that in the numerical implementation of rotation, the quaternions are used instead of Euler angles.Quaternions can avoid the defect of gimbal lock, which the Euler angle method suffers from [51].Moreover, they are easy to be integrated and reformulated into the rotation matrix.The quaternions have been extensively applied in DEM modelling [17,51,52].The total contact force F c k and torque T k of a particle are obtained by summing up all the contributions from its contact nodes where f c j is the contact force at the contact node j, N k is the total number of contact nodes of the k-th particle.r j = x j − X k is the relative position vector between the contact node x j and the centroid X k .The contact detection and force computation will be detailed in Section 2.3 and 2.4, respectively.
The governing equations of granular particles can be solved using explicit time integration schemes.As the particles are considered rigid, the positions and velocities at contact nodes can be determined by the position of centroid, rotated angle, and particle velocity as follows where x j and v j are the position and velocity of the j-th contact node of the k-th particle.R k is the rotation matrix, which is a function of the vector of the accumulated rotated angle Θ k .x j0 is the initial position of x j , and X k0 indicates the initial position of the centroid of the k-th particle.

Contact detection
Generally, two types of contact detection approaches are employed in DEM.One is the sphereto-sphere method, which features simple contact detection and is employed in many DEM schemes, including the common single-sphere and multi-spheres methods.The other approach is node-tosurface contact method, which is used in the super-quadrics DEM, GJK-DEM, and LS-DEM.Generally, the sphere-to-sphere method is straightforward and can be easily extended to threedimensions.However, it is less accurate when the contact involves flat planes or the particle surface is insufficiently represented.As the normal contact force in the sphere-to-sphere method is always aligned with the line connecting the contact particles, it has oscillations as the surface normal cannot be correctly considered.On the other hand, node-to-surface approaches usually involve higher computational cost and more complex formulation and implementation.The ability to detect contact between arbitrarily shaped particles is critical for the proposed method.To achieve simple yet reliable contact detection, a hybrid contact method is employed.As shown in Figure 2, in SMR-DEM the contact detection is based on the contact nodes, which carry nodal unit normals obtained through Eq. (1).For two particles with potential contact, they are first identified as the slave and master particle by their particle index.In this work, the particle with smaller index is treated as the slave particle.For a contact node x s on the slave particle, first we search for its nearest contact node x m on the master particle.Then the distance between the two contact nodes is defined as where n c is the unit normal of the contact.The determination of n c is of utmost importance in the hybrid contact method.As shown in Figure 3, considering the possible presence of sharp corners, there are in total four situations: (a) Neither of the two contact nodes are sharp corner.In this case, the contact normal n c is taken as the unit normal of the master surface, i.e. n c = n m ; (b) When the slave node is a sharp corner while the master node is not, the contact normal is taken as the unit normal of the master surface n c = n m ; (c) When the master node is a sharp corner and the slave node is not, the contact normal is taken as the unit normal of the slave surface n c = n s ; (d) When both the master and slave nodes are sharp corner, the contact normal is defined as Once the contact normal is determined, the distance between the master and slave contact node x m and x s can be computed using Eq. ( 18).The contact is activated when For situations (a), (b), and (c) where ∆ is the local average spacing of contact nodes.κ is a parameter set as 0.15 for the stability of simulations in this work.The special treatment of sharp corners ensures a smooth transition of the contact normals along surfaces.Once the contact is determined to be active, the overlap d p is calculated as For situations (a), (b), and (c) The overlap d p is used to calculate the contact forces between the two contact nodes, which will be detailed in Section 2.4.
The main processes of the contact detection are to find the nearest pairs of contact nodes, determine the contact normals, and compute the overlaps.Note that the contact detection in SMR-DEM is different from the sphere-to-sphere and the node-to-surface contact detection approaches.All contact detection in SMR-DEM is based on contact nodes without any surface mesh.This contact detection has similar accuracy as that of the node-to-surface approach, while its implementation is as simple as the shpere-to-sphere approach because it avoids complex geometrical operations in three-dimensions.Additionally, we introduces flexible treatment of sharp corners, which ensures the smooth transition of contact normals along contact surfaces.It is therefore termed as hybrid contact method.Figure 4 shows the actual surface representation in the hybrid contact method.A contact surface is implicitly described in SMR-DEM, as shown in Figure 4(c).This surface is slightly different from the actual surface and the segment/triangle discretization of the surface.If more and more contact nodes are used, the implicit surface automatically converges to the actual surface.The hybrid contact method is versatile as it can be applied to both flat and curved surfaces, as well as sharp corners.

Contact force model
To compute the interaction forces, various contact models are available in the literature.In this work, a contact model similar to that in Andrade et al. [53] is employed.The normal contact force acting on the slave contact node is calculated using a linear elastic contact model such that where k n is the normal contact stiffness, and β is the contact damping parameter.∆v n is the normal component of relative velocity between contact nodes x m and x s where v m and v s are the velocities of contact nodes x m and x s , respectively.For the calculation of tangential friction forces, an incremental formulation is used.The tangential force is initialized when the contact occurs and nullified when the slave contact node loses contact with the master particle.The tangential force at time step n can be expressed as where s is the value of the tangential force at the last time step n−1.k s is the tangential contact stiffness, and ∆d t is the incremental relative displacement in the tangential direction, which can be computed from The tangential contact force has an upper limit given by the Coulomb law of friction.Therefore, the final tangential contact force can be obtained as where µ is the frictional coefficient between the two contact particles.Finally, the total contact force f c s acting on the slave contact node is obtained as Similarly, the contact force on the master contact node is Applying the contact force of each contact node to Eqs. ( 14) and (15), the total contact force and contact moment acting on each particle are obtained.

GPU ACCELERATION AND IMPLEMENTATION
In SMR-DEM, the surface of particle is represented by a series of contact nodes, which are used for contact detection and force calculation.Like the multi-spheres DEM, a large number of contact nodes are required to obtain accurate surface representation and force calculation.To increase the numerical efficiency, the GPU parallelization is employed to accelerate the proposed SMR-DEM method.Graphics Processing Units (GPUs) have thousands of computing cores in one graphic card and is effective in cost and energy consumption.GPU acceleration has been applied in parallel computing of DEM [54][55][56] and other particle-based methods [57][58][59].The computations of a large number of particles can be distributed over many threads and performed simultaneously, which can greatly increase the efficiency.In this work, the GPU parallelization is based on the Compute Unified Device Architecture (CUDA).
Algorithm 1: The algorithm for parallel implementation in one computational step Prediction: integrate to n + 1/2 Step 1. Generate the neighbor list; Step 2. Detect contact states; Step 3. Compute contact forces; Step 4. Add contact forces to governing equations of particles; Step 5. Predict new velocity and position of particles by total acceleration: Step 6. Update velocity v n+1/2 i and position x n+1/2 i of contact nodes; Correction: integrate to n + 1 Step 7. Update the neighbor list; Step 8. Repeat steps 2 to 4; Step 9. Correct the new velocity and position of particles: Step 10.Update velocity v n+1 The most time-consuming part of SMR-DEM is contact detection and solving the motion equations.These computational tasks are executed by GPU in a massively parallel fashion.For contact detection, the Cell-linked list described by Domínguez [60] is employed in this work.The searching algorithm is originally proposed for generating neighbor list in smoothed particle hydrodynamic (SPH) method [57,59].The main idea of this method is to divide the domain into square cells, which helps to determine the nearest master-slave pair of contact nodes.The process of this neighbour search is almost the same as that in [58,59,61], except that the neighbour searching is only to find the nearest master contact node.Once the possible pair of contact nodes are determined, it is used for contact detection and force calculation.The total force and torque on each particle are obtained by summing up contributions from all its contact nodes.This summation is performed using the atomic operation provided by CUDA, which ensures that writing to the same memory address by multiple parallel threads can be performed without memory racing.
The motion equations of particles as discussed in Section 2.2 are solved using a predictorcorrector explicit integration scheme.The time integration is also performed by GPU, where the calculations on each particle is executed using one GPU thread.Once the translation and rotation of particles are obtained, the velocity and position of contact nodes can be updated simultaneously using CUDA kernels.The computations in one step are summarized in Algorithm 1.The implementation and GPU parallelization of the proposed SMR-DEM are based on the open-source solver LOQUAT [59], which is originally developed for GPU-accelerated SPH simulation.

NUMERICAL EXAMPLES
In this section, five numerical examples are presented to demonstrate the performance of SMR-DEM.Firstly, a cubic block sliding along a slope is simulated to validate the method in modeling frictional contact.The second case is the collision between two spheres to check the convergence property of the method.In the third one, a hopper flow is simulated and the numerical and experimental results are compared.The fourth example is a rotating drum filled with spherical and tetrahedral particles.In the last example, the collapse of a cylindrical column consisting of granules with irregular shapes is simulated.

A block slides on a rough slope
As shown in Figure 5, a cubic block sliding on a slope is considered as the first benchmark for SMR-DEM.The side length of the block is 1 m and the density is 2500 kg/m 3 .The slope angle of the inclined plane is θ = 30 • .The sliding block has an initial velocity of 6 m/s.Three different friction coefficients, i.e. 0.2, 0.4, and 0.6, are used in the simulation.The damping effect is not considered.In SMR-DEM, the surface of the cubic block is discretized into 1200 triangles and 1008 contact nodes.The surface of the slope is also represented by contact nodes with the same numerical resolution as that of the block.The normal and tangential contact stiffness are 2 × 10 5 N/m and 1.6 × 10 5 N/m, respectively.The damping effect is not considered.At the beginning, the block slides up the slope with a decreasing velocity due to the gravity and frictional force.When the block's velocity decreases to zero, it reaches the maximum height and then slides down the slope.The analytical solutions of the displacement s and velocity v are as follows: where v 0 and s 0 denote the initial velocity and maximum displacement, respectively.µ is the frictional coefficient.t 0 is the time when the block reaches the maximum displacement.
Figure 6 shows the time history of displacement and velocity of the block in the simulation with different frictional coefficients.There is a good agreement between the numerical and theoretical results.Furthermore, in the case with the highest frictional coefficient, the simulation correctly captures the static state of the block after reaching the largest displacement.This case demonstrates that the frictional contact can be modelled correctly using the hybrid contact method, and SMR-DEM has high accuracy in simulations with flat surfaces.

Collision between two rigid spheres
The proposed SMR-DEM employs contact nodes to discretize particle surfaces.As shown in Figure 4, in the contact detection an implicit discretized surface is used.This implicit surface is slightly different from the actual surface, which may induce error particularly on the single particle level.
In this section, the accuracy of SMR-DEM is assessed by studying the interaction between two rigid spheres.Initially, the two spheres are placed along the z axis, as shown in Figure 7.The particles have a diameter of 0.2 m and density of 2500 kg/m 3 .The distance between the mass centres of the two particles is 1.0 m.During the test, particle A is released and falls towards particle B under gravity, while particle B is kept fixed.Both particles are represented by contact nodes with the same average spacing.To check the effect of numerical resolution on the result of SMR-DEM, different numbers of contact nodes are used, ranging from 50 to 1500.In order to take the influence of random alignments into consideration, the orientation of both spheres are randomly initialized.The frictional and damping forces are neglected in this case.Figure 7 shows the motion of particle A with 1000 contact nodes.Once particle A is released, it falls towards particle B with an increasing velocity.As collision happens between two particles, particle A stops falling and rebounds to a certain height.Theoretically, the rebounding should be in the vertical direction and the particle should reach its original location.However, due to discretization errors, there is a rebounding angle δ as shown in Figure 7.The reasons of this discrepancy are twofold.First, the surface representation introduces errors that the direction of the normal contact force does not necessarily point towards the mass centre of the particle, which introduces a rotational velocity.Second, the contact force is not exactly along the z axis in the simulation due to the approximate surface representation.This results in a rotational velocity and a rebounding angle as show in Figure 7.
Three quantities are measured to check the accuracy of SMR-DEM, i.e. the rotation velocity, the rebounding angle, and the coefficient of normal restitution.The rotation velocity and rebounding angle are recorded after the collision finishes.The coefficient of normal restitution is calculated as the ratio between the velocity along the z axis after and before the interaction.For ideal collision between two spheres, the rotation velocity and rebounding angle should be zero, and the coefficient of normal restitution equals one.To check the convergence property of the method, different numbers of contact nodes are used to discretize the two spheres.For the case with a particular number of contact nodes, the orientations of the two spheres are randomly initialized and the falling tests are repeated for ten times.The results are then averaged.
The results of rotation velocity and rebounding angle with different numbers of contact node are given in Figure 8.A clear correlation between the number of contact node used for surface representation and the accuracy of the simulation is observed.The increased number of contact nodes results in smaller deviations.
For the coefficient of normal restitution, we calculated the relative error as e = u a z /u b z − 1 , where u a z and u b z are the velocity along the z axis after and before the collision, respectively.The result of the relative error is shown in Figure 9.It can be seen that the increase in number of  contact nodes gives rise to an improvement of the result.For simulations with more than 200 contact nodes, the deviation is less than 3%.
As observed from the results, there is a clear interconnection between accuracy and number of contact nodes in SMR-DEM, which indicates that it is convergent.It can be seen that the errors decrease as the number of contact nodes increases.Specifically, the errors decrease rapidly when the number of contact nodes is less than 200, whereas the decrease is much slower when more contact nodes are used.Despite the error in the single particle level, its influence on particle assemblies may be not significant.As argued by Kruggel-Emden et al. [34], the errors on the single particle level may be reduced in particle assemblies with a large number of contact.The capability of SMR-DEM in modeling particle assemblies is investigated in the subsequent sections.

Hopper experiment with spherical particles
A hopper experiment with spherical particles is studied numerically to check the performance of SMR-DEM in modelling particle assemblies.The laboratory test shown in Figure 10, which was carried out by Jin et al. [27], is considered as a benchmark for this simulation.In SMR-DEM, the surfaces of particles, the hopper, and the rectangular box are discretized using contact nodes obtained from triangulation.The parameters used in the simulation are listed in Table 1.The simulation parameters including particle radii, particle number, and frictional coefficients are identical to those in the experiment.Particles of two sizes are randomly mixed with a ratio of 1:1.A static state under gravity is achieved by putting all particles in the hopper and damping out the velocities.Figure 10 shows the obtained static state of the system.The profiles of hopper flows from experiment and numerical simulation at different time instances are compared in Figure 11.With the hopper being lifted at a speed of 4.2 mm/s, the particles gradually discharge from the hopper into the rectangular box and form a pile.The profiles of the pile and particles remaining in the hopper obtained using SMR-DEM are in good agreement with the experimental results.The discharge takes about 20 s both in simulation and experiment.
After the system reaches a stable state, the surface of the heap is extracted using Matlab's image processing toolbox.The slope of the left and right surfaces can be fitted using the least squares method, and then the angle of repose is calculated based on the fitted straight line.The results of angle of repose are listed in Table 2.It can be seen that the angle of repose measured from experiment varies from 22.34 • to 25.78 • , and the average value is 24.19 • .The angle of repose obtained through SMR-DEM simulation is 24.07 • , very close to the experimental measurement.
Figure 11: The profiles of the particle system from experiment [27] and numerical simulation at different time instances.

Rotating drum
In this section, the proposed SMR-DEM is applied to the classical rotating drum, which was extensively studied experimentally [62] and numerically [17,63].Two types of particles are considered in this work, i.e. spherical particles used to validate SMR-DEM in modelling dynamic flow, and tetrahedra to check the performance of SMR-DEM in simulating non-spherical shapes.The horizontal drum considered in the simulations has a radius of 500 mm and a width of 24 mm.Three sets of simulations are performed.In the first one, the drum is filled with 3000 uniform spherical particles of radius r s = 1.5 mm.The other two have the drum filled with 2600 and 3000 uniform tetrahedral particles, respectively.The tetrahedra have a side length l t = 4.9 mm.In the simulations with tetrahedral particles, the sharp corner treatment introduced in Section 2.1 is used with the criterion λ = 0.5.With this criterion, all contact nodes located at edges and corners are identified as sharp corner and specially treated using the method presented in Section 2.3.The configurations of the drum and particles are the same as those used in Wachs et al. [17].The parameters employed in the simulations are listed in Table 3.The same contact parameters for particle-particle and particle-drum interaction are assumed.To eliminate the influence of the side walls, the friction coefficient of the side walls is set as zero.The initial state is obtained by releasing the particles from random locations inside the drum.The particles are subjected to gravity and pack near the bottom of the drum.Once the particle assemblies reach a stable state, the drum is rotated at a specific speed.As shown in the literature, different flow patterns manifest progressively with the increase of rotating velocity [64].In this work, we focus on the rolling regime, in which a thin layer of particles flow down near the free surface while the remaining particles are conveyed upwards by the drum.The rotating rates in each set of simulations are ω = 5, 20, 42, 65 rpm.
Figure 12 presents the snapshots of materials at different time instances with ω = 65 rpm in the three sets of simulations.The profiles of final steady states agree well with the numerical results in [17] (Fig. 7 in their paper).Both the two types of particles move with the drum as a whole at the beginning of the simulation.As the slope angle of the free surface grows, a shear band develops under the free surface and a dynamic particle flow is initiated.The particles below the shear band move with the drum until they reach the critical height and start flowing down.The comparison between spherical and tetrahedral particles indicates that the particle shape has a significant influence on the flow dynamics.In the steady state with ω = 65 rpm, a flat free surface is observed for spherical particles, while an S-shape free surface is found for tetrahedral particles.Furthermore, it is obvious that the dilation caused by shear deformation is more apparent for the tetrahedral particles if one compares the results of spherical and tetrahedral particles.It is also observed that in the simulations with tetrahedral particles, the actual solid filling ratio have a significant influence on the flow dynamics in the rotating drum.With lower solid filling ratio (less particles), the S-shape free-surface is more evident and the dynamic angle of repose is higher.As shown in Figure 12, the velocity distribution obtained using SMR-DEM is reasonable and the flow patterns are in good agreement with those from Wachs et al. [17].To quantitatively compare our results and those from Wachs et al. [17], the dynamic angle of repose at different rotation rate is computed and shown in Figure 13.For the simulations with spherical particles, the dynamic angle of repose is calculated by fitting the flat free surface using the least squares method.As for the tetrahedral particles, the problem is more complex since the free surface is a S-shape curve.To address this issue, the central region of the free surface, which is between -0.5R drum and 0.5R drum , is chosen as the fitting area for the least squares method.This method was also adopted in Wachs et al. [17].As shown in Figure 13, the dynamic angle of repose for tetrahedral particles is much larger than that for spherical particles, even the same frictional coefficient is employed for the two types of particles.As for the simulations with different number of tetrahedral particles, it is found that the lower filling ratio leads to the increase of the angle of repose.Generally, our numerical results are in consistent with those from Wachs et al. [17] and the experimental data from Parker et al. [62].
To demonstrate the effect of sharp corner treatment, an additional simulation of 2600 tetrahedral particles without sharp corner treatment is performed.The particle distributions and velocities at two instances are shown in Figure 14.Compared with the results obtained with the sharp corner treatmen, there are mainly two differences.First, the deposition takes up more filling volume at the initial packing state.Second, without shape corner treatment the shape of the free-surface is less curved and more like that from the spherical particles.The reason of these discrepancies is that the normal contact constraint is too strong when unit normals at corner particles are taken as the average of the unit normals of the associated triangles.This leads to a increases in the implicit contact volume and a decrease of the sharpness of the tetrahedra, which result in the increase of packing volume and the smoothing of the free-surface, respectively.The comparison indicates that the sharp corner treatment is necessary for the simulation with irregular particles, particularly for arbitrarily shaped particles with many sharp corners and edges.

The collapse of a cylinder of rocks with irregular shapes
To further demonstrate the effectiveness of the proposed method in modelling assemblies of arbitrarily shaped particles, a rockfall analysis is carried out.Gao et al. [65] studied the effect of rock shape on the impact force of a falling rock assembly on structures.In their simulation, a maximum of 99 rocks with four shapes are considered using the multi-spheres method.According to their study, the particle shape has a significant influence on the falling pattern and impact force of rock assemblies.In this section, the collapse of a rock assembly is studied.The rock assembly consists of 10,000 rocks generated from a rock template [66] with nineteen natural rock shapes shown in Figure 15.The density of the rocks is assumed to be 2,000 kg/m 3 .The mass of each stone ranges between 0.84 and 7.88 kg.The discretization of the rock geometry and calculation of particle properties including normal vectors at contact nodes, mass, centroid and moment of inertia follow the particle representation approach described in Section 2.1.The total number of contact nodes used to represent the surfaces of the 10,000 rocks is 3,814,088.
In the numerical simulation, the rock assembly is generated within a vertical cylinder of diameter D = 3 m.To obtain the initial static state, each rock is generated from a random rock model in the template.Note that each rock is initialized with random position and orientation before being released.All generated rocks are subjected to gravity and damping.The final packing status is shown in Figure 16.The final rock column has a height of 4.7 m.The total mass of the rock assembly is 41,539 kg, the initial void ratio is about 0.6.The parameters used in the simulation are listed in Table 4.After the rock assembly achieves the equilibrium state, the cylinder wall is removed to allow the rocks to fall under gravity and spread on a flat surface.The flat surface is discretized using 1,000,000 contact nodes with a uniform spacing of 0.02 m.The contact parameters for particleparticle and particle-flat surface are assumed the same.The simulated profiles of the rock assembly at different time instances after the release are shown in Figure 17.It is observed that the large deformation of the rock assembly is well described by the simulation.At first, there is a fast flow region near the free surface of the rock assembly, resulting in a discontinuity in the collapse surface.As the collapse continues, the discontinuity disappears and a stone heap with a parabolic surface is formed when the system gradually comes to rest.The rockfall process lasts about 20 s.Although from literature we cannot find any exact numerical or experimental results with the same configuration and materials, the simulated process can be compared to the collapse of sand column (with similar aspect ratio) in Lube et al. [67].Qualitative agreement between the two processes can be observed.This test case demonstrates that SMR-DEM can give stable and reasonable results in modelling particles involving sharp corners and complex natural shapes.

Numerical efficiency
In SMR-DEM, each particle is represented using a large number of contact nodes, which may induce high computational cost.Table 5 summarizes the information of efficiency of some of the numerical examples.The numerical efficiency of the GJK-DEM from Wachs et al. [17], which was based on serial CPU computing, is also listed for comparison.In the simulations, physical durations and time steps are different; therefore, the total simulation time is of little use in efficiency analysis.Here we use the more objective frame per second (FPS) and execution time per particle and step (ETPS) to measure the numerical efficiency.FPS is the number of computational steps executed in one second of wall-clock time, and ETPS is the computational time used to model one particle in one step.First we focus on the SMR-DEM simulations of spherical and tetrahedral particles in the hopper and rotating drum cases.In these simulations, each particle is discretized using around 200 contact nodes.It is found that for these simulations the FPS is directly related to the particle number.The ETPS for these simulations are almost the same, regardless of the shape of the particle.This is expected, as in SMR-DEM the computational cost per particle is dependent on the number of contact nodes, not on its shape.As a matter of fact, as long as the numbers of contact nodes are the same, two particles with very different shape can have the same computational cost.This is different from the GJK-DEM, in which particle shape has significant influence on the numerical cost.For instance, in GJK-DEM, the cost of tetrahedral particles is about seven times higher than that of spherical particles.
This feature of SMR-DEM makes it very suitable for particles of complex shape, but for particles with simple shapes the computational cost is usually too high.For example, for simple spherical particles, the ETPS in SMR-DEM is between 3.35 and 3.56 µs, while in GJK-DEM the ETPS is 5.1 µs.The efficiency of SMR-DEM seems slightly better.However, the simulations of SMR-DEM employs GPU parallelization, whereas the GJK-DEM simulation is based on serial CPU computing, indicating that the numerical efficiency of SMR-DEM is much lower than that of GJK-DEM when simulating spherical particles.With tetrahedral particles, the situation improves and the computational cost of GJK-DEM is about 10 times higher than that of SMR-DEM.Note that tetrahedra cannot be considered as complex shape; therefore, we can expect that with increasing complexity the efficiency of SMR-DEM does not change much but the GJK-DEM has higher computational cost.This can be confirmed by the last example with irregular particles, where the ETPS is 5.74 µs.This ETPS is slightly higher than those of spherical and tetrahedral particles only because in the last numerical example the average contact nodes per particle is 381, larger than that in other simulations.

Strengths, limitations, and future developments
As shown in previous analysis, the SMR-DEM can model particles of arbitrarily irregular shapes.It employs the conventional surface mesh for shape representation, which makes the particle description straightforward.Its formulation and implementation are simple.The method is stable and convergent, and is particularly suitable for problems with complex irregular shapes.On the other hand, the proposed SMR-DEM is not suitable for particles of simple shape, e.g.sphere.Even for common non-spherical shapes such as cubes, tetrahedra, and cylinders, we believe there are other more efficient DEM methods.Besides, the present SMR-DEM discretizes all particles with similar numerical resolution.This approach is not efficient if the particle size varies greatly in a simulation.Future developments may include optimized treatment for simple particles and improved contact detection suitable for adaptive discretization.

CONCLUSIONS
A surface mesh represented discrete element method (SMR-DEM) is developed for arbitrarily shaped particles.In SMR-DEM, the particle surface is represented with contact nodes, which can be conveniently obtained from surface mesh.Based on surface mesh, the particle properties such as unit normals at contact nodes, total mass, centroid, and moment of inertia can be directly calculated.A hybrid contact detection method suitable for contact node representation is proposed, which enables easy contact detection.Furthermore, the GPU acceleration is employed to accelerate the SMR-DEM.
A couple of examples are performed to validate the proposed method.It is found that the method gives accurate results for particle level frictional contact, and it converges as the number of contact node per particle increases.When used to simulate shperical and tetrahedral particles, the SMR-DEM provides results that are well corroborated by experiments and numerical results from the literature.Particularly, SMR-DEM gives reasonable three-dimensional simulations for particle assembly of complex irregular shapes, which indicates that it can be employed in various practical applications with arbitrarily shaped particles.The computational cost of SMR-DEM is independent of the particle shape; therefore, it is efficient in modelling particle systems of complex irregular shapes but less efficient in handling particles of simple shapes.

Figure 1 :
Figure 1: The schematic diagram of particle surface representation of three typical shapes and one irregular shape (top: original shapes of particles; middle: surface meshes with triangles; bottom: distribution of contact nodes at particle surfaces).

Figure 2 :
Figure 2: Illustration of two contacting particles in two-dimensions.

Figure 3 :
Figure 3: The four cases for determining the contact normal n c .

Figure 4 :
Figure 4: The actual surface discretization of SMR-DEM: (a) the original surface and contact nodes; (b) segment/triangle discretization of the surface and the unit normals; (c) the implicit surface discretization in SMR-DEM.

i and position x n+1 i of contact nodes; Step 11 .
Add elapsed simulation time t = t + ∆t;A k = dV k /dt represents the acceleration of particle k.

Figure 5 :
Figure 5: The configuration of the slope and the cubic block.

Figure 6 :
Figure 6: The development of displacement and velocity with different frictional coefficients.

Figure 7 :
Figure 7: The predicted motion of particles with 1000 contact nodes, together with the paticle velocity.* Figure (c) and (d) are the profiles when collision and seperation happen.

Figure 8 :
Figure 8: The rotation velocity and rebound angle with different numbers of contact node.

Figure 9 :
Figure 9: Relative error of the coefficient of normal restitution with different numbers of contact node.

Table 1 . 4 ×
Parameters used in the simulation 10 −6 Number of contact nodes per particle 212 Number of contact nodes for the hopper 82,149 Number of contact nodes for the rectangular box 169,463

Figure 12 :
Figure 12: The snapshots of materials at different time instances with ω = 65 rpm.

Figure 13 :
Figure 13: The dynamic angle of repose at different rotation rate.

Figure 14 :
Figure 14: The snapshots of the flowing particles without sharp corner treatment.

Figure 15 :
Figure 15: The rock models used in the simulation [66].

Figure 16 :
Figure 16: The final packing status of the rock column in front and top views.

Table 4 .× 10 − 5
Parameters used in the rockfall simulation Number of contact nodes for the rocks 3,814,088 Number of contact nodes for the flat surface 1,000,000

Figure 17 :
Figure 17: The collapsing process of the rock assembly.

Table 3 .
Parameters used in the rotating drum simulation

Table 5 .
Summary of numerical efficiency of SMR-DEM.