*Article* **Algorithm for Virtual Aggregates' Reconstitution Based on Image Processing and Discrete-Element Modeling**

#### **Danhua Wang <sup>1</sup> , Xunhao Ding <sup>2</sup> , Tao Ma 2,\*, Weiguang Zhang <sup>2</sup> and Deyu Zhang <sup>3</sup>**


Received: 15 February 2018; Accepted: 3 May 2018; Published: 7 May 2018

**Abstract:** Based on the Aggregate Imaging Measurement System (AIMS) and the Particle Flow Code in Two Dimensions (PFC2D), an algorithm for modeling two-dimensional virtual aggregates was proposed in this study. To develop the virtual particles precisely, the realistic shapes of the aggregates were captured by the AIMS firstly. The shape images were then processed, and the morphological characteristics of aggregates were quantified by the angularity index. By dividing the particle irregular shape into many triangle areas and adjusting the positions of the generated balls via coordinate systems' conversion within PFC2D, the virtual particles could be reconstructed accurately. By calculating the mapping area, the gradations in two-dimensions could be determined. Controlled by two variables (µ\_1 and µ\_2), which were drawn from the uniform distribution (0, 1), the virtual particles forming the specimens could be developed with random sizes and angular shapes. In the end, the rebuilt model of the SMA-13 aggregate skeleton was verified by the virtual penetration tests. The results indicated that the proposed algorithm can not only model the realistic particle shape and gradations precisely, but also predict its mechanical behavior well.

**Keywords:** algorithm; aggregates; reconstitution; image processing; virtual modelling

### **1. Introduction**

As the basis materials within the asphalt mixture, coarse aggregates are widely used in asphalt pavements, which form the main skeleton of the total structure. It has been highlighted by many researchers that the aggregates' shape has a significant impact on controlling the mechanical properties of the asphalt mixtures [1–3]. Tests conducted by Campen and Smith [4] showed that the stability of hot-mix asphalt (HMA) concrete could be improved obviously by using crushed aggregates instead of rounded particles. Similar results were shown by Kandhal and Parker [5], that the HMA mixtures consisting of excessive flat and elongated aggregate particles tended to break down, affecting the mix durability. The fatigue performance of the asphalt mixture was also affected by the particles' roundness and angularity, concluded based on Cheung and Dawson's research [6]. In their studies, with the angularity increasing, the ultimate shear strength increased, and permanent deformation decreased, respectively. Thus, it is strongly recommended to use the complex angular particles to improve the grains' interlocking and packing effectively.

However, it is really difficult and time-consuming to reveal the shape influences of coarse aggregates just by laboratory tests due to the appropriate definition of the morphological characteristics and the following shape acquisition, which is the basis to guarantee the results [7,8]. The discrete-element method (DEM) was proposed by the Cundall [9] firstly, and then, a software program named Particle Flow Code (PFC) was developed to help the mechanical calculation of heterogeneous materials, which has been widely used in many fields [10–12]. As the basic calculation elements, the balls in the software (PFC2D, Itasca Consulting Group: Minneapolis, USA) are all round shaped without any angular shape, which is not consistent with the realistic particles in pavement construction. Therefore, prior to numerical simulation, it is important to rebuild the realistic particle shape, which will influence the simulation results obviously [13–15]. Many modeling methods have been proposed within PFC to characterize the shape features by some researchers. Zhang et al. [16] regarded the aggregate particles as hexahedrons, pentahedrons and tetrahedrons within PFC, respectively. Then, the equal-radius balls were arranged regularly within the particles and were defined as a rigid body. CPA (Computerized Particle Analysis) was utilized to help develop the virtual shape within PFC by Stahl and Konietzky [17]. The input index, known as the grain-shape coefficient, was used to redefine the particle contour. Then, numerical particles were reconstructed by filling the redefined area with balls. Huang [18] reconstructed 3D ballast particles with three images (top, side and front) of the particles with the help of the image processing technology. Similar methods can also be found in the works of Liu and You [19,20], Tutumluer et al. [21] and Dondi and Vignali [22–24], as well. Most of the existing modeling methods can only model the shape roughly with excessive simplification. Although some virtual models were rebuilt precisely by X-ray scanning, the angular shape could not be quantified precisely [25–29]. Liu scanned the particles and rebuilt their virtual shapes according to the acquired images [30,31]. In their studies, the stl-format file were developed, which can reflect the shapes of aggregates precisely. By the designed routine, the virtual particles in three dimensions were developed successfully with enough balls filled. However, their methods can only model the realistic shapes, but cannot evaluate the shape properties of mixtures by a specific index. Therefore, the simulation accuracy mainly depends on the samples' selection and cannot model the material heterogeneity well. Compared to X-ray scanning, the Aggregate Imaging Measurement System (AIMS) apparatus is a promising way for performing particle shape modeling and has been utilized gradually in material analyses [32]. It can capture the morphological characteristics of aggregates and quantifies the shape properties by some design indices such as the angularity index, texture index, size, and so on. However, this apparatus is mainly used for pavement skid resistance measurements nowadays [33–35]. More uses of it should be explored, especially combined with other technologies.

