# Pavement Analysis and Design by Multiphysics Reconstructing Algorithm for the Virtual Asphalt Mixture Based on the Discrete-Element Method

Authors: Danhua Wang; Xunhao Ding; Tao Ma; Weiguang Zhang; Deyu Zhang
Journal: Advances in Civil Engineering
 (2018)

---

## Abstract

Based on the Particle Flow Code in Two dimensions (PFC2D), an algorithm for modeling the two-dimensional virtual asphalt mixture was proposed in this study. By combining the AIMS scanning technology (Aggregate Imaging Measurement System) with the designed stochastic algorithm, the virtual coarse aggregates could be generated rapidly and precisely. Different from the conventional methods, the contour shapes of the coarse aggregates were rebuilt only to balance the shape modeling precision and simulation efficiency. Then by distributing the coarse particles within container, virtual skeletons were formed firstly. An innovate algorithm was proposed afterwards to distinguish the external and internal area of the coarse aggregates and then model the mastic part by filling the irregular hollow shape with uniformly arranged balls. By deleting the mastic balls randomly, the voids were reconstructed consistent with the actual ratio. In the end, the virtual uniaxial compressive tests of AC-16 were simulated within PFC2D and the dynamic modulus at different load frequencies was predicted. The results indicated that the proposed algorithm could not only model the asphalt mixture precisely but also characterized its mechanical behavior as well.


---

## Body

## 1. Introduction

Asphalt mixtures are three-phase structures consisting of aggregates, asphalt binder, and air voids. It has been highlighted by many researchers that the mechanical performance of asphalt mixture was influenced by its heterogeneous components significantly [1–7]. Tests conducted by Liu and Dai [8] focused on the gradation influence on the rut resistance of the asphalt mixtures. The results showed that based on the Superpave Gyratory Compactors (SGC), the optimum asphalt content is lower 0.2%∼0.5% than the one by the Marshall method. As the aggregate size increased, the mixtures could have larger resistance compared to the smaller particles. Similar results were also concluded based on the studies of Coleri et al. [9] and Gokhale et al., respectively [10]. Focusing on the aggregate degradation during compaction, several designed tests were carried out by Airey et al. [11] to reveal the inner mechanism. In their tests, two asphalt mixtures including a continuous and gap-graded gradation were compacted and compared. The results indicated that gap-graded asphalt mixture was degraded more easily compared to the continuously grade mixture. Isailović et al. [12] pointed out that asphalt ageing played an important role in reducing the effects of material recovery based on fatigue tests. In addition, the mixtures prepared with additionally 0.5% of bitumen by mass showed the best recovery characteristics compared to others. The void effects on fatigue life of asphalt mixture were analyzed by Hasan and Ahmad in early 1973 [13]. The voids within mixtures were varied by the changes of asphalt content and aggregate gradation. The fatigue response of different mixtures was investigated showing the important roles the voids played in. Apart from the fatigue performance, the rutting-resistance, moisture damage-resistance, and strength of asphalt mixture would also be influenced by air void, which was founded by some researchers [14, 15]. Current design methods for asphalt mixtures mostly focused on the volumetric properties (aggregate gradation, asphalt content, and void) and paid little attention on the morphological characteristics of the components which were really hard to investigate just by laboratory apparatus. It is believed that, apart from the volumetric properties, the shapes of the inner components, especially the aggregates, indeed have impact on the properties of asphalt mixtures as well [16–18].As an effective numerical technology, the Discrete-Element Method (DEM) was proposed by Cundall [19], and then a software named as Particle Flow Code (PFC) was developed to help the mechanical calculation of heterogeneous materials, which has been widely used in many fields [20, 21]. Not only the macromechanical behavior but also the microresponse of the asphalt mixtures could be predicted based on the PFC, and many virtual tests were conducted under complex conditions for various goals by some researchers [22–24]. However, prior to simulations, it is significantly important to develop the precise virtual models based on their volumetric and shape property which is a guarantee for the results validations [25]. Thus, many methods were proposed for heterogeneous component reconstitution within asphalt mixtures which fell into two categories in general. One is the random model and another is the image-based models. The random models tended to develop the virtual particles by some designed stochastic algorithm. By changing the control parameters related to shapes and sizes, a large number of the simplified virtual shapes could be developed rapidly. Such methods were proposed by Lu and Mcdowell [25], Das et al. [26], and Zhang et al. [27]. In Lu and Mcdowell’s studies, the virtual particles were modeled by the overlapping balls [25], and the modeling algorithm was optimized further by Das et al. [26]. Coarse aggregates were assumed to be hexahedrons, pentahedrons, and tetrahedrons by Zhang et al. [27]. By cutting the particle clumps randomly with the help of three variables, the irregular shape was rebuilt preliminarily which was roughly consistent with the realistic ones. The random models for virtual particles could be generated in larger amounts quickly without preparing specimens in laboratory. Although it was performed with better efficiency compared to the image-based models, the shape modeling was not precise enough with excessive simplifications. The image-based models are mostly generated with the help of the X-ray computed tomography (CT) scanning technology [28]. By processing the component images from the scanning, the heterogeneous materials could be reconstructed within PFC precisely. Such methods were introduced and utilized in some studies [29–33]. However, the use of the image-based models are limited by the test environments and conditions which is a barrier preventing the virtual modeling to be conducted efficiently. First, a scanning device is must which provides the fundamental images. Second, when rebuilding image-based models, test specimens should be developed for scanning in advance. Therefore, it is time-consumed and cannot meet the needs of large amount simulation tasks. Thus, further studies are needed to optimize the particle modeling which can be generated in quantities rapidly and maintain the realistic shapes at the same time.


## 2. Objective and Scope

The objective of this study is to rebuild the virtual asphalt mixtures rapidly and precisely based on the Discrete-Element Method (DEM). To achieve this, with the help of AIMS scanning technology, an algorithm was designed for modeling coarse aggregates, asphalt mastic, and voids. The relevant issues include the shape and size measurements through AIMS, irregular particle reconstitution derived from the standard virtual aggregates, asphalt mastic generation by irregular area judgments, random voids modeling, and the mechanical behavior prediction.


## 3. Shape Measurement through Scanning

