Algorithm for Virtual Aggregates’ Reconstitution Based on Image Processing and Discrete-Element Modeling

: 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 ﬁrstly. The shape images were then processed, and the morphological characteristics of aggregates were quantiﬁed 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 veriﬁed 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.


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][2][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 Appl. Sci 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. 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) 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 Ycoordinate 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. 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. 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): where 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 i , y i are the X-and Y-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.

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: 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.

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: 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.

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 1 and S 2 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 4 and d 5 were calculated as shown in Equation (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. 1 and 2 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. 1 , 2 , 3 , 4 and 5 were calculated as shown in Equation (2). When 2 ≤ 1 (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);  ②Determine the variable 1 . It is calculated as the following equation: When S 2 ≤ S 1 (Figure 3a), the algorithm for filling the triangle area with balls by PFC2D was conducted as following: • 1 Start the ball generation from Point A with the coordinate of (0, 0); • 2 Determine the variable x 1 . It is calculated as the following equation: 6 of 16 where x 1 is the X-coordinate of the next column of balls that will be generated; x 1 f is the X-coordinate of the former column of balls that has been generated; rad is the radius of the filled balls. • 3 Determine the variable h, calculated as shown in Equation (4): where h is the required height of the next column of balls; the others are the same as the former. • 4 Generate uniform balls with coordinates of (x 1 , 0), (x 1 , 2 × rad), (x 1 , 4 × rad), (x 1 , 6 × rad) ... until the Y-coordinate exceeds h; • 5 Start to generate the next column of balls by cycling from Steps 2-4 until the X-coordinate exceeds S 1 .
When S 2 > S 1 (Figure 3b), the algorithm is conducted as follows: • 1 Start the ball generation from Point A with the coordinate of (0,0); • 2 Determine the variable x 1 . It is calculated as Equation (3) (5) and (6): where h 1 is the required height of the next column of balls; the others are the same as the former.

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).
where x and y 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.  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 3 (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 3 (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 3 (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.

Mapping Area Calculation for Two-Dimensional Specimens
When it comes to the virtual aggregate specimens, several assumptions were made before the simulations as follows:  ①The differences of the densities among particles were not obvious, so the same density value of 1 can be assigned to all the virtual particles in the simulations;  ②The virtual specimens consist of particles and voids only without any water.
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.

Gradation
Passing Ratio (%) for Different Sieving Sizes (mm) 16  The particle numbers of different grades in the two-dimensional simulation were determined based on the mapping area, as the following equations shown: where M is the total mass of the specimen in three dimensions, g; is the density of the specimen, g cm 3 ⁄ ; 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 is the converted mass of the specimen in two dimensions, g; is the total volume of the specimen in three dimensions, cm 3 ; is the area of the specimen in two dimensions, cm 2 ; the others are the same as the former.
Thus, the required area for generating the total virtual particles could be calculated based on Equation (10).
where is the required area for generating the total virtual particles, cm 2 ; 1 is the density of all the coarse aggregates, which was assumed to be the same value in the simulations, g cm 3 ⁄ ; the others are the same as the former.

Mapping Area Calculation for Two-Dimensional Specimens
When it comes to the virtual aggregate specimens, several assumptions were made before the simulations as follows: • 1 The differences of the densities among particles were not obvious, so the same density value of ρ 1 can be assigned to all the virtual particles in the simulations; • 2 The virtual specimens consist of particles and voids only without any water.
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 particle numbers of different grades in the two-dimensional simulation were determined based on the mapping area, as the following equations shown: where M is the total mass of the specimen in three dimensions, g; ρ is the density of the specimen, g/cm 3 ; 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 converted mass of the specimen in two dimensions, g; V is the total volume of the specimen in three dimensions, cm 3 ; S is the area of the specimen in two dimensions, cm 2 ; the others are the same as the former. Thus, the required area for generating the total virtual particles could be calculated based on Equation (10). where S total is the required area for generating the total virtual particles, cm 2 ; ρ 1 is the density of all the coarse aggregates, which was assumed to be the same value in the simulations, g/cm 3 ; 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: 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: 9.5, 13.2, 16 (12) where

Generation Process for Virtual Particles with Random Shapes and Sizes
To distribute the virtual particles randomly within the specimen, two parameters (µ 1 and µ 2 ) were given to control the generation process better. When generating the next particle, µ 1 represented the probability of the size, while µ 2 determined its angular shape. Two random values drawn from the uniform distribution (0, 1) were assigned to µ 1 and µ 2 at first. Based on µ 1 and µ 2 at the current calculation step, a specific size and shape could be determined. For example, before generating the next virtual particle, µ 1 would be checked within PFC2D as follows. , only 16-mm particles would be generated; 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.

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. Enough particles were scanned and rebuilt within PFC2D. Thus, when µ 2 was 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 µ 1 and µ 2 , the virtual specimen could be developed with random particles obviously.

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.

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. 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.

Experiments
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 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.

Experiments
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].

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: 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 n is the overlap of the two contact entities along the normal direction; ∆U s is the overlap of the two contact entities along the shear direction.
Sliding model: 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 n and U n 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 s and ∆U s , 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 n 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 n and shear stiffness k s , respectively. The appropriate values of k n and k s 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. In the second step, the performance prediction of the SMA13 aggregate skeleton was done by the virtual penetration test. Moreover, the validation of the proposed algorithm was verified by comparing the prediction results with the actual ones. The virtual aggregate skeleton of SMA13 was formed and 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. simulation is close to the actual result. By the former process, the calibrated parameters can be obtained as shown in Table 6. In the second step, the performance prediction of the SMA13 aggregate skeleton was done by the virtual penetration test. Moreover, the validation of the proposed algorithm was verified by comparing the prediction results with the actual ones. The virtual aggregate skeleton of SMA13 was formed and 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.

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

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 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