The objective of this study is to propose an algorithm for the virtual aggregates' reconstitution based on the Discrete-Element Method (DEM). To achieve this, the AIMS apparatus with the image processing technology were utilized before the simulations. Then, an innovative method for developing the virtual particles and gradations was proposed precisely. In the end, based on the virtual penetration tests, the precision of the developed models was verified. The relevant issues include: capturing the shape image of particles by AIMS; processing the captured images and measuring their angularity index; modeling the virtual particles with different sizes and shapes; modeling the virtual gradations of specimens; developing the virtual penetration tests to predict and evaluate the performance of coarse aggregates.

### **2. Image Preprocessing and Shape Measurement**

Prior to the virtual particles' reconstitution, several steps were performed that could provide the proposed algorithm with the fundamental data. Firstly, enough aggregate particles were selected and were scanned by the AIMS to obtain their shape image and the angularity index. The particles were put in the turntable of the AIMS and were rotated slowly to go through the scanning area. A digital camera above the particles was utilized to capture the image of the maximum section from the top side. Then, the original images of the particles were converted to the binary ones within the MATLAB routine as shown in Figure 1. As shown, two forms of binary images were acquired via the image processing technology including the shape (Figure 1a) and its contour (Figure 1b).

**Figure 1.** Binary images of aggregate particles: (**a**) shape; (**b**) shape contour. **Figure 1.** Binary images of aggregate particles: (**a**) shape; (**b**) shape contour.

Due to the size differences between the binary images and realistic particles, the processed images should be scaled back based on the actual ones. Because the binary images consisted of pixels only, the sizes were represented by the pixel number at different directions in the images rather than the actual length. In this proposed method, the relationship between the actual length and "pixel length" were quantified by uniform standards as one pixel equal to 0.005 mm. After this standard was developed, based on the actual particle size, the required reduction scale of binary images could be calculated to match the actual size, as shown in Table 1. Furthermore, after the images' adjustment completed, the actual area of various particles was measured by the ImagePro software. The black area of the shape image in Figure 1a could be distinguished within the software, and the corresponding area was calculated and stored as shown in Table 1, as well. Due to the size differences between the binary images and realistic particles, the processed images should be scaled back based on the actual ones. Because the binary images consisted of pixels only, the sizes were represented by the pixel number at different directions in the images rather than the actual length. In this proposed method, the relationship between the actual length and "pixel length" were quantified by uniform standards as one pixel equal to 0.005 mm. After this standard was developed, based on the actual particle size, the required reduction scale of binary images could be calculated to match the actual size, as shown in Table 1. Furthermore, after the images' adjustment completed, the actual area of various particles was measured by the ImagePro software. The black area of the shape image in Figure 1a could be distinguished within the software, and the corresponding area was calculated and stored as shown in Table 1, as well.


**Table 1.** Summaries of the two-dimensional images. **Table 1.** Summaries of the two-dimensional images.

Although there was an original coordinate system after image processing, another two coordinate systems were introduced in this proposed method to help the virtual particles' reconstitution further, including the local coordinate system and global coordinate system. The local one will be introduced in the algorithm part, while the global one is developed as the following. It is necessary to develop a uniform global coordinate system for all particles because the original coordinate system in each image was not the same at all. The central position of each particle was defined firstly as shown in Equation (1): Although there was an original coordinate system after image processing, another two coordinate systems were introduced in this proposed method to help the virtual particles' reconstitution further, including the local coordinate system and global coordinate system. The local one will be introduced in the algorithm part, while the global one is developed as the following. It is necessary to develop a uniform global coordinate system for all particles because the original coordinate system in each image was not the same at all. The central position of each particle was defined firstly as shown in Equation (1):

$$\text{cen\\_x} = \frac{\sum\_{1}^{n} \text{x}\_{i}}{n} \text{ cen\\_y} = \frac{\sum\_{1}^{n} y\_i}{n} \tag{1}$$

where \_, \_ are the X- and Y-coordinate values of the particle central position, respectively, in the original coordinate system; *n* is the total number of pixels in Figure 1b; , are the X- and Ywhere *cen*\_*x*, *cen*\_*y* are the X- and Y-coordinate values of the particle central position, respectively, in the original coordinate system; *n* is the total number of pixels in Figure 1b; *x<sup>i</sup>* , *y<sup>i</sup>* are the X- and Y-coordinate values of the *i*-th pixel in Figure 1b in the original coordinate system.

coordinate values of the *i*-th pixel in Figure 1b in the original coordinate system. The original coordinate system of each image was converted to the uniform global coordinate system by setting the central positions as the new origin of the coordinates, and the adjustments of the original coordinate are summarized as shown in Table 2. It is noted that the data in Tables 1 and 2 are only a part of the total virtual particles limited by space. The data shown were utilized to illustrate the preparation procedure. The original coordinate system of each image was converted to the uniform global coordinate system by setting the central positions as the new origin of the coordinates, and the adjustments of the original coordinate are summarized as shown in Table 2. It is noted that the data in Tables 1 and 2 are only a part of the total virtual particles limited by space. The data shown were utilized to illustrate the preparation procedure.


**Table 2.** Adjustments of the image coordinates. **Table 2.** Adjustments of the image coordinates. **Sample Adjustments of X-Coordinate Adjustments of Y-Coordinate**

#### **3. Image Preprocessing and Shape Measurement 3. Image Preprocessing and Shape Measurement**

