# Automated Particle Shape Identification and Quantification for DEM Simulation of Rockfill Materials in Subgrade Construction

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

---

## Abstract

Rockfill materials, conducted by impermeable stone, are frequently used in subgrade construction projects. The irregularity and variability of particle shape are demonstrated to affect the mechanical properties of rockfill subgrade, such as void ratio and coordination number. This study first identifies the subgrade rockfill particle contour by machine learning algorithms, including AdaBoost, Cascade, and sliding windows. Then, the shape evaluation indexes of length flatness, edge angle, and roughness are quantified, and the statistical analysis of each index is presented. In addition, the discrete element method (DEM) simulation is implemented on the compaction of rockfill subgrade to explore the impact of roundness on characteristics of particles. Finally, the macroanalysis on the void ratio and cumulative settlement and the microanalysis on particle coordination number, rotation momentum, and displacement are studied. The results illustrate that roundness has a significant effect on the mechanical characteristics of subgrade rockfill materials. With the increase of rolling passes, the porosity of packing decreases, whereas the settlement increases gradually. The change rate starts fast and ends slowly.


---

## Body

## 1. Introduction

Rockfill materials have excellent engineering characteristics such as good compaction, strong water permeability, high filling density, high shear strength, small subsidence and deformation, high bearing capacity, and not easy to liquefy under the action of seismic load. Rockfill materials are applied in construction fields frequently for their excellent characteristics [1], meeting the requirements of in situ materials. In subgrade engineering, rockfill is mostly used, and it can make full use of local natural materials and adapt to different geological conditions.The rockfill material method is relatively simple and possesses good seismic performance. Also, the evaluation of seismic reliability is an essential factor that needs to be considered. Zhang et al. evaluated geometric anisotropy modeling and shear behavior of graded crushed rocks [2]. There is research showing the seismic reliability of rockfill in hydraulic engineering under stochastic earthquake excitation considering the strain-softening behavior of rockfill materials. A new and efficient methodology that combines generalized probability density evolution method (GPDEM) with a spectral representation-random function method is presented to assess the seismic reliability [3].To calculate the properties and state of rockfill materials more accurately, numerical simulation is a popular tool to make further research on rockfill materials [4]. For example, Li et al. researched by considering moisture and stress sensitivity of subgrade soils [5]. Hui et al. investigated the fracture evolution of asphalt mixture compared with acoustic emission [6]. The research found that particle breakage has some significant effects on the behavior of granular materials, and the shear strength and compressibility are also included. Raisianzadeh et al. [7] conducted a micromechanical study on particle breakage in 2D angular rockfill materials under biaxial compression loading using a combined DEM and extended finite element method (XFEM) approach. Song et al. [8] investigated the skeleton behaviors of open-graded friction courses using DEM. Zhou et al. [9] presented simulating tri-axial compression tests on numerical samples reveals the relationships among the micromechanical parameters and macroscopic response of rockfill materials.In engineering, the identification and quantification of particles help to realize the needs of material classification and risk prediction, intelligent construction and management, and computer deep learning. To explore the effect of particle shape, Li et al. simulated three dimensions of aggregate and asphalt mixture using parameterized shape and size gradation [10]. And Ma et al. used another discrete method by a microscopic modeling approach to simplify the discrete element’s shape by introducing irregular convex polyhedrons to reproduce the geometry-dependent behavior of rockfills [11]. In addition, in the research of Zhou et al., the deformable DEM was employed to study the effect of particle shape on the macroscopic mechanical behavior of rockfill materials [9], indicating that the particle shape has a significant impact on the macroscopic response of rockfill. On top of that, Nie et al. [12] explored the impacts of Fourier-based particle shape on the macro-mesoscopic shear responses of rockfill materials by numerical biaxial compression tests. However, these works on the impact of particle shape are limited, since most of them were conducted by numerical tests, which is time-consuming and expensive. Furthermore, the complete particle shape is not extracted and processed.The previous methods have achieved the identification of particles to a certain extent, but they have not studied the effects of particle size and shape in depth. In this study, the rockfill particle contour is firstly extracted based on machine learning algorithms such as AdaBoost, Cascade, and sliding windows. Then, the shape evaluation indexes are quantified and the statistical analysis of each index is presented. In addition, the DEM simulation is implemented on the deformation and compaction of rockfill. Based on the simulation, the macroanalysis on void ratio and cumulative settlement is given. Finally, the microanalysis on particle coordination number, rotation momentum, and displacement is studied. Overall, this study lays a foundation for further research on the effect of particle shape and size of rockfill.


## 2. Automated Shape Identification

In this study, low-precision photography equipment (such as mobile phones and digital cameras) is used to obtain multiple subgrade rockfill materials images to save time and cost, considering the characteristics of the large number and low-precision requirements of rockfill material samples, as shown in Figure1. The target of this study is to extract the irregular shape of subgrade rockfill materials.Figure 1
Schematic of subgrade rockfill particles.Due to the simplicity, high precision, and high computational efficiency of the algorithm, the improved Viola–Jones algorithm [13] is used to automatically identify subgrade rockfill particles.


### 2.1. Training Set Preparation

Firstly, the most common shape of subgrade rockfill particles is considered. All the images containing a particle have an aspect ratio of 1 : 1. About 1206 pictures of these rockfill particles are taken, and then, some of these particles are cropped and zoomed as a training dataset for the model. Each cropped image is a small square in the original image. Among them, the images containing the complete outline of the subgrade rockfill particle are marked as positive images, and the other images are marked as negative images. All the cropped images have a size of 32 × 32 pixels. Taking these images as input data, the computer can train the model and learn to obtain Haar-Like features [14]. Then, the feature results are stored in the classifier. The sample images are shown in Figure 2.Figure 2
(a) Samples of positive images; (b) samples of negative images.
(a)(b)


### 2.2. Scaling and Rotation of Training Images

Mostly, the aspect ratio of subgrade rockfill particles is 1 : 1 to 2 : 1. Therefore, for particles beyond this range of ratio, it is difficult to collect a large amount of image data by photographing. However, there exist some real rockfill particles with a high aspect ratio, so the classifier is expected to detect these particles as well. Therefore, a new method of acquiring a large number of images is needed. In this study, to quickly obtain a large number of slender rockfill particle images that can be used for training, a stretching method is used to digitally stretch the existing images with an aspect ratio of 1 : 1 [10]. By adjusting the size of the previous 1 : 1 training image (32 × 32 pixels), 1.5 : 1 training image (48 × 32 pixels), 2 : 1 training image (64 × 32 pixels), 3 : 1 training image (96 × 32 pixels), and 4 : 1 training images (128 × 32 pixels) can be obtained.The aspect ratio of most railway subgrade rockfill particles is 1 : 1 to 2 : 1. Therefore, for rockfill that is not in this ratio, it is difficult to collect a large amount of image data by shooting. However, these high aspect ratio rockfill particles are real, so the classifier also needs to detect this rockfill. Therefore, we need a new method of acquiring a large number of images. To quickly obtain a large number of slender rockfill particle images that can be used for training, this article digitally stretches the existing images with a rockfill aspect ratio of 1 : 1 [15] by adjusting the size of the previous 1 : 1 training image (32 × 32 pixels), 1.5 : 1 training image (48 × 32 pixels), 2 : 1 training image (64 × 32 pixels), 3 : 1 training image (96 × 32 pixels) and 4 : 1 training images (128 × 32 pixels), as shown in Figure 3.Figure 3
Schematic of stretching images: (a) 1.5 : 1 image; (b) 2 : 1 image; (c) 3 : 1 image; (d) 4 : 1 image.
(a)(b)(c)(d)In addition, to recognize slender subgrade rockfill particles in any direction, the original image is rotated counterclockwise in 15° increments, as shown in Figure4. For each rotated image, the sliding window scans the entire image in turn. The sliding window method may not be able to detect some subgrade rockfill particles due to the oblique placement of these particles. However, if the image is rotated, these particles will be converted to the vertical direction, which is easy to identify. After rotating the image by 15°, 30°, 45°, 60°, 75°, 90°, 105°, 120°, 135°, 150°, and 165°, 13 rotated images can be obtained. Figure 4 shows a simple example of 0°, 30°, and 60°.Figure 4
Schematic of rotation: (a) 0°; (b) 30°; (c) 60°.
(a)(b)(c)Then, all the results of these images are merged. In each image, if a subgrade rockfill particle is successfully identified, a bounding box is added. Then, all bounding boxes will be superimposed when merging. However, due to the rotation, the same rockfill particle will generate some different bounding boxes. Thus, the image is rotated to the 0° position, which is equivalent to only rotating the sliding window. Then, the bounding box with the smallest area is selected as the final bounding box of this rockfill particle, as shown in Figure5.Figure 5
(a) The merges results of all rotation; (b) the average window of particles.
(a)(b)


