# Automated Shape Analysis and DEM Study on Graded Crushed Stone

Authors: Hao Bai; Ruidong Li; Xiangyu Hu; Fei Chen; Zhiyong Liao
Journal: Advances in Materials Science and Engineering
 (2021)

---

## Abstract

Graded crushed stone (GCS), as a cheap and essential component, is of great importance in road construction. The irregularity and variability of particle shape is known to affect the packing characteristics of GCS, such as compactness and void ratio. In this study, the realistic particle outline is first automatically extracted based on digital image processing and deep learning algorithms. Then, the elongation (EI), roundness (Rd), and roughness (Rg) of GCS are quantified by shape evaluation algorithms. Moreover, based on the establishment of the GCS shape library, the gravity deposition with various elongations is simulated using the discrete element method to study the packing features of GCS. The elongation effects on the macroscopic and microscopic quantities are explored. Finally, the shear behavior of GCS is studied. The results illustrate that elongation has a significant effect on the packing of GCS.


---

## Body

## 1. Introduction

With the development of the highway industry and the deepening of engineering practice, the asphalt pavement structure used in China has been diversified from the single semirigid asphalt base, and the graded crushed stone flexible base asphalt pavement has become one of the main pavement structure forms in China. As a mixture composed of aggregates of different sizes, the graded crushed stone (GCS) is widely applied between the asphalt surface and semirigid base and can effectively weaken the cracks in the top of the semirigid base and the stress concentration effect. Additionally, GCS can also reduce the temperature shrinkage and dry shrinkage stress caused by the change of temperature and humidity of the semirigid base [1, 2].Due to the great significance of GCS to highway construction, it is particularly essential to study its mechanical properties, deformation, and failure mechanism. However, it is difficult to analyze its mechanical properties quantitatively and qualitatively by conventional methods because of its granular structural characteristics and nonlinear mechanical properties [3]. On the one hand, the conventional laboratory test for GCS is complex and costly, which makes it difficult to reveal the macro-meso mechanical characteristics of GCS under load. On the other hand, the granular material such as GCS has great variability and the data results of different specimens in the same group are relatively low in correlation, which makes experimental data unreliable [4]. Hence, it is very necessary to use the numerical test method to study the mechanical properties of GCS from the microscopic point of view, to make up for the lack of laboratory tests and provide a powerful means for the in-depth understanding of GCS [5, 6].As a basic road material, the shape of crushed stone particles has a significant effect on its mechanical properties. Scholars have done much research on the mechanical properties of GCS influenced by GCS shape. An [7] carried out a quantitative analysis on the particle shape of GCS to explore the relationship between the real particle shape and particle accumulation porosity. Le Pen et al. [8] set out methods for evaluating form and roundness (aspects of shape) and proposed a new measure for evaluating roundness, termed ellipses for crushed stone. Zhang et al. [9] simulated the geometric anisotropy of the GCS using randomly generated models. The results showed that the changes of both friction and anisotropy had a similar trend with the particle size enlarging in the dense assembly, which maintains stability when the size is less than 3.0 mm.Numerical simulation methods are often used in exploring the properties of GCS. Jiang et al. [10] used PFC2D to establish a numerical test method for California bearing ratio (CBR) and found that CBR values increase with the friction coefficient at the contact and shear modulus of the rocks, while the influence of Poisson’s ratio on the CBR values is insignificant. Ren et al. [11] proposed a random calculation model for GCS and studied the Direct Shear Tests (DST) of GCS by PFC2D. Xu et al. [12] proposed Plastic Deformation Index (PDI) as the evaluation indicator and founded that the PDI should be less than or equal to 1900 mm number for expressway and first-grade highway. Jiang et al. [13] established a numerical model of dynamic triaxial test of GCS based on PFC2D and founded that the size of the specimen has a slight influence on the simulated value of axial strain when specimen height is greater than 40 cm and diameter is greater than 20 cm.In this study, the particle outlines are firstly extracted based on deep learning algorithms. The U-NET model is trained for contour extraction. Then, the shape evaluation indexes including elongation, roundness, and roughness are quantified. In addition, the gravity deposition with a distinctEI value is simulated by the discrete element method to explore the elongation effects on the macroscopic and microscopic quantities. Finally, the shear behavior of GCS is studied. Overall, this study lays a foundation for future research about GCS.


## 2. Materials and Equipment

In construction projects, GCS are transported by large-scale transport equipment. To simulate the real conveyance cases, a conveyor belt is used to carry GCS obtained from the laboratory of civil engineering materials in Shenzhen University. The GCS is first screened and the stone with a particle size of 4.75 mm∼25 mm is selected for preparation. Then, in order to obtain the data set with images of GCS in batches, the equipment shown in Figure1(a) is used. The specific components are illustrated as follows:Component (1): A Fujinon X-T20 camera as shown in Figure1(b) with a resolution of 24 million pixels (6000 × 4000 pixels) and two lenses: Fujinon XF60 mm F2.4 R Macro and Fujinon XF35 mm F2 R WR. At least 90 pixels are used to represent the smallest diameter of the particle.Component (2): A conveyor belt as shown in Figure1(c) used to carry GCS to simulate the real transport of GCS.Component (3): Two lamps composed of 256 LEDs as shown in Figure1(d) provide controllable lighting conditions. The light is scattered uniformly so the lamps have good color rendering. The LED lamp arm is installed on the copy stand and is perpendicular to the plane of each LED lamp. In addition, the lamp is placed 45 cm above the transparent plastic plate. Besides, the side of the lamp is parallel to the side of the base of the copy stand. The color temperature is 5600 K at the maximum brightness value.Figure 1
Materials and equipment for GCS images capture. (a) The equipment to photograph GCS. (b) Fujinon camera. (c) Conveyor belt. (d) LED light.
(a)(b)(c)(d)


## 3. Automated Extraction of GCS Contours

### 3.1. Implementation of Algorithms

For image processing, convolutional neural network (CNN) and full convolutional network (FCN) are the main frames for image segmentation. However, CNN has two obvious drawbacks: (1) slow training process because of too much redundancy and (2) asynchronism between obtaining sensing field and positioning accuracy. Mostly, CNN is used for image classification at the image level. By contrast, FCN can classify images at pixel level with good training results. In order to avoid these defects as much as possible, in this study, according to Liang et al. [14] a lightweight U-NET [15] method of a fully convolutional neural network is adopted for image segmentation.


#### 3.1.1. Lightweight U-NET and Operation

In this study, a U-NET deep neural network is introduced. Firstly, the lightweight U-NET particle image is converted into a 3D matrixW × C × h as the input data, where W and h are the width and height of the input image, respectively; C is the number of color channels of input images, where C = 3 represents RGB image and C = 1 represents gray/binary image. The gray images are input data, and binary images are output data. Then, the network generates a mask by a matrix with a size of W × h × 1 to record the pixels on the projection of GCS. For clarity, three parts of the operation of lightweight U-NET are briefly explained: convolution, Maxpool, and deconvolution.The kernel of convolution used in this study is defined as a square matrix. And the overwritten elements in the input matrix are transformed into feature mapping elements (yij) by multiplying elements:(1)yij=AF∑m=i−si+s∑n=j−sj+sXmnWmn+b,s=kn2,where kn is the size of the kernel matrix; xmn is the element at position (m, n) of the input matrix; wmn is the element of the nuclear matrix overlapping with xmn; b is the bias term of the output feature graph, and AF(·) is the activation function. During the training process, both parameters wmn and b are the targets of training. According to the recommended values [16], the default value of kn is 3 and the default activation function is ReLu (AF(x) = Max(0, x)) in this study. The principle of convolution is shown in Figure 2.Figure 2
Convolution operation [14].The Maxpool and Upconv operations use different kernels and computations to downsample and upsample the input matrix, respectively. But they are both executed by sliding Windows, just like convolution operations. The Maxpool operation can be regarded as a special convolution operation, withkn = 2 and S = 2, and the function of formula yij which outputs the maximum element covered by the kernel. Deconvolution, as opposed to convolution, has different interpretations like upsampling or transposed convolution. In this study, we used a UpSampling2D built-in Keras to implement the upsampling operation.


#### 3.1.2. Model Optimization

Considering the limited size and computing power of the model, the size of the network was reduced by half; thus, the number of channels (c) of all operations was reduced to half of the original number except for the number of channels at the end and the beginning. At the same time, in order to prevent overfitting, part of neural network units are temporarily dropped from the network according to a certain probability, which is equivalent to finding a thinner network from the original network. Dropout [17] should be performed after the fourth Maxpool operation with the parameter keep_prob_(the probability of which part of the network would not be dropped) set to 0.6. It is a parameter of the dropout method. The overall structure of U-NET is shown in Figure 3.Figure 3
Improved overall U-NET architecture.The commonly used optimization method Adam of deep learning is used to dynamically adjust the learning rate of each parameter by using the first-order moment estimation and second-order moment estimation of the gradient. For the loss function, cross-entropy was used in this study.


