A Random Angular Bend Algorithm for Two- Dimensional Discrete Modeling of Granular Materials

Generation of particles with irregular shape and the overlap detection are crucial for numerical simulation of granular materials. This paper presents a systematic approach to develop a two-dimensional random particle model for numerical simulation of granular materials. Firstly, a random angular bend (RAB) algorithm is proposed and coded in Python to simulate the geometric model of individual particle with irregular shape. Three representative parameters are used to quantitatively control the shape feature of generated polygons in terms of three major aspects, respectively. Then, the generated geometrical models are implemented into particle flow code PFC2D to construct the clump library. The clumps are created via the mid-surface method. Besides, an overlap detection algorithm is developed to address the difficulties associated with spatial allocation of irregularly shaped particles. Finally, two application examples are adopted to validate the feasibility of the proposed algorithm in the numerical modeling of realistic granular materials. The study provides a solid foundation for the generation and simulation of the granular materials based on angular bend theory.


Introduction
Granular materials are conglomerations of discrete solids [1], which are widely used in various engineering fields, such as geotechnical engineering, mining engineering, and food and pharmaceutical industries. Due to the discontinuous nature and irregularity of particle shape, the study of properties of granular materials remains an open issue. In the study of granular material behavior, some experimental tests were conducted to analyze the effect of particle shape and size on the mechanical and physical properties of granular materials [2,3]. However, it is difficult to precisely control particle shape and size in these tests. Besides, conventional experimental tests cannot provide sufficient microscale information of samples, such as interparticle forces, particle velocities, etc.
Overcoming the disadvantage of model test in these respects, the numerical analysis method is commonly accepted by many researchers in the last decade. The discrete element method (DEM), originally developed by Cundall in 1971 [4,5], provides a convenient way to obtain particle information on multiple scales, and gradually becomes an effective tool to study the mechanical behaviors of granular materials [6][7][8][9]. To investigate the effect of particle shape, various shapes have been created to approximate the realistic shape of particles (such as circles and spheres [10][11][12][13], ellipses and ellipsoids [14,15], cylinders [16], polygons, and polyhedrons [17][18][19]). However, a considerable number of them are idealized geometric shape, which is unable to capture the behavior of actual or natural granular materials better.
To precisely generate arbitrary-shaped particle in two-dimensional simulation, many methods have been proposed in the literature. Typical approaches include R(θ) method [20], Voronoi grain-based method [21,22], digital image-based method [23], image-based clump library method [24,25], Fourier-based method [26][27][28], etc. However, the R(θ) method is only applicable for star-like particles [29,30]. The Voronoi grain-based method has limited capability in precisely controlling generated particle shape and size. During the modeling of granular materials based on digital image technology, it is difficult to obtain digital images of the whole research area or obtain the digital images of idealized section for specified research purpose. Image-based clump library method overcomes the disadvantage of above method and provides a new method to create virtual specimens possessing various shapes and particle gradations. However, tremendous efforts are required to develop the clump library (a clump is a rigid collection of n rigid spherical pebbles, which are unit-thickness disks in PFC 2D ). Among various Fourier-based approaches, the shape function R(θ) based method is widely used in shape analysis and generation. However, this method is only applicable for star-like particles. For nonstar-like particles, the radial line may have multiple intersections with the boundary of particle contour [30]. Compared with the methods above, the shape function ϕ(l) (Z-R shape function) based method which is developed by Zahn and Roskies (1972) has distinct advantages in precisely representing the realistic particles [31]. However, there is limited research on the random generation of irregular particles based on this theory.
Inspired by the theory of Z-R shape function, this paper presents a novel method to generate two-dimensional random particle model for the numerical simulation of granular materials. The flowchart of this paper is shown in Figure 1. Firstly, a random angular bend (RAB) algorithm is proposed and coded in Python [32] to simulate the geometric model of individual particle with irregular shape. Three representative parameters are used to quantitatively control the shape feature of generated polygons in term of three major aspects, namely, the form, roundness, and surface textures, respectively. The geometric models are imported into PFC 2D [33] to generate realistic granular media. Clumps consisting of pebbles are created based on the mid-surface method [34] and used to construct clump library. An overlap detection algorithm is developed to address difficulties associated with spatial allocation of irregularly shaped particles. Finally, two application examples are further employed to show the feasibility of proposed algorithm for the numerical modeling of realistic granular materials. considerable number of them are idealized geometric shape, which is unable to capture the behavior of actual or natural granular materials better.
To precisely generate arbitrary-shaped particle in two-dimensional simulation, many methods have been proposed in the literature. Typical approaches include R(θ) method [20], Voronoi grainbased method [21,22], digital image-based method [23], image-based clump library method [24,25], Fourier-based method [26][27][28], etc. However, the R(θ) method is only applicable for star-like particles [29,30]. The Voronoi grain-based method has limited capability in precisely controlling generated particle shape and size. During the modeling of granular materials based on digital image technology, it is difficult to obtain digital images of the whole research area or obtain the digital images of idealized section for specified research purpose. Image-based clump library method overcomes the disadvantage of above method and provides a new method to create virtual specimens possessing various shapes and particle gradations. However, tremendous efforts are required to develop the clump library (a clump is a rigid collection of n rigid spherical pebbles, which are unitthickness disks in PFC 2D ). Among various Fourier-based approaches, the shape function R(θ) based method is widely used in shape analysis and generation. However, this method is only applicable for star-like particles. For nonstar-like particles, the radial line may have multiple intersections with the boundary of particle contour [30]. Compared with the methods above, the shape function φ(l) (Z-R shape function) based method which is developed by Zahn and Roskies (1972) has distinct advantages in precisely representing the realistic particles [31]. However, there is limited research on the random generation of irregular particles based on this theory.
Inspired by the theory of Z-R shape function, this paper presents a novel method to generate two-dimensional random particle model for the numerical simulation of granular materials. The flowchart of this paper is shown in Figure 1. Firstly, a random angular bend (RAB) algorithm is proposed and coded in Python [32] to simulate the geometric model of individual particle with irregular shape. Three representative parameters are used to quantitatively control the shape feature of generated polygons in term of three major aspects, namely, the form, roundness, and surface textures, respectively. The geometric models are imported into PFC 2D [33] to generate realistic granular media. Clumps consisting of pebbles are created based on the mid-surface method [34] and used to construct clump library. An overlap detection algorithm is developed to address difficulties associated with spatial allocation of irregularly shaped particles. Finally, two application examples are further employed to show the feasibility of proposed algorithm for the numerical modeling of realistic granular materials.