### 2.3. AdaBoost and Cascade Classifier

All the labeled images are used to train the model to obtain the feature difference between the positive and negative images in a dataset. Firstly, AdaBoost and Cascade methods [16] are used. Then, through the sliding window method [17], the subgrade rockfill particles can be successfully identified.AdaBoost is a classification algorithm that can integrate some weak classifiers into a strong classifier. This strong classifier can be used in the Cascade method. The classifier is used to test an instance, and each classifierft will output the probability distribution of two classes of positive or negative. These probabilities are combined in proportion to the relevant multipliers. The detailed algorithm is shown in Figure 6.Figure 6
AdaBoost algorithm.The training process of the Cascade method is divided into multiple stages. At each stage, a strong classifier is trained through the AdaBoost algorithm. In addition, three parameters need to be set in advance: the number of stagess, the detection rate d, and the false positive rate f, where d represents the minimum percentage of correct detections and f represents the maximum percentage of negative images that are incorrectly identified as positive. In each stage, it is necessary to search for the smallest T in AdaBoost algorithm, so that the detection rate of the strong classifier C(x) trained by AdaBoost is greater than d and the false positive rate is less than f.In the hypothetical example of Figure7, 100 positive images and 100 negative images are used as input data. The original parameters of each stage are s = 3, d = 0.99, and f = 0.25. Then, in each stage, a strong classifier with appropriate T is trained through AdaBoost, with the value of T 3, 10, and 35, respectively. In the first stage, the model has 200 input images, and the classifier meets the requirements that the detection rate is greater than d and the false positive rate is less than f. Then, the second stage is similarly performed. In stage 3, the total detection rate D = d3 = 0.97 and the false positive rate F = f3 = 0.015 can be obtained.Figure 7
A hypothetical example of the Cascade algorithm.


### 2.4. Identifying Particles Based on the Sliding Window Method

Through the previously obtained classifier, images containing real particles can be obtained. The next task is to identify the complete contour of subgrade rockfill particles. Given a new image of rockfill particles, at each scanning position, the area in the sliding window is input to the detector using a sliding window algorithm. The detector will find and extract features in the window area of the image. Based on these features, the classifier finally determines whether the image is a positive image or a negative image. If it is a positive image, a bounding box is added to the scanning position to identify it as an area containing complete outline particles. The specific instructions are as follows.Since the smallest pixel of the particle is about 50 × 50, a 50 × 50 sliding window is first used to scan the entire image. In each window, the crop image is input to the classifier. Since the previous classifiers are trained based on 32 × 32 images, cropped images larger than 32 × 32 should be scaled to 32 × 32, and theC1(x) classifier in the Cascade algorithm is tested. If the output is positive, move to stage 2 of the Cascade method. Otherwise, the negative images are rejected, and the sliding window will move to the next position by 1 pixel. If the cropped image successfully passes all s classifiers, it will be a positive image containing the complete contour rockfill particles with a high probability. Then, the sliding window at this position is marked as the bounding box of the particle. Next, continue to move the sliding window until the entire image is scanned. The moving increment of the sliding window in the horizontal direction and the vertical direction is 1 pixel, and the moving path is shown in Figure 8(a).Figure 8
(a) The minimum sliding window; (b) the increment is 5 pixels.
(a)(b)When the 50 × 50 sliding window successfully scans the entire image, we add 5 pixels to the size of the sliding window in both horizontal and vertical directions, as shown in Figure8(b). The previous sliding process is repeated until the sliding window expands to the maximum size of 300 × 300 pixels. Since most of the cropped images will be classified as negative images and rejected, only a small number of images can successfully pass all s classifiers, so the algorithm has high computational efficiency.Note that when the sliding window is large enough, several movements of the sliding window may contain the same particle, resulting in multiple bounding boxes. In this case, we take the average value as a new single window, as shown in Figure9.Figure 9
(a) Several bounding box of a particle; (b) the average bounding box.
(a)(b)


## 2.1. Training Set Preparation

Firstly, the most common shape of subgrade rockfill particles is considered. All the images containing a particle have an aspect ratio of 1 : 1. About 1206 pictures of these rockfill particles are taken, and then, some of these particles are cropped and zoomed as a training dataset for the model. Each cropped image is a small square in the original image. Among them, the images containing the complete outline of the subgrade rockfill particle are marked as positive images, and the other images are marked as negative images. All the cropped images have a size of 32 × 32 pixels. Taking these images as input data, the computer can train the model and learn to obtain Haar-Like features [14]. Then, the feature results are stored in the classifier. The sample images are shown in Figure 2.Figure 2
(a) Samples of positive images; (b) samples of negative images.
(a)(b)


## 2.2. Scaling and Rotation of Training Images