The AIMS apparatus was utilized to record the particle shapes including the angularity index and size index. While scanning the samples, coarse aggregates should be put into slots of designed trays. With the trays rotated slowly, each particle would go through the scanning area, and a digital camera was set to capture their features. Shape properties of particles were calculated and stored and could be output for the virtual modeling within PFC2D. The related shape measurements are introduced as follows.


### 3.1. Angularity Index

Angularity index applies to coarse aggregate sizes and describes variations at the particle boundary that influence the overall shape. The angularity index quantifies changes along a particle boundary with higher gradient values indicating a more angular shape. Angularity index has a relative scale of 0 to 10000 with a perfect circle having a small nonzero value as shown in Figure1. The angularity index is analyzed by quantifying the change in the gradient on a particle boundary and is related to the sharpness of the corners of 2-dimensional images of aggregate particles as shown in the following equation:(1)GA=1n/3−1∑i=1n−3θi−θi+3,where θ is the angle of orientation of the edge points; n is the total number of points; and subscript i denotes the ith point on the edge of the particle.Figure 1
Angularity index for various shapes. (a) 0 ≤ low ≤ 2100; (b) 2100 < moderate ≤ 3975; (c) 3975 < high ≤ 5400; (d) 5400 < extreme ≤ 10000.
(a)
(b)
(c)
(d)


### 3.2. Size Index

The particle’s shortest, intermediate, and longest lengths are measured to describe the size characteristics as shown in Figure2. The value L is the particle’s longest length representing the largest size in all the directions. And the particle’s intermediate length w is the largest size in the planes perpendicular to the particle’s longest length. The particle’s shortest length t is the largest size of the planes perpendicular to the particle’s longest and intermediate lengths at the same time.Figure 2
Size index for aggregate particles. (a) Side view; (b) top view.
(a)
(b)


## 3.1. Angularity Index

Angularity index applies to coarse aggregate sizes and describes variations at the particle boundary that influence the overall shape. The angularity index quantifies changes along a particle boundary with higher gradient values indicating a more angular shape. Angularity index has a relative scale of 0 to 10000 with a perfect circle having a small nonzero value as shown in Figure1. The angularity index is analyzed by quantifying the change in the gradient on a particle boundary and is related to the sharpness of the corners of 2-dimensional images of aggregate particles as shown in the following equation:(1)GA=1n/3−1∑i=1n−3θi−θi+3,where θ is the angle of orientation of the edge points; n is the total number of points; and subscript i denotes the ith point on the edge of the particle.Figure 1
Angularity index for various shapes. (a) 0 ≤ low ≤ 2100; (b) 2100 < moderate ≤ 3975; (c) 3975 < high ≤ 5400; (d) 5400 < extreme ≤ 10000.
(a)
(b)
(c)
(d)


## 3.2. Size Index

The particle’s shortest, intermediate, and longest lengths are measured to describe the size characteristics as shown in Figure2. The value L is the particle’s longest length representing the largest size in all the directions. And the particle’s intermediate length w is the largest size in the planes perpendicular to the particle’s longest length. The particle’s shortest length t is the largest size of the planes perpendicular to the particle’s longest and intermediate lengths at the same time.Figure 2
Size index for aggregate particles. (a) Side view; (b) top view.
(a)
(b)


## 4. Algorithm for Developing Virtual Coarse Aggregates

### 4.1. Standard Particles Generations

The virtual coarse aggregates were rebuilt within PFC2D as shown in Figure3. It is noted that with the goal of developing virtual specimen rapidly and randomly, not all the particles in actual mixtures should be scanned which would consume excessive time. It is meaningful to conduct simulations by generating representative virtual particles with enough irregular shapes. Thus, as an effective mean, when conducting other simulations, the representative virtual particles can be generated rapidly without scanning the components again.Figure 3
Rebuilding virtual particles by filling contours with balls.
(a)
(b)As shown in (1), it is known that the angularity index quantifies changes along a particle boundary regardless of its size. Thus, shapes were modeled firstly without size variations. Enough coarse aggregates were selected firstly for rebuilding the standard particles with their intermediate lengths w (defined as the minimum sieve size that the particle is able to pass) converted to 19 mm compulsively. And the areas of all the adjusted particles were measured. Then, the particles of other sizes were derived from the standard ones by size adjustments. While generating the standard particles, the selected samples for particle reconstitution should have a large angularity index range from 1000 to 6500 which most realistic aggregates are distributed in. Then, the shape images were processed by filtering the redundant black pixels. Only some key black pixels in the shape contour were retained and then were imported to the PFC2D. The key black pixels could describe the boundary changes and capture the main properties of the irregular shapes. By generating balls in the positions of the retained pixels, the virtual particles could be developed as a rigid body through the input Fish command “clump,” as shown in Figure 3. 27 particles with various angular shapes were rebuilt by the proposed methods. More angular shapes could also be modeled in the same way depending on the simulation needs but were not included here. Since it is very time-consumed to rebuild each particle in reality, the representative shape modeling, as a more efficient mean, could save time and develop virtual mixtures rapidly. In this study, the following proposed algorithm was introduced emphatically.Assumptions were made that there was no fracture or break of virtual particles; thus, the inner balls of the virtual particles could be released which could conduct the simulations more effectively. The particle surfaces would bear the contact force and were simulated by the contour balls as shown in Figure3. Compared to the solid virtual particles, the mass of the proposed methods was smaller than the reality and should be calibrated for all virtual particles as in the following equation.(2)p2=p1×m×π×r2n×π×0.5/k2=p1×4×m×rk2n,where p2 is the calibrated density (g/cm3), p1 is the initial density (g/cm3), m is the number of retained pixels, n is the number of all the pixels in Figure 1(a), k is the scaling ratio in image processing, and r is the radius of the ball located in the contour (mm).


### 4.2. Particle Size Adjustments