Z-R Shape Function
Shape function ϕ(l) was introduced by Zahn and Roskies (1972) to represent the shape of closed polygonal curve in term of arc length l and cumulative angular bends ϕ. To describe the shape of a polygon, which is made up of n vertices P 0 , P 1 , . . . , P n with angular bends ∆ϕ 0 , ∆ϕ 1 , . . . , ∆ϕ n and edge lengths ∆l 0 , ∆l 1 , . . . , ∆l n (as shown in Figure 2), the first thing that need to be done is to trace clockwise around the outline of the polygon from some initial starting point to collect Cartesian (x, y) coordinates [35]. Once the coordinates of all the vertices are determined, they can be used to calculate edge lengths ∆l i and angular bend ∆ϕ i . The clockwise arc lengths l along the curve from starting point P 0 up to the ith vertex P i is given as Equation (1) Next step is to determine the angle bend ∆ϕ i between adjacent polygon segments. Taking the adjacent polygon segments as a vector, the magnitude and orientation of the angle can be determined based on vector operation. For three adjacent point, P i-1 (x i-1 ,y i-1 ), P i (x i , y i ) and P i+1 (x i+1 ,y i+1 ), the magnitude of ∆ϕ i can be calculated using the dot products of vector where ∆ϕ i is the absolute value of ∆ϕ i . The polarity of ∆ϕ i depends upon the sign of the cross products of the two vectors, which is given as the determinant of a matrix S i [36].
Once the sign of this cross product is determined, the value of ∆ϕ i can be obtained by applying the following rules.
where |S i | is the absolute value of S i . It should be noted that a cross product of 0 means that vector → P i−1 P i and → P i P i+1 . are collinear. The cumulative angular function ϕ(l), also known as Z-R shape function, is the summation of all the angular bends from vertex P 0 up to the vertex P i (corresponding to arc length l). This function is given by ϕ(l) = i k=1 ∆ϕ k . If the outline is traced clockwise and is not spiral, the value of ϕ(L all ) equal to -2π, where L all is the total length of the polygon. Throughout the above procedures described in Equations (1)-(5), the coordinates of points on the polygon are converted from their Cartesian (x,y) form to the ϕ(l) form of Z-R shape function [37]. Shape analysis can be conducted based on the result of Z-R shape function [31] (Figure 2). Materials 2019, 12, x FOR PEER REVIEW 4 of 17 Figure 2. Description of Particle shape using the result of the Z-R shape function.