Mostly, the aspect ratio of subgrade rockfill particles is 1 : 1 to 2 : 1. Therefore, for particles beyond this range of ratio, it is difficult to collect a large amount of image data by photographing. However, there exist some real rockfill particles with a high aspect ratio, so the classifier is expected to detect these particles as well. Therefore, a new method of acquiring a large number of images is needed. In this study, to quickly obtain a large number of slender rockfill particle images that can be used for training, a stretching method is used to digitally stretch the existing images with an aspect ratio of 1 : 1 [10]. By adjusting the size of the previous 1 : 1 training image (32 × 32 pixels), 1.5 : 1 training image (48 × 32 pixels), 2 : 1 training image (64 × 32 pixels), 3 : 1 training image (96 × 32 pixels), and 4 : 1 training images (128 × 32 pixels) can be obtained.The aspect ratio of most railway subgrade rockfill particles is 1 : 1 to 2 : 1. Therefore, for rockfill that is not in this ratio, it is difficult to collect a large amount of image data by shooting. However, these high aspect ratio rockfill particles are real, so the classifier also needs to detect this rockfill. Therefore, we need a new method of acquiring a large number of images. To quickly obtain a large number of slender rockfill particle images that can be used for training, this article digitally stretches the existing images with a rockfill aspect ratio of 1 : 1 [15] by adjusting the size of the previous 1 : 1 training image (32 × 32 pixels), 1.5 : 1 training image (48 × 32 pixels), 2 : 1 training image (64 × 32 pixels), 3 : 1 training image (96 × 32 pixels) and 4 : 1 training images (128 × 32 pixels), as shown in Figure 3.Figure 3
Schematic of stretching images: (a) 1.5 : 1 image; (b) 2 : 1 image; (c) 3 : 1 image; (d) 4 : 1 image.
(a)(b)(c)(d)In addition, to recognize slender subgrade rockfill particles in any direction, the original image is rotated counterclockwise in 15° increments, as shown in Figure4. For each rotated image, the sliding window scans the entire image in turn. The sliding window method may not be able to detect some subgrade rockfill particles due to the oblique placement of these particles. However, if the image is rotated, these particles will be converted to the vertical direction, which is easy to identify. After rotating the image by 15°, 30°, 45°, 60°, 75°, 90°, 105°, 120°, 135°, 150°, and 165°, 13 rotated images can be obtained. Figure 4 shows a simple example of 0°, 30°, and 60°.Figure 4
Schematic of rotation: (a) 0°; (b) 30°; (c) 60°.
(a)(b)(c)Then, all the results of these images are merged. In each image, if a subgrade rockfill particle is successfully identified, a bounding box is added. Then, all bounding boxes will be superimposed when merging. However, due to the rotation, the same rockfill particle will generate some different bounding boxes. Thus, the image is rotated to the 0° position, which is equivalent to only rotating the sliding window. Then, the bounding box with the smallest area is selected as the final bounding box of this rockfill particle, as shown in Figure5.Figure 5
(a) The merges results of all rotation; (b) the average window of particles.
(a)(b)


## 2.3. AdaBoost and Cascade Classifier

All the labeled images are used to train the model to obtain the feature difference between the positive and negative images in a dataset. Firstly, AdaBoost and Cascade methods [16] are used. Then, through the sliding window method [17], the subgrade rockfill particles can be successfully identified.AdaBoost is a classification algorithm that can integrate some weak classifiers into a strong classifier. This strong classifier can be used in the Cascade method. The classifier is used to test an instance, and each classifierft will output the probability distribution of two classes of positive or negative. These probabilities are combined in proportion to the relevant multipliers. The detailed algorithm is shown in Figure 6.Figure 6
AdaBoost algorithm.The training process of the Cascade method is divided into multiple stages. At each stage, a strong classifier is trained through the AdaBoost algorithm. In addition, three parameters need to be set in advance: the number of stagess, the detection rate d, and the false positive rate f, where d represents the minimum percentage of correct detections and f represents the maximum percentage of negative images that are incorrectly identified as positive. In each stage, it is necessary to search for the smallest T in AdaBoost algorithm, so that the detection rate of the strong classifier C(x) trained by AdaBoost is greater than d and the false positive rate is less than f.In the hypothetical example of Figure7, 100 positive images and 100 negative images are used as input data. The original parameters of each stage are s = 3, d = 0.99, and f = 0.25. Then, in each stage, a strong classifier with appropriate T is trained through AdaBoost, with the value of T 3, 10, and 35, respectively. In the first stage, the model has 200 input images, and the classifier meets the requirements that the detection rate is greater than d and the false positive rate is less than f. Then, the second stage is similarly performed. In stage 3, the total detection rate D = d3 = 0.97 and the false positive rate F = f3 = 0.015 can be obtained.Figure 7
A hypothetical example of the Cascade algorithm.


## 2.4. Identifying Particles Based on the Sliding Window Method

Through the previously obtained classifier, images containing real particles can be obtained. The next task is to identify the complete contour of subgrade rockfill particles. Given a new image of rockfill particles, at each scanning position, the area in the sliding window is input to the detector using a sliding window algorithm. The detector will find and extract features in the window area of the image. Based on these features, the classifier finally determines whether the image is a positive image or a negative image. If it is a positive image, a bounding box is added to the scanning position to identify it as an area containing complete outline particles. The specific instructions are as follows.Since the smallest pixel of the particle is about 50 × 50, a 50 × 50 sliding window is first used to scan the entire image. In each window, the crop image is input to the classifier. Since the previous classifiers are trained based on 32 × 32 images, cropped images larger than 32 × 32 should be scaled to 32 × 32, and theC1(x) classifier in the Cascade algorithm is tested. If the output is positive, move to stage 2 of the Cascade method. Otherwise, the negative images are rejected, and the sliding window will move to the next position by 1 pixel. If the cropped image successfully passes all s classifiers, it will be a positive image containing the complete contour rockfill particles with a high probability. Then, the sliding window at this position is marked as the bounding box of the particle. Next, continue to move the sliding window until the entire image is scanned. The moving increment of the sliding window in the horizontal direction and the vertical direction is 1 pixel, and the moving path is shown in Figure 8(a).Figure 8
(a) The minimum sliding window; (b) the increment is 5 pixels.
(a)(b)When the 50 × 50 sliding window successfully scans the entire image, we add 5 pixels to the size of the sliding window in both horizontal and vertical directions, as shown in Figure8(b). The previous sliding process is repeated until the sliding window expands to the maximum size of 300 × 300 pixels. Since most of the cropped images will be classified as negative images and rejected, only a small number of images can successfully pass all s classifiers, so the algorithm has high computational efficiency.Note that when the sliding window is large enough, several movements of the sliding window may contain the same particle, resulting in multiple bounding boxes. In this case, we take the average value as a new single window, as shown in Figure9.Figure 9
(a) Several bounding box of a particle; (b) the average bounding box.
(a)(b)


## 3. Quantitative Shape Analysis

### 3.1. Calculation of Length Flatness of Particles

Common first-scale morphological indicators of particles include flatness, form factor, sphericity, and others. They are mostly used to measure the similarity between the overall contour of particles and the contour of the standard graph of the same diameter (round, square particles, etc.). Also, they are the indicators that have the greatest influence on the overall morphology and mechanical behavior of particles. Among them, the length flatness is widely used in engineering because of its simple measurement and simple concept. Although the existing definitions of length flatness indexes are very diversified, some basic geometric properties of particles are often adopted by most indexes and become common factors of the index system (as shown in Table1). Therefore, the calculation of these basic indexes is the premise of the calculation of the length flatness index.Table 1
Basic particle morphology indicators that are widely used.
ParameterSymbolGeometric meaningParticle sizeArThe area bounded by a particle contourParticle circumferencePerParticle contour perimeterFeret’s diameterFera“Ruler length” of corresponding direction angleaLength scaleSmajorMaximum Feret’s diameterVertical scaleSmajor90Smajor90Feret’s diameter orthogonal to the direction corresponding to the long scale (short scale)Short scaleSminorMinimum Feret’s diameterFeret’s diameterFera of a particle describes the vernier caliper length of a particle at a certain direction angle, which is the basis of the index of particle length flatness. For the directional angle ai of a particular particle Op1, the Feret diameter Ferai corresponding to it can be given by the following equation:(1)Op1′=i1′j1′i2′j2′i3′j3′⋯⋯ik′jk′=i1j1i2j2i3j3⋯⋯ikjkcos−aisin−ai−sin−aicos−ai,(2)Ferai=maxi1′,i2′…ik′−mini1′,i2′…ik′,where the function of equation (1) is to rotate the particle so that the direction corresponding to ai is parallel to the X-axis.The corresponding Feret diameter can be easily obtained with known direction angleai using equation (2). Therefore, when calculating the maximum (minimum) Feret diameter of particles, the most important problem is to find their corresponding direction angle ai. The rotation method [18] is commonly used to solve this problem in existing research studies. The essence of the rotation method is to enumerate the direction angles between 0 and 180° with a certain angle step length and obtain the approximate maximum (minimum) Feret diameter through comparison by calculating the Feret diameter corresponding to each direction angle. This process is represented geometrically by rotating at angular intervals and measuring the particle profile. In general, the rotation method is simple and intuitive, with high calculation accuracy when the rotation step size is small. It is not affected by the concave part of particles. Therefore, the rotation method can well accomplish the calculation of the particle Feret diameter. Since the rotation method has a certain error when the rotation step is large, the calculation time is relatively great.The calculation of the length flatness of particles can be divided into two basic ideas: the ratio of direct use of maximum (Smajor) and minimum (Smajor) Feret diameter and the aspect ratio of particle envelopment box. Among them, the definition of particle envelopment box is still controversial in academic circles. A more widely used definition is to take Smajor as the length of the enclosing box and define length flatness EI1=Sminor90/Smajor. However, this method will give incorrect results when calculating the particles of the rectangle (or approximately rectangle). The other definition builds the bounding box indicator EI2=Sminor/Sminor90 based on Smajor. This method tends to get a smaller bounding box, and it cannot obtain the correct results when calculating rhombic particles. In fact, since EI1 and EI2 contain Smajor and Smajor, they are also statistically strongly related to another commonly used metric, EI3=Sminor/Sminor.In this study, referring to existing studies [19–21], EI3 is selected as the index to evaluate the length and flatness of particles. A total of 991 real images of subgrade rockfill particles were extracted, and the calculation results of EI3 can be obtained, as shown in Figure 10. From the figure, it can be seen that the results of EI3 approximately follow a normal distribution, with a mean value of 0.73 and a width from 0.4 to 1.2.Figure 10
Statistical results of the length flatness index calculated by the rotation method.


