Multisphere Representation of Convex Polyhedral Particles for DEM Simulation

The representation of particles of complex shapes is one of the key challenges of numerical simulations based on the discrete element method (DEM). A novel algorithm has been developed by the authors to accurately represent 2D arbitrary particles for DEM modelling. In this paper, the algorithm is extended from 2D to 3D to model convex polyhedral particles based on multisphere methods, which includes three steps: the placement of spheres at the corners, along the edges, and on the facets in sequence. To give a good representation of a polyhedral particle, the spheres are placed tangent to the particle surface in each step. All spheres placed in the three steps are clumped together into a clump in DEM. In addition, the mass properties of the clump are determined based on the corresponding polyhedral particle to obtain accurate simulation results. Finally, an example is used to validate the robust and automatic performance of the algorithm in generating a sphere clump model for an assembly of polyhedral particles. A current FORTRAN version of the algorithm is available by contacting the authors.


Introduction
Discrete element method (DEM), initially proposed by CUNDALL and STRACK in 1979 [1], provides a powerful tool to perform numerical simulations of discrete particle systems. In most discrete element software, complex particles are simplified as discs (in 2D) or spheres (in 3D) for easy implementation of contact detection and timesaving. However, it is not accurate to represent a particle with only one disc or sphere, because one disc or sphere on its own cannot capture all the characteristics of the complex geometrical shape of most granular particles. erefore, DEM modelling of granular materials in which the particles are represented by discs or spheres suffers from limited capability in capturing essential aspects of interparticle interactions and corresponding mechanical behaviors [2,3].
To obtain more realistic behavior of particles, some particle shape approximations have been introduced into DEM towards realistic particle shape modelling [4]. Basically, there are two different methods: new primitive methods and multisphere methods. e new primitive methods mainly lie in the direct implementation of nonspherical primitives such as ellipsoids [5][6][7], polyhedra [2,8,9], and spherocylinders [10,11] in DEM. It is a common drawback for new primitive methods that some complicated algorithms need to be developed for contact detection and determination of contact forces, leading to the deterioration of computational efficiency in the simulations. On the other hand, multisphere methods are based on that a nonspherical particle is represented by a clump of rigid spheres which may vary in size and even overlap each other [12]. A fairly comprehensive review of some available multisphere methods can be founded in scientific literatures [13][14][15][16][17][18]. Currently, multisphere methods have been applied widely in DEM simulations and implemented into some DEM-packages such as PFC3D [19] and YADE [20] because contact detection and force calculation are still based on some simple algorithms valid for spheres and are, therefore, very efficient and robust.
Polyhedral particles are widely encountered in the field of granular materials such as railway ballast [21,22], concrete aggregate [23], and soil-rock mixture (S-RM) [24]. In those numerical simulations, they are often simplified as convex and irregular polyhedra whose surfaces consist of several triangle facets. DEM offers a means of obtaining micromechanical insight into the behaviors of granular materials mentioned above. However, how to accurately represent polyhedral particles with clumps of spheres is still one of the key challenges for researchers in DEM simulations. Although the available multisphere methods have been proposed to represent such polyhedral particles with clumps of spheres, they cannot give a good approximation of the particle shape. In this work, a novel algorithm developed by ZHANG et al. [25] is extended from 2D to 3D to accurately represent convex and irregular polyhedral particles based on multisphere methods for DEM simulations.