Generation of Irregularly Shaped Polygon
Inspired by the theory of Z-R shape function, the RAB algorithm is proposed in this section and coded in Python. This method is the inverse operation and modification of the theory of Z-R shape function. As shown in Figure 3, the algorithm can be described as follows. 1. Define a random point O as the center and a point P1 (x1, y1) as initial starting point of the polygon to be generated. 2. Determine the position of the sequential points to be generated. To determine the position of sequential points, a series of angular turns need to be conducted. The common approach to solve this problem is to treat the polygon segments as vectors. The counterclockwise angle that vector P i P i+1 ⃗ makes with the positive x-axis is βi. For i=1, βi is a random number varying uniformly between 0 and 2π. For i >1, the value of this angel can be calculated by Equation (6) βi = βi-1 + Δφi (i >1) (6) where Δφi is the angle bend at point Pi, which is treated as a uniformly distributed random variable. In the procedures of generating randomly shaped polygon, the orientation of the angle should be taken into account and are randomly assigned. For i >1, the coordinates of Pi in a Cartesian coordinate system can be represented by Equation (7) where lstep is constant step length. 3. To ensure the points on the polygon are generated counterclockwise, the point Pi should be at the left side of line OPi-1 or on the line (for example, point P5 in Figure 3). This can be validated using the cross product of vector OP i ⃗ and OP i-1 ⃗ .The calculation of the cross product is conducted through Equation (3). If this condition is satisfied, proceed with step 4; otherwise, repeat from step 2. 4. Calculate DOPi -the distance between O and the generated point Pi(xi, yi). There are two parameters required to be inputted in advance in this part. They are the minimum distance Dmin and the maximum distance Dmax between the center O and the generated points. If the distance is within the specified range (Dmin ≤ DOPi ≤ Dmax) then the point is stored as a new vertex. Otherwise repeat from step 2. 5. If a point cannot be generated within specified times of trials, then remove the vertex stored in the previous step and repeat from step 2. 6. The summation of all the angular bends from the initial starting point and the point generated in the previous step is φ (φ=∑Δφi). In case of the generation of spiral polygon (a spiral polygon is a simple polygon whose boundary chain contains exactly one concave subchain [38]), φ should equal to 2π [39]. Consequently, repeat steps 2, 3, and 4 until the value of φ is larger than or equal to 2π.

Generation of Irregularly Shaped Polygon
Inspired by the theory of Z-R shape function, the RAB algorithm is proposed in this section and coded in Python. This method is the inverse operation and modification of the theory of Z-R shape function. As shown in Figure 3, the algorithm can be described as follows.

1.
Define a random point O as the center and a point P 1 (x 1 , y 1 ) as initial starting point of the polygon to be generated.

2.
Determine the position of the sequential points to be generated. To determine the position of sequential points, a series of angular turns need to be conducted. The common approach to solve this problem is to treat the polygon segments as vectors. The counterclockwise angle that vector → P i P i+1 makes with the positive x-axis is β i . For i=1, β i is a random number varying uniformly between 0 and 2π. For i >1, the value of this angel can be calculated by Equation (6) where ∆ϕ i is the angle bend at point P i , which is treated as a uniformly distributed random variable. In the procedures of generating randomly shaped polygon, the orientation of the angle should be taken into account and are randomly assigned. For i >1, the coordinates of P i in a Cartesian coordinate system can be represented by Equation (7) x where l step is constant step length.

3.
To ensure the points on the polygon are generated counterclockwise, the point P i should be at the left side of line OP i-1 or on the line (for example, point P 5 in Figure 3). This can be validated using the cross product of vector The calculation of the cross product is conducted through Equation (3). If this condition is satisfied, proceed with step 4; otherwise, repeat from step 2.

4.
Calculate D OPi -the distance between O and the generated point P i (x i , y i ). There are two parameters required to be inputted in advance in this part. They are the minimum distance D min and the maximum distance D max between the center O and the generated points. If the distance is within the specified range (D min ≤ D OPi ≤ D max ) then the point is stored as a new vertex. Otherwise repeat from step 2.

5.
If a point cannot be generated within specified times of trials, then remove the vertex stored in the previous step and repeat from step 2. The summation of all the angular bends from the initial starting point and the point generated in the previous step is ϕ (ϕ= ∆ϕ i ). In case of the generation of spiral polygon (a spiral polygon is a simple polygon whose boundary chain contains exactly one concave subchain [38]), ϕ should equal to 2π [39]. Consequently, repeat steps 2, 3, and 4 until the value of ϕ is larger than or equal to 2π.
If the value of ϕ satisfies the convergence criteria in step 6, then delete the last point and connect the stored vertexes counterclockwise. Thus, individual polygon with irregular shape is successfully generated. If necessary, the above steps are repeated until the specified amount of polygon are obtained.
It should be noted that the distance between point P n and P 1 may not equal to l step . However, this situation does not have a significant influence on the polygon contour. If the value of φ satisfies the convergence criteria in step 6, then delete the last point and connect the stored vertexes counterclockwise. Thus, individual polygon with irregular shape is successfully generated. If necessary, the above steps are repeated until the specified amount of polygon are obtained. It should be noted that the distance between point Pₙ and P₁ may not equal to lstep. However, this situation does not have a significant influence on the polygon contour.

Shape Descriptors
Many attempts have been made to quantitatively analyze the geometrical characteristics of particles in the literature. Based on the dimensional scales, the morphological characteristics of these particles are expressed in terms of three major aspects: form (overall shape), roundness (angularity), and surface texture (roughness) [26,40,41].
Form, the property of the first level, represents spatial irregularities of particle shape in the large dimension. The common shape include circle, ellipse, rectangle, etc. Shape descriptors, such as elongation index (EI), aspect ratio (AR), and flatness index (FI) are used to quantify the basic shape characteristic of particle. The EI can be determined as follows Equation 8 [22] EI = S / L (8) Where S and L equal to the width and length of the smallest rectangular box containing the particle, respectively ( Figure 4a). The elongation index varies from 0 to 1, where the higher value of EI indicates a lower degree of elongation. For a circle, the EI equals to 1. Particles with different elongation index are shown in Figure 5.