### 3.2. Calculation of Edge Angles and Roundness of Particles

Barrett [22] expounded the advantages of the independence of grain morphology indexes among scales. However, in the existing research, the edge angle index is often difficult to maintain its independence from the first- and third-scale indexes. For example, RI calculation formula contains the basic indexes of the first scale, such as particle circumference and area, which are highly dependent on the morphological characteristics of the first scale. The AT index cannot guarantee its independence from the third-scale morphology because the grain angular and grain roughness effects are considered simultaneously. To solve this problem, Liu et al. [23] proposed an edge angle index AI based on a convex hull. The index is defined as the ratio of the equivalent elliptic perimeter of the particle to the circumference of the convex hull of the particle, where the equivalent ellipse of a particle is an ellipse with the same length flatness and area as the particle, and its semimajor axis (a) and semiminor axis (b) are given by the following equation:(3)a=ArπEI3,b=aEI3.AI indicator has the following characteristics:(1)
For smooth particles, their circumference should be close to its equivalent ellipse circumference. The more trenchant the edges and corners particles are, the longer the perimeter will be, leading toAI value declining. Therefore, AI indicator can be used to reflect the degree of angular.(2)
By convex hull, the influence of particle contour roughness is excluded and the independence with the third-scale index is maintained.(3)
By comparing the equivalent ellipse with the same length flatness, the influence of length flatness is normalized and the dependence with the first-scale index is reduced.The characteristics (2) and (3) guarantee the independence of this index, which is an advantage compared with other indexes. However, when the convex hull is used to exclude the roughness of the contour, the angular fluctuation of the particle will also be excluded, so that the concave part of the particle cannot be considered in this index. Therefore, this index has been modified in this study. Firstly, B-spline resampling technique was used to process the original particle outline, eliminating particle roughness and retaining edges and corners. Then, the improved edge angleAI2 is obtained by directly comparing the smoothed particle circumference pers with the equivalent ellipse circumference Pere:(4)AI2=PerePers.A total of 991 real subgrade rockfill particles is used to calculate the contour edge angle values, and the statistical curve is drawn in Figure11. It can be seen that the edge angle approximately follows a negative skewness distribution.Figure 11
Statistical results of edge angles and improved posterior edge angles.


### 3.3. Calculation of Particle Contour Roughness

Particle roughness describes the degree of fine fluctuation of the particle surface. In this section, the conventional particle roughness based on polar radius is improved. By using B-spline resampling to calculate the particle smooth and normal vector, an index for quantitative analysis of particle contour roughness was proposed. Based on the above three-segment intersection determination algorithm, a roughness index with scale invariance is proposed in this study. The calculation is as follows:(1)
B-spline resampling is carried out for the particle contour, and the normal vector of each sampling point is calculated(2)
A line segment is made at the sampling point along the normal vector, so that its length is equal to the radiusReq of the equal-area circle (to ensure that the line segment has an intersection)(3)
The line segment is judged to be intersecting with the original particle contour, and the distancedisti from the nearest intersection point to the sampling point is calculated(4)
Repeat steps (2) and (3) for all sampling pointsBased on the above steps, roughness indexRg can be obtained as follows:(5)Rg=∑i=1ndistinReq,where Req ensures the scale invariability of the index; that is, Rg remains unchanged when the particle contour is equilateral magnified.The angle values of edges of 991 real subgrade rockfill particle contours are calculated, and the statistical curves are plotted in Figure12, which is approximately a positive skewness distribution.Figure 12
Statistical results of roughness index.


## 3.1. Calculation of Length Flatness of Particles

Common first-scale morphological indicators of particles include flatness, form factor, sphericity, and others. They are mostly used to measure the similarity between the overall contour of particles and the contour of the standard graph of the same diameter (round, square particles, etc.). Also, they are the indicators that have the greatest influence on the overall morphology and mechanical behavior of particles. Among them, the length flatness is widely used in engineering because of its simple measurement and simple concept. Although the existing definitions of length flatness indexes are very diversified, some basic geometric properties of particles are often adopted by most indexes and become common factors of the index system (as shown in Table1). Therefore, the calculation of these basic indexes is the premise of the calculation of the length flatness index.Table 1
Basic particle morphology indicators that are widely used.
ParameterSymbolGeometric meaningParticle sizeArThe area bounded by a particle contourParticle circumferencePerParticle contour perimeterFeret’s diameterFera“Ruler length” of corresponding direction angleaLength scaleSmajorMaximum Feret’s diameterVertical scaleSmajor90Smajor90Feret’s diameter orthogonal to the direction corresponding to the long scale (short scale)Short scaleSminorMinimum Feret’s diameterFeret’s diameterFera of a particle describes the vernier caliper length of a particle at a certain direction angle, which is the basis of the index of particle length flatness. For the directional angle ai of a particular particle Op1, the Feret diameter Ferai corresponding to it can be given by the following equation:(1)Op1′=i1′j1′i2′j2′i3′j3′⋯⋯ik′jk′=i1j1i2j2i3j3⋯⋯ikjkcos−aisin−ai−sin−aicos−ai,(2)Ferai=maxi1′,i2′…ik′−mini1′,i2′…ik′,where the function of equation (1) is to rotate the particle so that the direction corresponding to ai is parallel to the X-axis.The corresponding Feret diameter can be easily obtained with known direction angleai using equation (2). Therefore, when calculating the maximum (minimum) Feret diameter of particles, the most important problem is to find their corresponding direction angle ai. The rotation method [18] is commonly used to solve this problem in existing research studies. The essence of the rotation method is to enumerate the direction angles between 0 and 180° with a certain angle step length and obtain the approximate maximum (minimum) Feret diameter through comparison by calculating the Feret diameter corresponding to each direction angle. This process is represented geometrically by rotating at angular intervals and measuring the particle profile. In general, the rotation method is simple and intuitive, with high calculation accuracy when the rotation step size is small. It is not affected by the concave part of particles. Therefore, the rotation method can well accomplish the calculation of the particle Feret diameter. Since the rotation method has a certain error when the rotation step is large, the calculation time is relatively great.The calculation of the length flatness of particles can be divided into two basic ideas: the ratio of direct use of maximum (Smajor) and minimum (Smajor) Feret diameter and the aspect ratio of particle envelopment box. Among them, the definition of particle envelopment box is still controversial in academic circles. A more widely used definition is to take Smajor as the length of the enclosing box and define length flatness EI1=Sminor90/Smajor. However, this method will give incorrect results when calculating the particles of the rectangle (or approximately rectangle). The other definition builds the bounding box indicator EI2=Sminor/Sminor90 based on Smajor. This method tends to get a smaller bounding box, and it cannot obtain the correct results when calculating rhombic particles. In fact, since EI1 and EI2 contain Smajor and Smajor, they are also statistically strongly related to another commonly used metric, EI3=Sminor/Sminor.In this study, referring to existing studies [19–21], EI3 is selected as the index to evaluate the length and flatness of particles. A total of 991 real images of subgrade rockfill particles were extracted, and the calculation results of EI3 can be obtained, as shown in Figure 10. From the figure, it can be seen that the results of EI3 approximately follow a normal distribution, with a mean value of 0.73 and a width from 0.4 to 1.2.Figure 10
Statistical results of the length flatness index calculated by the rotation method.