Each standard particle was coded and saved as an executable file. When generations started, corresponding code files were evoked and executed by PFC2D. To model the size variations, a designed controls parameter known asSr were embedded in the code files, as shown in Table 1. As shown, the size ratio is defined as the intermediate length ratio between the other grade aggregates and the standard ones. And the size ratios of each grade aggregates were coded with the help of the input “urand” value which was drawn from the uniform distribution U0,1 within PFC2D. Before the code files were evoked, the Sr was determined randomly by the system leading to the size generations. Moreover, the areas of different grades were adjusted together by multiplying the standard particle area with their size ratio squared.Table 1
Size adjustments for simulations.
Size (mm)
Size ratio
Controls parameter
Size (mm)
Size ratio
Controls parameter
19
1
1
13.2–16
13.2
19
,
16
19
13.2
19
+
16
19
−
13.2
19
×
urand
19–25
1
,
25
19
1
+
25
19
−
1
×
urand
9.5–13.2
9.5
19
,
13.2
19
9.5
19
+
13.2
19
−
9.5
19
×
urand
16–19
16
19
,
1
16
19
+
1
−
16
19
×
urand
4.75–9.5
4.75
19
,
9.5
19
4.75
19
+
9.5
19
−
4.75
19
×
urandFigure4 shows the virtual size reconstitutions based on standard particles. As shown, the standard particles with 19 mm intermediate lengths were developed firstly, and then the others were derived from the standard ones with same angular shapes. By this way, 135 virtual particles with 27 angular shapes of 5 different grades were developed rapidly.Figure 4
Virtual size reconstitutions based on standard particles.


### 4.3. Coarse Particles Distribution within Mixture

The virtual asphalt mixtures were developed in 3 steps: the coarse aggregate generation, the asphalt mastic filling, and the voids modeling. It is unpractical to rebuild all the particles’ shapes especially the numerous finer aggregates, which will exceed the computer capacity limitations. On the other side, the shape influences of the finer aggregates are not obvious to some extent compared to the coarse parts. So the major assumptions are made that the aggregates finer than 4.75 mm are regarded as the part of the asphalt mastic and are modeled by round balls directly within PFC2D.Coarse aggregates were distributed randomly within the container area. The gradations of particles and the asphalt content were calculated and converted to mapping area in two dimensions. The mapping area of different grades particles were taken as control values which determined the end time of coarse aggregate generations. The coarse aggregates distributions were conducted as follows:(1)
Determine the size of the virtual container and develop corresponding walls to model the boundary conditions.(2)
Start the coarse particles generation with the size ranging from 16 mm to 19 mm. The virtual particles are generated one by one and the shape is selected randomly by executing the particle code files.(3)
Develop three variables known asx1, y1, and rot which is calculated as shown in (3). These three variables are all drawn from the uniform distribution with different ranges. Before a virtual particle is generated, x1, y1, and rot embedded in the code files are determined randomly. And then the local coordinate system (coordinate system of particles) is converted to a global coordinate system (coordinate system of container). By using the random values of x1, y1, and rot, the coarse particles can be generated at a random position within the container area.(3)x1=U0,xw,y1=U0,yh,rot=U0,2π,where x1 and y1 are the random x- and y-coordinates of the initial positions, respectively; rot is the random angle of the generated particles. xw and yh are the width and height of the virtual containers.(4)
Check the overlap of the generated particles. Since the particles are generated within container area one by one. So when a particle has been generated at a random position, it is necessary to check if it is overlapped with the others. If there is no overlap, the particle can be generated, else, changex1, y1, and rot until the particle can be generated correctly;(5)
Start to generate the next particle by cycling from step 2 to 4 until reach the required mapping area of 16–19 mm particles.By the steps from 1 to 5, the virtual particles with the size ranging from 16 mm to 19 mm could be developed successfully. Then similar processes were carried out for the particles of other grades as shown in Figure5. To distribute the different grades particles correctly, the particles with largest sizes should be generated first and then the smaller ones. For example, start the particle generation of 16–19 mm firstly, and then followed by 13.2–16 mm, 9.5–13.2 mm, and 4.75–9.5 mm in turns.Figure 5
Coarse particles distributions within the virtual container.


## 4.1. Standard Particles Generations

The virtual coarse aggregates were rebuilt within PFC2D as shown in Figure3. It is noted that with the goal of developing virtual specimen rapidly and randomly, not all the particles in actual mixtures should be scanned which would consume excessive time. It is meaningful to conduct simulations by generating representative virtual particles with enough irregular shapes. Thus, as an effective mean, when conducting other simulations, the representative virtual particles can be generated rapidly without scanning the components again.Figure 3
Rebuilding virtual particles by filling contours with balls.
(a)
(b)As shown in (1), it is known that the angularity index quantifies changes along a particle boundary regardless of its size. Thus, shapes were modeled firstly without size variations. Enough coarse aggregates were selected firstly for rebuilding the standard particles with their intermediate lengths w (defined as the minimum sieve size that the particle is able to pass) converted to 19 mm compulsively. And the areas of all the adjusted particles were measured. Then, the particles of other sizes were derived from the standard ones by size adjustments. While generating the standard particles, the selected samples for particle reconstitution should have a large angularity index range from 1000 to 6500 which most realistic aggregates are distributed in. Then, the shape images were processed by filtering the redundant black pixels. Only some key black pixels in the shape contour were retained and then were imported to the PFC2D. The key black pixels could describe the boundary changes and capture the main properties of the irregular shapes. By generating balls in the positions of the retained pixels, the virtual particles could be developed as a rigid body through the input Fish command “clump,” as shown in Figure 3. 27 particles with various angular shapes were rebuilt by the proposed methods. More angular shapes could also be modeled in the same way depending on the simulation needs but were not included here. Since it is very time-consumed to rebuild each particle in reality, the representative shape modeling, as a more efficient mean, could save time and develop virtual mixtures rapidly. In this study, the following proposed algorithm was introduced emphatically.Assumptions were made that there was no fracture or break of virtual particles; thus, the inner balls of the virtual particles could be released which could conduct the simulations more effectively. The particle surfaces would bear the contact force and were simulated by the contour balls as shown in Figure3. Compared to the solid virtual particles, the mass of the proposed methods was smaller than the reality and should be calibrated for all virtual particles as in the following equation.(2)p2=p1×m×π×r2n×π×0.5/k2=p1×4×m×rk2n,where p2 is the calibrated density (g/cm3), p1 is the initial density (g/cm3), m is the number of retained pixels, n is the number of all the pixels in Figure 1(a), k is the scaling ratio in image processing, and r is the radius of the ball located in the contour (mm).


## 4.2. Particle Size Adjustments