#### *3.1. Triangle Area Divisions 3.1. Triangle Area Divisions*

As the independent element for calculation, the ball within PFC2D should be generated with a specific coordinate value. However, it is very hard to fill the particle shape area with balls directly because of the shape complexity, which prevents precise virtual particle modelling. The proposed algorithm converted the irregular particle shape into several triangles to make it easier to fill the required area with balls. The algorithm is illustrated as follows: As the independent element for calculation, the ball within PFC2D should be generated with a specific coordinate value. However, it is very hard to fill the particle shape area with balls directly because of the shape complexity, which prevents precise virtual particle modelling. The proposed algorithm converted the irregular particle shape into several triangles to make it easier to fill the required area with balls. The algorithm is illustrated as follows:

After the shape contour (Figure 1b) was scaled back, its irregular shape was divided into 36 triangle areas firstly, as shown in Figure 2. Measuring lines were drawn every 10 degrees by a user-defined routine based on MATLAB. They were drawn from the central position to the shape contour. Each triangle could represent a part of the shape area approximately. Then, the lengths of 36 lines were measured respectively by the routine. It is noted that the measuring lines in Figure 2 are not true lines, but a routine developed within MATLAB to measure the distance from the central position to the contour along the designed directions. The modelling precision and efficiency are definitely influenced by the number of divided triangles. As shown in Figure 2, with more triangles divided, the boundary of the virtual particle can have better accordance with the realistic one, but needs more complex routines for length measurements, which means more time consumption. To balance the simulation precision and efficiency, 36 triangle divisions were determined finally and were thought as the best after many tries. After the shape contour (Figure 1b) was scaled back, its irregular shape was divided into 36 triangle areas firstly, as shown in Figure 2. Measuring lines were drawn every 10 degrees by a user-defined routine based on MATLAB. They were drawn from the central position to the shape contour. Each triangle could represent a part of the shape area approximately. Then, the lengths of 36 lines were measured respectively by the routine. It is noted that the measuring lines in Figure 2 are not true lines, but a routine developed within MATLAB to measure the distance from the central position to the contour along the designed directions. The modelling precision and efficiency are definitely influenced by the number of divided triangles. As shown in Figure 2, with more triangles divided, the boundary of the virtual particle can have better accordance with the realistic one, but needs more complex routines for length measurements, which means more time consumption. To balance the simulation precision and efficiency, 36 triangle divisions were determined finally and were thought as the best after many tries.

**Figure Figure 2. 2.** Triangle area divisions for measurement Triangle area divisions for measurement..

#### *3.2. Filling Area Judgments 3.2. Filling Area Judgments*

By labelling the lines with numbers increasing counter-clockwise as shown in Figure 2, the divided triangles could fall into two categories. One is that the top side of the triangle is longer than the bottom side; another is that the top one is shorter than the bottom one, as shown in Figure 3. The top side is always labelled with a larger number compared to the bottom side in one triangle. *S*<sup>1</sup> and *S*<sup>2</sup> are known quantities based on the former measurement in Figure 3, and "*a*" represents the included angle between two adjacent measuring lines with a magnitude of 10 degrees. *d*1, *d*2, *d*3, *d*<sup>4</sup> and *d*<sup>5</sup> were calculated as shown in Equation (2). By labelling the lines with numbers increasing counter-clockwise as shown in Figure 2, the divided triangles could fall into two categories. One is that the top side of the triangle is longer than the bottom side; another is that the top one is shorter than the bottom one, as shown in Figure 3. The top side is always labelled with a larger number compared to the bottom side in one triangle. <sup>1</sup> and <sup>2</sup> are known quantities based on the former measurement in Figure 3, and "*a*" represents the included angle between two adjacent measuring lines with a magnitude of 10 degrees. <sup>1</sup> , <sup>2</sup> , <sup>3</sup> , <sup>4</sup> and <sup>5</sup> were calculated as shown in Equation (2). <sup>1</sup> = <sup>2</sup> ×

$$\begin{aligned} d\_1 &= S\_2 \times \cos a \\ d\_2 &= S\_1 - d\_1 \\ d\_3 &= S\_2 \times \sin a \\ d\_4 &= S\_2 \times \cos a \\ d\_5 &= d\_4 - S\_1 \end{aligned}$$

(2)

(2)

**Figure 3.** Algorithm for filling the triangle area with discrete balls: (**a**) <sup>2</sup> ≤ 1; (**b**) <sup>2</sup> > 1. **Figure 3.** Algorithm for filling the triangle area with discrete balls: (**a**) *S*<sup>2</sup> ≤ *S*<sup>1</sup> ; (**b**) *S*<sup>2</sup> > *S*<sup>1</sup> .

(Figure 3a), the algorithm for filling the triangle area with balls by PFC2D was

conducted as following: ①Start the ball generation from Point A with the coordinate of (0, 0); When *S*<sup>2</sup> ≤ *S*<sup>1</sup> (Figure 3a), the algorithm for filling the triangle area with balls by PFC2D was conducted as following:

 ②Determine the variable <sup>1</sup> . It is calculated as the following equation: • <sup>1</sup> Start the ball generation from Point A with the coordinate of (0, 0);

When <sup>2</sup> ≤ <sup>1</sup>