## 3.2. Calculation of Edge Angles and Roundness of Particles

Barrett [22] expounded the advantages of the independence of grain morphology indexes among scales. However, in the existing research, the edge angle index is often difficult to maintain its independence from the first- and third-scale indexes. For example, RI calculation formula contains the basic indexes of the first scale, such as particle circumference and area, which are highly dependent on the morphological characteristics of the first scale. The AT index cannot guarantee its independence from the third-scale morphology because the grain angular and grain roughness effects are considered simultaneously. To solve this problem, Liu et al. [23] proposed an edge angle index AI based on a convex hull. The index is defined as the ratio of the equivalent elliptic perimeter of the particle to the circumference of the convex hull of the particle, where the equivalent ellipse of a particle is an ellipse with the same length flatness and area as the particle, and its semimajor axis (a) and semiminor axis (b) are given by the following equation:(3)a=ArπEI3,b=aEI3.AI indicator has the following characteristics:(1)
For smooth particles, their circumference should be close to its equivalent ellipse circumference. The more trenchant the edges and corners particles are, the longer the perimeter will be, leading toAI value declining. Therefore, AI indicator can be used to reflect the degree of angular.(2)
By convex hull, the influence of particle contour roughness is excluded and the independence with the third-scale index is maintained.(3)
By comparing the equivalent ellipse with the same length flatness, the influence of length flatness is normalized and the dependence with the first-scale index is reduced.The characteristics (2) and (3) guarantee the independence of this index, which is an advantage compared with other indexes. However, when the convex hull is used to exclude the roughness of the contour, the angular fluctuation of the particle will also be excluded, so that the concave part of the particle cannot be considered in this index. Therefore, this index has been modified in this study. Firstly, B-spline resampling technique was used to process the original particle outline, eliminating particle roughness and retaining edges and corners. Then, the improved edge angleAI2 is obtained by directly comparing the smoothed particle circumference pers with the equivalent ellipse circumference Pere:(4)AI2=PerePers.A total of 991 real subgrade rockfill particles is used to calculate the contour edge angle values, and the statistical curve is drawn in Figure11. It can be seen that the edge angle approximately follows a negative skewness distribution.Figure 11
Statistical results of edge angles and improved posterior edge angles.


## 3.3. Calculation of Particle Contour Roughness

Particle roughness describes the degree of fine fluctuation of the particle surface. In this section, the conventional particle roughness based on polar radius is improved. By using B-spline resampling to calculate the particle smooth and normal vector, an index for quantitative analysis of particle contour roughness was proposed. Based on the above three-segment intersection determination algorithm, a roughness index with scale invariance is proposed in this study. The calculation is as follows:(1)
B-spline resampling is carried out for the particle contour, and the normal vector of each sampling point is calculated(2)
A line segment is made at the sampling point along the normal vector, so that its length is equal to the radiusReq of the equal-area circle (to ensure that the line segment has an intersection)(3)
The line segment is judged to be intersecting with the original particle contour, and the distancedisti from the nearest intersection point to the sampling point is calculated(4)
Repeat steps (2) and (3) for all sampling pointsBased on the above steps, roughness indexRg can be obtained as follows:(5)Rg=∑i=1ndistinReq,where Req ensures the scale invariability of the index; that is, Rg remains unchanged when the particle contour is equilateral magnified.The angle values of edges of 991 real subgrade rockfill particle contours are calculated, and the statistical curves are plotted in Figure12, which is approximately a positive skewness distribution.Figure 12
Statistical results of roughness index.


## 4. Application to DEM Simulation

### 4.1. Simulation of Earth Compaction of Rockfill Subgrade

Affected by the actual construction scale of the subgrade project, particle gradation and the irregular shape particles are adopted in the simulation, and the adoption of the 3D model will greatly increase the workload of the computer. Therefore, the two-dimensional model can also be used to simulate and analyze the particle motion characteristics of rockfill in the compaction process from a microscopic perspective.In this study, discrete element software PFC2D was used to simulate the compaction process of rockfill in two dimensions. To ensure that the roller vibrator rolls rather than sliding during the movement, the relationship between the angular velocityw and the horizontal velocity v of the vibrator is shown in the following equation:(6)v=w×R,where, v,w, and R are the horizontal velocity, angular velocity, and radius of the vibration wheel, respectively.The simulated material is rockfill for subgrade construction, and the linear elastic contact model is adopted in the discrete element model, which is widely used in similar studies [24]. The mesoscopic parameters used in the simulation refer to the research results of Liu et al. [24], and the specific parameters are shown in Table 2.Table 2
Microscale parameters used in the DEM simulations.
Microscopic parametersThe normal stiffness,kn (N/m)The tangential stiffness,ks (N/m)The coefficient of friction, fricLocal dampNormal/tangential dampingdpn, dptDensity (kg/m3)Value (particle)1.5 × 1072.9 × 1070.50.90.52600Value (wheel)1.0 × 1091.0 × 1090.30.00.57800Referring to the research of Goldenberg and Goldhirsch [25], Ma et al. [11], Zhou et al. [9], Alaei and Mahboubi [26] and Nie et al. [12], in this study, 5 groups of two-dimensional particles with different edge angles were generated, with 8 samples in each group, as shown in Figure 13. To exclude the influence of elongation on the simulation results, the elongation of particles was controlled between 0.97 and 1.00. Note that particles of roundness 1 are standard circles and are not shown in Figure 13. In addition, the effect of particle roughness was not considered in this study, that is, D8 = 0.0.Figure 13
Two-dimensional particle models with different roundnesses.In PFC2D, the clump function can realize the fitting of complex particle contour, with the fitting accuracy control parameters distance (range 0–180) and ratio (range 0.0–1.0). The larger the distance, the smaller the ratio, and the fitted results are closer to the real model. However, more computational memory and running time are consumed. In this study, the distance was set as 162, and the ratio was set as 0.25. Figure14 shows a clump model filled and generated under this setting. It can be seen that the set of two parameters is sufficient to reflect the external contour of the particle.Figure 14
A DEM-simulated virtual particle from a Fourier-based particle contour.Coarse particle packing with a diameter greater than 10 mm was used in this simulation to improve the calculation efficiency, and the distribution of the size of the packing was calculated. The results show that the content percentage of diameter ranging from 10 to 20 mm, 20 to 30 mm, 30 to 40 mm, 40 to 50 mm, and 50 to 60 mm is 17%, 20%, 25%, 20%, and 18%, respectively.The packing model was set in the range of height (H) × width (W) = 0.8 m × 8 m. The initial model is generated by gravity sedimentation. The process is as follows:(1)
In the 2H × W range, contactless particles slightly larger than the total volume of the target were randomly generated according to the particle size range and content, as shown in Figure 15(a).(2)
Under the action of gravity, particles fall freely, when the particle system reaches a quasi-static state and stops running, as shown in Figure15(b).(3)
In order to ensure that fillers with different roundnesses are compacted at the same thickness, redundant particles above the thickness line are removed to simulate the leveling process. The final leveling filler is shown in Figure15(c).(4)
A rigid body circular element is generated to simulate the roller. The movement parameters and attribute parameters of the roller are determined in Tables2 and 3. The chord wave vibration load is applied to the filler during the movement of the roller, as shown in Figure 15(d). And the simulation is stopped after repeated rolling 12 times.Figure 15
PFC simulation of vibratory compaction, including dumping, spreading, and compaction. (a) The contactless model; (b) sedimentation accumulation; (c) leveling; (d) vibratory compaction.
(a)(b)(c)(d)Table 3
External load parameters.
Rolling wheel parametersQuality of grinding wheel,G (kg)Grinding wheel diameter,D (mm)Grinding wheel width,b (mm)Vibration frequency,f (Hz)Exciting force,F0 (kN)Speed,v (m/s)Rolling passThe actual parameter1 × 1041500200040300112Simulation parameters5 × 1031500100040150112