Each standard particle was coded and saved as an executable file. When generations started, corresponding code files were evoked and executed by PFC2D. To model the size variations, a designed controls parameter known asSr were embedded in the code files, as shown in Table 1. As shown, the size ratio is defined as the intermediate length ratio between the other grade aggregates and the standard ones. And the size ratios of each grade aggregates were coded with the help of the input “urand” value which was drawn from the uniform distribution U0,1 within PFC2D. Before the code files were evoked, the Sr was determined randomly by the system leading to the size generations. Moreover, the areas of different grades were adjusted together by multiplying the standard particle area with their size ratio squared.Table 1
Size adjustments for simulations.
Size (mm)
Size ratio
Controls parameter
Size (mm)
Size ratio
Controls parameter
19
1
1
13.2–16
13.2
19
,
16
19
13.2
19
+
16
19
−
13.2
19
×
urand
19–25
1
,
25
19
1
+
25
19
−
1
×
urand
9.5–13.2
9.5
19
,
13.2
19
9.5
19
+
13.2
19
−
9.5
19
×
urand
16–19
16
19
,
1
16
19
+
1
−
16
19
×
urand
4.75–9.5
4.75
19
,
9.5
19
4.75
19
+
9.5
19
−
4.75
19
×
urandFigure4 shows the virtual size reconstitutions based on standard particles. As shown, the standard particles with 19 mm intermediate lengths were developed firstly, and then the others were derived from the standard ones with same angular shapes. By this way, 135 virtual particles with 27 angular shapes of 5 different grades were developed rapidly.Figure 4
Virtual size reconstitutions based on standard particles.


## 4.3. Coarse Particles Distribution within Mixture

The virtual asphalt mixtures were developed in 3 steps: the coarse aggregate generation, the asphalt mastic filling, and the voids modeling. It is unpractical to rebuild all the particles’ shapes especially the numerous finer aggregates, which will exceed the computer capacity limitations. On the other side, the shape influences of the finer aggregates are not obvious to some extent compared to the coarse parts. So the major assumptions are made that the aggregates finer than 4.75 mm are regarded as the part of the asphalt mastic and are modeled by round balls directly within PFC2D.Coarse aggregates were distributed randomly within the container area. The gradations of particles and the asphalt content were calculated and converted to mapping area in two dimensions. The mapping area of different grades particles were taken as control values which determined the end time of coarse aggregate generations. The coarse aggregates distributions were conducted as follows:(1)
Determine the size of the virtual container and develop corresponding walls to model the boundary conditions.(2)
Start the coarse particles generation with the size ranging from 16 mm to 19 mm. The virtual particles are generated one by one and the shape is selected randomly by executing the particle code files.(3)
Develop three variables known asx1, y1, and rot which is calculated as shown in (3). These three variables are all drawn from the uniform distribution with different ranges. Before a virtual particle is generated, x1, y1, and rot embedded in the code files are determined randomly. And then the local coordinate system (coordinate system of particles) is converted to a global coordinate system (coordinate system of container). By using the random values of x1, y1, and rot, the coarse particles can be generated at a random position within the container area.(3)x1=U0,xw,y1=U0,yh,rot=U0,2π,where x1 and y1 are the random x- and y-coordinates of the initial positions, respectively; rot is the random angle of the generated particles. xw and yh are the width and height of the virtual containers.(4)
Check the overlap of the generated particles. Since the particles are generated within container area one by one. So when a particle has been generated at a random position, it is necessary to check if it is overlapped with the others. If there is no overlap, the particle can be generated, else, changex1, y1, and rot until the particle can be generated correctly;(5)
Start to generate the next particle by cycling from step 2 to 4 until reach the required mapping area of 16–19 mm particles.By the steps from 1 to 5, the virtual particles with the size ranging from 16 mm to 19 mm could be developed successfully. Then similar processes were carried out for the particles of other grades as shown in Figure5. To distribute the different grades particles correctly, the particles with largest sizes should be generated first and then the smaller ones. For example, start the particle generation of 16–19 mm firstly, and then followed by 13.2–16 mm, 9.5–13.2 mm, and 4.75–9.5 mm in turns.Figure 5
Coarse particles distributions within the virtual container.


## 5. Algorithm for Developing Asphalt Mixture

### 5.1. Asphalt Mastic Modeling

The asphalt mastic was modeled by the uniform arranged balls with a diameter of 0.001 m. As common practices, the conventional methods tended to fill the particle external area with finer balls (mastic part) by judging if there are balls already [27]. Because the virtual particles in conventional methods were filled with balls (solid shape), the external and internal areas of the particle could be distinguished easily by balls judgments. However, it is hard to fill in the mastic balls directly when it comes to the coarse aggregate skeleton developed in this study, as shown in Figure 5. The differences were not obvious between the external and internal particle areas which were all irregular in shape. To fill the irregular hollow shape (particle external area) with balls and ignore the particle internal area, the algorithm for the asphalt mastic was introduced as followed.(1)
Measure the size of total specimens to get all the coordinate values which are prepared for filling judgments.(2)
Develop two parameters known askright and kleft to record the position properties. They are all set zero initially. Then check all the coordinate values from left to right and bottom to top in turns. The positions will fall into four categories as shown in Figures 6–9. Before a group of coordinate values is checked, two parallel control lines are utilized with a 0.001 m interval (diameter of the mastic ball). Check the two sides of the position and find the nearest entities (walls and particle balls without mastic balls generated already) between two control lines as follows. The control lines were not real lines plotted in the images but a designed routine coded in MATLAB to help the distance judgments.(1)
If there is a nearest particle ball in the right side, convertkright from 0 to 1, else, kright is still equal to 0.(2)
If there is a nearest particle ball in left side, convertkleft from 0 to 1, else, kleft is still equal to 0.Figure 6
Distinguishing the position of mastic balls between the particle clump and the wall on the left side.Figure 7
Distinguishing the position of mastic balls between the particle clump and the wall on the right side.Figure 8
Distinguishing the position of mastic balls between two particle clumps (inside or outside the aggregates).Figure 9
Distinguishing the position of mastic balls without any entities nearby in two sides.Then the external and internal area of the particles can be distinguished by the combination ofkright and kleft. If kright and kleft are equal to 1 and 0, respectively, the position is between the particle clump and the left wall as shown in Figure 6. Then the coordinate value is effective, and a mastic ball with the diameter of 0.001 m is generated here.Ifkright and kleft are equal to 0 and 1, respectively, the position is between the particle clump and the right wall as shown in Figure 7. Then the coordinate value is effective, and a mastic ball with the diameter of 0.001 m is generated here as well.Ifkright and kleft are equal to 1 and 1, respectively, two conditions should be analyzed, as shown in Figure 8. As shown, the points A and B are between two particle clumps with same kright and kleft. The point A is inside the aggregates while B is outside the model between two different particles. In this condition, the nearest particle balls on both sides were analyzed further. With the help of the command pointer, the clump id of all the balls can be read and recognized. So comparing the clump id of the nearest balls, it is easy to determine whether the nearest balls belong to a same clump. If the clump id of two balls is equal, they come from the same particle clump. Otherwise, they belong to a different particle clump. Thus, if the nearest balls on left and right sides belong to a same particle clump, as the point A shown in Figure 8, the coordinate value is noneffective without any mastic balls generated here. When the nearest balls belong to different clumps, a mastic ball should be generated.A specific condition should also be included as shown in Figure9 due to the space between coarse aggregates. kright and kleft are 0 and 0, respectively, when the space is enough among particles. When comes to this condition, mastic balls should also be generated. Following the proposed rules and coded it within PFC2D, the virtual asphalt mastic could be rebuilt successfully as shown in Figure 10. As shown, the uniform arranged balls were filled in the external particle area to model the mastic while the coarse aggregates kept same shapes and positions all the time.Figure 10
Filling the irregular area with balls to model the asphalt mastic.