### 3.2. Training of the Model

#### 3.2.1. Data Set Preparation

For the collected pictures, the “labelme” image marking tool was used to manually mark the GCS. A total of 400 pictures is collected as the training set. Meanwhile, 50 pictures are collected as the validation set.


#### 3.2.2. Data Enhancement

Image enhancement is a widely used training sample expansion method. In order to improve the diversity of training samples and enhance the antinoise performance of the model, image enhancement technology was used to expand the obtained samples. A total of 4 categories and 8 image enhancement technologies were used, including (1) image rotation; (2) image flipping; (3) random brightness transformation; and (4) random pepper and salt noise, as shown in Figures4.Figure 4
Data enhancement techniques: (a) initial image; (b–d) rotated image; (e) flipped image; (f-g) brightness transformed image; (h-i) image with salt and pepper noise.
(a)(b)(c)(d)(e)(f)(g)(h)(i)


#### 3.2.3. Recognition Effect

Parameters are adjusted manually. Through multiple experiments of different parameter setting, the best parameter is chosen. By adjusting the parameters, the learning rate was set to 1 × 10−4, the scale of the U-NET network was reduced, and the effect of the data set was enhanced. The performance of U-NET model is validated by the validation data set. After each iteration of 30 times and 300 rounds of epochs training, the accuracy of the validation set reached a value of 96%, which shows a strong generalization ability of this U-NET model. The loss function reached a value of about 0.01. The curves of accuracy and loss value of the trained model are shown in Figures 5 and 6, respectively.Figure 5
Accuracy of the validation set.Figure 6
Loss value of loss function.The particle identification effect is shown in Figure7.Figure 7
(a) Grayscale. (b) Recognition effect.
(a)(b)


### 3.3. Segmentation of GCS

Generally, images obtained through the trained network still have some defects. The red marked area in Figure8 shows the adhesion between particles, which affects further morphological analysis and calculation.Figure 8
Problems existing in neural network output images.The binary image output by the neural network often contains some defects, and postprocessing technology is needed to improve the accuracy. Among all the defects, the two most obvious problems are as follows: (1) there are holes in the particles. (2) There are interconnected particles. To solve the problem of holes, a relatively mature pore filling algorithm can be used [18]. For the segmentation of connected particles, algorithms based on image morphology, such as the watershed method [19] and the corrosion-flooding method [20] are mostly adopted. Among them, the classic corrosion-flooding method used by Liu et al. [20] has a simple concept and has a good effect on the segmentation of objects with small contact sizes. The main steps are shown in Figure 9.(1)
The corrosion algorithm is used to separate the connected particles into several independent parts, and each part is labeled as the expanded seed of the flooding method.(2)
The seed is taken as the base point and expanded outward in the way of coloring. When the pixel meets the boundary of other seeds or particles, freeze the pixel and stop the expansion of this part.(3)
When the expansion of all seed pixels stops, the segmented particles are obtained.Figure 9
Main steps of the corrosion-flooding process [14].The radiusrero of the corrosion mask is a decisive parameter in the corrosion operation, and a fixed value is often used in the existing research. Liu et al. [20] suggested the value of rero as follows: (2)rero=0.3S0.


## 3.1. Implementation of Algorithms

For image processing, convolutional neural network (CNN) and full convolutional network (FCN) are the main frames for image segmentation. However, CNN has two obvious drawbacks: (1) slow training process because of too much redundancy and (2) asynchronism between obtaining sensing field and positioning accuracy. Mostly, CNN is used for image classification at the image level. By contrast, FCN can classify images at pixel level with good training results. In order to avoid these defects as much as possible, in this study, according to Liang et al. [14] a lightweight U-NET [15] method of a fully convolutional neural network is adopted for image segmentation.


### 3.1.1. Lightweight U-NET and Operation

In this study, a U-NET deep neural network is introduced. Firstly, the lightweight U-NET particle image is converted into a 3D matrixW × C × h as the input data, where W and h are the width and height of the input image, respectively; C is the number of color channels of input images, where C = 3 represents RGB image and C = 1 represents gray/binary image. The gray images are input data, and binary images are output data. Then, the network generates a mask by a matrix with a size of W × h × 1 to record the pixels on the projection of GCS. For clarity, three parts of the operation of lightweight U-NET are briefly explained: convolution, Maxpool, and deconvolution.The kernel of convolution used in this study is defined as a square matrix. And the overwritten elements in the input matrix are transformed into feature mapping elements (yij) by multiplying elements:(1)yij=AF∑m=i−si+s∑n=j−sj+sXmnWmn+b,s=kn2,where kn is the size of the kernel matrix; xmn is the element at position (m, n) of the input matrix; wmn is the element of the nuclear matrix overlapping with xmn; b is the bias term of the output feature graph, and AF(·) is the activation function. During the training process, both parameters wmn and b are the targets of training. According to the recommended values [16], the default value of kn is 3 and the default activation function is ReLu (AF(x) = Max(0, x)) in this study. The principle of convolution is shown in Figure 2.Figure 2
Convolution operation [14].The Maxpool and Upconv operations use different kernels and computations to downsample and upsample the input matrix, respectively. But they are both executed by sliding Windows, just like convolution operations. The Maxpool operation can be regarded as a special convolution operation, withkn = 2 and S = 2, and the function of formula yij which outputs the maximum element covered by the kernel. Deconvolution, as opposed to convolution, has different interpretations like upsampling or transposed convolution. In this study, we used a UpSampling2D built-in Keras to implement the upsampling operation.


### 3.1.2. Model Optimization

Considering the limited size and computing power of the model, the size of the network was reduced by half; thus, the number of channels (c) of all operations was reduced to half of the original number except for the number of channels at the end and the beginning. At the same time, in order to prevent overfitting, part of neural network units are temporarily dropped from the network according to a certain probability, which is equivalent to finding a thinner network from the original network. Dropout [17] should be performed after the fourth Maxpool operation with the parameter keep_prob_(the probability of which part of the network would not be dropped) set to 0.6. It is a parameter of the dropout method. The overall structure of U-NET is shown in Figure 3.Figure 3
Improved overall U-NET architecture.The commonly used optimization method Adam of deep learning is used to dynamically adjust the learning rate of each parameter by using the first-order moment estimation and second-order moment estimation of the gradient. For the loss function, cross-entropy was used in this study.


## 3.1.1. Lightweight U-NET and Operation

In this study, a U-NET deep neural network is introduced. Firstly, the lightweight U-NET particle image is converted into a 3D matrixW × C × h as the input data, where W and h are the width and height of the input image, respectively; C is the number of color channels of input images, where C = 3 represents RGB image and C = 1 represents gray/binary image. The gray images are input data, and binary images are output data. Then, the network generates a mask by a matrix with a size of W × h × 1 to record the pixels on the projection of GCS. For clarity, three parts of the operation of lightweight U-NET are briefly explained: convolution, Maxpool, and deconvolution.The kernel of convolution used in this study is defined as a square matrix. And the overwritten elements in the input matrix are transformed into feature mapping elements (yij) by multiplying elements:(1)yij=AF∑m=i−si+s∑n=j−sj+sXmnWmn+b,s=kn2,where kn is the size of the kernel matrix; xmn is the element at position (m, n) of the input matrix; wmn is the element of the nuclear matrix overlapping with xmn; b is the bias term of the output feature graph, and AF(·) is the activation function. During the training process, both parameters wmn and b are the targets of training. According to the recommended values [16], the default value of kn is 3 and the default activation function is ReLu (AF(x) = Max(0, x)) in this study. The principle of convolution is shown in Figure 2.Figure 2
Convolution operation [14].The Maxpool and Upconv operations use different kernels and computations to downsample and upsample the input matrix, respectively. But they are both executed by sliding Windows, just like convolution operations. The Maxpool operation can be regarded as a special convolution operation, withkn = 2 and S = 2, and the function of formula yij which outputs the maximum element covered by the kernel. Deconvolution, as opposed to convolution, has different interpretations like upsampling or transposed convolution. In this study, we used a UpSampling2D built-in Keras to implement the upsampling operation.


## 3.1.2. Model Optimization