### 4.2. The Macroresult

Figure16(a) shows the porosity of fillers with different roundnesses varying with the number of rolls. Obviously, with the increase in the number of rolling times, the pores first decreased sharply and then remained unchanged. But the correlation between roundness and vibratory compaction was not strong. The initial porosity was extracted (shown in Figure 16(b)), and it was found that with the increase of roundness, the initial porosity first decreased and then increased, which was consistent with the previous research results related to particle roundness [27]. Figure 16(c) shows the normalized porosity rate ∆n obtained after conversion using equation (7). It is found that the small roundness of the particle indicates irregular shape and great drop in the particle porosity compared to the initial state. The normalized porosity rate is obtained as(7)Δn=n0−nin0×100%,where n0 and ni represent the porosity corresponding to the initial porosity and the porosity corresponding to the ith pass after rolling, respectively.Figure 16
The relationship between the compactor pass and pore change of fillings with different roundnesses. (a) Porosity; (b) initial porosity; (c) normalized void change rate.
(a)(b)(c)Figure17 shows the cumulative settlement amount of the pavement corresponding to each roll pass of the filler, which shows a trend similar to the normalized pore change rate. It shows a trend of rapid settlement in the early stage and a trend of slow settlement in the later stage. This trend is consistent with the actual change in the monitored value in construction [24]. The indirect results showed that the particles sank under the vibratory compaction condition, which made the fillers squeeze and approach each other and rearrange the particles. This results in a smaller porosity, showing the changing trend as shown in Figure 16(a). In addition, a small value of roundness of the particles represents the great settlement of the fillers. Therefore, in practical engineering, the irregular shape of the packing requires more rolling times.Figure 17
Cumulative settlement-compactor pass responses for fillings with different roundnesses.


### 4.3. The Micromechanical Results

#### 4.3.1. Coordination Number

Coordination number is an important index to characterize the internal structure of particle systems. For the irregular particle model composed of Clump, there are two different definitions of coordination number: one is the number of particles in contact with the specified particles, called the coordination number of particles; and the other is the number of contacts belonging to a given particle, called the contact coordination number.Thornton [28] pointed out that there are particles in granular materials that do not contact with other particles or only contact with one particle, and such particles have no contribution to the overall stability of the material. Therefore, the mechanical average particle coordination number is defined as follows:(8)CN=2C−N1N−N0−N1,where C is the total number of particle contact, N is the total number of particles, and N0 and N1 are the numbers of particles in contact with zero or one particle in the particle system. Figure 18 shows the relationship between the mechanical particle coordination number and the rolling pass number of fillers with different roundnesses. For each roundness case, the particle coordination number showed a slight upward trend as the number of rolling passes increased; that is, the number of particles in contact with a certain particle increased. This was related to the overall looseness of the packing and reflected that the packing compactness increased as the number of rolling passes increased. For a certain roundness filler, the increase in the coordination number will lead to the enhancement of interlock between the particles, thus improving the antirotation ability of the particles. In addition, large roundness values of particles represent a large coordination number of the particles, indicating that the rounder particles are more likely to contact multiple particles.Figure 18
Mean particle coordination number-compactor pass responses for fillings with different roundnesses.Figure19 shows the relationship between average contact coordination number and rolling pass number of fillers with different roundnesses. Similar to the variation trend of the particle coordination number, the contact coordination number of a certainly shaped particle increased only slightly on the whole. However, for the particles with different roundnesses, smaller roundness corresponds to a larger contact coordination number of the particles, indicating that the interlock phenomenon is more obvious.Figure 19
Mean contact coordination number-compactor pass responses for fillings with different roundnesses.


#### 4.3.2. Particle Rotation Momentum

Rotation is an important relative motion between particles in contact. The rotation between particles is related to the rotation behavior of individual particles along the reference axis. Figures20(a) and 20(b) show the cumulative rotational momentum of fillers with different roundnesses along with clockwise and counterclockwise directions as a result of rolling. It can be seen from Figure 20 that the rotation momentum of the packing in the clockwise and counterclockwise directions is basically the same. Under a certain roundness, the rotation value of the particles increases with the increase in the number of rolling times, and the increase tends to slow down. This is because the increase in the particle coordination number and contact coordination number enhances the interlock between particles and hinders the rotation of particles. In addition, the smaller roundness of the particle represents the smaller rotational momentum of the particle, indicating that the more irregular the particle is, the stronger the antirotation ability is. The presumed reason is that the particles with lower roundness have a larger contact coordination number, which increases the interlock effect between the particles and thus hinders the rotation of the particles.Figure 20
Particle rotation magnitude-compactor pass responses for fillings with different roundnesses. (a) Clockwise momentum; (b) rotational momentum in reverse time.
(a)(b)Figure21(b) shows the vertical displacement of the packing. It showed a similar trend to the cumulative settlement; that is, with the increase in the number of rolling times, the particles gradually sank and the downward trend gradually slowed down. In addition, for fillers with different roundnesses, the smaller the particle roundness is, the greater the vertical displacement is. Estrada et al. [29] pointed out that when the particles in contact have relative movement, irregular particles with stronger antirotation ability will be forced to slide along the contact surface, so as to regulate the overall deformation of the particle system. Therefore, more irregular particles are more prone to relative sliding between particles, resulting in larger vertical displacement of particles.Figure 21
Displacement components of fillings in the vibration rolling process: vertical displacement.


## 4.1. Simulation of Earth Compaction of Rockfill Subgrade