### 5.2. Escaped Balls Deletion and Voids Modeling

Prior to modeling the voids, enough calculation steps should be done within PFC2D to make the ball system stable. The mastic balls would move continuously until reach the equilibrium state. However, due to the initial overlaps among balls (Figure10), some of the mastic balls would go through the particle boundary and then escaped. So after the ball system reached the equilibrium state, an initial upward velocity was assigned to all balls and then simulated for a very short time. The balls without any contact force during the short-time simulation were identified as the escaped balls and should be deleted as shown in Figure 11.Figure 11
Escaped balls deletion before random voids modelling.The void content in three dimensions was converted to mapping area in two dimensions which was the same with the coarse aggregates and asphalt mastic. Since the diameters of the mastic balls have been known, the number of the void balls could be determined. By deleting the mastic balls randomly until the number of the void balls reached requirements, the virtual void could be modeled well as shown in Figure12(d).Figure 12
Modeling process of virtual asphalt mixtures: (a) aggregate skeleton; (b) uniform balls filling; (c) virtual asphalt mixture; (d) voids modeling; (e) asphalt mastic modelling.
(a)
(b)
(c)
(d)
(e)


### 5.3. Example of Generating a Virtual Specimen

The asphalt mixture of AC-16 was developed for simulation. The gradation of the AC-16 is shown in Table2. Based on the gradation, the mapping areas of three-phase structures including the coarse aggregates, asphalt mastic, and the void content were calculated using (4) and (5). And the results of the AC-16 were summarized in Table 3.(4)M=ρ×π×D24×h,Where M is the total mass of the specimen in three dimensions, g; ρ is the density of the specimen, g/cm3; D is the diameter of the specimen, cm; and h is the height of the specimen, cm.Table 2
Gradations for AC-16.
Gradation
Passing ratio (%)for different sieving size (mm)
19
16
13.2
9.5
4.75
2.36
1.18
0.6
0.3
0.15
0.075
AC-16
100
95
84
70
48
34
24.5
17.5
12.5
9.5
6.0Table 3
Mapping area calculations for three-phase structures.
Container area (m2)
Total mapping area for void (m2)
Total mapping area for coarse aggregates (m2)
Total mapping area for asphalt mastic (m2)
0.015
0.0006
0.0071
0.0073When it comes to the two-dimensional specimen, the mass can be converted as following:(5)m=MV×S=V×ρV×D×h=Dhρ,where m is the converted mass of specimen in two dimensions (g), V is the total volume of the specimen in three dimensions (cm3), S is the area of the specimen in two dimensions (cm2), and the others are the same as those of three dimensions.The void content was 4% and the asphalt content was 4.5% in mixtures. Assumptions were made that the air void is distributed uniformly in mixtures so the void contents in two dimensions and three dimensions can be set as a same value of 4%. Moreover, the densities of the coarse aggregates and asphalt mastic were 2.7 g/cm3 and 2.0 g/cm3, respectively. By the proposed algorithm, the virtual container with a height of 150 mm and a width of 100 mm were developed firstly, then followed by the coarse aggregate, mastic and void generation, respectively, within PFC2D as shown in Figure 12.To verify the precision of the inner components modeling, the final generated area of coarse aggregates, asphalt mastic, and voids were measured afterwards by the user-defined routine within PFC2D. The results were summarized in Tables4 and 5. As shown, the error is 5.1%, 4.7%, 2.1%, and 1.7% for the 16 mm, 13.2 mm, 9.5 mm, and 4.75 mm aggregates, respectively. The modeling error decreased as the size decreased. This is due to the overflow of the lastly generated coarse particles of each grade and it is inevitable. When comes to the asphalt mastic, the modeling is significantly precise with 0.014% error only. The sole difference between the realistic and virtual mastic content was caused by the rounding error in mapping area calculations.Table 4
Summaries of the coarse aggregates.
Size (mm)
16
13.2
9.5
4.75
Required mapping area (m2)
6.83e − 4
1.50e − 3
1.91e − 3
3.01e − 3
Generated area (m2)
7.18e − 4
1.57e − 3
1.95e − 3
3.06e − 3
Error (%)
5.1
4.7
2.1
1.7Table 5
Summaries of the asphalt mastic and voids.
A
B
C
D
E
F
Error (%)
7.3e − 3
7.301e − 3
9300
10076
764
12
0.014
A: required mapping area for asphalt mastic (m2); B: generated mapping area for asphalt mastic (m2); C: required number of balls within asphalt mastic; D: actual number of balls within asphalt mastic before void modeling; E: number of the deleted balls for modeling voids; F: number of the deleted balls which have been escaped.


## 5.1. Asphalt Mastic Modeling

