# Multisphere Representation of Convex Polyhedral Particles for DEM Simulation

Authors: Qiang Zhang; Chaojun Jia; Jun Yu; Jiawen Zhou
Journal: Advances in Civil Engineering
 (2021)

---

## Abstract

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.


---

## Body

## 1. 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 time-saving. 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. Therefore, 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. The new primitive methods mainly lie in the direct implementation of nonspherical primitives such as ellipsoids [5–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–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.


## 2. Representation of Convex Polyhedral Particles

### 2.1. Algorithm

The 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. Then, 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(d)–1(f).Figure 1
Schematic operations of the algorithm to represent a star-like particle: (a) placement of discs at corners; (b) placement of discs along edges; (c) placement of discs in the domain; (d), (e), and (f) the first, thirteenth, and nineteenth substeps in the third steps.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 Figure2. 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). The 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.Figure 2
A flow chart of the algorithm for clump generation for a polyhedral particle:(a) the tetrahedral particle; (b) placement of spheres at corners; (c) placement of spheres along edges; (d), (e), (f), (g), (h), and (i) placement of spheres on a facet.In the first step, some spheres of the identical radius (rs′) 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 xs,ys,zs can be determined by the following equations:(1)nx1xs−xp1+ny1ys−yp2+nz1zs−zp3=rs′,nx2xs−xp2+ny1ys−yp2+nz1zs−zp3=rs′,nx2xs−xp2+ny1ys−yp2+nz1zs−zp3=rs′,in which nxi,nyi,nzi and xpi,ypi,zpi 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 radiusrs″ are placed one by one along each edge, as shown in Figure 2(c). This 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)). The 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 xs,ys,zs can be obtained according to the following equations:(2)xs−xs12+ys−ys12+ys−ys12=r″+r″2,nx1xs−xp1+ny1ys−yp1+nz1zs−zp1=r″,nx2xs−xp2+ny2ys−yp2+nz2zs−zp2=r″,where nxi,nyi,nzi and xpi,ypi,zpi 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 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 Figures2(d)–2(l), spheres of radius rs″′ between rminandrmax are placed on each facet by a series of substeps based on a local Delaunay tessellation in the last step, where rmin and rmax are the minimum and maximum radii of the spheres to be placed. If rmin is equal to rmax, the radii of the placed sphere are identical. And the smaller the value of rmin, the more spheres to be placed on each facet, which gives a more accurate representation of a polyhedral particle. This 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 xs1,ys1,zs1, xs2,ys2,zs2, and (rs1, rs2), the central position of a new sphere to be placed based on the front element xs,ys,zs can be determined according to the following equations:(3)xs−xs12+ys−ys12+ys−ys12=r‴+rs12,xs−xs22+ys−ys22+zs−zs22=r‴+rs22,nx1xs−xp1+ny1ys−yp1+nz1zs−zp1=r‴,where nx1,ny1,nz1 and xp1,yp1,zp1 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. The 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). Then, 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). The 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 Figure2(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. Therefore, 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:(4)Lij−ri−rj≤2rmax,where ri and rj are the radii of two spheres of a front element, respectively, and Lij 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 ln 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:(5)l≥ln−h,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. Therefore, 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).


### 2.2. 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,rs′ and rs″ denote the radius of a sphere to be placed at a corner and along an edge, respectively. In general, there is no relation between rs′ and rs″. For simplicity, a relation of rs′=rs″ 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 Ns, denoting the number of spheres placed along the shortest edge, is proposed to determine rs′ and rs″ indirectly. For example, for a polyhedral particle, the length of its shortest edge (Lmin) are determined by comparing those of all its edges. Then, rs′ and rs″ are determined based on the relation of rs′=rs″=Lmin/2Ns when Ns is given by users.In the third step,rs″′ denotes the radius of a sphere to be placed on a facet, which ranges randomly and uniformly in an interval of rmin,rmax. Two parameters (fr and rmin/rmax) are proposed to determine the interval of rmin,rmax. fr is defined as a ratio of rmin to rs′ or rs″, and rmin/rmax is defined as a ratio of rmin to rmax. The 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 includingNs, fr, rmin/rmax, 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 Ns. As shown in Figure 3, two clumps of different numbers of spheres for a polyhedral particle are generated by the algorithm when Ns 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.Figure 3
Representation of a polyhedral particle by a clump of spheres by the proposed algorithm: (a) the polyhedral particle; (b) Cluster-1; (c) Cluster-2.
(a)(b)(c)Table 1
Parameters of the algorithm to represent a polyhedral particle by a clump of spheres.
ClumpsNsfrrmax/rminhClump-1150.51.56Clump-2100.51.56In 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 Figures3 and 4 that the clump generated by the proposed algorithm gives a better approximation to the particle shape.Figure 4
Representation of a polyhedral particle by a clump of spheres by the algorithm in PFC5.0: (a) and (b) the polyhedral particle and the clump; (c) the clump.
(a)(b)(c)


### 2.3. 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. The 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 (I3×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 includingm, x, and I3×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.


## 2.1. Algorithm