Affected by the actual construction scale of the subgrade project, particle gradation and the irregular shape particles are adopted in the simulation, and the adoption of the 3D model will greatly increase the workload of the computer. Therefore, the two-dimensional model can also be used to simulate and analyze the particle motion characteristics of rockfill in the compaction process from a microscopic perspective.In this study, discrete element software PFC2D was used to simulate the compaction process of rockfill in two dimensions. To ensure that the roller vibrator rolls rather than sliding during the movement, the relationship between the angular velocityw and the horizontal velocity v of the vibrator is shown in the following equation:(6)v=w×R,where, v,w, and R are the horizontal velocity, angular velocity, and radius of the vibration wheel, respectively.The simulated material is rockfill for subgrade construction, and the linear elastic contact model is adopted in the discrete element model, which is widely used in similar studies [24]. The mesoscopic parameters used in the simulation refer to the research results of Liu et al. [24], and the specific parameters are shown in Table 2.Table 2
Microscale parameters used in the DEM simulations.
Microscopic parametersThe normal stiffness,kn (N/m)The tangential stiffness,ks (N/m)The coefficient of friction, fricLocal dampNormal/tangential dampingdpn, dptDensity (kg/m3)Value (particle)1.5 × 1072.9 × 1070.50.90.52600Value (wheel)1.0 × 1091.0 × 1090.30.00.57800Referring to the research of Goldenberg and Goldhirsch [25], Ma et al. [11], Zhou et al. [9], Alaei and Mahboubi [26] and Nie et al. [12], in this study, 5 groups of two-dimensional particles with different edge angles were generated, with 8 samples in each group, as shown in Figure 13. To exclude the influence of elongation on the simulation results, the elongation of particles was controlled between 0.97 and 1.00. Note that particles of roundness 1 are standard circles and are not shown in Figure 13. In addition, the effect of particle roughness was not considered in this study, that is, D8 = 0.0.Figure 13
Two-dimensional particle models with different roundnesses.In PFC2D, the clump function can realize the fitting of complex particle contour, with the fitting accuracy control parameters distance (range 0–180) and ratio (range 0.0–1.0). The larger the distance, the smaller the ratio, and the fitted results are closer to the real model. However, more computational memory and running time are consumed. In this study, the distance was set as 162, and the ratio was set as 0.25. Figure14 shows a clump model filled and generated under this setting. It can be seen that the set of two parameters is sufficient to reflect the external contour of the particle.Figure 14
A DEM-simulated virtual particle from a Fourier-based particle contour.Coarse particle packing with a diameter greater than 10 mm was used in this simulation to improve the calculation efficiency, and the distribution of the size of the packing was calculated. The results show that the content percentage of diameter ranging from 10 to 20 mm, 20 to 30 mm, 30 to 40 mm, 40 to 50 mm, and 50 to 60 mm is 17%, 20%, 25%, 20%, and 18%, respectively.The packing model was set in the range of height (H) × width (W) = 0.8 m × 8 m. The initial model is generated by gravity sedimentation. The process is as follows:(1)
In the 2H × W range, contactless particles slightly larger than the total volume of the target were randomly generated according to the particle size range and content, as shown in Figure 15(a).(2)
Under the action of gravity, particles fall freely, when the particle system reaches a quasi-static state and stops running, as shown in Figure15(b).(3)
In order to ensure that fillers with different roundnesses are compacted at the same thickness, redundant particles above the thickness line are removed to simulate the leveling process. The final leveling filler is shown in Figure15(c).(4)
A rigid body circular element is generated to simulate the roller. The movement parameters and attribute parameters of the roller are determined in Tables2 and 3. The chord wave vibration load is applied to the filler during the movement of the roller, as shown in Figure 15(d). And the simulation is stopped after repeated rolling 12 times.Figure 15
PFC simulation of vibratory compaction, including dumping, spreading, and compaction. (a) The contactless model; (b) sedimentation accumulation; (c) leveling; (d) vibratory compaction.
(a)(b)(c)(d)Table 3
External load parameters.
Rolling wheel parametersQuality of grinding wheel,G (kg)Grinding wheel diameter,D (mm)Grinding wheel width,b (mm)Vibration frequency,f (Hz)Exciting force,F0 (kN)Speed,v (m/s)Rolling passThe actual parameter1 × 1041500200040300112Simulation parameters5 × 1031500100040150112


## 4.2. The Macroresult

Figure16(a) shows the porosity of fillers with different roundnesses varying with the number of rolls. Obviously, with the increase in the number of rolling times, the pores first decreased sharply and then remained unchanged. But the correlation between roundness and vibratory compaction was not strong. The initial porosity was extracted (shown in Figure 16(b)), and it was found that with the increase of roundness, the initial porosity first decreased and then increased, which was consistent with the previous research results related to particle roundness [27]. Figure 16(c) shows the normalized porosity rate ∆n obtained after conversion using equation (7). It is found that the small roundness of the particle indicates irregular shape and great drop in the particle porosity compared to the initial state. The normalized porosity rate is obtained as(7)Δn=n0−nin0×100%,where n0 and ni represent the porosity corresponding to the initial porosity and the porosity corresponding to the ith pass after rolling, respectively.Figure 16
The relationship between the compactor pass and pore change of fillings with different roundnesses. (a) Porosity; (b) initial porosity; (c) normalized void change rate.
(a)(b)(c)Figure17 shows the cumulative settlement amount of the pavement corresponding to each roll pass of the filler, which shows a trend similar to the normalized pore change rate. It shows a trend of rapid settlement in the early stage and a trend of slow settlement in the later stage. This trend is consistent with the actual change in the monitored value in construction [24]. The indirect results showed that the particles sank under the vibratory compaction condition, which made the fillers squeeze and approach each other and rearrange the particles. This results in a smaller porosity, showing the changing trend as shown in Figure 16(a). In addition, a small value of roundness of the particles represents the great settlement of the fillers. Therefore, in practical engineering, the irregular shape of the packing requires more rolling times.Figure 17
Cumulative settlement-compactor pass responses for fillings with different roundnesses.


## 4.3. The Micromechanical Results

### 4.3.1. Coordination Number

Coordination number is an important index to characterize the internal structure of particle systems. For the irregular particle model composed of Clump, there are two different definitions of coordination number: one is the number of particles in contact with the specified particles, called the coordination number of particles; and the other is the number of contacts belonging to a given particle, called the contact coordination number.Thornton [28] pointed out that there are particles in granular materials that do not contact with other particles or only contact with one particle, and such particles have no contribution to the overall stability of the material. Therefore, the mechanical average particle coordination number is defined as follows:(8)CN=2C−N1N−N0−N1,where C is the total number of particle contact, N is the total number of particles, and N0 and N1 are the numbers of particles in contact with zero or one particle in the particle system. Figure 18 shows the relationship between the mechanical particle coordination number and the rolling pass number of fillers with different roundnesses. For each roundness case, the particle coordination number showed a slight upward trend as the number of rolling passes increased; that is, the number of particles in contact with a certain particle increased. This was related to the overall looseness of the packing and reflected that the packing compactness increased as the number of rolling passes increased. For a certain roundness filler, the increase in the coordination number will lead to the enhancement of interlock between the particles, thus improving the antirotation ability of the particles. In addition, large roundness values of particles represent a large coordination number of the particles, indicating that the rounder particles are more likely to contact multiple particles.Figure 18
Mean particle coordination number-compactor pass responses for fillings with different roundnesses.Figure19 shows the relationship between average contact coordination number and rolling pass number of fillers with different roundnesses. Similar to the variation trend of the particle coordination number, the contact coordination number of a certainly shaped particle increased only slightly on the whole. However, for the particles with different roundnesses, smaller roundness corresponds to a larger contact coordination number of the particles, indicating that the interlock phenomenon is more obvious.Figure 19
Mean contact coordination number-compactor pass responses for fillings with different roundnesses.