• <sup>2</sup> Determine the variable *x*1. It is calculated as the following equation:

$$\mathbf{x}\_1 = \mathbf{x}\_{1f} + \mathbf{2} \times \mathbf{r}ad \tag{3}$$

where *x*<sup>1</sup> is the X-coordinate of the next column of balls that will be generated; *x*<sup>1</sup> *<sup>f</sup>* is the X-coordinate of the former column of balls that has been generated; *rad* is the radius of the filled balls.

• <sup>3</sup> Determine the variable *h*, calculated as shown in Equation (4):

$$h = \begin{cases} \begin{array}{ll} \mathfrak{x}\_1 \times \text{tan}\mathfrak{a} & \mathfrak{x}\_1 \le d\_1 \\ \frac{S\_1 - \mathfrak{x}\_1}{d\_2} \times d\_3 & d\_1 < \mathfrak{x}\_1 \le \mathfrak{S}\_1 \end{array} \tag{4}$$

where *h* is the required height of the next column of balls; the others are the same as the former.


When *S*<sup>2</sup> > *S*<sup>1</sup> (Figure 3b), the algorithm is conducted as follows:


$$h\_1 = \mathbf{x}\_1 \times \mathbf{t} \mathbf{n} a \qquad \mathbf{x}\_1 \le \mathbf{S}\_1 \tag{5}$$

where *h*<sup>1</sup> is the required height of the next column of balls; the others are the same as the former.

$$h\_i = \begin{cases} \frac{X\_1 \times \tan a}{\frac{X\_1 - S\_1}{d\_\mathcal{g}} \times S\_2 \times \sin a} & i = 1\\ \frac{X\_1 - S\_1}{d\_\mathcal{g}} \times S\_2 \times \sin a & i = 2 \end{cases} \tag{6}$$

where *h*<sup>1</sup> is the required height of the next column of balls; *h*<sup>2</sup> is the lower limit of the column height when *x*<sup>1</sup> exceeds *S*1; the others are the same as the former;


### *3.3. Coordinate Systems' Conversion*

By the proposed algorithm, the various triangle areas within the particle images could be filled with uniformly-arranged balls successfully. However, when combining all of the triangles to form the whole virtual particle, it is necessary to calibrate the coordinate value from the local coordinate system to the global one. The relationships between the two coordinate systems are shown in Figure 4. The Xand Y-coordinate of each filled ball should be adjusted based on Equation (7).

$$\begin{cases} \mathbf{x} = \mathbf{x}' \times \cos b + y' \times \sin b \\ y = -\mathbf{x}' \times \sin b + y' \times \cos b \end{cases} \tag{7}$$

where *x* 0 and *y* 0 are the X-, Y-coordinates in the local coordinate system of the triangles, respectively; *x* and *y* are the X-, Y-coordinates in the global coordinate system of the particles, respectively; *b* is the angle between the bottom side of the triangle (*S*1) and the *X*-axis.

*Appl. Sci.* **2018**, *2*, 15 FOR PEER REVIEW 7 of 16

**Figure 4.** Conversions between local and global coordinate systems. **Figure 4.** Conversions between local and global coordinate systems. **Figure 4.** Conversions between local and global coordinate systems.

It is noted that two different coordinate systems were used in the modeling process. The local