Representation of Convex
Polyhedral Particles 2.1. Algorithm. e algorithm developed by Zhang et al. [25] in 2D is briefly introduced firstly in this section. Figure 1 shows schematic operations of the algorithm for the representation of a star-like particle by a clump of discs. As shown in Figure 1, the algorithm includes three steps: placement of discs at the corners (Figure 1(a)), placement of discs along the edges (Figure 1(b)), and placement of discs in the domain (Figure 1(c)). All discs placed in the three steps can be combined as a clump in DEM simulation. To obtain a good representation of the particle, some discs are placed tangent to the interior side of the particle boundary in the first two steps, as shown in Figure 1(b). To be specific, some identical discs are placed at each corner, whose centers are located on the interior angle bisector of each corner in the first step, as shown in Figure 1(a); some identical discs tangent to each other are placed along the interior side of each edge in the second step, as shown in Figure 1(b). In the last step, a series of substeps are performed to place some discs in the domain step by step based on a local Delaunay tessellation, as shown in Figures 1(d)-1(f ). In each substep, a set of the centers of some placed discs are selected to create a local Delaunay mesh based on Delaunay triangulation first. en, the segments are extracted from the local Delaunay mesh. Afterward, two discs connected by each of the segments are called a front element, and some front elements are selected to place new disks, as shown in Figures 1 Similarly, the algorithm is extended from 2D to 3D for an accurate representation of convex polyhedral particles by a clump of spheres, as shown in Figure 2. In 3D, it requires that the surface of an irregular polyhedral particle is represented by a triangle mesh where a vertex (denoted by Vi) corresponds to a corner (denoted by Ci) surrounded by the facets containing the vertex, as shown in Figure 2(a). e algorithm in 3D also includes three steps: placement of spheres at the corners (Figure 2(b)), placement of spheres along the edges (Figure 2(c)), and placement of spheres on the facets (Figure 2(i)). All spheres placed in the three steps are clumped together into a clump in DEM simulations. For the sake of clarity, a tetrahedral particle is used to demonstrate each step of the algorithm, as shown in Figure 2.
In the first step, some spheres of the identical radius (r s ′ ) are placed at each corner based on the condition that a sphere placed at a corner needs to be tangent to any three of the facets surrounding the corner, whose central position (x s , y s , z s ) can be determined by the following equations: in which (n xi , n yi , n zi ) and (x pi , y pi , z pi ) are, respectively, the unit normal vector of a facet and the position of a point on the facet, and subscripts 1, 2, and 3 denote any three of the facets surrounding a corner. As shown in Figure 2(b), three spheres are placed at the corners C1, C2, and C3 according to the condition of placement of spheres at the corners.
In the second step, some spheres of identical radius r s ″ are placed one by one along each edge, as shown in Figure 2(c). is process is similar to that of the placement of discs along the edges. Taking an edge V1V2, for example, the process of placement of spheres along this edge is given as follows: the first sphere (the blue one shown in Figure 2(c)) is placed according to the condition of being tangent to the two facets (V1V2V3 and V1V4V2) connected by the edge V1V2 and the sphere placed at the corner C1 in the first step (the red one shown in Figure 2(c)) and then the second one is placed next to the first placed one according to the similar condition (the yellow one shown in Figure 2(c)). e process of placement of spheres along the edge V1V2 is completed until the placed sphere overlaps with the sphere located at the corner C2 in the first step, and it is replaced by a new one (the green one shown in Figure 2(c)) which is tangent to the two facets, the last-placed sphere, and the sphere placed at the corner C2 in the first step. For the spheres placed along the edge V1V2 (the blue and yellow ones in Figure 2(c)), their central positions (x s , y s , z s ) can be obtained according to the following equations: where (n xi , n yi , n zi ) and (x pi , y pi , z pi ) are, respectively, the unit normal vector of a facet and the position of a point on the facet, and subscripts 1 and 2 denote the two facets containing the same edge. Similarly, spheres are placed along 2 Advances in Civil Engineering other edges such as V2V3 and V1V3 by repeating the above operations. As shown in Figure 2(c), twenty-four spheres are, respectively, placed along each of the edges (V1V2, V2V3, and V1V3) in this example. As shown in Figures 2(d)-2(l), spheres of radius r ″ ′ s between [r min and r max ] are placed on each facet by a series of substeps based on a local Delaunay tessellation in the last step, where r min and r max are the minimum and maximum radii of the spheres to be placed. If r min is equal to r max , the radii of the placed sphere are identical. And the smaller the value of r min , the more spheres to be placed on each facet, which gives a more accurate representation of a polyhedral particle. is process is similar to that of placement of discs in the domain. Taking a facet V1V2V3, for example, the process of the placement of spheres on the facet is given as follows: a Delaunay mesh is created based on a set of points made up of the centers of the spheres including the ones placed at the corners V1, V2, and V3 in the first step and placed along the three edges of the facet in the second step, as shown in Figure 2(d); the segments are extracted from the Delaunay mesh, and similarly, two spheres connected by each of them are taken as a front element, as shown in Figure 2(d); some rational front elements are selected based on a constraint (it will be presented later), and they are used to place new spheres on the facet according to the condition of that a newly placed sphere is tangent to both spheres of the corresponding front element and the facet V1V2V3, as shown in Figure 2(e). If the central positions and radii of two spheres of a front element are, respectively, represented by (x s1 , y s1 , z s1 ), (x s2 , y s2 , z s2 ) , and (r s1 , r s2 ), the central position of a new sphere to be placed based on the front element (x s , y s , z s ) can be determined according to the following equations: where (n x1 , n y1 , n z1 ) and (x p1 , y p1 , z p1 ) are, respectively, a unit normal vector of a facet and the position of a point on the facet. A new sphere is placed on the facet if it does not overlap with the existing spheres placed on the facet and located inside the polyhedral particle. e first substep is finished when all rational segments are used to attempt to place new spheres on the facet, as shown in Figure 2(e). en, the centers of the newly placed spheres are also added into the set of points, and the next substep will be performed by repeating the operations, as shown in Figures 2(i)-2(k). e process of placement of spheres on the facet V1V2V3 is completed until no new spheres can be placed on the facet in the current substep, as shown in Figure 2(l).  As shown in Figure 2(f ), some unexpected segments extracted from a Delaunay mesh may exist. However, a new sphere cannot be successfully placed on the facet based on a front element corresponding to a long segment. erefore, a constraint for a front element is suggested to save the calculating time. A front element is called a rational one and used to place a new sphere on the facet only if the length of its corresponding segment satisfies the following relation: where r i and r j are the radii of two spheres of a front element, respectively, and L ij is the length of the segment corresponding to the front element. Another constraint proposed by Liu et al. [26] is considered in the third step to save the calculating time, which focuses on ones which are selected from the existing spheres for the generation of a Delaunay mesh in the next substep. Similarly, a sequence number indicating the sequence of substeps is assigned to spheres placed on the facet in each substep. For example, the sequence number of all spheres placed on a facet in the nth substep is marked as n, and the spheres have the same sequence and are marked by the same color if they are placed in the same substep, as shown in Figures 2(h) and 2(i). Suppose that l n represents the sequence number of the spheres placed in the current substep; then, the center of an existing sphere placed in previous substeps would be regarded as a member of a set of points for the next generation of a Delaunay mesh if the sequence number of the sphere l satisfies as follows: where h called an advancing-front thickness in Ref. [26] is an integer number greater than 0. If h is set to a larger number, it will require more calculating time to complete the processes of placements of spheres on a facet because all the placed spheres are selected for the generation of a Delaunay mesh in the next substep. On the other hand, if h is set to 1, only the spheres placed in the current substep are selected for the generation of a Delaunay mesh in the next substep. In this case, the algorithm may be not convergent, especially for complex polyhedral particles. erefore, it is a good strategy that the spheres placed in the last several substeps are only selected for generation of a Delaunay mesh in the next substep, and the generated Delaunay mesh is called a local Delaunay mesh, as shown in Figures 2(g) and 2(h).