### 4.3.2. Particle Rotation Momentum

Rotation is an important relative motion between particles in contact. The rotation between particles is related to the rotation behavior of individual particles along the reference axis. Figures20(a) and 20(b) show the cumulative rotational momentum of fillers with different roundnesses along with clockwise and counterclockwise directions as a result of rolling. It can be seen from Figure 20 that the rotation momentum of the packing in the clockwise and counterclockwise directions is basically the same. Under a certain roundness, the rotation value of the particles increases with the increase in the number of rolling times, and the increase tends to slow down. This is because the increase in the particle coordination number and contact coordination number enhances the interlock between particles and hinders the rotation of particles. In addition, the smaller roundness of the particle represents the smaller rotational momentum of the particle, indicating that the more irregular the particle is, the stronger the antirotation ability is. The presumed reason is that the particles with lower roundness have a larger contact coordination number, which increases the interlock effect between the particles and thus hinders the rotation of the particles.Figure 20
Particle rotation magnitude-compactor pass responses for fillings with different roundnesses. (a) Clockwise momentum; (b) rotational momentum in reverse time.
(a)(b)Figure21(b) shows the vertical displacement of the packing. It showed a similar trend to the cumulative settlement; that is, with the increase in the number of rolling times, the particles gradually sank and the downward trend gradually slowed down. In addition, for fillers with different roundnesses, the smaller the particle roundness is, the greater the vertical displacement is. Estrada et al. [29] pointed out that when the particles in contact have relative movement, irregular particles with stronger antirotation ability will be forced to slide along the contact surface, so as to regulate the overall deformation of the particle system. Therefore, more irregular particles are more prone to relative sliding between particles, resulting in larger vertical displacement of particles.Figure 21
Displacement components of fillings in the vibration rolling process: vertical displacement.


## 4.3.1. Coordination Number

Coordination number is an important index to characterize the internal structure of particle systems. For the irregular particle model composed of Clump, there are two different definitions of coordination number: one is the number of particles in contact with the specified particles, called the coordination number of particles; and the other is the number of contacts belonging to a given particle, called the contact coordination number.Thornton [28] pointed out that there are particles in granular materials that do not contact with other particles or only contact with one particle, and such particles have no contribution to the overall stability of the material. Therefore, the mechanical average particle coordination number is defined as follows:(8)CN=2C−N1N−N0−N1,where C is the total number of particle contact, N is the total number of particles, and N0 and N1 are the numbers of particles in contact with zero or one particle in the particle system. Figure 18 shows the relationship between the mechanical particle coordination number and the rolling pass number of fillers with different roundnesses. For each roundness case, the particle coordination number showed a slight upward trend as the number of rolling passes increased; that is, the number of particles in contact with a certain particle increased. This was related to the overall looseness of the packing and reflected that the packing compactness increased as the number of rolling passes increased. For a certain roundness filler, the increase in the coordination number will lead to the enhancement of interlock between the particles, thus improving the antirotation ability of the particles. In addition, large roundness values of particles represent a large coordination number of the particles, indicating that the rounder particles are more likely to contact multiple particles.Figure 18
Mean particle coordination number-compactor pass responses for fillings with different roundnesses.Figure19 shows the relationship between average contact coordination number and rolling pass number of fillers with different roundnesses. Similar to the variation trend of the particle coordination number, the contact coordination number of a certainly shaped particle increased only slightly on the whole. However, for the particles with different roundnesses, smaller roundness corresponds to a larger contact coordination number of the particles, indicating that the interlock phenomenon is more obvious.Figure 19
Mean contact coordination number-compactor pass responses for fillings with different roundnesses.


## 4.3.2. Particle Rotation Momentum

Rotation is an important relative motion between particles in contact. The rotation between particles is related to the rotation behavior of individual particles along the reference axis. Figures20(a) and 20(b) show the cumulative rotational momentum of fillers with different roundnesses along with clockwise and counterclockwise directions as a result of rolling. It can be seen from Figure 20 that the rotation momentum of the packing in the clockwise and counterclockwise directions is basically the same. Under a certain roundness, the rotation value of the particles increases with the increase in the number of rolling times, and the increase tends to slow down. This is because the increase in the particle coordination number and contact coordination number enhances the interlock between particles and hinders the rotation of particles. In addition, the smaller roundness of the particle represents the smaller rotational momentum of the particle, indicating that the more irregular the particle is, the stronger the antirotation ability is. The presumed reason is that the particles with lower roundness have a larger contact coordination number, which increases the interlock effect between the particles and thus hinders the rotation of the particles.Figure 20
Particle rotation magnitude-compactor pass responses for fillings with different roundnesses. (a) Clockwise momentum; (b) rotational momentum in reverse time.
(a)(b)Figure21(b) shows the vertical displacement of the packing. It showed a similar trend to the cumulative settlement; that is, with the increase in the number of rolling times, the particles gradually sank and the downward trend gradually slowed down. In addition, for fillers with different roundnesses, the smaller the particle roundness is, the greater the vertical displacement is. Estrada et al. [29] pointed out that when the particles in contact have relative movement, irregular particles with stronger antirotation ability will be forced to slide along the contact surface, so as to regulate the overall deformation of the particle system. Therefore, more irregular particles are more prone to relative sliding between particles, resulting in larger vertical displacement of particles.Figure 21
Displacement components of fillings in the vibration rolling process: vertical displacement.


## 5. Conclusion

Rockfill material is of great significance in the subgrade construction industry. The mechanical properties can be affected by the irregularity and variability of particle shape， which can be explored by discrete element simulation. In this study, the subgrade rockfill particle outlines are firstly extracted based on machine learning algorithms, including AdaBoost, Cascade, and sliding windows. Then, the shape evaluation indexes are quantified and the statistical analysis of each index is presented. In addition, the DEM simulation is implemented on the compaction of rockfill subgrade. Based on the simulation, the macroanalysis on the void ratio and cumulative settlement is given. Finally, the microanalysis on the particle coordination number, rotation momentum, and displacement is studied. The specific results are as follows:(1)
1206 pictures of the subgrade rockfill particles are taken as the training data. By training a strong classifier through AdaBoost in each stage of Cascade, only a small amount of pictures can pass alls stages in the Cascade algorithm. These pictures contain a complete outline of particles. Then, using sliding window and adding boundaries of particles, the final window of a particle can be obtained indicating that the particle is well extracted.(2)
The quantitative shape analysis is presented. For length flatness, 991 real particles show an approximately normal distribution, whereas the edge angle and roughness show a negative skewness distribution and positive skewness distribution, respectively.(3)
With the increase in the number of rolling passes, the porosity of the packing decreases gradually and the settlement increases gradually, and the change rate is fast at the beginning and slow in the end. In addition, smaller roundness of the particles represents a more obvious porosity drop and the greater settlement of the filler compared with the initial state.(4)
Under the rolling load, the particle coordination number and contact coordination number of the filler increased slightly. The smaller the roundness is, the larger the contact coordination number is, which increases the interlock effect and hinders the rotation of the particles. The horizontal displacement of particles shows a certain fluctuation, whereas the vertical displacement increases gradually. In addition, the smaller the roundness is, the larger the vertical displacement is and the smaller the horizontal displacement is.



---

*Source: 5043729-2022-03-04.xml*