The 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. Then, 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(d)–1(f).Figure 1
Schematic operations of the algorithm to represent a star-like particle: (a) placement of discs at corners; (b) placement of discs along edges; (c) placement of discs in the domain; (d), (e), and (f) the first, thirteenth, and nineteenth substeps in the third steps.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 Figure2. 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). The 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.Figure 2
A flow chart of the algorithm for clump generation for a polyhedral particle:(a) the tetrahedral particle; (b) placement of spheres at corners; (c) placement of spheres along edges; (d), (e), (f), (g), (h), and (i) placement of spheres on a facet.In the first step, some spheres of the identical radius (rs′) 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 xs,ys,zs can be determined by the following equations:(1)nx1xs−xp1+ny1ys−yp2+nz1zs−zp3=rs′,nx2xs−xp2+ny1ys−yp2+nz1zs−zp3=rs′,nx2xs−xp2+ny1ys−yp2+nz1zs−zp3=rs′,in which nxi,nyi,nzi and xpi,ypi,zpi 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 radiusrs″ are placed one by one along each edge, as shown in Figure 2(c). This 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)). The 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 xs,ys,zs can be obtained according to the following equations:(2)xs−xs12+ys−ys12+ys−ys12=r″+r″2,nx1xs−xp1+ny1ys−yp1+nz1zs−zp1=r″,nx2xs−xp2+ny2ys−yp2+nz2zs−zp2=r″,where nxi,nyi,nzi and xpi,ypi,zpi 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 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 Figures2(d)–2(l), spheres of radius rs″′ between rminandrmax are placed on each facet by a series of substeps based on a local Delaunay tessellation in the last step, where rmin and rmax are the minimum and maximum radii of the spheres to be placed. If rmin is equal to rmax, the radii of the placed sphere are identical. And the smaller the value of rmin, the more spheres to be placed on each facet, which gives a more accurate representation of a polyhedral particle. This 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 xs1,ys1,zs1, xs2,ys2,zs2, and (rs1, rs2), the central position of a new sphere to be placed based on the front element xs,ys,zs can be determined according to the following equations:(3)xs−xs12+ys−ys12+ys−ys12=r‴+rs12,xs−xs22+ys−ys22+zs−zs22=r‴+rs22,nx1xs−xp1+ny1ys−yp1+nz1zs−zp1=r‴,where nx1,ny1,nz1 and xp1,yp1,zp1 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. The 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). Then, 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). The 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 Figure2(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. Therefore, 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:(4)Lij−ri−rj≤2rmax,where ri and rj are the radii of two spheres of a front element, respectively, and Lij 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 ln 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:(5)l≥ln−h,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. Therefore, 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).


## 2.2. 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,rs′ and rs″ denote the radius of a sphere to be placed at a corner and along an edge, respectively. In general, there is no relation between rs′ and rs″. For simplicity, a relation of rs′=rs″ 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 Ns, denoting the number of spheres placed along the shortest edge, is proposed to determine rs′ and rs″ indirectly. For example, for a polyhedral particle, the length of its shortest edge (Lmin) are determined by comparing those of all its edges. Then, rs′ and rs″ are determined based on the relation of rs′=rs″=Lmin/2Ns when Ns is given by users.In the third step,rs″′ denotes the radius of a sphere to be placed on a facet, which ranges randomly and uniformly in an interval of rmin,rmax. Two parameters (fr and rmin/rmax) are proposed to determine the interval of rmin,rmax. fr is defined as a ratio of rmin to rs′ or rs″, and rmin/rmax is defined as a ratio of rmin to rmax. The 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 includingNs, fr, rmin/rmax, 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 Ns. As shown in Figure 3, two clumps of different numbers of spheres for a polyhedral particle are generated by the algorithm when Ns 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.Figure 3
Representation of a polyhedral particle by a clump of spheres by the proposed algorithm: (a) the polyhedral particle; (b) Cluster-1; (c) Cluster-2.
(a)(b)(c)Table 1
Parameters of the algorithm to represent a polyhedral particle by a clump of spheres.
ClumpsNsfrrmax/rminhClump-1150.51.56Clump-2100.51.56In 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 Figures3 and 4 that the clump generated by the proposed algorithm gives a better approximation to the particle shape.Figure 4
Representation of a polyhedral particle by a clump of spheres by the algorithm in PFC5.0: (a) and (b) the polyhedral particle and the clump; (c) the clump.
(a)(b)(c)


## 2.3. 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. The 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 (I3×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 includingm, x, and I3×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.


## 3. 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 Figure3. 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(b), the assembly includes 366 polyhedral 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.Figure 5
A flowchart of clump generation for an assembly of convex polyhedral particles by the proposed algorithm.Figure 6
Generation of a sphere clump model for an assembly of polyhedral particles: (a) construction of a polyhedral particle based on a sphere; (b) an assembly of polyhedral particles constructed based on different spheres; (c) the corresponding sphere clump model imported in PFC3D; (d) the sphere clump model after settlement under gravity.
(a)(b)(c)(d)


## 4. 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)
The 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.



---

*Source: 8846004-2021-02-02.xml*