coordinate values represented the relative positions of balls in the virtual particles, while the global coordinate values represented relative positions of balls in the virtual container. When developing the single particle, the local coordinate systems were utilized to fill the irregular area with arranged balls. After that, the positions of all the virtual particles were converted to the global one so that the numerous virtual particles could be distributed within the designed area to form the whole mixtures. After the coordinates were calibrated, the virtual particles could be rebuilt as shown in Figure 5. As shown, the inner area of the virtual particles was filled with balls gradually after each triangle area was processed based on the proposed algorithm. However, although the irregular shapes of the particles were divided into several triangles, the side of <sup>3</sup> (Figure 3) could only match the actual shape contour similarly. The shape modeling was not precise enough. Therefore, further improvements were made by modeling the shape contour firstly and then filling the inner area of the shape contour with balls by the former algorithm, as shown in Figure 5. The shape contour in Figure 1b was processed by a user-defined routine to read and store the pixel locations of the images. The black pixels were distinguished. Then, a threshold value of 0.005 was selected to filter the redundant black pixels such that the adjacent black pixels within a distance of 0.005 could be ignored and the others retained. By importing the locations of the retained pixels into PFC2D, the balls could be generated in the shape contour, as shown in Figure 5. Then, the whole particle was filled with balls by triangle area filling in the end. Based on the scanning images and with the help of the proposed algorithm, aggregate particles varying in size and shape were reconstructed as shown in Figure 6. It is noted that two different coordinate systems were used in the modeling process. The local coordinate values represented the relative positions of balls in the virtual particles, while the global coordinate values represented relative positions of balls in the virtual container. When developing the single particle, the local coordinate systems were utilized to fill the irregular area with arranged balls. After that, the positions of all the virtual particles were converted to the global one so that the numerous virtual particles could be distributed within the designed area to form the whole mixtures. After the coordinates were calibrated, the virtual particles could be rebuilt as shown in Figure 5. As shown, the inner area of the virtual particles was filled with balls gradually after each triangle area was processed based on the proposed algorithm. However, although the irregular shapes of the particles were divided into several triangles, the side of *S*<sup>3</sup> (Figure 3) could only match the actual shape contour similarly. The shape modeling was not precise enough. Therefore, further improvements were made by modeling the shape contour firstly and then filling the inner area of the shape contour with balls by the former algorithm, as shown in Figure 5. The shape contour in Figure 1b was processed by a user-defined routine to read and store the pixel locations of the images. The black pixels were distinguished. Then, a threshold value of 0.005 was selected to filter the redundant black pixels such that the adjacent black pixels within a distance of 0.005 could be ignored and the others retained. By importing the locations of the retained pixels into PFC2D, the balls could be generated in the shape contour, as shown in Figure 5. Then, the whole particle was filled with balls by triangle area filling in the end. Based on the scanning images and with the help of the proposed algorithm, aggregate particles varying in size and shape were reconstructed as shown in Figure 6. It is noted that two different coordinate systems were used in the modeling process. The local coordinate values represented the relative positions of balls in the virtual particles, while the global coordinate values represented relative positions of balls in the virtual container. When developing the single particle, the local coordinate systems were utilized to fill the irregular area with arranged balls. After that, the positions of all the virtual particles were converted to the global one so that the numerous virtual particles could be distributed within the designed area to form the whole mixtures. After the coordinates were calibrated, the virtual particles could be rebuilt as shown in Figure 5. As shown, the inner area of the virtual particles was filled with balls gradually after each triangle area was processed based on the proposed algorithm. However, although the irregular shapes of the particles were divided into several triangles, the side of <sup>3</sup> (Figure 3) could only match the actual shape contour similarly. The shape modeling was not precise enough. Therefore, further improvements were made by modeling the shape contour firstly and then filling the inner area of the shape contour with balls by the former algorithm, as shown in Figure 5. The shape contour in Figure 1b was processed by a user-defined routine to read and store the pixel locations of the images. The black pixels were distinguished. Then, a threshold value of 0.005 was selected to filter the redundant black pixels such that the adjacent black pixels within a distance of 0.005 could be ignored and the others retained. By importing the locations of the retained pixels into PFC2D, the balls could be generated in the shape contour, as shown in Figure 5. Then, the whole particle was filled with balls by triangle area filling in the end. Based on the scanning images and with the help of the proposed algorithm, aggregate particles varying in size and shape were reconstructed as shown in Figure 6.

**Figure 5.** Procedure of the virtual particles' modeling. **Figure 5.** Procedure of the virtual particles' modeling.

*Appl. Sci.* **2018**, *2*, 15 FOR PEER REVIEW 8 of 16

**Figure 6.** Various virtual particles of different morphological characteristics. **Figure 6.** Various virtual particles of different morphological characteristics.

#### **4. Algorithm for Developing Virtual Specimens 4. Algorithm for Developing Virtual Specimens**

#### *4.1. Mapping Area Calculation for Two-Dimensional Specimens 4.1. Mapping Area Calculation for Two-Dimensional Specimens*

When it comes to the virtual aggregate specimens, several assumptions were made before the simulations as follows: When it comes to the virtual aggregate specimens, several assumptions were made before the simulations as follows:


The virtual specimen of SMA-13 in the penetration test was then developed, and the gradations of SMA-13 are shown in Table 3. To improve the simulation efficiency, only the coarse aggregates were developed despite the fine particles. The virtual specimen of SMA-13 in the penetration test was then developed, and the gradations of SMA-13 are shown in Table 3. To improve the simulation efficiency, only the coarse aggregates were developed despite the fine particles.


**Table 3.** Gradation for SMA-13. **Table 3.** Gradation for SMA-13.

The particle numbers of different grades in the two-dimensional simulation were determined based on the mapping area, as the following equations shown: The particle numbers of different grades in the two-dimensional simulation were determined based on the mapping area, as the following equations shown:

2

$$M = \rho \times \frac{\pi \times D^2}{4} \times h \tag{8}$$

g cm<sup>3</sup> ⁄ ; *D* is the diameter of the specimen, cm; *h* is the height of the specimen, cm. When it comes to the two-dimensional specimen, the mass can be converted as follows: where *M* is the total mass of the specimen in three dimensions, g; *ρ* is the density of the specimen, g/cm<sup>3</sup> ; *D* is the diameter of the specimen, cm; *h* is the height of the specimen, cm.

 = × = × × × ℎ = ℎ (9) When it comes to the two-dimensional specimen, the mass can be converted as follows:

others are the same as the former.

$$m = \frac{M}{V} \times S = \frac{V \times \rho}{V} \times D \times h = Dh\rho \tag{9}$$

are the same as the former. Thus, the required area for generating the total virtual particles could be calculated based on Equation (10). where *m* is the converted mass of the specimen in two dimensions, g; *V* is the total volume of the specimen in three dimensions, cm<sup>3</sup> ; *S* is the area of the specimen in two dimensions, cm<sup>2</sup> ; the others are the same as the former.

 = 1 = ℎ 1 (10) Thus, the required area for generating the total virtual particles could be calculated based on Equation (10).

$$S\_{total} = \frac{m}{\rho\_1} = \frac{Dh\rho}{\rho\_1} \tag{10}$$