Shape Descriptors
Many attempts have been made to quantitatively analyze the geometrical characteristics of particles in the literature. Based on the dimensional scales, the morphological characteristics of these particles are expressed in terms of three major aspects: form (overall shape), roundness (angularity), and surface texture (roughness) [26,40,41].
Form, the property of the first level, represents spatial irregularities of particle shape in the large dimension. The common shape include circle, ellipse, rectangle, etc. Shape descriptors, such as elongation index (EI), aspect ratio (AR), and flatness index (FI) are used to quantify the basic shape characteristic of particle. The EI can be determined as follows Equation (8) where S and L equal to the width and length of the smallest rectangular box containing the particle, respectively (Figure 4a). The elongation index varies from 0 to 1, where the higher value of EI indicates a lower degree of elongation. For a circle, the EI equals to 1. Particles with different elongation index are shown in Figure 5.
elongation index (EI), aspect ratio (AR), and flatness index (FI) are used to quantify the basic shape characteristic of particle. The EI can be determined as follows Equation 8 [22] EI = S / L (8) Where S and L equal to the width and length of the smallest rectangular box containing the particle, respectively (Figure 4a). The elongation index varies from 0 to 1, where the higher value of EI indicates a lower degree of elongation. For a circle, the EI equals to 1. Particles with different elongation index are shown in Figure 5.  Roundness-the property of the second level-describes the variations of particle shape in the medium dimension, which reflects the average sharpness at corners. It is independent of the overall shape. Available shape descriptors including angularity index, roundness index (RI) and others are of this expression. The degree of roundness is defined and calculated via the following Equation 9 [42].
Where A is the projected area and P is the perimeter of the particle contour (see Figure 4b). Figure 6 displays six particles generated using the proposed algorithm. It is obvious that as RI increase from 0.58 to 0.92, high rounded particles can be obtained ( Figure 6). Roundness-the property of the second level-describes the variations of particle shape in the medium dimension, which reflects the average sharpness at corners. It is independent of the overall shape. Available shape descriptors including angularity index, roundness index (RI) and others are of this expression. The degree of roundness is defined and calculated via the following Equation (9) [42].
where A is the projected area and P is the perimeter of the particle contour (see Figure 4b). Figure 6 displays six particles generated using the proposed algorithm. It is obvious that as RI increase from 0.58 to 0.92, high rounded particles can be obtained ( Figure 6).
Where A is the projected area and P is the perimeter of the particle contour (see Figure 4b). Figure 6 displays six particles generated using the proposed algorithm. It is obvious that as RI increase from 0.58 to 0.92, high rounded particles can be obtained ( Figure 6).  Surface texture, the property of the third level, refers to the surface features which are small-scale relative to the size of the object [41]. Roughness index and regularity index are generally used to describe the surface texture features. The regularity index (RE) is defined as Equation (10) [22] RE = log P P − P conv (10) where P is the perimeter and P conv is the convex perimeter of the particle (see Figure 4b). The value of regularity greater than 3 or 4 indicated the particle surface is very smooth. The example particles shown in Figure 7 illustrate this point. Surface texture, the property of the third level, refers to the surface features which are smallscale relative to the size of the object [41]. Roughness index and regularity index are generally used to describe the surface texture features. The regularity index (RE) is defined as Equation 10 [22] RE= log P P-P conv (10) where P is the perimeter and Pconv is the convex perimeter of the particle (see Figure 4b). The value of regularity greater than 3 or 4 indicated the particle surface is very smooth. The example particles shown in Figure 7 illustrate this point. From Figure 5-7, it is obvious that three shape descriptors mentioned above-elongation index, roundness index and regularity index-have a good performance in obtaining the multiaspect characteristics of particle shape (form, roundness, and surface texture of a particle outline). So these three shape descriptors are chosen for the quantitative shape analysis in the following section.

Relationship between Input Parameters and Shape Descriptors
There are several parameters required to be inputted manually in RAB algorithm. They are the From Figures 5-7, it is obvious that three shape descriptors mentioned above-elongation index, roundness index and regularity index-have a good performance in obtaining the multiaspect Materials 2019, 12, 2169 8 of 17 characteristics of particle shape (form, roundness, and surface texture of a particle outline). So these three shape descriptors are chosen for the quantitative shape analysis in the following section.