Considering the limited size and computing power of the model, the size of the network was reduced by half; thus, the number of channels (c) of all operations was reduced to half of the original number except for the number of channels at the end and the beginning. At the same time, in order to prevent overfitting, part of neural network units are temporarily dropped from the network according to a certain probability, which is equivalent to finding a thinner network from the original network. Dropout [17] should be performed after the fourth Maxpool operation with the parameter keep_prob_(the probability of which part of the network would not be dropped) set to 0.6. It is a parameter of the dropout method. The overall structure of U-NET is shown in Figure 3.Figure 3
Improved overall U-NET architecture.The commonly used optimization method Adam of deep learning is used to dynamically adjust the learning rate of each parameter by using the first-order moment estimation and second-order moment estimation of the gradient. For the loss function, cross-entropy was used in this study.


## 3.2. Training of the Model

### 3.2.1. Data Set Preparation

For the collected pictures, the “labelme” image marking tool was used to manually mark the GCS. A total of 400 pictures is collected as the training set. Meanwhile, 50 pictures are collected as the validation set.


### 3.2.2. Data Enhancement

Image enhancement is a widely used training sample expansion method. In order to improve the diversity of training samples and enhance the antinoise performance of the model, image enhancement technology was used to expand the obtained samples. A total of 4 categories and 8 image enhancement technologies were used, including (1) image rotation; (2) image flipping; (3) random brightness transformation; and (4) random pepper and salt noise, as shown in Figures4.Figure 4
Data enhancement techniques: (a) initial image; (b–d) rotated image; (e) flipped image; (f-g) brightness transformed image; (h-i) image with salt and pepper noise.
(a)(b)(c)(d)(e)(f)(g)(h)(i)


### 3.2.3. Recognition Effect

Parameters are adjusted manually. Through multiple experiments of different parameter setting, the best parameter is chosen. By adjusting the parameters, the learning rate was set to 1 × 10−4, the scale of the U-NET network was reduced, and the effect of the data set was enhanced. The performance of U-NET model is validated by the validation data set. After each iteration of 30 times and 300 rounds of epochs training, the accuracy of the validation set reached a value of 96%, which shows a strong generalization ability of this U-NET model. The loss function reached a value of about 0.01. The curves of accuracy and loss value of the trained model are shown in Figures 5 and 6, respectively.Figure 5
Accuracy of the validation set.Figure 6
Loss value of loss function.The particle identification effect is shown in Figure7.Figure 7
(a) Grayscale. (b) Recognition effect.
(a)(b)


## 3.2.1. Data Set Preparation

For the collected pictures, the “labelme” image marking tool was used to manually mark the GCS. A total of 400 pictures is collected as the training set. Meanwhile, 50 pictures are collected as the validation set.


## 3.2.2. Data Enhancement

Image enhancement is a widely used training sample expansion method. In order to improve the diversity of training samples and enhance the antinoise performance of the model, image enhancement technology was used to expand the obtained samples. A total of 4 categories and 8 image enhancement technologies were used, including (1) image rotation; (2) image flipping; (3) random brightness transformation; and (4) random pepper and salt noise, as shown in Figures4.Figure 4
Data enhancement techniques: (a) initial image; (b–d) rotated image; (e) flipped image; (f-g) brightness transformed image; (h-i) image with salt and pepper noise.
(a)(b)(c)(d)(e)(f)(g)(h)(i)


## 3.2.3. Recognition Effect

Parameters are adjusted manually. Through multiple experiments of different parameter setting, the best parameter is chosen. By adjusting the parameters, the learning rate was set to 1 × 10−4, the scale of the U-NET network was reduced, and the effect of the data set was enhanced. The performance of U-NET model is validated by the validation data set. After each iteration of 30 times and 300 rounds of epochs training, the accuracy of the validation set reached a value of 96%, which shows a strong generalization ability of this U-NET model. The loss function reached a value of about 0.01. The curves of accuracy and loss value of the trained model are shown in Figures 5 and 6, respectively.Figure 5
Accuracy of the validation set.Figure 6
Loss value of loss function.The particle identification effect is shown in Figure7.Figure 7
(a) Grayscale. (b) Recognition effect.
(a)(b)


## 3.3. Segmentation of GCS

Generally, images obtained through the trained network still have some defects. The red marked area in Figure8 shows the adhesion between particles, which affects further morphological analysis and calculation.Figure 8
Problems existing in neural network output images.The binary image output by the neural network often contains some defects, and postprocessing technology is needed to improve the accuracy. Among all the defects, the two most obvious problems are as follows: (1) there are holes in the particles. (2) There are interconnected particles. To solve the problem of holes, a relatively mature pore filling algorithm can be used [18]. For the segmentation of connected particles, algorithms based on image morphology, such as the watershed method [19] and the corrosion-flooding method [20] are mostly adopted. Among them, the classic corrosion-flooding method used by Liu et al. [20] has a simple concept and has a good effect on the segmentation of objects with small contact sizes. The main steps are shown in Figure 9.(1)
The corrosion algorithm is used to separate the connected particles into several independent parts, and each part is labeled as the expanded seed of the flooding method.(2)
The seed is taken as the base point and expanded outward in the way of coloring. When the pixel meets the boundary of other seeds or particles, freeze the pixel and stop the expansion of this part.(3)
When the expansion of all seed pixels stops, the segmented particles are obtained.Figure 9
Main steps of the corrosion-flooding process [14].The radiusrero of the corrosion mask is a decisive parameter in the corrosion operation, and a fixed value is often used in the existing research. Liu et al. [20] suggested the value of rero as follows: (2)rero=0.3S0.


## 4. Quantitatively Analysis of GCS Shapes

### 4.1. Methods for Quantifying the Particle Shape Indexes

#### 4.1.1. Quantification of Elongation

To reflect the particle contour, elongation (EI) is often used:(3)EI=DminorDmajor,where Dminor is the width along the minor principal axis and Dmajor is the length along the major principal axis. Apparently, a smaller EI value represents the greater slenderness of the particle. Through a second-order tensor matrix based on the normal vector of particle contour, the two principal axes can be determined by (4)Ωij=1Lp∑k=1nlkTikTjk,C=Ω11Ω12Ω21Ω22=cosθ−sinθsinθcosθλa00λbcosθsinθ−sinθcosθ,where Tik is the component of the unit normal vector Tk and the kth micro arc in i direction and Tjk is the counterpart in j direction; Lk is the length of the kth micro arc; LP is the circumference of the particle contour, which is the sum of lengths of n micro arcs; and λa and λb are eigenvalues of matrix C, λa≥λb. The value of λb is the minimum percentage of the arc length in particle contour circumference along the major principal axis direction of Ωij. The value of λa is the maximum percentage along the minor principal axis direction of Ωij. After having two principle axis directions, Dminor and Dmajor can be calculated in a rectangular boundary containing the particle contour (Figure 10).Figure 10
Rectangular box surrounding the particle contour.


#### 4.1.2. Quantification of Roundness

Roundness (Rd) reflects the sharpness and smoothness of the particle corners. Rd can be calculated by the ratio between the mean curvature of particle corner and the radius of the maximum inscribed circle of particle contour, as shown below:(5)Rd=1N∑i=1NriRinsc,where N is the number of corners; ri is the curvature radius of the ith corner, and Rinsc is the radius of the maximum inscribed circle of the particle contour.Wadell [21] presented that a point belongs to the corner region when the curvature radius of the point is less than Rinsc. Thus, the curvature radius of all contour points should be calculated to classify the corner region. For Fourier-presented contour curve, the curvature radius of the contour point can be calculated as (6)rθ=rθ2+drθ/dθ23/2rθ2+2drθ/dθ2−rθd2rθ/dθ2,where r(θ) is the curvature radius.The contour pointP(θ, r(θ)) can be seen as a corner point when r(θ) < Rinsc. Thus, the corner region is the contour of all corner points (Figure 11(a)). The ODEC algorithm [22] can be used to calculate the radius of the inscribed circle of corner points, as shown in Figure 11(b). To be specific, the inner normal vector Ni⟶ of the ith corner point is first calculated. Then, the inscribed circle along the direction of Ni⟶ is determined, by starting from the ith corner point with a radius of ∆r. After that, whether the inscribed circle is tangent to any other contour point is determined. If not, the radius ∆r is increased and the last two steps are repeated; otherwise, terminate the process and take ∆r as the radius of the inscribed circle of the ith corner point.Figure 11
(a) Corner region of the particle contour, (b) determining the inscribed circle of a corner point, and (c) optimal inscribed circles of all corners.
(a)(b)(c)The optimal inscribed circle of the corner can be determined by having the inscribed circles of all corner points. To value the fit performance of the inscribed circle of a corner point relative to other corner points,efit is defined:(7)efit=∑i=1ncpxi−xc2+yi−yc2−rcncp,where ncp is the number of corner points; (xi, yi) are the coordinates of the ith corner point; (xc, yc) are the center coordinates of the inscribed circle; and rc is the radius of inscribed circle radius. Particularly, the minimum value of efit represents the optimal inscribed circle of the corner. Through the optimal inscribed circles of all corners shown in Figure 11(c), the curvature radius of each corner and Rd can be calculated.