where *Stotal* is the required area for generating the total virtual particles, cm<sup>2</sup> ; *ρ*<sup>1</sup> is the density of all the coarse aggregates, which was assumed to be the same value in the simulations, g/cm<sup>3</sup> ; the others are the same as the former.

Based on the assumption that there is no water within the specimen, the following relation can be given:

$$
\rho = \frac{\rho\_1 \times V \times (1 - K)}{V} = \rho\_1 \times (1 - K) \tag{11}
$$

where *K* is the void content of the total specimen in three dimensions, %; the others are the same as the former.

The mapping area of different grades of particles can be developed in the end as follows:

$$\begin{aligned} \mathbf{S}\_{\text{total}} &= \mathbf{D}h \times (1 - \mathbf{K}) \\ \mathbf{S}\_{i} = \mathbf{S}\_{\text{total}} \times \mathbf{P}\_{i} &= \mathbf{D}h \times (1 - \mathbf{K}) \times \mathbf{P}\_{i} \quad i = 4.75, \ 9.5, \ 13.2, \ 16 \end{aligned} \tag{12}$$

where *S*4.75, *S*9.5, *S*13.2, *S*<sup>16</sup> are the mapping areas of the 4.75-, 9.5-, 13.2-, 16-mm virtual particles, respectively, cm<sup>2</sup> ; *P*4.75, *P*9.5, *P*13.2, *P*<sup>16</sup> are the proportions of the 4.75-, 9.5-, 13.2-, 16-mm virtual particles to the total mass, respectively; the others are the same as the former.

### *4.2. Generation Process for Virtual Particles with Random Shapes and Sizes*

To distribute the virtual particles randomly within the specimen, two parameters (*µ*<sup>1</sup> and *µ*2) were given to control the generation process better. When generating the next particle, *µ*<sup>1</sup> represented the probability of the size, while *µ*<sup>2</sup> determined its angular shape. Two random values drawn from the uniform distribution (0, 1) were assigned to *µ*<sup>1</sup> and *µ*<sup>2</sup> at first. Based on *µ*<sup>1</sup> and *µ*<sup>2</sup> at the current calculation step, a specific size and shape could be determined. For example, before generating the next virtual particle, *µ*<sup>1</sup> would be checked within PFC2D as follows.


After the particle size was determined, *µ*2. would be checked to develop the various angular shapes of the particles based on the angularity index distribution, as shown in Figure 7. For example, for the 4.75-mm particles, the process was conducted as following.


*Appl. Sci.* **2018**, *2*, 15 FOR PEER REVIEW 10 of 16

**Figure 7.** Angularity index distribution of different sizes of particles. **Figure 7.** Angularity index distribution of different sizes of particles.

*D*, *h* and *K* were set as 15 cm, 12 cm and 35% during the modeling process, which was in accordance with the laboratory tests. The calculated mapping area is an important parameter to control the particles' generation. The virtual aggregates were generated three times when forming the numerical models, as shown in Figure 8. When simulated, only one third of the mapping area was taken to determine the required number of different grades of aggregates. The virtual particles Enough particles were scanned and rebuilt within PFC2D. Thus, when *µ*2was checked and the angularity index range was determined, corresponding virtual particles with the required shape would be generated by the user-defined routine. By the control of *µ*<sup>1</sup> and *µ*2, the virtual specimen could be developed with random particles obviously. 0.00% 2000 2500 3000 3500 4000 4500 5000 5500 **angularity index**

#### were generated one by one above the container and then fell down under gravity. Each time the generation completed, the specimen would be compacted by a virtual wall developed within PFC2D. *4.3. Example of Generating a Virtual Specimen* **Figure 7.** Angularity index distribution of different sizes of particles.

20.00%

**accumulative proportion**

*4.3. Example of Generating a Virtual Specimen*

The compacted wall would move downwards until reaching the designed height so as to achieving the required void content. By three cyclical steps, the virtual specimen could be developed as shown in Figure 8. *D*, *h* and *K* were set as 15 cm, 12 cm and 35% during the modeling process, which was in accordance with the laboratory tests. The calculated mapping area is an important parameter to control the particles' generation. The virtual aggregates were generated three times when forming the numerical models, as shown in Figure 8. When simulated, only one third of the mapping area was taken to determine the required number of different grades of aggregates. The virtual particles were generated one by one above the container and then fell down under gravity. Each time the generation completed, the specimen would be compacted by a virtual wall developed within PFC2D. The compacted wall would move downwards until reaching the designed height so as to achieving the required void content. By three cyclical steps, the virtual specimen could be developed as shown in Figure 8. *4.3. Example of Generating a Virtual Specimen D*, *h* and *K* were set as 15 cm, 12 cm and 35% during the modeling process, which was in accordance with the laboratory tests. The calculated mapping area is an important parameter to control the particles' generation. The virtual aggregates were generated three times when forming the numerical models, as shown in Figure 8. When simulated, only one third of the mapping area was taken to determine the required number of different grades of aggregates. The virtual particles were generated one by one above the container and then fell down under gravity. Each time the generation completed, the specimen would be compacted by a virtual wall developed within PFC2D. The compacted wall would move downwards until reaching the designed height so as to achieving the required void content. By three cyclical steps, the virtual specimen could be developed as shown in Figure 8.