Parameters and Comparation.
As presented in the last section, the spheres are placed based on different conditions in the three steps, and their central positions can be obtained easily by solving the corresponding equations such as equations (1)-(3) by the method presented in Ref. [27]. Determination of the radii of the spheres placed in each step and the parameters of the algorithm are introduced as follows.
In the first two steps, r s ′ and r s ″ denote the radius of a sphere to be placed at a corner and along an edge, respectively. In general, there is no relation between r s ′ and r s ″ . For simplicity, a relation of r s ′ � r s ″ is adopted here, which means that the radii of all spheres placed at the corners and along the edges are identical. In addition, a parameter N s , denoting the number of spheres placed along the shortest edge, is proposed to determine r s ′ and r s ″ indirectly. For example, for a polyhedral particle, the length of its shortest edge (L min ) are determined by comparing those of all its edges. en, r s ′ and r s ″ are determined based on the relation of r s ′ � r s ″ � L min /2N s when N s is given by users. In the third step, r ″ ′ s denotes the radius of a sphere to be placed on a facet, which ranges randomly and uniformly in an interval of [r min , r max ]. Two parameters (f r and r min /r max ) are proposed to determine the interval of [r min , r max ]. f r is defined as a ratio of r min to r s ′ or r s ″ , and r min /r max is defined as a ratio of r min to r max . e relation between the radii of the spheres placed in the three steps are established by these two parameters, which is easy and convenient in use. In addition, the advancing-front thickness (h) is a key parameter for the efficiency and convergence of the algorithm, and it is suggested based on our experience that h is between 5 and 7.
As introduced previously, the algorithm only includes four input parameters including N s , f r , r min /r max , and h for clump generation of a polyhedral particle. And it is also very convenient for users to generate the clumps of different numbers of spheres by adjusting N s . As shown in Figure 3, two clumps of different numbers of spheres for a polyhedral particle are generated by the algorithm when N s is, respectively, 5 and 10. It can be clearly seen from Figure 3 that the clumps generated by the algorithm can give a good approximation to the particle shape. Table 1 lists all parameters of the algorithm for the polyhedral particle.
In order to show the validity of the algorithm for the representation of convex polyhedral particles, the comparation between the proposed algorithm and the algorithm in PFC 5.0 is conducted. By contrast, it can be seen from Figures 3 and 4 that the clump generated by the proposed algorithm gives a better approximation to the particle shape.