#### 4.1.3. Quantification of Roughness

Roughness is a shape index that evaluates the concavity of the particle contour and it can show the shape difference between the real contour and the smoothed contour of particles. Fourier transform can be used to calculate the smoothed contour of a particle. The larger Fourier series (N) corresponds to a higher similarity between the real contour and the smoothed contour [23]. In this study, the Fourier series N is set to 16 to calculate the smoothed contour. The deviation between the smoothed contour and the real contour is shown in Figure 12.Figure 12
Smoothed contour and real contour of a particle.Next, the deviation regions and the geometric types are extracted and determined. For nonintersecting parts between smoothed contour and real contour, the deviation region can be calculated as a quadrilateral area. For intersecting parts, the deviation region is simplified to two triangles. At theith microunit, the deviation distance (△di) is the ratio of the deviation area to the corresponding length. Based on the average deviation distance of all microunits, Rg can be expressed as (8)△di=Aili,Rg=△d¯r¯=πAp⋅1Lp∑i=1Nd△di×li,where △d¯ is the average deviation distance, Nd is the number of all microunits, Ap is the area of the particle contour, Ai is the deviation area, li is the corresponding arc length, and Lp is the perimeter of the particle contour.


### 4.2. Statistics of Shape Indexes

Through the quantification algorithms introduced above, the elongation (EI), roundness (Ed), and roughness (Rg) of 1200 GCS were calculated. The fitting results are shown in Table 1.Table 1
Fitting function and related parameters.
Shape indexFitting functionMean valueGoodness of fit R2RangeEINormal distribution0.680.98From 0.25 to 1RdNormal distribution0.330.98From 0.2 to 0.6RgNormal distribution0.00390.97From 0.002 to 0.006Normal distribution tests were implemented on the probability density of the three shape indexes. The results showed that all three indexes satisfied normal distribution appropriately. Within the results, the width of the distribution of elongation is the maximum and the width of the distribution of roundness is the minimum, with ranges from 0.4 to 1 and from 0.2 to 0.5, respectively. Furthermore, the mean values of the three shape indexesEI, Rd, and Rg are 0.68, 0.33, and 0.0039, respectively.


### 4.3. Establish Shape Library

To facilitate retrieval and sample preparation, the shape indicators of the GCS are stored together with the coordinates, with the storage format [number, outline, slender length, elongation, roundness, roughness].The steps for sample preparation from the particle library are as follows:(1)
The required shape parameter range is input(2)
According to the range of these values, the two-dimensional contour that meets the requirements can be searched automatically(3)
The visualized two-dimensional contour is selected and is output to the “dxf” fileAfter establishing the particle library, it will be very convenient to select the two-dimensional contours of GCS that meet specific requirements, so as to facilitate the discrete element modeling of GCS. Moreover, it lays the foundation for further study on the relationship between the geometric information and mechanical properties of GCS.


## 4.1. Methods for Quantifying the Particle Shape Indexes

### 4.1.1. Quantification of Elongation

To reflect the particle contour, elongation (EI) is often used:(3)EI=DminorDmajor,where Dminor is the width along the minor principal axis and Dmajor is the length along the major principal axis. Apparently, a smaller EI value represents the greater slenderness of the particle. Through a second-order tensor matrix based on the normal vector of particle contour, the two principal axes can be determined by (4)Ωij=1Lp∑k=1nlkTikTjk,C=Ω11Ω12Ω21Ω22=cosθ−sinθsinθcosθλa00λbcosθsinθ−sinθcosθ,where Tik is the component of the unit normal vector Tk and the kth micro arc in i direction and Tjk is the counterpart in j direction; Lk is the length of the kth micro arc; LP is the circumference of the particle contour, which is the sum of lengths of n micro arcs; and λa and λb are eigenvalues of matrix C, λa≥λb. The value of λb is the minimum percentage of the arc length in particle contour circumference along the major principal axis direction of Ωij. The value of λa is the maximum percentage along the minor principal axis direction of Ωij. After having two principle axis directions, Dminor and Dmajor can be calculated in a rectangular boundary containing the particle contour (Figure 10).Figure 10
Rectangular box surrounding the particle contour.


### 4.1.2. Quantification of Roundness

Roundness (Rd) reflects the sharpness and smoothness of the particle corners. Rd can be calculated by the ratio between the mean curvature of particle corner and the radius of the maximum inscribed circle of particle contour, as shown below:(5)Rd=1N∑i=1NriRinsc,where N is the number of corners; ri is the curvature radius of the ith corner, and Rinsc is the radius of the maximum inscribed circle of the particle contour.Wadell [21] presented that a point belongs to the corner region when the curvature radius of the point is less than Rinsc. Thus, the curvature radius of all contour points should be calculated to classify the corner region. For Fourier-presented contour curve, the curvature radius of the contour point can be calculated as (6)rθ=rθ2+drθ/dθ23/2rθ2+2drθ/dθ2−rθd2rθ/dθ2,where r(θ) is the curvature radius.The contour pointP(θ, r(θ)) can be seen as a corner point when r(θ) < Rinsc. Thus, the corner region is the contour of all corner points (Figure 11(a)). The ODEC algorithm [22] can be used to calculate the radius of the inscribed circle of corner points, as shown in Figure 11(b). To be specific, the inner normal vector Ni⟶ of the ith corner point is first calculated. Then, the inscribed circle along the direction of Ni⟶ is determined, by starting from the ith corner point with a radius of ∆r. After that, whether the inscribed circle is tangent to any other contour point is determined. If not, the radius ∆r is increased and the last two steps are repeated; otherwise, terminate the process and take ∆r as the radius of the inscribed circle of the ith corner point.Figure 11
(a) Corner region of the particle contour, (b) determining the inscribed circle of a corner point, and (c) optimal inscribed circles of all corners.
(a)(b)(c)The optimal inscribed circle of the corner can be determined by having the inscribed circles of all corner points. To value the fit performance of the inscribed circle of a corner point relative to other corner points,efit is defined:(7)efit=∑i=1ncpxi−xc2+yi−yc2−rcncp,where ncp is the number of corner points; (xi, yi) are the coordinates of the ith corner point; (xc, yc) are the center coordinates of the inscribed circle; and rc is the radius of inscribed circle radius. Particularly, the minimum value of efit represents the optimal inscribed circle of the corner. Through the optimal inscribed circles of all corners shown in Figure 11(c), the curvature radius of each corner and Rd can be calculated.


### 4.1.3. Quantification of Roughness

Roughness is a shape index that evaluates the concavity of the particle contour and it can show the shape difference between the real contour and the smoothed contour of particles. Fourier transform can be used to calculate the smoothed contour of a particle. The larger Fourier series (N) corresponds to a higher similarity between the real contour and the smoothed contour [23]. In this study, the Fourier series N is set to 16 to calculate the smoothed contour. The deviation between the smoothed contour and the real contour is shown in Figure 12.Figure 12
Smoothed contour and real contour of a particle.Next, the deviation regions and the geometric types are extracted and determined. For nonintersecting parts between smoothed contour and real contour, the deviation region can be calculated as a quadrilateral area. For intersecting parts, the deviation region is simplified to two triangles. At theith microunit, the deviation distance (△di) is the ratio of the deviation area to the corresponding length. Based on the average deviation distance of all microunits, Rg can be expressed as (8)△di=Aili,Rg=△d¯r¯=πAp⋅1Lp∑i=1Nd△di×li,where △d¯ is the average deviation distance, Nd is the number of all microunits, Ap is the area of the particle contour, Ai is the deviation area, li is the corresponding arc length, and Lp is the perimeter of the particle contour.


## 4.1.1. Quantification of Elongation