The asphalt mastic was modeled by the uniform arranged balls with a diameter of 0.001 m. As common practices, the conventional methods tended to fill the particle external area with finer balls (mastic part) by judging if there are balls already [27]. Because the virtual particles in conventional methods were filled with balls (solid shape), the external and internal areas of the particle could be distinguished easily by balls judgments. However, it is hard to fill in the mastic balls directly when it comes to the coarse aggregate skeleton developed in this study, as shown in Figure 5. The differences were not obvious between the external and internal particle areas which were all irregular in shape. To fill the irregular hollow shape (particle external area) with balls and ignore the particle internal area, the algorithm for the asphalt mastic was introduced as followed.(1)
Measure the size of total specimens to get all the coordinate values which are prepared for filling judgments.(2)
Develop two parameters known askright and kleft to record the position properties. They are all set zero initially. Then check all the coordinate values from left to right and bottom to top in turns. The positions will fall into four categories as shown in Figures 6–9. Before a group of coordinate values is checked, two parallel control lines are utilized with a 0.001 m interval (diameter of the mastic ball). Check the two sides of the position and find the nearest entities (walls and particle balls without mastic balls generated already) between two control lines as follows. The control lines were not real lines plotted in the images but a designed routine coded in MATLAB to help the distance judgments.(1)
If there is a nearest particle ball in the right side, convertkright from 0 to 1, else, kright is still equal to 0.(2)
If there is a nearest particle ball in left side, convertkleft from 0 to 1, else, kleft is still equal to 0.Figure 6
Distinguishing the position of mastic balls between the particle clump and the wall on the left side.Figure 7
Distinguishing the position of mastic balls between the particle clump and the wall on the right side.Figure 8
Distinguishing the position of mastic balls between two particle clumps (inside or outside the aggregates).Figure 9
Distinguishing the position of mastic balls without any entities nearby in two sides.Then the external and internal area of the particles can be distinguished by the combination ofkright and kleft. If kright and kleft are equal to 1 and 0, respectively, the position is between the particle clump and the left wall as shown in Figure 6. Then the coordinate value is effective, and a mastic ball with the diameter of 0.001 m is generated here.Ifkright and kleft are equal to 0 and 1, respectively, the position is between the particle clump and the right wall as shown in Figure 7. Then the coordinate value is effective, and a mastic ball with the diameter of 0.001 m is generated here as well.Ifkright and kleft are equal to 1 and 1, respectively, two conditions should be analyzed, as shown in Figure 8. As shown, the points A and B are between two particle clumps with same kright and kleft. The point A is inside the aggregates while B is outside the model between two different particles. In this condition, the nearest particle balls on both sides were analyzed further. With the help of the command pointer, the clump id of all the balls can be read and recognized. So comparing the clump id of the nearest balls, it is easy to determine whether the nearest balls belong to a same clump. If the clump id of two balls is equal, they come from the same particle clump. Otherwise, they belong to a different particle clump. Thus, if the nearest balls on left and right sides belong to a same particle clump, as the point A shown in Figure 8, the coordinate value is noneffective without any mastic balls generated here. When the nearest balls belong to different clumps, a mastic ball should be generated.A specific condition should also be included as shown in Figure9 due to the space between coarse aggregates. kright and kleft are 0 and 0, respectively, when the space is enough among particles. When comes to this condition, mastic balls should also be generated. Following the proposed rules and coded it within PFC2D, the virtual asphalt mastic could be rebuilt successfully as shown in Figure 10. As shown, the uniform arranged balls were filled in the external particle area to model the mastic while the coarse aggregates kept same shapes and positions all the time.Figure 10
Filling the irregular area with balls to model the asphalt mastic.


## 5.2. Escaped Balls Deletion and Voids Modeling

Prior to modeling the voids, enough calculation steps should be done within PFC2D to make the ball system stable. The mastic balls would move continuously until reach the equilibrium state. However, due to the initial overlaps among balls (Figure10), some of the mastic balls would go through the particle boundary and then escaped. So after the ball system reached the equilibrium state, an initial upward velocity was assigned to all balls and then simulated for a very short time. The balls without any contact force during the short-time simulation were identified as the escaped balls and should be deleted as shown in Figure 11.Figure 11
Escaped balls deletion before random voids modelling.The void content in three dimensions was converted to mapping area in two dimensions which was the same with the coarse aggregates and asphalt mastic. Since the diameters of the mastic balls have been known, the number of the void balls could be determined. By deleting the mastic balls randomly until the number of the void balls reached requirements, the virtual void could be modeled well as shown in Figure12(d).Figure 12
Modeling process of virtual asphalt mixtures: (a) aggregate skeleton; (b) uniform balls filling; (c) virtual asphalt mixture; (d) voids modeling; (e) asphalt mastic modelling.
(a)
(b)
(c)
(d)
(e)


## 5.3. Example of Generating a Virtual Specimen

The asphalt mixture of AC-16 was developed for simulation. The gradation of the AC-16 is shown in Table2. Based on the gradation, the mapping areas of three-phase structures including the coarse aggregates, asphalt mastic, and the void content were calculated using (4) and (5). And the results of the AC-16 were summarized in Table 3.(4)M=ρ×π×D24×h,Where M is the total mass of the specimen in three dimensions, g; ρ is the density of the specimen, g/cm3; D is the diameter of the specimen, cm; and h is the height of the specimen, cm.Table 2
Gradations for AC-16.
Gradation
Passing ratio (%)for different sieving size (mm)
19
16
13.2
9.5
4.75
2.36
1.18
0.6
0.3
0.15
0.075
AC-16
100
95
84
70
48
34
24.5
17.5
12.5
9.5
6.0Table 3
Mapping area calculations for three-phase structures.
Container area (m2)
Total mapping area for void (m2)
Total mapping area for coarse aggregates (m2)
Total mapping area for asphalt mastic (m2)
0.015
0.0006
0.0071
0.0073When it comes to the two-dimensional specimen, the mass can be converted as following:(5)m=MV×S=V×ρV×D×h=Dhρ,where m is the converted mass of specimen in two dimensions (g), V is the total volume of the specimen in three dimensions (cm3), S is the area of the specimen in two dimensions (cm2), and the others are the same as those of three dimensions.The void content was 4% and the asphalt content was 4.5% in mixtures. Assumptions were made that the air void is distributed uniformly in mixtures so the void contents in two dimensions and three dimensions can be set as a same value of 4%. Moreover, the densities of the coarse aggregates and asphalt mastic were 2.7 g/cm3 and 2.0 g/cm3, respectively. By the proposed algorithm, the virtual container with a height of 150 mm and a width of 100 mm were developed firstly, then followed by the coarse aggregate, mastic and void generation, respectively, within PFC2D as shown in Figure 12.To verify the precision of the inner components modeling, the final generated area of coarse aggregates, asphalt mastic, and voids were measured afterwards by the user-defined routine within PFC2D. The results were summarized in Tables4 and 5. As shown, the error is 5.1%, 4.7%, 2.1%, and 1.7% for the 16 mm, 13.2 mm, 9.5 mm, and 4.75 mm aggregates, respectively. The modeling error decreased as the size decreased. This is due to the overflow of the lastly generated coarse particles of each grade and it is inevitable. When comes to the asphalt mastic, the modeling is significantly precise with 0.014% error only. The sole difference between the realistic and virtual mastic content was caused by the rounding error in mapping area calculations.Table 4
Summaries of the coarse aggregates.
Size (mm)
16
13.2
9.5
4.75
Required mapping area (m2)
6.83e − 4
1.50e − 3
1.91e − 3
3.01e − 3
Generated area (m2)
7.18e − 4
1.57e − 3
1.95e − 3
3.06e − 3
Error (%)
5.1
4.7
2.1
1.7Table 5
Summaries of the asphalt mastic and voids.
A
B
C
D
E
F
Error (%)
7.3e − 3
7.301e − 3
9300
10076
764
12
0.014
A: required mapping area for asphalt mastic (m2); B: generated mapping area for asphalt mastic (m2); C: required number of balls within asphalt mastic; D: actual number of balls within asphalt mastic before void modeling; E: number of the deleted balls for modeling voids; F: number of the deleted balls which have been escaped.