**Figure 8.** *Cont.*

*Appl. Sci.* **2018**, *2*, 15 FOR PEER REVIEW 11 of 16

**Figure 8.** Modeling procedure of the virtual specimens: (**a**) first step of generation; (**b**) second step of generation; (**c**) third step of generation. **Figure 8.** Modeling procedure of the virtual specimens: (**a**) first step of generation; (**b**) second step of generation; (**c**) third step of generation.

The generation routine of each grade of aggregate would be carried out until the total area of the

generated particles exceeded the requirements. After the modeling completed, the total area of generated particles with different sizes was measured within PFC2D, as shown in Table 4. As shown, the generated virtual specimen was in good accordance with the actual gradations. The error was 5.51%, 3.86% and 3.70% for the 13.2-, 9.5- and 4.75-mm coarse aggregates, respectively. The error was attributed to the last generated particles, which exceeded the required mapping area, and it can be ignored. Thus, it is concluded that the proposed algorithm can well characterize the composition features of the aggregate specimen with precise particle shape reconstitutions and consistent gradation modeling. **Table 4.** Summaries of the virtual modeling. Sieving size (mm) 13.2 9.5 4.75 The generation routine of each grade of aggregate would be carried out until the total area of the generated particles exceeded the requirements. After the modeling completed, the total area of generated particles with different sizes was measured within PFC2D, as shown in Table 4. As shown, the generated virtual specimen was in good accordance with the actual gradations. The error was 5.51%, 3.86% and 3.70% for the 13.2-, 9.5- and 4.75-mm coarse aggregates, respectively. The error was attributed to the last generated particles, which exceeded the required mapping area, and it can be ignored. Thus, it is concluded that the proposed algorithm can well characterize the composition features of the aggregate specimen with precise particle shape reconstitutions and consistent gradation modeling.