To reflect the particle contour, elongation (EI) is often used:(3)EI=DminorDmajor,where Dminor is the width along the minor principal axis and Dmajor is the length along the major principal axis. Apparently, a smaller EI value represents the greater slenderness of the particle. Through a second-order tensor matrix based on the normal vector of particle contour, the two principal axes can be determined by (4)Ωij=1Lp∑k=1nlkTikTjk,C=Ω11Ω12Ω21Ω22=cosθ−sinθsinθcosθλa00λbcosθsinθ−sinθcosθ,where Tik is the component of the unit normal vector Tk and the kth micro arc in i direction and Tjk is the counterpart in j direction; Lk is the length of the kth micro arc; LP is the circumference of the particle contour, which is the sum of lengths of n micro arcs; and λa and λb are eigenvalues of matrix C, λa≥λb. The value of λb is the minimum percentage of the arc length in particle contour circumference along the major principal axis direction of Ωij. The value of λa is the maximum percentage along the minor principal axis direction of Ωij. After having two principle axis directions, Dminor and Dmajor can be calculated in a rectangular boundary containing the particle contour (Figure 10).Figure 10
Rectangular box surrounding the particle contour.


## 4.1.2. Quantification of Roundness

Roundness (Rd) reflects the sharpness and smoothness of the particle corners. Rd can be calculated by the ratio between the mean curvature of particle corner and the radius of the maximum inscribed circle of particle contour, as shown below:(5)Rd=1N∑i=1NriRinsc,where N is the number of corners; ri is the curvature radius of the ith corner, and Rinsc is the radius of the maximum inscribed circle of the particle contour.Wadell [21] presented that a point belongs to the corner region when the curvature radius of the point is less than Rinsc. Thus, the curvature radius of all contour points should be calculated to classify the corner region. For Fourier-presented contour curve, the curvature radius of the contour point can be calculated as (6)rθ=rθ2+drθ/dθ23/2rθ2+2drθ/dθ2−rθd2rθ/dθ2,where r(θ) is the curvature radius.The contour pointP(θ, r(θ)) can be seen as a corner point when r(θ) < Rinsc. Thus, the corner region is the contour of all corner points (Figure 11(a)). The ODEC algorithm [22] can be used to calculate the radius of the inscribed circle of corner points, as shown in Figure 11(b). To be specific, the inner normal vector Ni⟶ of the ith corner point is first calculated. Then, the inscribed circle along the direction of Ni⟶ is determined, by starting from the ith corner point with a radius of ∆r. After that, whether the inscribed circle is tangent to any other contour point is determined. If not, the radius ∆r is increased and the last two steps are repeated; otherwise, terminate the process and take ∆r as the radius of the inscribed circle of the ith corner point.Figure 11
(a) Corner region of the particle contour, (b) determining the inscribed circle of a corner point, and (c) optimal inscribed circles of all corners.
(a)(b)(c)The optimal inscribed circle of the corner can be determined by having the inscribed circles of all corner points. To value the fit performance of the inscribed circle of a corner point relative to other corner points,efit is defined:(7)efit=∑i=1ncpxi−xc2+yi−yc2−rcncp,where ncp is the number of corner points; (xi, yi) are the coordinates of the ith corner point; (xc, yc) are the center coordinates of the inscribed circle; and rc is the radius of inscribed circle radius. Particularly, the minimum value of efit represents the optimal inscribed circle of the corner. Through the optimal inscribed circles of all corners shown in Figure 11(c), the curvature radius of each corner and Rd can be calculated.


## 4.1.3. Quantification of Roughness

Roughness is a shape index that evaluates the concavity of the particle contour and it can show the shape difference between the real contour and the smoothed contour of particles. Fourier transform can be used to calculate the smoothed contour of a particle. The larger Fourier series (N) corresponds to a higher similarity between the real contour and the smoothed contour [23]. In this study, the Fourier series N is set to 16 to calculate the smoothed contour. The deviation between the smoothed contour and the real contour is shown in Figure 12.Figure 12
Smoothed contour and real contour of a particle.Next, the deviation regions and the geometric types are extracted and determined. For nonintersecting parts between smoothed contour and real contour, the deviation region can be calculated as a quadrilateral area. For intersecting parts, the deviation region is simplified to two triangles. At theith microunit, the deviation distance (△di) is the ratio of the deviation area to the corresponding length. Based on the average deviation distance of all microunits, Rg can be expressed as (8)△di=Aili,Rg=△d¯r¯=πAp⋅1Lp∑i=1Nd△di×li,where △d¯ is the average deviation distance, Nd is the number of all microunits, Ap is the area of the particle contour, Ai is the deviation area, li is the corresponding arc length, and Lp is the perimeter of the particle contour.


## 4.2. Statistics of Shape Indexes

Through the quantification algorithms introduced above, the elongation (EI), roundness (Ed), and roughness (Rg) of 1200 GCS were calculated. The fitting results are shown in Table 1.Table 1
Fitting function and related parameters.
Shape indexFitting functionMean valueGoodness of fit R2RangeEINormal distribution0.680.98From 0.25 to 1RdNormal distribution0.330.98From 0.2 to 0.6RgNormal distribution0.00390.97From 0.002 to 0.006Normal distribution tests were implemented on the probability density of the three shape indexes. The results showed that all three indexes satisfied normal distribution appropriately. Within the results, the width of the distribution of elongation is the maximum and the width of the distribution of roundness is the minimum, with ranges from 0.4 to 1 and from 0.2 to 0.5, respectively. Furthermore, the mean values of the three shape indexesEI, Rd, and Rg are 0.68, 0.33, and 0.0039, respectively.


## 4.3. Establish Shape Library

To facilitate retrieval and sample preparation, the shape indicators of the GCS are stored together with the coordinates, with the storage format [number, outline, slender length, elongation, roundness, roughness].The steps for sample preparation from the particle library are as follows:(1)
The required shape parameter range is input(2)
According to the range of these values, the two-dimensional contour that meets the requirements can be searched automatically(3)
The visualized two-dimensional contour is selected and is output to the “dxf” fileAfter establishing the particle library, it will be very convenient to select the two-dimensional contours of GCS that meet specific requirements, so as to facilitate the discrete element modeling of GCS. Moreover, it lays the foundation for further study on the relationship between the geometric information and mechanical properties of GCS.


## 5. Application to DEM Simulation of GCS Packing

### 5.1. Simulation of Gravity Deposition

The mechanical parameters of GCS like shear strength and deformation modulus can be affected by the packing characteristics. The shape indexEI is set as the DEM [24] simulation variable to explore the influence of particle shape on packing properties, as shown in Figure 13(a). Considering the elongation variable, the GCS were selected and divided into 7 groups where EI is set from 0.4 to 1.0. In this study, the particles in DEM are not generated randomly; on the contrary, real GCS were used as samples to enhance the reliability of the simulations.Figure 13
(a) GCS particles with different elongation; (b) stratified gravity deposition model.
(a)(b)Stratified gravity deposition simulation of the natural packing process of particles is the common way to generate particle models. In the model, DEM particles fall naturally with gravity. After completing the previous layer deposition, a new layer of particles is generated and the deposition process is repeated until the whole simulation is completed, as shown in Figure13(b) [25]. The gravity deposition process is simulated by PFC2D, and the interparticle contact and the wall-particle contact are both linear contacts. The parameters of the model are shown in Table 2.Table 2
Particle number and model parameters.
ParameterValueParticle number1000 × 3Particle density (kg/m3)2650Interparticle frictional coefficient0.5Wall-particle frictional coefficient0Effective modulus (Pa)1 × 108Stiffness ratio4/3Damping factor0.7


### 5.2. Analysis of Macro-Meso Parameter

Macro-meso parameters reflect the particle packing properties, such as void ratio, compactness, and coordination number. In this section, the effect onEI is considered.


#### 5.2.1. Macro Parameter

The degree of density of particles can be described by a macro parameter, void ratio (e), which is affected by particle shape.The change curves of the void ratio with respect toEI are shown in Figure 14. Void ratio e decreases with the growth of EI, indicating that the particle packing gets closer with EI increasing. Then, the void ratio reaches a bottom value of 0.187 and starts to soar at the point where EI is 0.9.Figure 14
Change curves of the void ratio with respect to elongation.


#### 5.2.2. Micro Parameter

The degree of particle contact can be measured by a micro parameter, the mean coordination number (Cmean), which can reflect the internal structure of particle packing [26]:(9)Cmean=1Np∑i=1NpCi,where Np is the number of particles and Ci is the contact particle number of the ith particle.In the change curve of the mean coordination number with respect toEI shown in Figure 15, a reflection point also appears at the point where EI is 0.9. Similar to the curve of void ratio, Cmean exhibits a decreasing trend with EI increasing. After reaching the bottom of the curve, Cmean stays stable and has a slight increase.Figure 15
Change curve of the mean coordination number with respect to elongation.


#### 5.2.3. Fabric Anisotropy