Relationship between Input Parameters and Shape Descriptors
There are several parameters required to be inputted manually in RAB algorithm. They are the minimum distance D min and the maximum distance D max between the center O and the generated points, the range of angle bend ∆ϕ i and the step length l step (see Figure 3). In this paper, the minimum value of ∆ϕ i is set to 0 and the step length is treated as a constant. To better control the shape of generated particle, based on the geometric relationship between these parameters, these parameters are summarized and classed into three representative parameters, which are R dist -the ratio between D min and D max , R l/d -the ratio between l step and D max and the maximum value of angle bend ∆ϕ max (see Equations (11)-(13)).
To explore the influence of these three parameters on shape descriptors (EI, RI and RE), each of them is considered as a random variable, within a specified range. Based on the range of R dist , R l/d , ∆ϕ max , which are provided in Table 1, 1000 irregularly shaped polygons are generated for each group. The values of shape descriptors are obtained using Equations (8)-(10), respectively. Figure 8 presents the effect of the three representative parameters on particle shape. The influence can be determined through the dispersion degree of the data points, which is obtained by range analysis. The greater dispersion degree indicates the input parameter has marginal control effect on corresponding shape descriptors. By comparison of the graphs in the same column, combined with the dispersion degree of data point, the effect of three representative parameters on shape descriptors are investigated through qualitative analysis. Besides, some fitting curves are almost horizontal in Figure 8, indicating that the control effect of corresponding parameters on a certain shape descriptors is rather limited or can be neglected.  Figure 8a illustrates the influence of parameters ∆ϕ max on the generated particle shapes. It is obvious that ∆ϕ max is well correlated with the particle roundness and regularity. Nevertheless, an increase of ∆ϕ max has limited effect on the particle elongation. However, with the increase of R dist , the variation of shape descriptors is different with the above situation (see Figure 8b). Based on the results of shape analysis, the EI of generated shape is almost always greater than that of corresponding R dist . Thus, in the process of generating random particles, R dist can be used to control the minimum value of elongation, while R dist has limited influence on RI and RE. Finally, the discrete degree of the data points shown in Figure 8c reflects the influence of R l/d . It can be observed that this parameter have negligible influences on elongation and roundness. Nevertheless, increase of R l/d leads to a significant increase of particle regularity. variation of shape descriptors is different with the above situation (see Figure 8b). Based on the results of shape analysis, the EI of generated shape is almost always greater than that of corresponding Rdist. Thus, in the process of generating random particles, Rdist can be used to control the minimum value of elongation, while Rdist has limited influence on RI and RE. Finally, the discrete degree of the data points shown in Figure 8c reflects the influence of Rl/d. It can be observed that this parameter have negligible influences on elongation and roundness. Nevertheless, increase of Rl/d leads to a significant increase of particle regularity.  Figure 8. Influence of input parameters on shape descriptor: (a) Influence of Δφmax on particle shape; (b) Influence of Rdist on particle shape; (c) Influence of Rl/d on particle shape. Regularity R l/d Figure 8. Influence of input parameters on shape descriptor: (a) Influence of ∆ϕ max on particle shape; (b) Influence of R dist on particle shape; (c) Influence of R l/d on particle shape.
Based on the analysis above, the three shape descriptors are controlled by single or multiple input parameters with an obvious effect. The relationship between input parameter and shape descriptor is obtained by using the curve fitting method (Figure 9). The correlation coefficient R 2 is adapted to evaluate fitting precision [43,44]. Among these shape descriptors, the minimum and mean value of elongation (EI min and EI mean ) are determined by R dist (Figure 9a), which can be expressed as EI mean = 0.329R dist + 0.702 (15) Besides, as shown in Figure 9b, the ∆ϕ max has a statistically significant relationship with the mean value of roundness (RI mean ). The value is calculated based on the following fitting formula.
Different from EI and RI, the mean value of regularity (RE mean ) is controlled by both ∆ϕ max and R l/d (Figure 9c). The predictive formula could be obtained by the curve fitting method with R 2 =0.989 and expressed as follows RE mean = 3.43 + 23.99R l/d − 3.578∆ϕ max − 6.985 R 2 l/d − 11.57 R l/d · ∆ϕ max + 1.504∆ϕ 2 max (17) Thus, in the procedure of generating random particles, the shape of particles can be quantitatively controlled based on Equations (14)- (17).
Besides, as shown in Figure 9b, the Δφmax has a statistically significant relationship with the mean value of roundness (RImean). The value is calculated based on the following fitting formula. RImean = 0.353Δφmax -0.524 +0.285 (16) Different from EI and RI, the mean value of regularity (REmean) is controlled by both Δφmax and Rl/d (Figure 9c). The predictive formula could be obtained by the curve fitting method with R 2 =0.989 and expressed as follows REmean=3.43+23.99Rl/d-3.578Δφmax-6.985 R 2 l/d -11.57 Rl/d · Δφmax+ 1.504Δφ 2 max (17) Thus, in the procedure of generating random particles, the shape of particles can be quantitatively controlled based on Equations (14-17). Figure 9. Relationship between input parameters and shape descriptor: (a) Relationship between input parameters and EImean, (b) relationship between input parameters and RImean, and (c) relationship between input parameters and REmean.  Figure 9. Relationship between input parameters and shape descriptor: (a) Relationship between input parameters and EI mean , (b) relationship between input parameters and RI mean , and (c) relationship between input parameters and RE mean.