Mass Properties of a Clump.
A polyhedral particle is represented by a clump consisting of multiple spheres as a rigid body in the discrete element modelling. e shape representation of the clump is only used for contact judgment in DEM simulations, while its motion is dominated by the equations of translational and rotational motions, which are determined by its mass properties including the total mass (m), center of mass (x), and inertia tensor (I 3×3 ). To obtain accurate simulation results, the mass properties of the clump should give a precise approximation to that of the polyhedral particle. In this work, the mass properties of a clump are determined based on the corresponding polyhedral particle rather than the spheres, while the spheres are only used for shape representation. By separating the mass properties and shape representation, the clump can give an accurate representation of a convex polyhedral particle in terms of particle shape and mass properties.
For a polyhedral particle with uniform mass distribution, its basic mass properties including m, x, and I 3×3 can be determined easily by the method proposed in Refs. [28,29] and are regarded as those of the clump to represent the polyhedral particle.

Examples
In the last section, an example is used to show the capability of the algorithm to generate a clump of multispheres for a polyhedral particle, as shown in Figure 3. In this section, another example is used to validate the robust and automatic performance of the algorithm in generating the clumps for an assembly of polyhedral particles. Figure 5 shows a flowchart of clump generation for an assembly of polyhedral particles. In this example, each convex polyhedral particle in the assembly is constructed based on the convex hull of several points randomly selected on the surface of a sphere, as shown in Figure 6(a). As shown in Figure 6(    Input an assembly of convex polyhedral particles

Loop each polyhedral particle
Loop all corners and place the spheres at each corner base on equation (1) Loop all edges and place the spheres along each edge base on equation (2) Loop all facets and place the spheres on each facet base on equation (3) to equation (5) Export the model into AutoCAD No Yes Figure 5: A flowchart of clump generation for an assembly of convex polyhedral particles by the proposed algorithm. 6 Advances in Civil Engineering particles constructed based on different spheres. Based on the assembly of polyhedral particles, the clump for each polyhedral particle is generated by the algorithm based on the parameters listed in Table 1 for Clump-2. For DEM simulation later, an interface is also implemented into the algorithm to facilitate users to export the clumps into PFC3D [14]. Figure 6(c) shows the clumps imported in PFC3D, and Figure 6(d) shows the clumps after settlement under gravity in PFC3D.

Conclusions
(1) A novel algorithm was extended from 2D to 3D to accurately represent convex polyhedral particles based on the multisphere method for DEM modelling in this work. In 3D, the algorithm includes three steps, and all spheres placed in the three steps are clumped together into a clump in DEM simulations. Clumps can give a better approximation to the particle shape because the spheres are placed tangent to the particle surface in each step. On the other hand, the mass properties of clumps are determined based on the corresponding polyhedral particles to obtain accurate simulation results. (2) e algorithm includes four input parameters, where the advancing-front thickness (h) is a key parameter for the efficiency and convergence of the algorithm. Based on our experience, it is suggested that h is between 5 and 7 for complex polyhedral particles.
(3) A current FORTRAN version of the algorithm is available upon request for academic use, and an interface is provided to import the sphere clump model into PFC3D. It is quite easy and convenient for users to generate sphere clump models by the algorithm for DEM modelling.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest. Advances in Civil Engineering 7