The fabric anisotropy is the anisotropy between the spatial distribution of particle contacts and the length of branch vectors. The branch vector is the vector connecting the mass centers of two particles. Fabric anisotropy coefficients can be used to quantify the fabric anisotropy properties, including coefficient of contact (ac), coefficient of normal contact forces (an), coefficient of tangential contact forces (at), and coefficient of branch vectors (ad). The anisotropy coefficients can be calculated as (10)a∗=32aij∗aij∗,where a∗ is the anisotropy coefficients (ac, an, at, and ad) and aij∗ is the tensors (aijc, aijn, aijt, and aijd) indicating the directional distribution of ac, an, at, and ad. The tensors can be calculated as(11)aijc=152Φijc′,aijn=15Fijn′2fn¯,aijt=15Fijt′3fn¯,aijd=15Dijn′3d0¯,where fn¯ is the average normal contact force, d0¯ is the average length of the branch vector, Φijc′ is the partial tensor of Φijc, Fijn′ is the partial tensor of Fijn, Fijt′ is the partial tensor of Fijt, and Dijn′ is the partial tensor of Dijn. Φijc, Fijn, Fijt, and Dijn can be calculated as(12)Φijc=1Nc∑c∈Ncnicnjc,Fijn=1Nc∑c∈Ncfcnnicnjc1+aklcnkcnlc,Fijt=1Nc∑c∈Ncfctnicnjc1+aklcnkcnlc,Dijn=1Nc∑c∈Ncdcnnicnjc1+aklcnkcnlc,where nic and njc are the ith and jth components of the unit normal vector of contact c, respectively, and fcn, fct, and dcn are the normal contact force, tangential contact force and normal branch vector of contact c, respectively.Figure16 shows the change curves of the fabric anisotropy coefficients with respect to EI. Among all coefficients, ac, an, and at slightly decrease with the increase of EI, that is, the anisotropy degree of the normal contact force and tangential contact force decreases with EI collapsing. On the contrary, ad shows a positive correlation with EI, indicating that the anisotropy degree of branch vector declines with EI increasing. Furthermore, ad slightly stabilizes to 0 with EI increasing from 0.4 to 1.0, showing that the anisotropy of branch vector almost disappears. Thus, the degree of fabric anisotropy decreases with EI decreasing. Note that when EI ≥ 0.8, the fabric anisotropy coefficients stay relatively stable; thus, the degree of fabric anisotropy is not affected by EI.Figure 16
Change curves of the fabric anisotropy coefficients with respect toEI.


## 5.1. Simulation of Gravity Deposition

The mechanical parameters of GCS like shear strength and deformation modulus can be affected by the packing characteristics. The shape indexEI is set as the DEM [24] simulation variable to explore the influence of particle shape on packing properties, as shown in Figure 13(a). Considering the elongation variable, the GCS were selected and divided into 7 groups where EI is set from 0.4 to 1.0. In this study, the particles in DEM are not generated randomly; on the contrary, real GCS were used as samples to enhance the reliability of the simulations.Figure 13
(a) GCS particles with different elongation; (b) stratified gravity deposition model.
(a)(b)Stratified gravity deposition simulation of the natural packing process of particles is the common way to generate particle models. In the model, DEM particles fall naturally with gravity. After completing the previous layer deposition, a new layer of particles is generated and the deposition process is repeated until the whole simulation is completed, as shown in Figure13(b) [25]. The gravity deposition process is simulated by PFC2D, and the interparticle contact and the wall-particle contact are both linear contacts. The parameters of the model are shown in Table 2.Table 2
Particle number and model parameters.
ParameterValueParticle number1000 × 3Particle density (kg/m3)2650Interparticle frictional coefficient0.5Wall-particle frictional coefficient0Effective modulus (Pa)1 × 108Stiffness ratio4/3Damping factor0.7


## 5.2. Analysis of Macro-Meso Parameter

Macro-meso parameters reflect the particle packing properties, such as void ratio, compactness, and coordination number. In this section, the effect onEI is considered.


### 5.2.1. Macro Parameter

The degree of density of particles can be described by a macro parameter, void ratio (e), which is affected by particle shape.The change curves of the void ratio with respect toEI are shown in Figure 14. Void ratio e decreases with the growth of EI, indicating that the particle packing gets closer with EI increasing. Then, the void ratio reaches a bottom value of 0.187 and starts to soar at the point where EI is 0.9.Figure 14
Change curves of the void ratio with respect to elongation.


### 5.2.2. Micro Parameter

The degree of particle contact can be measured by a micro parameter, the mean coordination number (Cmean), which can reflect the internal structure of particle packing [26]:(9)Cmean=1Np∑i=1NpCi,where Np is the number of particles and Ci is the contact particle number of the ith particle.In the change curve of the mean coordination number with respect toEI shown in Figure 15, a reflection point also appears at the point where EI is 0.9. Similar to the curve of void ratio, Cmean exhibits a decreasing trend with EI increasing. After reaching the bottom of the curve, Cmean stays stable and has a slight increase.Figure 15
Change curve of the mean coordination number with respect to elongation.


### 5.2.3. Fabric Anisotropy

The fabric anisotropy is the anisotropy between the spatial distribution of particle contacts and the length of branch vectors. The branch vector is the vector connecting the mass centers of two particles. Fabric anisotropy coefficients can be used to quantify the fabric anisotropy properties, including coefficient of contact (ac), coefficient of normal contact forces (an), coefficient of tangential contact forces (at), and coefficient of branch vectors (ad). The anisotropy coefficients can be calculated as (10)a∗=32aij∗aij∗,where a∗ is the anisotropy coefficients (ac, an, at, and ad) and aij∗ is the tensors (aijc, aijn, aijt, and aijd) indicating the directional distribution of ac, an, at, and ad. The tensors can be calculated as(11)aijc=152Φijc′,aijn=15Fijn′2fn¯,aijt=15Fijt′3fn¯,aijd=15Dijn′3d0¯,where fn¯ is the average normal contact force, d0¯ is the average length of the branch vector, Φijc′ is the partial tensor of Φijc, Fijn′ is the partial tensor of Fijn, Fijt′ is the partial tensor of Fijt, and Dijn′ is the partial tensor of Dijn. Φijc, Fijn, Fijt, and Dijn can be calculated as(12)Φijc=1Nc∑c∈Ncnicnjc,Fijn=1Nc∑c∈Ncfcnnicnjc1+aklcnkcnlc,Fijt=1Nc∑c∈Ncfctnicnjc1+aklcnkcnlc,Dijn=1Nc∑c∈Ncdcnnicnjc1+aklcnkcnlc,where nic and njc are the ith and jth components of the unit normal vector of contact c, respectively, and fcn, fct, and dcn are the normal contact force, tangential contact force and normal branch vector of contact c, respectively.Figure16 shows the change curves of the fabric anisotropy coefficients with respect to EI. Among all coefficients, ac, an, and at slightly decrease with the increase of EI, that is, the anisotropy degree of the normal contact force and tangential contact force decreases with EI collapsing. On the contrary, ad shows a positive correlation with EI, indicating that the anisotropy degree of branch vector declines with EI increasing. Furthermore, ad slightly stabilizes to 0 with EI increasing from 0.4 to 1.0, showing that the anisotropy of branch vector almost disappears. Thus, the degree of fabric anisotropy decreases with EI decreasing. Note that when EI ≥ 0.8, the fabric anisotropy coefficients stay relatively stable; thus, the degree of fabric anisotropy is not affected by EI.Figure 16
Change curves of the fabric anisotropy coefficients with respect toEI.


## 5.2.1. Macro Parameter

The degree of density of particles can be described by a macro parameter, void ratio (e), which is affected by particle shape.The change curves of the void ratio with respect toEI are shown in Figure 14. Void ratio e decreases with the growth of EI, indicating that the particle packing gets closer with EI increasing. Then, the void ratio reaches a bottom value of 0.187 and starts to soar at the point where EI is 0.9.Figure 14
Change curves of the void ratio with respect to elongation.


## 5.2.2. Micro Parameter

The degree of particle contact can be measured by a micro parameter, the mean coordination number (Cmean), which can reflect the internal structure of particle packing [26]:(9)Cmean=1Np∑i=1NpCi,where Np is the number of particles and Ci is the contact particle number of the ith particle.In the change curve of the mean coordination number with respect toEI shown in Figure 15, a reflection point also appears at the point where EI is 0.9. Similar to the curve of void ratio, Cmean exhibits a decreasing trend with EI increasing. After reaching the bottom of the curve, Cmean stays stable and has a slight increase.Figure 15
Change curve of the mean coordination number with respect to elongation.


## 5.2.3. Fabric Anisotropy