Overlap Detection Algorithm
After a series of arbitrarily shaped polygons are introduced, the generated geometric models are import into the Particle Flow Code PFC 2D (Version 5.0., Itasca Consulting Group Inc.: Minneapolis, MN, USA) [33] to construct clump library. PFC 2D is a two-dimensional DEM software. The particles are treated as disks or clumps, and can both translate and rotate. The movement of particle obeys Newton's laws of motion. Based on the mid-surface method, the clumps are created to represent realistic particle shapes [34] and appropriately packed in the model domain. To facilitate random and quick allocation of irregularly shaped clumps, an overlap detection algorithm is proposed in this section. This algorithm can be summarized as follows.

1.
Construct a clump library based on the randomly generated polygons.

2.
Select a clump from the library and specify the position randomly in the model domain. It is obvious that the shape of a clump is determined by the position and radii of the comprising pebbles. Based on the pebble information, the minimum distance between the clump and the boundary of model can be determined by calculating the distance (D i ) between each pebble and the model boundary (see Figure 10).

3.
Detect overlap between the clump and model boundary. If the minimum distance between them is greater than or equal to zero, it suggest that the clump does not intersect with the boundaries of model. If this condition is satisfied, store the clump information; otherwise, repeat from step 2.

4.
Put next clump in the model domain and calculate the minimum distance between this clump and the packed clumps (or the model boundary) based on the position and radii of all the pebbles. If this clump do not overlap with the packed clumps or exceed the model boundary, store the clump and pebble information; otherwise, repeat this step. 5.
Repeat step 4 until a given condition is satisfied (e.g., the number of all the clumps, the area of all the clumps).
The most important part of the algorithm is the calculation of the minimum distance between clumps or the distance between clump and the model boundary. This is an iterative procedure in which the minimum distance is determined by calculating the distance between pebbles in new selected clump and pebbles in the packed clumps (and the distance between pebbles in new selected clump and model boundary). As shown in Figure 10, the distance between clumps C1 and C2 is defined by Equation (18): where Dist_ij is the distance between pebble peb_i in clump C1 and pebble peb_j in C2. In a similar way, the distance between clump and model boundary can be calculated. Thus, in the procedures of clump packing, the distance between each clumps in the model domain can be precisely controlled.
Newton's laws of motion. Based on the mid-surface method, the clumps are created to represent realistic particle shapes [34] and appropriately packed in the model domain. To facilitate random and quick allocation of irregularly shaped clumps, an overlap detection algorithm is proposed in this section. This algorithm can be summarized as follows.
1. Construct a clump library based on the randomly generated polygons. 2. Select a clump from the library and specify the position randomly in the model domain. It is obvious that the shape of a clump is determined by the position and radii of the comprising pebbles. Based on the pebble information, the minimum distance between the clump and the boundary of model can be determined by calculating the distance (Di) between each pebble and the model boundary (see Figure 10). 3. Detect overlap between the clump and model boundary. If the minimum distance between them is greater than or equal to zero, it suggest that the clump does not intersect with the boundaries of model. If this condition is satisfied, store the clump information; otherwise, repeat from step 2. 4. Put next clump in the model domain and calculate the minimum distance between this clump and the packed clumps (or the model boundary) based on the position and radii of all the pebbles. If this clump do not overlap with the packed clumps or exceed the model boundary, store the clump and pebble information; otherwise, repeat this step. 5. Repeat step 4 until a given condition is satisfied (e.g., the number of all the clumps, the area of all the clumps).
The most important part of the algorithm is the calculation of the minimum distance between clumps or the distance between clump and the model boundary. This is an iterative procedure in which the minimum distance is determined by calculating the distance between pebbles in new selected clump and pebbles in the packed clumps (and the distance between pebbles in new selected clump and model boundary). As shown in Figure 10, the distance between clumps C1 and C2 is defined by Equation (18): Dist(C1,C2) = Minimum(Dist_11, Dist_12, Dist_13…, Dist_ij) (18) where Dist_ij is the distance between pebble peb_i in clump C1 and pebble peb_j in C2. In a similar way, the distance between clump and model boundary can be calculated. Thus, in the procedures of clump packing, the distance between each clumps in the model domain can be precisely controlled.

Application of Proposed Algorithm
To validate the feasibility of the proposed RAB and overlap detection algorithms in modeling of granular materials, two group of granular mixture samples for application examples are generated in the following sections.

Modeling of Granular Mixture Materials
Granular mixtures are inhomogeneous multiphase materials, which are consisted of coarse particles and fine particles [45]. In this numerical model, the coarse particles are simulated using clumps of overlap pebbles and fine particles by single disk. The procedures of modeling granular mixture are as follows (Figure 11).

