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- 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 ectively model particle assemblies of various regular and complex irregular shapes. Keywords: Discrete element method (DEM), Granular materials, Arbitrarily shaped particle, Hybrid contact method, GPU acceleration 1. 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, ﬂow dynamics, shear stress, contact forces, and energy dissipation. It has achieved notable success in a wide range of applications, e.g. granular ﬂows [2–4], natural hazards [5, 6], industrial processes [7–10], and geotechnical problems [11– 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 parti- cles 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 signiﬁcant Corresponding author. E-mail: pengchong@boku.ac.at Preprint submitted to Powder Technology March 23, 2020This is a preprint of the paper published in Powder Technology Accepted citation:Zhan, L., Peng, C., Zhang, B. and Wu, W., 2021. A surface mesh represented discrete element method (SMR-DEM) for particles of arbitrary shape. Powder Technology, 377, pp.760-779. DOI: https://doi.org/10.1016/j.powtec.2020.09.046 inﬂuence 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 segre- gation [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 con- tact detection. To date, many approaches have been developed to model non-spherical particles in DEM. An incomplete list includes: polygon/polyhedron [20–22], multi-spheres [23, 24], super- quadrics [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 eciency were extensively discussed. Many methods are being developed in parallel but no method is to be favoured in terms of accuracy and eciency. 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 eciency, we believe that versatility and simplicity are also impor- tant 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 diculties in modelling particles with ﬂat surfaces. The Gilbert-Johanson-Keethi algorithm-based DEM (GJK-DEM) can model convex shapes [17] and by gluing multiple parti- cles 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 chal- lenge 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–45] to ﬁll 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 ecient, or at least can be conve- niently and eciently parallelized; (2) its formulations and implementation are simple in three- dimensions; (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, polyhe- dral, and arbitrarily irregular, using the same algorithm; (5) the data preparation should be sim- ple. 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 addi- tional 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 2 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, in- cluding surface discretization, motion equations, and contact computation. Numerical implemen- tation 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 eciency and limitations is provided in Section 5. Finally, some conclusions are given in Section 6. 2. SURFACE MESH REPRESENTED DEM 2.1. 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. 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). There are mainly two beneﬁts 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 di- rectly utilizing the triangular meshes, additional algorithms for describing particle shapes, e.g. the ﬁlling algorithms in the multi-sphere method [42–45], super-quadric equations [28], or level-set functions [30, 31], are no longer needed. This greatly simpliﬁes the formulation, implementation, 3 and data preparation. Second, with triangular surface mesh, arbitrarily complex shapes can be represented without diculty. All regular and arbitrarily irregular particle shapes can be described using the same uniﬁed 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, reﬁnement, 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 Ai;i=1;:::; Nbe the Ntriangles on the surface of particle P. For each triangle, the outward unit normal nican be deﬁned as ni=ˆni=jˆnij, where the unnormalized vector ˆniis computed as ˆni=(biai)(ciai), with ai,biandcibeing the vertices of triangle Aiin a counter-clockwise order. Consequently, the unnormalized vector at contact node jis obtained by averaging the normal vectors of associated triangles ˆnj=1 NjNjX i=1ni (1) where Njis the number of triangles associated with the contact node. The unit normal vector at contact node jis obtained as nj=ˆnj=jˆnjj. At each contact node, a nodal shape factor jis deﬁned as j=max i ninj  ;i=1;2;;Nj (2) which quantiﬁes the maximum di erence 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 di erent 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 Mof a host particle can be expressed as M=Z (x)d (3) 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 M= 3Z rxd = 3Z SxndS= 3NX i=1Z AixnidS (4) where Nis the number of triangles on the surface of particle. As xniis constant across Ai, and jˆnijequals two times the area of Ai, the mass of a host particle can be computed as 4 M= 6NX i=1aiˆni (5) where aiis the position of the ﬁrst vertex of triangle Ai. Similarly, the centroid Xof a particle can be calculated from X=1 VZ xd = MZ xd (6) where Vis the particle volume. Applying the divergence theorem once more, the above equation can be rewritten as Xl= MZ S1 2x2 lnldS= 2MNX i=1Z Aix2 lnildS (7) where lis the index of spatial coordinates, Xl,xl,nldenote the l-th component of the vectors X, x, and n, respectively. nilis 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 ﬁnal expression of the coordinate of centroid is written as Xl= 48MNX i=1ˆnil (ail+bil)2+(ail+cil)2+(bil+cil)2 (8) where ail,bilandcilare the l-th coordinate component of vertices ai,biandci, 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 Ntetrahedra by con- necting 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 J=NX i=1Z i[(rr)Ir r]d (9) where iis the integral domain of the i-th tetrahedron, and r=xXis the relative position vector between a point xin the i-th tetrahedron and the center of mass Xof 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. 5 2.2. 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 MkdVk dt=Fk (10) dXk dt=Vk (11) where Mk,Vk, and Xkdenote the mass, translational velocity and position of the centroid of the k-th particle, respectively. Fkis 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 Jkd k dt=Tk kJk k (12) d k dt= k (13) where Jk, kandkare the inertial matrix, angular velocity, and rotated angle of particle k, respectively. Tkis 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 su ers 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 Fc kand torque Tkof a particle are obtained by summing up all the contributions from its contact nodes Fc k=NkX jfc j (14) Tk=NkX jrjfc j (15) where fc jis the contact force at the contact node j,Nkis the total number of contact nodes of the k-th particle. rj=xjXkis the relative position vector between the contact node xjand the centroid Xk. 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 6 xj=Xk+Rk(xj0Xk0) (16) vj=Vk+ k(xjXk) (17) where xjandvjare the position and velocity of the j-th contact node of the k-th particle. Rkis the rotation matrix, which is a function of the vector of the accumulated rotated angle k.xj0is the initial position of xj, and Xk0indicates the initial position of the centroid of the k-th particle. 2.3. Contact detection Generally, two types of contact detection approaches are employed in DEM. One is the sphere- to-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-to- surface 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 three- dimensions. However, it is less accurate when the contact involves ﬂat planes or the particle surface is insuciently 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. Figure 2: Illustration of two contacting particles in two-dimensions. 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 ﬁrst identiﬁed 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 xson the slave particle, ﬁrst we search for its nearest contact node xmon the master particle. Then the distance between the two contact nodes is deﬁned as 7 d=(xmxs)nc (18) where ncis the unit normal of the contact. Figure 3: The four cases for determining the contact normal nc. The determination of ncis 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 ncis taken as the unit normal of the master surface, i.e. nc=nm; (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 nc=nm; (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 nc=ns; (d) When both the master and slave nodes are sharp corner, the contact normal is deﬁned as follows nc=xsxm kxsxmk(19) Once the contact normal is determined, the distance between the master and slave contact node xmandxscan be computed using Eq. (18). The contact is activated when (d>0 For situations (a); (b);and (c) d>For situation (d)(20) 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 8 the contact normals along surfaces. Once the contact is determined to be active, the overlap dpis calculated as (dp=d For situations (a); (b);and (c) dp= +dFor situation (d)(21) The overlap dpis 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 ﬁnd the nearest pairs of contact nodes, deter- mine the contact normals, and compute the overlaps. Note that the contact detection in SMR-DEM is di erent 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 ﬂexible 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: 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. 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 di erent 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 ﬂat and curved surfaces, as well as sharp corners. 9 2.4. 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 fn s=kn(dp+ v n)nc (22) where knis the normal contact sti ness, and is the contact damping parameter. vnis the normal component of relative velocity between contact nodes xmandxs vn=(vmvs)nc (23) where vmandvsare the velocities of contact nodes xmandxs, respectively. For the calculation of tangential friction forces, an incremental formulation is used. The tan- gential force is initialized when the contact occurs and nulliﬁed when the slave contact node loses contact with the master particle. The tangential force at time step ncan be expressed as ft(n) s=ft(n1) s+ksdt (24) where ft(n1) s is the value of the tangential force at the last time step n1.ksis the tangential contact sti ness, and dtis the incremental relative displacement in the tangential direction, which can be computed from dt= (vmvs)nc(vmvs)nc t (25) The tangential contact force has an upper limit given by the Coulomb law of friction. There- fore, the ﬁnal tangential contact force can be obtained as ft s=8>>>><>>>>: fn s ft(n) s ft(n) s ft(n) s  fn s ft(n) s ft(n) s < fn s (26) whereis the frictional coecient between the two contact particles. Finally, the total contact force fc sacting on the slave contact node is obtained as fc s=fn s+ft s (27) Similarly, the contact force on the master contact node is fc m=fc s. 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. 3. 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 10 of contact nodes are required to obtain accurate surface representation and force calculation. To increase the numerical eciency, 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 e ective in cost and energy consumption. GPU acceleration has been applied in parallel computing of DEM [54–56] and other particle-based methods [57–59]. The computations of a large number of particles can be distributed over many threads and performed simultaneously, which can greatly increase the eciency. In this work, the GPU parallelization is based on the Compute Uniﬁed 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: Vn+1=2 k=Vn k+1 2tAn k;Xn+1=2 k=Xn k+1 2tVn+1=2 k Step 6. Update velocity vn+1=2 iand position xn+1=2 iof 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: Vn+1 k=Vn k+ tAn+1=2 k;Xn+1 k=Xn k+ tVn+1=2 k Step 10. Update velocity vn+1 iand position xn+1 iof contact nodes; Step 11. Add elapsed simulation time t=t+ t; Ak=dVk=dtrepresents the acceleration of particle k. The most time-consuming part of SMR-DEM is contact detection and solving the motion equa- tions. 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 search- ing algorithm is originally proposed for generating neighbor list in smoothed particle hydrody- namic (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 ﬁnd the nearest master contact node. Once the possible pair of contact nodes are deter- mined, 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 predictor- corrector 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 11 using CUDA kernels. The computations in one step are summarized in Algorithm 1. The im- plementation 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. 4. NUMERICAL EXAMPLES In this section, ﬁve 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 conver- gence property of the method. In the third one, a hopper ﬂow is simulated and the numerical and experimental results are compared. The fourth example is a rotating drum ﬁlled with spherical and tetrahedral particles. In the last example, the collapse of a cylindrical column consisting of granules with irregular shapes is simulated. 4.1. A block slides on a rough slope As shown in Figure 5, a cubic block sliding on a slope is considered as the ﬁrst benchmark for SMR-DEM. The side length of the block is 1 m and the density is 2500 kg=m3. The slope angle of the inclined plane is =30. The sliding block has an initial velocity of 6 m/s. Three di erent friction coecients, i.e. 0:2, 0:4, and 0:6, are used in the simulation. The damping e ect is not considered. Figure 5: The conﬁguration of the slope and the cubic block. 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 sti ness are 2 105N/m and 1:6105N/m, respectively. The damping e ect 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 sand velocity vare as follows: s=(v0t1 2g(sin+cos)t2;tt0 s01 2g(sincos)t2;t>t0(28) 12 v=(v0g(sin+cos)t;tt0 g(sincos)t; t>t0(29) where v0ands0denote the initial velocity and maximum displacement, respectively. is the frictional coecient. t0is 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 di erent frictional coecients. There is a good agreement between the numerical and theoretical results. Furthermore, in the case with the highest frictional coecient, 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 ﬂat surfaces. Figure 6: The development of displacement and velocity with di erent frictional coecients. 4.2. 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 di erent 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 zaxis, as shown in Figure 7. The particles have a diameter of 0.2 m and density of 2500 kg=m3. 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 ﬁxed. Both particles are represented by contact nodes with the same average spacing. To check the e ect of numerical resolution on the result of SMR-DEM, di erent numbers of contact nodes are used, ranging from 50 to 1500. In order to take the inﬂuence of random alignments into consideration, the orientation of both spheres are randomly initialized. The frictional and damping forces are neglected in this case. 13 Figure 7: The predicted motion of particles with 1000 contact nodes, together with the paticle velocity. * Figure (c) and (d) are the proﬁles when collision and seperation happen. 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 zaxis 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 coecient of normal restitution. The rotation velocity and rebound- ing angle are recorded after the collision ﬁnishes. The coecient of normal restitution is calculated as the ratio between the velocity along the zaxis after and before the interaction. For ideal col- lision between two spheres, the rotation velocity and rebounding angle should be zero, and the coecient of normal restitution equals one. To check the convergence property of the method, di erent 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 di erent 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 coecient of normal restitution, we calculated the relative error as e=kua z=ub z1k, where ua zandub zare the velocity along the zaxis 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 14 Figure 8: The rotation velocity and rebound angle with di erent numbers of contact node. Figure 9: Relative error of the coecient of normal restitution with di erent numbers of contact node. 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. Speciﬁcally, 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 inﬂuence on particle assemblies may be not signiﬁcant. 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. 4.3. 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 15 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 coecients 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. Table 1. Parameters used in the simulation Parameter Value Particle radii (mm) 4, 4.5 Particle number 2020 Density (kg=m3) wall 1180 particle 1410 Friction coecient particle-particle 0.27 particle-wall 0.31 Normal sti ness (N/m) 1:0106 Tangential sti ness (N/m) 6:0105 Damping coecient (s) 0.001 Time step (s) 5:4106 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 10: Initial state of the hopper system. Left: experiment [27]; right: SMR-DEM simulation. The proﬁles of hopper ﬂows from experiment and numerical simulation at di erent time in- stances are compared in Figure 11. With the hopper being lifted at a speed of 4.2 mm/s, the parti- cles gradually discharge from the hopper into the rectangular box and form a pile. The proﬁles 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. 16 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 ﬁtted using the least squares method, and then the angle of repose is calculated based on the ﬁtted 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 proﬁles of the particle system from experiment [27] and numerical simulation at di erent time instances. 17 Table 2. Angle of repose from experiment [27] and SMR-DEM simulation (). No. Left Right Average 1 25.30 24.31 24.81 2 24.79 23.27 24.03 3 24.99 23.82 24.41 4 25.78 22.34 24.06 5 23.92 23.34 23.63 Simulation 24.69 23.45 24.07 4.4. 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 con- sidered in this work, i.e. spherical particles used to validate SMR-DEM in modelling dynamic ﬂow, 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 ﬁrst one, the drum is ﬁlled with 3000 uniform spherical particles of radius rs=1:5 mm. The other two have the drum ﬁlled with 2600 and 3000 uniform tetrahedral particles, respectively. The tetrahedra have a side length lt=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 identiﬁed as sharp corner and specially treated using the method presented in Section 2.3. The conﬁgurations 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 inﬂuence of the side walls, the friction coecient of the side walls is set as zero. Table 3. Parameters used in the rotating drum simulation Parameter Value Particle density (kg=m3) 2500 Friction coecient 0.55 Normal sti ness (N/m) 5105 Tangential sti ness (N/m) 3105 Damping coecient (s) 0.001 Time step (s) 1:35105 Number of contact nodes per particle Sphere 199 Tetrahedron 202 Number of contact nodes for the drum 115,269 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 18 assemblies reach a stable state, the drum is rotated at a speciﬁc speed. As shown in the literature, di erent ﬂow 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 ﬂow 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 di erent time instances with !=65 rpm in the three sets of simulations. The proﬁles of ﬁnal 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 ﬂow is initiated. The particles below the shear band move with the drum until they reach the critical height and start ﬂowing down. The comparison between spherical and tetrahedral particles indicates that the particle shape has a signiﬁcant inﬂuence on the ﬂow dynamics. In the steady state with !=65 rpm, a ﬂat 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 ﬁlling ratio have a signiﬁcant inﬂuence on the ﬂow dynamics in the rotating drum. With lower solid ﬁlling 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 ﬂow patterns are in good agreement with those from Wachs et al. [17]. Figure 12: The snapshots of materials at di erent time instances with !=65 rpm. 19 Figure 13: The dynamic angle of repose at di erent rotation rate. To quantitatively compare our results and those from Wachs et al. [17], the dynamic angle of repose at di erent rotation rate is computed and shown in Figure 13. For the simulations with spherical particles, the dynamic angle of repose is calculated by ﬁtting the ﬂat 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 drumand 0.5R drum, is chosen as the ﬁtting 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 coecient is employed for the two types of particles. As for the simulations with di erent number of tetrahedral particles, it is found that the lower ﬁlling 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 e ect of sharp corner treatment, an additional simulation of 2600 tetrahe- dral 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 di erences. First, the deposition takes up more ﬁlling 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. 20 Figure 14: The snapshots of the ﬂowing particles without sharp corner treatment. 4.5. The collapse of a cylinder of rocks with irregular shapes To further demonstrate the e ectiveness of the proposed method in modelling assemblies of arbitrarily shaped particles, a rockfall analysis is carried out. Gao et al. [65] studied the e ect 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 signiﬁcant inﬂuence 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=m3. 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 diam- eterD=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 ﬁnal packing status is shown in Figure 16. The ﬁnal 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. 21 Figure 15: The rock models used in the simulation [66]. Figure 16: The ﬁnal packing status of the rock column in front and top views. Table 4. Parameters used in the rockfall simulation Parameter Value Particle density (kg=m3) 2,000 Friction coecient 0.5 Normal sti ness (N/m) 8105 Tangential sti ness (N/m) 6105 Damping coecient (s) 0.02 Time step (s) 2:25105 Number of contact nodes for the rocks 3,814,088 Number of contact nodes for the ﬂat surface 1,000,000 22 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 ﬂat surface. The ﬂat surface is discretized using 1,000,000 contact nodes with a uniform spacing of 0.02 m. The contact parameters for particle- particle and particle-ﬂat surface are assumed the same. The simulated proﬁles of the rock assembly at di erent 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 ﬁrst, there is a fast ﬂow 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 ﬁnd any exact numerical or experimental results with the same conﬁguration 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. Figure 17: The collapsing process of the rock assembly. 23 5. DISCUSSION 5.1. Numerical eciency 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 eciency of some of the numerical examples. The numerical eciency of the GJK-DEM from Wachs et al. [17], which was based on serial CPU computing, is also listed for comparison. Table 5. Summary of numerical eciency of SMR-DEM. Case Particle shape Particle number Contact node number FPS ETPS ( s) Hopper ﬂow Spherical 2,020 679,852 139.0 3.56 Rotating drumSpherical 3,000 712,269 98.5 3.38 Tetrahedral (2600) 2,600 640,469 111.4 3.45 Tetrahedral (3000) 3,000 721,269 99.5 3.35 Spherical [17] 3,000 - 64.9 5.1 Tetrahedral [17] 3,000 - 9.4 35.3 Rock collapse Irregular 10,000 4,814,088 17.4 5.74 In the simulations, physical durations and time steps are di erent; therefore, the total simula- tion time is of little use in eciency analysis. Here we use the more objective frame per second (FPS) and execution time per particle and step (ETPS) to measure the numerical eciency. 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 di erent shape can have the same computational cost. This is di erent from the GJK-DEM, in which particle shape has signiﬁcant inﬂuence 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 par- ticles 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 eciency 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 eciency 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 increas- ing complexity the eciency of SMR-DEM does not change much but the GJK-DEM has higher 24 computational cost. This can be conﬁrmed 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. 5.2. 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 ecient DEM methods. Besides, the present SMR-DEM discretizes all particles with similar numerical resolution. This approach is not ecient 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. 6. CONCLUSIONS A surface mesh represented discrete element method (SMR-DEM) is developed for arbitrar- ily 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 pro- posed, 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 ecient in modelling particle systems of complex irregular shapes but less ecient in handling particles of simple shapes. ACKNOWLEDGMENT The authors would like to thank the funding from the National Key Research and Development Program of China (Project No. 2017YFC0404802) and the National Natural Science Foundation of China (Grant No. 51709230 and 51979143). The work in this paper was also partly funded by the EU Horizon 2020 RISE projects HERCULES (778360) and FRAMED (734485), and the Nazarbayev University Research Fund (SOE2017001). 25