The fabric anisotropy is the anisotropy between the spatial distribution of particle contacts and the length of branch vectors. The branch vector is the vector connecting the mass centers of two particles. Fabric anisotropy coefficients can be used to quantify the fabric anisotropy properties, including coefficient of contact (ac), coefficient of normal contact forces (an), coefficient of tangential contact forces (at), and coefficient of branch vectors (ad). The anisotropy coefficients can be calculated as (10)a∗=32aij∗aij∗,where a∗ is the anisotropy coefficients (ac, an, at, and ad) and aij∗ is the tensors (aijc, aijn, aijt, and aijd) indicating the directional distribution of ac, an, at, and ad. The tensors can be calculated as(11)aijc=152Φijc′,aijn=15Fijn′2fn¯,aijt=15Fijt′3fn¯,aijd=15Dijn′3d0¯,where fn¯ is the average normal contact force, d0¯ is the average length of the branch vector, Φijc′ is the partial tensor of Φijc, Fijn′ is the partial tensor of Fijn, Fijt′ is the partial tensor of Fijt, and Dijn′ is the partial tensor of Dijn. Φijc, Fijn, Fijt, and Dijn can be calculated as(12)Φijc=1Nc∑c∈Ncnicnjc,Fijn=1Nc∑c∈Ncfcnnicnjc1+aklcnkcnlc,Fijt=1Nc∑c∈Ncfctnicnjc1+aklcnkcnlc,Dijn=1Nc∑c∈Ncdcnnicnjc1+aklcnkcnlc,where nic and njc are the ith and jth components of the unit normal vector of contact c, respectively, and fcn, fct, and dcn are the normal contact force, tangential contact force and normal branch vector of contact c, respectively.Figure16 shows the change curves of the fabric anisotropy coefficients with respect to EI. Among all coefficients, ac, an, and at slightly decrease with the increase of EI, that is, the anisotropy degree of the normal contact force and tangential contact force decreases with EI collapsing. On the contrary, ad shows a positive correlation with EI, indicating that the anisotropy degree of branch vector declines with EI increasing. Furthermore, ad slightly stabilizes to 0 with EI increasing from 0.4 to 1.0, showing that the anisotropy of branch vector almost disappears. Thus, the degree of fabric anisotropy decreases with EI decreasing. Note that when EI ≥ 0.8, the fabric anisotropy coefficients stay relatively stable; thus, the degree of fabric anisotropy is not affected by EI.Figure 16
Change curves of the fabric anisotropy coefficients with respect toEI.


## 6. Application to DEM Simulation of GCS Shear Behavior

### 6.1. Simulation of Biaxial Compression Test

First, five samples with different elongation (EI = 0.4, 0.5, 0.6, 0.7, 0.9) were constructed. Each sample has an initial size of 9 m ∗ 18 m and contains approximately 5000 GCS particles. The linear elastic model is used as the contact law between particles. The samples are compressed by shifting the four rigid boundaries, and the boundary has a force at 100 kPa (stress boundary condition controlled by numerical servo). During the shearing process, in order to make the state of samples quasistatic (to make the loading process of the sample quasistatic loading), the rate of shear strain is supposed to be relatively small. Thus, the inertial parameters are introduced as follows:(13)Iintertia=ε˙1dσ0/ρ<10−3,where ε˙1 is the loading strain rate, d is mean particle size, ρ is material density, σ0 is the consolidation stress. Previous studies have shown that the value of the Inertia parameter should be less than 2.5 × 10−3. Thus, the shifting speed of the rigid wall is fixed to 0.05 m/s. The inertia parameter is less than 10−4 throughout the tests’ negative correlation.


### 6.2. Analysis of Macro-Meso Parameter

#### 6.2.1. Analysis of Stress Ratio

The mean stressp and deviator stress q are calculated to explore the sample shear strength. In the two-dimensional biaxial shear test, the valid average stress p' and deviator stress q are shown as follows:(14)p′=σ1+σ22,(15)q=σ1−σ2.For differentEI values of samples, Figure 17 shows the variation of the stress ratio q/p' with the axial strain ε1. In the initial stage, as the axial strain increases, the stress ratio of all samples rises significantly. After increasing to a peak value, the stress ratio slightly descends and tends to stabilize.Figure 17
The evolution of stress ratio versus axial strain.


#### 6.2.2. Analysis of Volumetric Deformation

In this section, the evolution of volumetric deformation is analyzed. Regarding the boundaries of the simulation being rigid walls, the axial strain and volume strain can be calculated as(16)ε1=h0−hh0,εv=v0−vv0,where h0 and h are the original height and the current height of the sample, respectively; v0 and v are the original volume and the current volume of the sample, respectively; negative volume strain indicates volume expansion. In order to study the critical state characteristics of the stabilized sample, all samples were sheared to the axial strain ε1 = 30%. Under this deformation, the typical conditions of the critical state are basically satisfied (i.e., the stress and volume remain constant with strain).As shown in Figure18, in the initial stage, the volume strain descends with the ascent of the axial strain. After the axial strain reaches about 1%, the volume strain reverses. As the axial strain increases to a range of 20% to 30% which is a critical state, the volume deformation of the sample remains constant, which indicates that all samples show shear-induced dilatancy deformation with strain softening. Moreover, as the sample gets a steady state, the volume strain shows a negative correlation with the EI value. Therefore, a sample with a smaller EI value has stronger dilatancy.Figure 18
The evolution of volume strain versus axial strain.


#### 6.2.3. Analysis of Mean Coordination Number

In this section, we focus on the mean coordination number (MCN), which can reflect the microscopic packing structure of the simulated GCS. The evolution of MCN versus axial strain is illustrated in Figure19 for samples with various values. In the beginning stage of shear, the MCN of samples decreases significantly with axial strain. Then, MCN declines with a lower speed when the axial strain is greater than 2%. When the axial strain is greater than 10%, the MCN is roughly steady. Then, as the EI value rises from 0.4 to 0.9, the MCN collapses from 5.0 to 4.0; thus, a negative correlation with the EI value can be found.Figure 19
The evolution of mean coordination number versus axial strain.


#### 6.2.4. Analysis of Sliding Contact Percentage

The Mohr–Coulomb law conducts the percentage of sliding contact. The sliding coefficient is shown as follows:(17)SC=ftcμfnc,where ftc is the tangential contact force of contact c, fnc is the normal contact force of contact c, and μ is the friction factor. When SC > 0.9999, it is considered that a contact slip has occurred. The percentage of sliding contact in the particle system is as follows:(18)SP=NSCNC×100%,where NSC indicates the number of sliding contacts in the particle system and NC is the total number of contacts. Figure 20 shows the variation of SP of each sample with axial strain. The variation of SP with an axial strain of all samples is basically the same. As ε1 increases, SP rapidly increases to a peak value and then gradually decreases to a stable value. Figure 20 also indicates that the contact slip rate descends with EI climbing. As the EI value rises from 0.4 to 0.9, the maximum value and the stability value of SP descends from 0.36 to about 0.16, and from 11% to approximately 5%, respectively.Figure 20
The evolution of the percentage of sliding contact versus axial strain.


## 6.1. Simulation of Biaxial Compression Test

First, five samples with different elongation (EI = 0.4, 0.5, 0.6, 0.7, 0.9) were constructed. Each sample has an initial size of 9 m ∗ 18 m and contains approximately 5000 GCS particles. The linear elastic model is used as the contact law between particles. The samples are compressed by shifting the four rigid boundaries, and the boundary has a force at 100 kPa (stress boundary condition controlled by numerical servo). During the shearing process, in order to make the state of samples quasistatic (to make the loading process of the sample quasistatic loading), the rate of shear strain is supposed to be relatively small. Thus, the inertial parameters are introduced as follows:(13)Iintertia=ε˙1dσ0/ρ<10−3,where ε˙1 is the loading strain rate, d is mean particle size, ρ is material density, σ0 is the consolidation stress. Previous studies have shown that the value of the Inertia parameter should be less than 2.5 × 10−3. Thus, the shifting speed of the rigid wall is fixed to 0.05 m/s. The inertia parameter is less than 10−4 throughout the tests’ negative correlation.


## 6.2. Analysis of Macro-Meso Parameter

### 6.2.1. Analysis of Stress Ratio

The mean stressp and deviator stress q are calculated to explore the sample shear strength. In the two-dimensional biaxial shear test, the valid average stress p' and deviator stress q are shown as follows:(14)p′=σ1+σ22,(15)q=σ1−σ2.For differentEI values of samples, Figure 17 shows the variation of the stress ratio q/p' with the axial strain ε1. In the initial stage, as the axial strain increases, the stress ratio of all samples rises significantly. After increasing to a peak value, the stress ratio slightly descends and tends to stabilize.Figure 17
The evolution of stress ratio versus axial strain.