1.
Code the RAB algorithm in Python and input parameters R dist , R l/d , and ∆ϕ max to simulate the geometric model of individual coarse particle with irregular shape.

2.
Import the generated geometric model into PFC 2D to construct clump library. The clumps are created based on mid-surface method. These clumps are used to represent realistic coarse particle.

3.
Put coarse particles into the model domain based on the overlap detection algorithm which is proposed above (Figure 11a).

4.
Generate sample that is composed of single disk to simulate fine particles in the same model domain (Figure 11b). The size of particle and the porosity of fine particles specimen should be in a given range.

5.
Remove fine particles that overlap with coarse particles (Figure 11c,d). The generated granular mixture sample is shown in Figure 11e.

Application of Proposed Algorithm
To validate the feasibility of the proposed RAB and overlap detection algorithms in modeling of granular materials, two group of granular mixture samples for application examples are generated in the following sections.

Modeling of Granular Mixture Materials
Granular mixtures are inhomogeneous multiphase materials, which are consisted of coarse particles and fine particles [45] . In this numerical model, the coarse particles are simulated using clumps of overlap pebbles and fine particles by single disk. The procedures of modeling granular mixture are as follows (Figure 11). 1. Code the RAB algorithm in Python and input parameters Rdist, Rl/d, and Δφmax to simulate the geometric model of individual coarse particle with irregular shape. 2. Import the generated geometric model into PFC 2D to construct clump library. The clumps are created based on mid-surface method. These clumps are used to represent realistic coarse particle. 3. Put coarse particles into the model domain based on the overlap detection algorithm which is proposed above (Figure 11a). 4. Generate sample that is composed of single disk to simulate fine particles in the same model domain (Figure 11b). The size of particle and the porosity of fine particles specimen should be in a given range. 5. Remove fine particles that overlap with coarse particles (Figure 11c,d). The generated granular mixture sample is shown in Figure 11e.

Granular Mixture with Predetermined Coarse Particle Shape Features
The relationship between input parameters (Rdist, Rl/d, and Δφmax) and shape descriptors are discussed in Section 3.2, as well as the input parameters used to control particle shape features (such as the mean value of some shape descriptors). This section illustrates the capabilities of the proposed method in generating sample with predetermined shape feature. In this example, the mean and minimum elongation are equal to 0.81 and 1/3, respectively. The value of mean roundness (RImean) and mean regularity (REmean) are set to 0.7 and 3.0, respectively. Based on Equations (14)(15)(16)(17), the three representative parameters (Rdist, Rl/d, and Δφmax) are back-calculated. The values are 1/3, 0.094, and 0.735, respectively. Besides, the equivalent diameters of the coarse particles varied uniformly from 0.01 to 1.0 cm and the coarse particles content is specified as 40%. It should be noted that because the clumps are used to represent realistic coarse particles shapes. Therefore, the area of all the coarse particles is obtained by calculating the area of clumps not that of randomly generated polygons. A

Granular Mixture with Predetermined Coarse Particle Shape Features
The relationship between input parameters (R dist , R l/d , and ∆ϕ max ) and shape descriptors are discussed in Section 3.2, as well as the input parameters used to control particle shape features (such as the mean value of some shape descriptors). This section illustrates the capabilities of the proposed method in generating sample with predetermined shape feature. In this example, the mean and minimum elongation are equal to 0.81 and 1/3, respectively. The value of mean roundness (RI mean ) and mean regularity (RE mean ) are set to 0.7 and 3.0, respectively. Based on Equations (14)- (17), the three representative parameters (R dist , R l/d , and ∆ϕ max ) are back-calculated. The values are 1/3, 0.094, and 0.735, respectively. Besides, the equivalent diameters of the coarse particles varied uniformly from 0.01 to 1.0 cm and the coarse particles content is specified as 40%. It should be noted that because the clumps are used to represent realistic coarse particles shapes. Therefore, the area of all the coarse particles is obtained by calculating the area of clumps not that of randomly generated polygons. A granular mixture sample with predetermined shape feature and the derived distributions of the three shape descriptors (elongation, roundness, and regularity) are shown in Figure 12. It can be observed that the shape descriptors have a good match with the target value, which validates the feasibility of the RAB algorithm and three representative parameters in precisely controlling the generation of granular materials. granular mixture sample with predetermined shape feature and the derived distributions of the three shape descriptors (elongation, roundness, and regularity) are shown in Figure 12. It can be observed that the shape descriptors have a good match with the target value, which validates the feasibility of the RAB algorithm and three representative parameters in precisely controlling the generation of granular materials.

