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

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 ImagingMeasurement 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.(en 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. (e results indicated that the proposed algorithm could not only model the asphalt mixture precisely but also characterized its mechanical behavior as well.


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][2][3][4][5][6][7].Tests conducted by Liu and Dai [8] focused on the gradation influence on the rut resistance of the asphalt mixtures.e 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.e 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.e void effects on fatigue life of asphalt mixture were analyzed by Hasan and Ahmad in early 1973 [13].e voids within mixtures were varied by the changes of asphalt content and aggregate gradation.e 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][17][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][23][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].
us, 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.e 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.
e 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.e imagebased 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][30][31][32][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.erefore, it is time-consumed and cannot meet the needs of large amount simulation tasks.
us, 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.

Objective and Scope
e 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.e 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.

Shape Measurement through Scanning
e 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.
e 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.e 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 Figure 1. e 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: 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.

Size
Index.e particle's shortest, intermediate, and longest lengths are measured to describe the size characteristics as shown in Figure 2. e 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.e 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.

Algorithm for Developing Virtual
Coarse Aggregates

Standard Particles Generations.
e virtual coarse aggregates were rebuilt within PFC2D as shown in Figure 3.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.us, as an e ective mean, when conducting other simulations, the representative virtual particles can be generated rapidly without scanning the components again.
As shown in (1), it is known that the angularity index quanti es changes along a particle boundary regardless of its size.us, shapes were modeled rstly without size variations.Enough coarse aggregates were selected rstly for rebuilding the standard particles with their intermediate lengths w (de ned 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.en, 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.en, the shape images were processed by ltering the redundant black pixels.Only some key black pixels in the shape contour were retained and then were imported to the PFC2D.e 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 e cient 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 e ectively.e particle surfaces would bear the contact force and were simulated by the contour balls as shown in Figure 3. 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.
where p 2 is the calibrated density (g/cm 3 ), p 1 is the initial density (g/cm 3 ), 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).

Particle Size Adjustments.
Each standard particle was coded and saved as an executable le.When generations started, corresponding code les were evoked and executed by PFC2D.To model the size variations, a designed controls parameter known as S r were embedded in the code les, as shown in Table 1.As shown, the size ratio is de ned as the intermediate length ratio between the other grade aggregates Advances in Civil Engineering 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 U(0, 1 within PFC2D.Before the code les were evoked, the S r was determined randomly by the system leading to the size generations.Moreover, the areas of di erent grades were adjusted together by multiplying the standard particle area with their size ratio squared.
Figure 4 shows the virtual size reconstitutions based on standard particles.As shown, the standard particles with 19 mm intermediate lengths were developed rstly, 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 di erent grades were developed rapidly.

Coarse Particles Distribution within Mixture.
e virtual asphalt mixtures were developed in 3 steps: the coarse aggregate generation, the asphalt mastic lling, and the voids modeling.It is unpractical to rebuild all the particles' shapes especially the numerous ner aggregates, which will exceed the computer capacity limitations.On the other side, the shape in uences of the ner aggregates are not obvious to some extent compared to the coarse parts.So the major assumptions are made that the aggregates ner 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.e gradations of particles and the asphalt content were calculated and converted to mapping area in two dimensions.
e mapping area of di erent grades particles were taken as control values which determined the end time of coarse aggregate generations.
e 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. e virtual particles are generated one by one and the shape is selected randomly by executing the particle code les.(3) Develop three variables known as x 1 , y 1 , and rot which is calculated as shown in (3).ese three variables are all drawn from the uniform distribution with di erent ranges.Before a virtual particle is generated, x 1 , y 1 , and rot embedded in the code les 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 x 1 , y 1 , and rot, the coarse particles can be generated at a random position within the container area.
where x 1 and y 1 are the random x-and y-coordinates of the initial positions, respectively; rot is the random angle of the generated particles.x w and y h 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, change x 1 , y 1 , 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.en similar processes were carried out for the particles of other grades as shown in Figure 5.To distribute the di erent grades particles correctly, the particles with largest sizes should be generated rst and then the smaller ones.For example, start the particle generation of 16-19 mm rstly, and then followed by 13.2-16 mm, 9.5-13.2mm, and 4.75-9.5 mm in turns.

Asphalt Mastic Modeling.
e 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. e 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 as k right and k left to record the position properties.ey are all set zero initially.en check all the coordinate values from left to right and bottom to top in turns.e 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.e control lines were not real lines plotted in the images but a designed routine coded in MATLAB to help the distance judgments.en the external and internal area of the particles can be distinguished by the combination of k right and k left .If k right and k left are equal to 1 and 0, respectively, the position is between the particle clump and the left wall as shown in Figure 6.en the coordinate value is effective, and a mastic ball with the diameter of 0.001 m is generated here.
If k right and k left are equal to 0 and 1, respectively, the position is between the particle clump and the right wall as shown in Figure 7. en the coordinate value is effective, and a mastic ball with the diameter of 0.001 m is generated here as well.
If k right and k left 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 k right and k left .e 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.us, 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 Figure 9 due to the space between coarse aggregates.k right and k left 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.

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.e mastic balls would move continuously until reach the equilibrium state.However, due to the initial overlaps among balls (Figure 10), 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.e 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.
e 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 Figure 12(d).

Example of Generating a Virtual Specimen.
e asphalt mixture of AC-16 was developed for simulation.e gradation of the AC-16 is shown in Table 2. 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.
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; and h is the height of the specimen, cm.When it comes to the two-dimensional specimen, the mass can be converted as following: where m is the converted mass of 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 ), and the others are the same as those of three dimensions.e 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  6 Advances in Civil Engineering value of 4%.Moreover, the densities of the coarse aggregates and asphalt mastic were 2.7 g/cm 3 and 2.0 g/cm 3 , 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 userdefined routine within PFC2D.e results were summarized in Tables 4 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.e modeling error decreased as the size decreased.
is 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.e sole difference between the realistic and virtual mastic content was caused by the rounding error in mapping area calculations.

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

Simulations for Dynamic Modulus Prediction.
e virtual asphalt mixtures were developed as shown in Figure 12(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.
e 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 k n and shear stiffness k s .e 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.en 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 k n , shear stiffness k s , 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 di erent contacts were shown in Tables 7  and 8.
Figure 13 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.e 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 e ective 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.
is 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 di erence 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 veri ed the correctness of the proposed algorithm.

Conclusions
is paper proposed an algorithm for the reconstitution of the virtual asphalt mixture based on the DEM methods.By lling 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 veri ed by the virtual uniaxial compressive test.
e 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) e virtual asphalt mixtures were rebuilt in three steps.Firstly, the virtual coarse aggregates were distributed randomly to form the skeletons.en an innovate algorithm was proposed to distinguish the external and internal area of the coarse particles.At last, by lling the designed area with uniform arranged balls and deleting mastic balls randomly, the asphalt mastic and void could be generated successfully.
e results shows that the proposed algorithm can well model the inner components consistent to the actual gradations, asphalt content, and voids.
(3) e validation of the developed models was veri ed by the virtual uniaxial compressive test.e error of    Advances in Civil Engineering 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.

Figure 4 :
Figure 4: Virtual size reconstitutions based on standard particles.

( 1 )
If there is a nearest particle ball in the right side, convert k right from 0 to 1, else, k right is still equal to 0.(2) If there is a nearest particle ball in left side, convert k left from 0 to 1, else, k left is still equal to 0.

Figure 9 :
Figure 9: Distinguishing the position of mastic balls without any entities nearby in two sides.

Figure 10 :
Figure 10: Filling the irregular area with balls to model the asphalt mastic.

Figure 13 :
Figure 13: Dynamic modulus prediction based on the proposed models.

Figure 14 :
Figure 14: Phase angle prediction based on the proposed models.

Table 1 :
Size adjustments for simulations.

Table 3 :
Mapping area calculations for three-phase structures.

Table 4 :
Summaries of the coarse aggregates.

Table 5 :
Summaries of the asphalt mastic and voids.
A: required mapping area for asphalt mastic (m 2 ); B: generated mapping area for asphalt mastic (m 2 ); 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.

Table 6 :
Microparameters of selected constitutive models for particles and walls.