The penetration test using a Universal Testing Machine (UTM: Fuel Instruments & Engineers

on the Chinese test standard (JTG E40-2007: T0134-1993) [36]. A cylinder container with a diameter

#### Pvt. Ltd. Maharashtra, India) was designed to measure the penetration force with the rod depth based **5. Performance Prediction of the Rebuilt Models**

#### of 150 mm and a height of 120 mm was adopted during the tests. Then, three specimens with only *5.1. Experiments*

the particles of 4.5, 9.5 and 13.2 mm were developed and were compacted based on the Chinese test standard (JTG E40-2007: T0131-2007) [36]. When conducting the penetration test, a rod would penetrate into the specimens slowly at a speed of 1.27 mm/min. Then, the force-displacement curves were recorded by the detectors. To help the reconstitution of the virtual specimens, several data should also be measured before or after the laboratory tests, including the particle density, the density and the void contents of the compacted specimens, as shown in Table 5. The particle density is an input variable in DEM simulations to model the entities' total mass, while the density and void contents of the compacted specimens were utilized to calculate the mapping area of coarse aggregates according to the proposed algorithm. The measurements for the particle density, the density and void contents of the compacted The penetration test using a Universal Testing Machine (UTM: Fuel Instruments & Engineers Pvt. Ltd. Maharashtra, India) was designed to measure the penetration force with the rod depth based on the Chinese test standard (JTG E40-2007: T0134-1993) [36]. A cylinder container with a diameter of 150 mm and a height of 120 mm was adopted during the tests. Then, three specimens with only the particles of 4.5, 9.5 and 13.2 mm were developed and were compacted based on the Chinese test standard (JTG E40-2007: T0131-2007) [36]. When conducting the penetration test, a rod would penetrate into the specimens slowly at a speed of 1.27 mm/min. Then, the force-displacement curves were recorded by the detectors.

To help the reconstitution of the virtual specimens, several data should also be measured before or after the laboratory tests, including the particle density, the density and the void contents of the compacted specimens, as shown in Table 5. The particle density is an input variable in DEM simulations to model the entities' total mass, while the density and void contents of the compacted specimens were utilized to calculate the mapping area of coarse aggregates according to the proposed algorithm. The measurements for the particle density, the density and void contents of the compacted specimens were carried out based on the Chinese test standard, respectively (JTG E42-2005: T0308-2005, T0309-2005) [37].


**Table 5.** Materials for laboratory penetration tests.

#### *5.2. Simulations for Performance Prediction*

Prior to simulations, several micromechanical models were adopted to characterize the contact behavior between coarse aggregates within PFC2D; including the contact-stiffness model and sliding model [38], as the following equation shows.

Contact-stiffness model:

$$\begin{aligned} F\_i^n &= \mathcal{K}^n \mathcal{U}^n\\ \Delta F\_i^s &= -k^s \Delta \mathcal{U}^s \end{aligned} \tag{13}$$

where *F n i* is the normal force at the i-th contact point; ∆*F s i* is the shear contact force increment at the *i*-th contact point; *K n* is the normal stiffness at the contact; *k s* is the shear stiffness at the contact; *U<sup>n</sup>* is the overlap of the two contact entities along the normal direction; ∆*U<sup>s</sup>* is the overlap of the two contact entities along the shear direction.

Sliding model:

$$F\_{\text{max}}^{\text{s}} = \mu |F\_i^{\text{n}}| \tag{14}$$

where *F s max* is the maximum allowable shear contact force; *F n i* is the normal contact force between two entities; *µ* is the friction coefficient in PFC2D.

As shown in Equations (13) and (14), the contact force between two entities (balls) in the normal/shear direction was determined by the contact normal/shear stiffness and their overlap in corresponding directions. As for the normal contact force, the formula is always effective in the simulation, such that the intensity of *F n i* depends on the value of *K <sup>n</sup>* and *U<sup>n</sup>* directly. However, the shear contact force was controlled and calculated based on two different modes for different periods. At the beginning, the shear contact force increment depends on *k <sup>s</sup>* and ∆*U<sup>s</sup>* , when the adjacent entities contact each other. As the accumulative shear force reaches the limiting condition, as shown in Equation (14), the calculation process is converted immediately, in which the intensity of the shear contact force (not the shear force increment) is determined by the friction coefficient and absolute value of the normal contact force currently. As shown, it is significantly important to determine the values of *K <sup>n</sup>* and *k s* for modeling the micro-contact process in the simulation. According to the PFC2D manual, these two key micro-parameters are calculated by the input value of the balls' normal stiffness *k<sup>n</sup>* and shear stiffness *ks* , respectively. The appropriate values of *k<sup>n</sup>* and *k<sup>s</sup>* can be obtained through a calibration process during the simulation based on the laboratory results. When calibrating the micro-parameters, the size influences should be taken into consideration. This means that the micro-parameters differed from each other for different sizes of particles and should be determined based on the size, respectively.

Thus, the simulations were conducted in two steps. In the first step, the virtual penetration tests of three different compositions were simulated. The three virtual specimens were developed by the proposed algorithm with 4.75-, 9.5- and 13.2-mm particles only, respectively. Then, the calibrated process for the best micro-parameters was done as follows: (1) conduct the laboratory penetration test to get the actual force-displacement curve of the three specimens; (2) conduct the simulations to get the virtual force-displacement curve of the three specimens; (3) obtain the best group of micro-parameters by adjusting their values continuously until the force-displacement curve of the simulation is close to the actual result. By the former process, the calibrated parameters can be obtained as shown in Table 6.


9.5 2.5e6 2.5e6 0.55

wall 1e10 1e10 0.6

ball

**Table 6.** Calibrated micro-parameters for balls and walls within Particle Flow Code in Two Dimensions (PFC2D). Dimensions (PFC2D).

*Appl. Sci.* **2018**, *2*, 15 FOR PEER REVIEW 13 of 16

simulation is close to the actual result. By the former process, the calibrated parameters can be

**Table 6.** Calibrated micro-parameters for balls and walls within Particle Flow Code in Two


tested with the rod moving downward at a speed of 1.27 mm/min. Each particle was assigned the best micro-parameters based on Table 6. The prediction results and the actual force-displacement curve are illustrated in Figure 9. Since assumptions were made that there was no breakage or fraction during modeling, only 10-mm penetration depths were simulated. When the depth was less than 10 mm, there was no obvious breakage of the particles in laboratory tests, which can be used for nonlinear mechanical behavior predictions. It is obviously seen that the developed virtual models indeed have a good agreement with the results of the laboratory tests. The proposed algorithm was validated, which can not only model the aggregate shape precisely, but also characterize its mechanical performance well with the calibrated micro-parameters. force-displacement curve are illustrated in Figure 9. Since assumptions were made that there was no breakage or fraction during modeling, only 10-mm penetration depths were simulated. When the depth was less than 10 mm, there was no obvious breakage of the particles in laboratory tests, which can be used for nonlinear mechanical behavior predictions. It is obviously seen that the developed virtual models indeed have a good agreement with the results of the laboratory tests. The proposed algorithm was validated, which can not only model the aggregate shape precisely, but also characterize its mechanical performance well with the calibrated micro-parameters.

**Figure 9.** Prediction results of SMA-13.

**Figure 9.** Prediction results of SMA-13.

#### **6. Conclusions**

**6. Conclusions** This paper proposed an algorithm for the reconstitution of the virtual aggregates and specimens based on the DEM methods. By the triangle area's division, the irregular shape of the coarse This paper proposed an algorithm for the reconstitution of the virtual aggregates and specimens based on the DEM methods. By the triangle area's division, the irregular shape of the coarse aggregates can be formed well, and their mechanical performance was validated through the designed simulations. The main conclusions drawn from the study are as follows:

(1) The scanning images from the AIMS are significant and guarantee the precise modelling of the coarse aggregates. The realistic shapes of the aggregates were captured and quantified through the inner software within AIMS. By generating the balls in the shape contours and then filling the inner area with balls by the proposed triangle-dividing algorithm, the virtual particles could be rebuilt, which had the same shapes as the realistic ones;


**Author Contributions:** D.W. and T.M. conceived and designed the experiments; X.D. performed the numerical experiments; W.Z. performed the laboratory experiments; D.Z. contributed analysis tools; D.W. and X.D. wrote the paper.

**Acknowledgments:** The study was financially supported by the Science Foundations of Nanjing Institute of Technology, China (YKJ201423), the National Natural Science Foundation of China (No. 51378006) and the Natural Science Foundation of Jiangsu Province, China (BK20161421).

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**


© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

**(MSCR) Test**