Granular Mixture with Different Coarse Particle Distances
The proposed overlap detection algorithm has an advantage in calculation of the minimum distance between coarse particles or the distance between coarse particle and the model boundary.
To illustrate the capability of this algorithm, four granular mixture samples with different minimum coarse particle distances are generated ( Figure 13). For these four samples, the minimum distance (dmin) between coarse particles are 0, 0.05, 0.10, and 0.20 cm, respectively. The values of Δφmax, Rdist, and Rl/d are set to π/4, 1/3, and 0.1, respectively. The value range of equivalent diameters are 0.3 to 1.5 cm. The content of the coarse particles is equal to 30%. The distributions of minimum distance between clumps are shown in Figure 14. All these results show that distance meet the requirements of the target value.

Granular Mixture with Different Coarse Particle Distances
The proposed overlap detection algorithm has an advantage in calculation of the minimum distance between coarse particles or the distance between coarse particle and the model boundary.
To illustrate the capability of this algorithm, four granular mixture samples with different minimum coarse particle distances are generated ( Figure 13). For these four samples, the minimum distance (d min ) between coarse particles are 0, 0.05, 0.10, and 0.20 cm, respectively. The values of ∆ϕ max , R dist , and R l/d are set to π/4, 1/3, and 0.1, respectively. The value range of equivalent diameters are 0.3 to 1.5 cm. The content of the coarse particles is equal to 30%. The distributions of minimum distance between clumps are shown in Figure 14. All these results show that distance meet the requirements of the target value.

Conclusion
This paper presents a systematic method to generate arbitrary two-dimensional particle for the numerical simulation of granular materials. A random angular bend (RAB) algorithm is proposed and coded in Python to simulate the geometric model of individual particle with irregular shape. Three representative parameters (Rdist, Rl/d, and Δφmax) are used to quantitatively control the shape feature of generated polygons in terms of three major aspects (form, roundness, and surface texture). Besides, a novel overlap detection algorithm is developed to address difficulties associated with spatial allocation of irregularly shaped particles.

Conclusion
This paper presents a systematic method to generate arbitrary two-dimensional particle for the numerical simulation of granular materials. A random angular bend (RAB) algorithm is proposed and coded in Python to simulate the geometric model of individual particle with irregular shape. Three representative parameters (Rdist, Rl/d, and Δφmax) are used to quantitatively control the shape feature of generated polygons in terms of three major aspects (form, roundness, and surface texture). Besides, a novel overlap detection algorithm is developed to address difficulties associated with spatial allocation of irregularly shaped particles.

Conclusions
This paper presents a systematic method to generate arbitrary two-dimensional particle for the numerical simulation of granular materials. A random angular bend (RAB) algorithm is proposed and coded in Python to simulate the geometric model of individual particle with irregular shape. Three representative parameters (R dist , R l/d , and ∆ϕ max ) are used to quantitatively control the shape feature of generated polygons in terms of three major aspects (form, roundness, and surface texture). Besides, a novel overlap detection algorithm is developed to address difficulties associated with spatial allocation of irregularly shaped particles. Finally, two examples are further employed to validate the feasibility of the proposed algorithm for the numerical modeling of realistic granular materials. The main conclusions are summarized as follows.

1.
The proposed RAB algorithm shows a good performance in the generation of randomly shaped polygons (especially for nonstar-like particle reconstruction based on image information).

2.
Three representative parameters (R dist , R l/d , and ∆ϕ max ) in RAB algorithm could quantitatively control the shape feature of generated polygons by control three shape descriptors (elongation index, roundness index, and regularity). Besides, these three parameters have definitude physical meaning, which makes the RAB algorithm easier to understand. 3.
The proposed overlap detection algorithm is able to allocate particle to the model domain easily by calculating the minimum distance between coarse particles or the distance between coarse particle and the model boundary.
It should be noted that the RAB algorithm is only suitable for the generation of two-dimensional particle. For real 3D particle modeling, this algorithm needs to be extended to a general 3D case in future works. The study provide a foundation for the construction of granular materials based on angular bend theory.

Conflicts of Interest:
The authors declare no conflicts of interest.

List of symbols
l Arc length from point P 0 up to the vertex P i ϕ Cumulative angular bends P i The ith vertex on the polygon ∆l i Edge length ∆ϕ i Angular bend at vertex P i L all The total length of the polygon.
β i The counterclockwise angle that vector → P i P i+1 makes with the positive x-axis l step Constant step length EI Elongation index S The width of the smallest rectangular box containing the particle L The length of the smallest rectangular box containing the particle RI Roundness index A The projected area of the particle P The perimeter of the particle RE Regularity index P conv The convex perimeter of the particle D OPi The distance between O and the generated point P i D min The minimum distance between the center O and the generated points D max The maximum distance between the center O and the generated points R dist The ratio between D min and D max R l/d The ratio between l steplength and D max ∆ϕ max The maximum value of angle bend EI min The minimum value of elongation EI mean The mean value of elongation RI mean The mean value of roundness RE mean The mean value of regularity D i The distance between pebble and the model boundary Dist_ij The distance between pebble peb_i in clump C1 and pebble peb_j in C2 d min The minimum distance between coarse particles