## 6. Performance Prediction of the Rebuilt Models

### 6.1. Experiments

Based on the Chinese test standards [34], uniaxial compressive test was conducted to evaluate the dynamic modulus of asphalt mixture in laboratory. The gradations of mixtures were shown in Table 2. And the mixtures were formed in a cylinder container with a height of 150 mm and a width of 100 mm. A Superpave gyratory compactor was utilized to compact the asphalt mixture to a targeted void level of 4% at the 4.5% asphalt content. The specimens were tested at 0°C and 138 kPa confining pressure and the loading frequencies are 0.1, 1, 5, 10, and 25 Hz, respectively.


### 6.2. Simulations for Dynamic Modulus Prediction

The virtual asphalt mixtures were developed as shown in Figure12(c). When simulated, the top and bottom walls were modeled as a sine load while the walls on left and right sides kept a constant confining stress based on the numerical servomechanism within PFC2D. During the simulation, the axial deviatoric stress and axial strain were recorded together for the dynamic modulus calculations.According to the PFC2D manuals [35], several constitutive models were used for characterizing the mechanical behavior of the heterogeneous materials including the contact-stiffness models, sliding models, and Burger’s models. The contact force in normal and shear direction between two entities was determined by the key microparameters in contact-stiffness models, known as the normal stiffness kn and shear stiffness ks. The sliding models were used for sliding movements by the frictional coefficient while Burger’s models characterized the viscoelasticity of asphalt part. To reduce the error of the selected microparameters, the parameter calibrations were conducted separately for the coarse aggregates and asphalt mastic. According to the Chinese test standards [36], the penetration tests of the 4.75 mm, 9.5 mm, 13.2 mm, and 16 mm particles were carried out to get the force-displacement curves. Then virtual simulations for corresponding penetration tests were developed within PFC2D and were conducted for virtual results. By adjusting the microparameters of different size particles continuously until the simulation curves match the realistic, the best group of normal stiffness kn , shear stiffness ks, and frictional coefficient could be determined as shown in Table 6. Two kinds of contact points were defined by Burger’s models. One is the contact point within asphalt, and another is between the asphalt mastic and aggregates. Based on the previous studies [37], the microparameters of Burger’s models at different contacts were shown in Tables 7 and 8.Table 6
Microparameters of selected constitutive models for particles and walls.
Normal stiffness,kn (N/m)
Shear stiffness,ks (N/m)
Friction coefficient
4.75
2.2e6
1.5e6
0.6
9.5
2.5e6
2.5e6
0.55
13.2
1.2e6
1.2e6
0.5
16
1.8e6
1.8e6
0.4
Wall
1e10
1e10
0.35Table 7
Microparameters at contacts within asphalt mastic.
K
mn (Pa·m)
K
kn (Pa·m)
C
mn (Pa·m·s)
C
kn (Pa·m·s)
K
ms (Pa·m)
K
ks (Pa·m)
C
ms (Pa·m·s)
C
ks (Pa·m·s)
1.14e3
7.92e2
1.95e6
5.58e4
3.79e2
2.64e1
6.49e5
1.86e4Table 8
Microparameters at contacts between asphalt mastic and aggregates.
K
mn
′ (Pa·m)
K
kn
′ (Pa·m)
C
mn
′ (Pa·m·s)
C
kn
′ (Pa·m·s)
K
ms
′ (Pa·m)
K
ks
′ (Pa·m)
C
ms
′ (Pa·m·s)
C
ks
′ (Pa·m·s)
2.27e3
1.58e3
3.89e6
1.12e5
1.51e3
5.28e5
1.30e6
3.72e4Figure13 shows the simulation results compared to the laboratory tests at the test temperature of 0°C. As shown, the virtual results have a good accordance with the truth. The prediction error of dynamic modulus is 12.7%, 6.7%, 6.6%, 3.8%, and 3.0% for the 25 Hz, 10 Hz, 5 Hz, 1 Hz, and 0.1 Hz, respectively. As shown in Figure 14, the predictions of phase angle is also effective with the minor error of 12.9%, 10.6%, 7.8%, 5.4% and 4.8% for the 25 Hz, 10 Hz, 5 Hz, 1 Hz, and 0.1 Hz, respectively. As the loading frequency decreased, the error decreased obviously. This is because the microparameters for simulations were all calibrated and determined by tests under the static loading rather than the cyclic loading. So when the load frequency is low in simulations, the cyclic loading could be taken as the static loading to some extent. When the load frequency increased, the difference between the static and cyclic loading grew leading to the error shown in the Figures 13 and 14. In summaries, the developed models can well predict the dynamic modulus of asphalt mixture, especially at low frequency, which verified the correctness of the proposed algorithm.Figure 13
Dynamic modulus prediction based on the proposed models.Figure 14
Phase angle prediction based on the proposed models.


## 6.1. Experiments