### 6.2.2. Analysis of Volumetric Deformation

In this section, the evolution of volumetric deformation is analyzed. Regarding the boundaries of the simulation being rigid walls, the axial strain and volume strain can be calculated as(16)ε1=h0−hh0,εv=v0−vv0,where h0 and h are the original height and the current height of the sample, respectively; v0 and v are the original volume and the current volume of the sample, respectively; negative volume strain indicates volume expansion. In order to study the critical state characteristics of the stabilized sample, all samples were sheared to the axial strain ε1 = 30%. Under this deformation, the typical conditions of the critical state are basically satisfied (i.e., the stress and volume remain constant with strain).As shown in Figure18, in the initial stage, the volume strain descends with the ascent of the axial strain. After the axial strain reaches about 1%, the volume strain reverses. As the axial strain increases to a range of 20% to 30% which is a critical state, the volume deformation of the sample remains constant, which indicates that all samples show shear-induced dilatancy deformation with strain softening. Moreover, as the sample gets a steady state, the volume strain shows a negative correlation with the EI value. Therefore, a sample with a smaller EI value has stronger dilatancy.Figure 18
The evolution of volume strain versus axial strain.


### 6.2.3. Analysis of Mean Coordination Number

In this section, we focus on the mean coordination number (MCN), which can reflect the microscopic packing structure of the simulated GCS. The evolution of MCN versus axial strain is illustrated in Figure19 for samples with various values. In the beginning stage of shear, the MCN of samples decreases significantly with axial strain. Then, MCN declines with a lower speed when the axial strain is greater than 2%. When the axial strain is greater than 10%, the MCN is roughly steady. Then, as the EI value rises from 0.4 to 0.9, the MCN collapses from 5.0 to 4.0; thus, a negative correlation with the EI value can be found.Figure 19
The evolution of mean coordination number versus axial strain.


### 6.2.4. Analysis of Sliding Contact Percentage

The Mohr–Coulomb law conducts the percentage of sliding contact. The sliding coefficient is shown as follows:(17)SC=ftcμfnc,where ftc is the tangential contact force of contact c, fnc is the normal contact force of contact c, and μ is the friction factor. When SC > 0.9999, it is considered that a contact slip has occurred. The percentage of sliding contact in the particle system is as follows:(18)SP=NSCNC×100%,where NSC indicates the number of sliding contacts in the particle system and NC is the total number of contacts. Figure 20 shows the variation of SP of each sample with axial strain. The variation of SP with an axial strain of all samples is basically the same. As ε1 increases, SP rapidly increases to a peak value and then gradually decreases to a stable value. Figure 20 also indicates that the contact slip rate descends with EI climbing. As the EI value rises from 0.4 to 0.9, the maximum value and the stability value of SP descends from 0.36 to about 0.16, and from 11% to approximately 5%, respectively.Figure 20
The evolution of the percentage of sliding contact versus axial strain.


## 6.2.1. Analysis of Stress Ratio

The mean stressp and deviator stress q are calculated to explore the sample shear strength. In the two-dimensional biaxial shear test, the valid average stress p' and deviator stress q are shown as follows:(14)p′=σ1+σ22,(15)q=σ1−σ2.For differentEI values of samples, Figure 17 shows the variation of the stress ratio q/p' with the axial strain ε1. In the initial stage, as the axial strain increases, the stress ratio of all samples rises significantly. After increasing to a peak value, the stress ratio slightly descends and tends to stabilize.Figure 17
The evolution of stress ratio versus axial strain.


## 6.2.2. Analysis of Volumetric Deformation

In this section, the evolution of volumetric deformation is analyzed. Regarding the boundaries of the simulation being rigid walls, the axial strain and volume strain can be calculated as(16)ε1=h0−hh0,εv=v0−vv0,where h0 and h are the original height and the current height of the sample, respectively; v0 and v are the original volume and the current volume of the sample, respectively; negative volume strain indicates volume expansion. In order to study the critical state characteristics of the stabilized sample, all samples were sheared to the axial strain ε1 = 30%. Under this deformation, the typical conditions of the critical state are basically satisfied (i.e., the stress and volume remain constant with strain).As shown in Figure18, in the initial stage, the volume strain descends with the ascent of the axial strain. After the axial strain reaches about 1%, the volume strain reverses. As the axial strain increases to a range of 20% to 30% which is a critical state, the volume deformation of the sample remains constant, which indicates that all samples show shear-induced dilatancy deformation with strain softening. Moreover, as the sample gets a steady state, the volume strain shows a negative correlation with the EI value. Therefore, a sample with a smaller EI value has stronger dilatancy.Figure 18
The evolution of volume strain versus axial strain.


## 6.2.3. Analysis of Mean Coordination Number

In this section, we focus on the mean coordination number (MCN), which can reflect the microscopic packing structure of the simulated GCS. The evolution of MCN versus axial strain is illustrated in Figure19 for samples with various values. In the beginning stage of shear, the MCN of samples decreases significantly with axial strain. Then, MCN declines with a lower speed when the axial strain is greater than 2%. When the axial strain is greater than 10%, the MCN is roughly steady. Then, as the EI value rises from 0.4 to 0.9, the MCN collapses from 5.0 to 4.0; thus, a negative correlation with the EI value can be found.Figure 19
The evolution of mean coordination number versus axial strain.


## 6.2.4. Analysis of Sliding Contact Percentage

The Mohr–Coulomb law conducts the percentage of sliding contact. The sliding coefficient is shown as follows:(17)SC=ftcμfnc,where ftc is the tangential contact force of contact c, fnc is the normal contact force of contact c, and μ is the friction factor. When SC > 0.9999, it is considered that a contact slip has occurred. The percentage of sliding contact in the particle system is as follows:(18)SP=NSCNC×100%,where NSC indicates the number of sliding contacts in the particle system and NC is the total number of contacts. Figure 20 shows the variation of SP of each sample with axial strain. The variation of SP with an axial strain of all samples is basically the same. As ε1 increases, SP rapidly increases to a peak value and then gradually decreases to a stable value. Figure 20 also indicates that the contact slip rate descends with EI climbing. As the EI value rises from 0.4 to 0.9, the maximum value and the stability value of SP descends from 0.36 to about 0.16, and from 11% to approximately 5%, respectively.Figure 20
The evolution of the percentage of sliding contact versus axial strain.


## 7. Conclusion

Graded crushed stone (GCS) is of great significance in road construction as an essential material used in the base or transition layer. The packing properties like compactness and void ratio can be affected by the irregularity and variability of particle shape. In this study, the particle outlines are firstly extracted based on deep learning algorithms. Then, the shape evaluation indexes are quantified. In addition, the gravity deposition with distinctEI value is simulated by the discrete element method. Finally, the shear behavior of GCS is studied. The specific results are as follows:(1)
The realistic particle contour is first automatically extracted based on digital image processing and deep learning algorithms. The lightweight U-NET model is trained and the results show the accuracy of the training set reached a value of 97% and the loss value is about 0.01.(2)
The elongation (EI), roundness (Rd), and roughness (Rg) of GCS are quantified by shape evaluation algorithms.(3)
Based on the establishment of the GCS shape library, the gravity deposition with various elongation is simulated using the discrete element method to study the packing properties of GCS. The elongation effects on the macroscopic and microscopic quantities are explored. The results illustrate that the void ratio decreases with increasing elongation before a reflection point whereEI value is 0.9. Similar to the void ratio, the mean coordinate number Cmean exhibits a decreasing trend with EI increasing. After reaching the bottom of the curve, Cmean stays stable and has a slight increase.(4)
The shear behavior of GCS is studied by simulation of the biaxial compression test. The stress ratio of all samples increases rapidly, and then reaches the peak stage, after which the stress ratio gradually decreases and tends to stabilize. The volume strain decreases with the increase of the axial strain, and when the axial strain reaches about 1%, the volume strain reverses. Then, in the range of 20% to 30%, the volume strain stays constant, indicating a shear-induced dilatancy deformation accompanied by strain softening. The mean coordinate number of all samples increases with axial strain at the beginning; then, after the value of axial strain reaching 2%, the MCN reaches its peak and decreases to be stable when the axial strain is greater than 10%. The variation of the percentage of sliding contactSP with the axial strain of all samples is basically the same. As ε1 increases, SP rapidly increases to a peak value and then gradually decreases to a stable value.



---

*Source: 3463363-2021-10-07.xml*