Based on the Chinese test standards [34], uniaxial compressive test was conducted to evaluate the dynamic modulus of asphalt mixture in laboratory. The gradations of mixtures were shown in Table 2. And the mixtures were formed in a cylinder container with a height of 150 mm and a width of 100 mm. A Superpave gyratory compactor was utilized to compact the asphalt mixture to a targeted void level of 4% at the 4.5% asphalt content. The specimens were tested at 0°C and 138 kPa confining pressure and the loading frequencies are 0.1, 1, 5, 10, and 25 Hz, respectively.


## 6.2. Simulations for Dynamic Modulus Prediction

The virtual asphalt mixtures were developed as shown in Figure12(c). When simulated, the top and bottom walls were modeled as a sine load while the walls on left and right sides kept a constant confining stress based on the numerical servomechanism within PFC2D. During the simulation, the axial deviatoric stress and axial strain were recorded together for the dynamic modulus calculations.According to the PFC2D manuals [35], several constitutive models were used for characterizing the mechanical behavior of the heterogeneous materials including the contact-stiffness models, sliding models, and Burger’s models. The contact force in normal and shear direction between two entities was determined by the key microparameters in contact-stiffness models, known as the normal stiffness kn and shear stiffness ks. The sliding models were used for sliding movements by the frictional coefficient while Burger’s models characterized the viscoelasticity of asphalt part. To reduce the error of the selected microparameters, the parameter calibrations were conducted separately for the coarse aggregates and asphalt mastic. According to the Chinese test standards [36], the penetration tests of the 4.75 mm, 9.5 mm, 13.2 mm, and 16 mm particles were carried out to get the force-displacement curves. Then virtual simulations for corresponding penetration tests were developed within PFC2D and were conducted for virtual results. By adjusting the microparameters of different size particles continuously until the simulation curves match the realistic, the best group of normal stiffness kn , shear stiffness ks, and frictional coefficient could be determined as shown in Table 6. Two kinds of contact points were defined by Burger’s models. One is the contact point within asphalt, and another is between the asphalt mastic and aggregates. Based on the previous studies [37], the microparameters of Burger’s models at different contacts were shown in Tables 7 and 8.Table 6
Microparameters of selected constitutive models for particles and walls.
Normal stiffness,kn (N/m)
Shear stiffness,ks (N/m)
Friction coefficient
4.75
2.2e6
1.5e6
0.6
9.5
2.5e6
2.5e6
0.55
13.2
1.2e6
1.2e6
0.5
16
1.8e6
1.8e6
0.4
Wall
1e10
1e10
0.35Table 7
Microparameters at contacts within asphalt mastic.
K
mn (Pa·m)
K
kn (Pa·m)
C
mn (Pa·m·s)
C
kn (Pa·m·s)
K
ms (Pa·m)
K
ks (Pa·m)
C
ms (Pa·m·s)
C
ks (Pa·m·s)
1.14e3
7.92e2
1.95e6
5.58e4
3.79e2
2.64e1
6.49e5
1.86e4Table 8
Microparameters at contacts between asphalt mastic and aggregates.
K
mn
′ (Pa·m)
K
kn
′ (Pa·m)
C
mn
′ (Pa·m·s)
C
kn
′ (Pa·m·s)
K
ms
′ (Pa·m)
K
ks
′ (Pa·m)
C
ms
′ (Pa·m·s)
C
ks
′ (Pa·m·s)
2.27e3
1.58e3
3.89e6
1.12e5
1.51e3
5.28e5
1.30e6
3.72e4Figure13 shows the simulation results compared to the laboratory tests at the test temperature of 0°C. As shown, the virtual results have a good accordance with the truth. The prediction error of dynamic modulus is 12.7%, 6.7%, 6.6%, 3.8%, and 3.0% for the 25 Hz, 10 Hz, 5 Hz, 1 Hz, and 0.1 Hz, respectively. As shown in Figure 14, the predictions of phase angle is also effective with the minor error of 12.9%, 10.6%, 7.8%, 5.4% and 4.8% for the 25 Hz, 10 Hz, 5 Hz, 1 Hz, and 0.1 Hz, respectively. As the loading frequency decreased, the error decreased obviously. This is because the microparameters for simulations were all calibrated and determined by tests under the static loading rather than the cyclic loading. So when the load frequency is low in simulations, the cyclic loading could be taken as the static loading to some extent. When the load frequency increased, the difference between the static and cyclic loading grew leading to the error shown in the Figures 13 and 14. In summaries, the developed models can well predict the dynamic modulus of asphalt mixture, especially at low frequency, which verified the correctness of the proposed algorithm.Figure 13
Dynamic modulus prediction based on the proposed models.Figure 14
Phase angle prediction based on the proposed models.


## 7. Conclusions

This paper proposed an algorithm for the reconstitution of the virtual asphalt mixture based on the DEM methods. By filling the shape contours with balls, the virtual coarse aggregates were developed precisely. Based on the algorithm for distinguishing the external and internal area of irregular particles, the asphalt mastic was generated properly by the designed rules. In the end, the validation of the proposed algorithm was verified by the virtual uniaxial compressive test. The main conclusions drawn from the study are as follows:(1)
Virtual shapes were rebuilt precisely based on the AIMS scanning. By the random control parameters, numerous virtual particles derived from the standard ones were developed rapidly with various shapes and sizes. By combining the scanning technology with the stochastic algorithm, this method can well meet the DEM simulation needs.(2)
The virtual asphalt mixtures were rebuilt in three steps. Firstly, the virtual coarse aggregates were distributed randomly to form the skeletons. Then an innovate algorithm was proposed to distinguish the external and internal area of the coarse particles. At last, by filling the designed area with uniform arranged balls and deleting mastic balls randomly, the asphalt mastic and void could be generated successfully. The results shows that the proposed algorithm can well model the inner components consistent to the actual gradations, asphalt content, and voids.(3)
The validation of the developed models was verified by the virtual uniaxial compressive test. The error of dynamic modulus is 12.7%, 6.7%, 6.6%, 3.8%, 3.0%, and the error of phase angle is 12.9%, 10.6%, 7.8%, 5.4%, and 4.8% for the 25 Hz, 10 Hz, 15 Hz, 1 Hz, and 0.1 Hz, respectively. It is indicated that the dynamic modulus of the asphalt mixtures can be well predicted based on the proposed algorithm, especially at low frequency.



---

*Source: 5180562-2018-05-07.xml*