Effects of Aggregate Mesostructure on Permanent Deformation of Asphalt Mixture Using Three-Dimensional Discrete Element Modeling

This paper investigated the effects of aggregate mesostructures on permanent deformation behavior of an asphalt mixture using the three-dimensional (3D) discrete element method (DEM). A 3D discrete element (DE) model of an asphalt mixture composed of coarse aggregates, asphalt mastic, and air voids was developed. Mesomechanical models representing the interactions among the components of asphalt mixture were assigned. Based on the mesomechanical modeling, the uniaxial static load creep tests were simulated using the prepared models, and effects of aggregate angularity, orientation, surface texture, and distribution on the permanent deformation behavior of the asphalt mixtures were analyzed. It was proven that good aggregate angularity had a positive effect on the permanent deformation performance of the asphalt mixtures, especially when approximate cubic aggregates were used. Aggregate packing was more stable when the aggregate orientations tended to be horizontal, which improved the permanent deformation performance of the asphalt mixture. The influence of orientations of 4.75 mm size aggregates on the permanent deformation behavior of the asphalt mixture was significant. Use of aggregates with good surface texture benefitted the permanent deformation performance of the asphalt mixture. Additionally, the non-uniform distribution of aggregates had a negative impact on the permanent deformation performance of the asphalt mixtures, especially when aggregates were distributed non-uniformly in the vertical direction.


Introduction
With the increase of traffic and heavy vehicles, rutting has become one of the main distresses of asphalt pavement, which seriously affects driving safety and comfort. The permanent deformation behavior of asphalt mixtures has become an important issue for road researchers. Within asphalt mixtures, aggregates account for more than 90% by weight and 80% by volume. Aggregates, as the main component of asphalt mixtures, have a significant impact on the permanent deformation behavior of asphalt mixtures. Previous studies on the effects of aggregates on the permanent deformation behavior of asphalt mixtures have been mainly related to aggregate gradations. It is difficult to investigate the effects of the mesostructural characteristics of aggregates on permanent deformation behavior of asphalt mixtures, due to the limitations of available laboratory tests.
In previous studies, X-ray computed tomography (CT) techniques and image processing technology have been used for aggregate mesostructural analysis. Masad proposed computer-automated image analysis procedures to quantify the internal structure of asphalt concrete in terms of aggregate orientation, aggregate contacts, and air void distribution, and developed a finite-element model of

Materials
According to the Marshall mix design method based on the Chinese specifications [26], a dense-graded asphalt mixture, AC-20, was prepared in the laboratory. Aggregate gradation with a nominal maximum aggregate size of 19 mm is shown in Figure 1. Based on the Marshall mix design method, the volumetric properties of the AC-20 asphalt mixture were determined as shown in Table 1. Aggregates smaller than 2.36 mm, mineral filler, and asphalt binder were mixed as asphalt mastic with an asphalt content of 11.5%. The air void content of the asphalt mastic was supposed to be zero to maximize its flowability.

Materials
According to the Marshall mix design method based on the Chinese specifications [26], a densegraded asphalt mixture, AC-20, was prepared in the laboratory. Aggregate gradation with a nominal maximum aggregate size of 19 mm is shown in Figure 1. Based on the Marshall mix design method, the volumetric properties of the AC-20 asphalt mixture were determined as shown in Table 1. Aggregates smaller than 2.36 mm, mineral filler, and asphalt binder were mixed as asphalt mastic with an asphalt content of 11.5%. The air void content of the asphalt mastic was supposed to be zero to maximize its flowability. .

Laboratory Tests
The uniaxial static creep test is one of the simplest and most practical test methods, and it is an effective way to investigate the permanent deformation behaviors of asphalt mixtures and mastics at high temperatures. In this test, an instantaneous load is applied to a cylinder specimen in the axial direction and the load is kept constant, and the creep curve of the asphalt mixtures or mastics can be obtained. Uniaxial static creep tests were conducted in this paper to evaluate the permanent deformation behaviors of asphalt mixtures and mastics at the temperature of 60 °C using a universal test machine (UTM, IPC Global, Melbourne, Australia) [27]. Cylindrical specimens of asphalt mixtures with diameter of 100 mm and height of 150 mm were prepared using a gyratory compactor for the uniaxial static creep test. Cylindrical asphalt mastic specimens with diameter and height of 100 mm were prepared by vibration. The applied axial stresses for the asphalt mixture and asphalt mastic were set at 0.7 MPa and 0.07 MPa, respectively, in this study. The asphalt mastic specimens and the testing process are shown in Figure 2. Lab test results were used for parameters of DE modeling and to be compared with simulation test results, and are presented in the following sections. The experimental program for the asphalt mixture and mastics is summarized in Table 2.

Laboratory Tests
The uniaxial static creep test is one of the simplest and most practical test methods, and it is an effective way to investigate the permanent deformation behaviors of asphalt mixtures and mastics at high temperatures. In this test, an instantaneous load is applied to a cylinder specimen in the axial direction and the load is kept constant, and the creep curve of the asphalt mixtures or mastics can be obtained. Uniaxial static creep tests were conducted in this paper to evaluate the permanent deformation behaviors of asphalt mixtures and mastics at the temperature of 60 • C using a universal test machine (UTM, IPC Global, Melbourne, Australia) [27]. Cylindrical specimens of asphalt mixtures with diameter of 100 mm and height of 150 mm were prepared using a gyratory compactor for the uniaxial static creep test. Cylindrical asphalt mastic specimens with diameter and height of 100 mm were prepared by vibration. The applied axial stresses for the asphalt mixture and asphalt mastic were set at 0.7 MPa and 0.07 MPa, respectively, in this study. The asphalt mastic specimens and the testing process are shown in Figure 2. Lab test results were used for parameters of DE modeling and to be compared with simulation test results, and are presented in the following sections. The experimental program for the asphalt mixture and mastics is summarized in Table 2.

Discrete Element Modeling
An asphalt mixture, as a multiphase composite, is composed of aggregates, asphalt binder, and air voids. If the particle size of the fine aggregates and mineral filler is small, it will lead to a large increase of discrete elements in the DE model, which is bound to significantly reduce the computational efficiency if the fine aggregate and mineral powder are considered fully in DE models of asphalt mixtures. Therefore, the asphalt mixture was simplified into coarse aggregates (bigger than 2.36 mm), asphalt mastics (asphalt binder, fine aggregates smaller than 2.36 mm and mineral filler), and air voids to improve the calculation efficiency. The asphalt mastic was considered as a homogeneous material within the DE model.
During modeling, the spatial range of the mixture model was first constructed by "wall". The quantity of coarse aggregates in each grade was calculated according to aggregate gradation, asphalt content and air void content. The coarse aggregate balls were then delivered into the spatial range of the mixture model constructed by "wall". Force was generated between overlapped balls, as there were overlaps between the coarse aggregate balls during delivery. This force was eliminated by the "cycle" command. The coarse aggregate balls within the spatial range of the mixture model are shown in Figure 3a.

Discrete Element Modeling
An asphalt mixture, as a multiphase composite, is composed of aggregates, asphalt binder, and air voids. If the particle size of the fine aggregates and mineral filler is small, it will lead to a large increase of discrete elements in the DE model, which is bound to significantly reduce the computational efficiency if the fine aggregate and mineral powder are considered fully in DE models of asphalt mixtures. Therefore, the asphalt mixture was simplified into coarse aggregates (bigger than 2.36 mm), asphalt mastics (asphalt binder, fine aggregates smaller than 2.36 mm and mineral filler), and air voids to improve the calculation efficiency. The asphalt mastic was considered as a homogeneous material within the DE model.
During modeling, the spatial range of the mixture model was first constructed by "wall". The quantity of coarse aggregates in each grade was calculated according to aggregate gradation, asphalt content and air void content. The coarse aggregate balls were then delivered into the spatial range of the mixture model constructed by "wall". Force was generated between overlapped balls, as there were overlaps between the coarse aggregate balls during delivery. This force was eliminated by the "cycle" command. The coarse aggregate balls within the spatial range of the mixture model are shown in Figure 3a.
A regular array of discrete elements was then filled into the mixture space as the base of coarse aggregates and asphalt mastic, as shown in Figure 3b. The radius of the discrete element was set to 1 mm, considering the calculation efficiency. quantity of coarse aggregates in each grade was calculated according to aggregate gradation, asphalt content and air void content. The coarse aggregate balls were then delivered into the spatial range of the mixture model constructed by "wall". Force was generated between overlapped balls, as there were overlaps between the coarse aggregate balls during delivery. This force was eliminated by the "cycle" command. The coarse aggregate balls within the spatial range of the mixture model are shown in Figure 3a. It has been proven that the geometry of coarse aggregates has a significant effect on the permanent deformation behavior of asphalt mixtures [28,29]. Therefore, the coarse aggregate geometry will directly determine the accuracy of the simulation results. In this paper, coarse aggregate was simplified as an irregular polyhedron. Random planes were used to cut a cube or a sphere to generate irregular polyhedral aggregates using a user-defined program. By traversing the regularly packed discrete elements, the positional relationship between the regular packing discrete elements and the irregular polyhedral aggregates could be evaluated. Discrete elements belonging to irregular polyhedral aggregate were regarded as aggregate elements, and were set as a clump, as seen in Figure 4. The original coarse aggregate balls were then deleted. Discrete elements outside the irregular polyhedron aggregate were considered as asphalt mastic. The initially developed 3D DE model of asphalt mixture is shown in Figure 5a.  A regular array of discrete elements was then filled into the mixture space as the base of coarse aggregates and asphalt mastic, as shown in Figure 3b. The radius of the discrete element was set to 1 mm, considering the calculation efficiency.
It has been proven that the geometry of coarse aggregates has a significant effect on the permanent deformation behavior of asphalt mixtures [28,29]. Therefore, the coarse aggregate geometry will directly determine the accuracy of the simulation results. In this paper, coarse aggregate was simplified as an irregular polyhedron. Random planes were used to cut a cube or a sphere to generate irregular polyhedral aggregates using a user-defined program. By traversing the regularly packed discrete elements, the positional relationship between the regular packing discrete elements and the irregular polyhedral aggregates could be evaluated. Discrete elements belonging to irregular polyhedral aggregate were regarded as aggregate elements, and were set as a clump, as seen in Figure 4. The original coarse aggregate balls were then deleted. Discrete elements outside the irregular polyhedron aggregate were considered as asphalt mastic. The initially developed 3D DE model of asphalt mixture is shown in Figure 5a. Considering the complexity of air void distribution within an asphalt mixture specimen, the air voids were assumed to be randomly distributed [30]. Asphalt mastic elements in the model were traversed and randomly deleted and regarded as air voids. The distribution of air voids is shown in Figure 5b.  A regular array of discrete elements was then filled into the mixture space as the base of coarse aggregates and asphalt mastic, as shown in Figure 3b. The radius of the discrete element was set to 1 mm, considering the calculation efficiency.
It has been proven that the geometry of coarse aggregates has a significant effect on the permanent deformation behavior of asphalt mixtures [28,29]. Therefore, the coarse aggregate geometry will directly determine the accuracy of the simulation results. In this paper, coarse aggregate was simplified as an irregular polyhedron. Random planes were used to cut a cube or a sphere to generate irregular polyhedral aggregates using a user-defined program. By traversing the regularly packed discrete elements, the positional relationship between the regular packing discrete elements and the irregular polyhedral aggregates could be evaluated. Discrete elements belonging to irregular polyhedral aggregate were regarded as aggregate elements, and were set as a clump, as seen in Figure 4. The original coarse aggregate balls were then deleted. Discrete elements outside the irregular polyhedron aggregate were considered as asphalt mastic. The initially developed 3D DE model of asphalt mixture is shown in Figure 5a. Considering the complexity of air void distribution within an asphalt mixture specimen, the air voids were assumed to be randomly distributed [30]. Asphalt mastic elements in the model were traversed and randomly deleted and regarded as air voids. The distribution of air voids is shown in Figure 5b. Considering the complexity of air void distribution within an asphalt mixture specimen, the air voids were assumed to be randomly distributed [30]. Asphalt mastic elements in the model were traversed and randomly deleted and regarded as air voids. The distribution of air voids is shown in Figure 5b.

Mesomechanical Models and Parameters
There are four types of contact within asphalt mixtures, including contacts between aggregate elements, contacts between adjacent aggregates, contacts between asphalt mastic elements, and contact between asphalt mastic and aggregates, as seen in Figure 6. In PFC3D, mesomechanical models are used to describe the contact behaviors between different components within asphalt mixtures. In this study, the mesomechanical models used in the model of asphalt mixtures included the stiffness model, slipping model, bonding model, and Burger's model. Due to the high stiffness of an aggregate, it can be approximately regarded as an elastic material. In this study, the stiffness model and slipping model were used to characterize the mesomechanical behavior between adjacent aggregates. As coarse aggregates within the model of asphalt mixtures were set as clumps, it was unnecessary to assign the mesomechanical model within aggregates. The mesomechanical parameters of the stiffness model could be obtained from the macro properties of the aggregates, as shown in Equations (1) and (2) [15,31,32]. The macro parameters for the aggregates are shown in Table 3 [9,12,33].
where E is the apparent Young's modulus of the aggregates, kn and ks are the stiffness in the normal and shear direction, respectively, R is the discrete element radius, υ' is the aggregate Poisson's ratio, μc is the friction coefficient between aggregates, and μa is the friction coefficient of aggregates. Asphalt mixtures show a macro viscoelastic behavior due to the viscoelastic characteristic of the asphalt mastic. Therefore, the mesomechanical model of the asphalt mastic directly affects the macro properties of asphalt mixtures. The meso Burger's model in PFC3D was well able to describe the mechanical properties of viscoelastic materials and was used to characterize the viscoelastic properties of the asphalt mastic in this study, as shown in Figure 7. It has been proven that there is a conversion between the parameters of the meso Burger's model and the macro Burger's model [15,31], as shown in Equations (3) and (4). Parameters of meso Burger's model could then be obtained from the macro Burger's model parameters. The macro Burger's model parameters of asphalt mastic were obtained by uniaxial static creep test, as shown in Table 3. Due to the high stiffness of an aggregate, it can be approximately regarded as an elastic material. In this study, the stiffness model and slipping model were used to characterize the mesomechanical behavior between adjacent aggregates. As coarse aggregates within the model of asphalt mixtures were set as clumps, it was unnecessary to assign the mesomechanical model within aggregates. The mesomechanical parameters of the stiffness model could be obtained from the macro properties of the aggregates, as shown in Equations (1) and (2) [15,31,32]. The macro parameters for the aggregates are shown in Table 3 [9,12,33].
where E is the apparent Young's modulus of the aggregates, k n and k s are the stiffness in the normal and shear direction, respectively, R is the discrete element radius, υ is the aggregate Poisson's ratio, µ c is the friction coefficient between aggregates, and µ a is the friction coefficient of aggregates. Table 3. Macro parameters for asphalt mastic and aggregates. properties of asphalt mixtures. The meso Burger's model in PFC3D was well able to describe the mechanical properties of viscoelastic materials and was used to characterize the viscoelastic properties of the asphalt mastic in this study, as shown in Figure 7. It has been proven that there is a conversion between the parameters of the meso Burger's model and the macro Burger's model [15,31], as shown in Equations (3) and (4). Parameters of meso Burger's model could then be obtained from the macro Burger's model parameters. The macro Burger's model parameters of asphalt mastic were obtained by uniaxial static creep test, as shown in Table 3.
where E 1 , η 1 , E 2 , and η 2 , are parameters of the macro Burger's model; K mn , C mn , K kn , and C kn are parameters of the meso Burger's model in the normal direction; K ms , C ms , K ks , and C ks are parameters of the meso Burger's model in the shear direction; L is the distance between adjacent discrete elements; and υ is Poisson's ratio of the asphalt mastic.
where E1, η1, E2, and η2, are parameters of the macro Burger's model; Kmn, Cmn, Kkn, and Ckn are parameters of the meso Burger's model in the normal direction; Kms, Cms, Kks, and Cks are parameters of the meso Burger's model in the shear direction; L is the distance between adjacent discrete elements; and υ is Poisson's ratio of the asphalt mastic.
The contact behavior between aggregates and the asphalt mastic can be characterized by the equivalent meso Burger's model, as shown in Figure 8. The equivalent meso model parameters can be obtained from the macro parameters of the aggregates and asphalt mastic, as expressed in Equations (5) and (6) [15,31]: where K'mm, C'mm, K'kn, and C'kn are the parameters of the equivalent meso Burger's model between the The contact behavior between aggregates and the asphalt mastic can be characterized by the equivalent meso Burger's model, as shown in Figure 8. The equivalent meso model parameters can be obtained from the macro parameters of the aggregates and asphalt mastic, as expressed in Equations (5) and (6) [15,31]: where K mm , C mm , K kn , and C kn are the parameters of the equivalent meso Burger's model between the aggregates and mastic in the normal direction and K ms , C ms , K ks , and C ks are the parameters of the equivalent meso Burger's model between the aggregates and the mastic in the shear direction.
where K'mm, C'mm, K'kn, and C'kn are the parameters of the equivalent meso Burger's model between the aggregates and mastic in the normal direction and K'ms, C'ms, K'ks, and C'ks are the parameters of the equivalent meso Burger's model between the aggregates and the mastic in the shear direction.

Simulation of Uniaxial Creep Test
Based on the above-mentioned DE modeling procedures for asphalt mixtures, the simulated uniaxial creep test was carried out using PFC3D. The results were compared with laboratory test results under the same conditions. Figure 9 shows the axial strain curve obtained from the simulation and lab test, and it can be seen that the simulation curve was very close to the curve of the laboratory test. This indicates that the simulation of uniaxial creep test conducted using DEM could precisely estimate the permanent deformation behavior of asphalt mixtures.

Simulation of Uniaxial Creep Test
Based on the above-mentioned DE modeling procedures for asphalt mixtures, the simulated uniaxial creep test was carried out using PFC3D. The results were compared with laboratory test results under the same conditions. Figure 9 shows the axial strain curve obtained from the simulation and lab test, and it can be seen that the simulation curve was very close to the curve of the laboratory test. This indicates that the simulation of uniaxial creep test conducted using DEM could precisely estimate the permanent deformation behavior of asphalt mixtures.

Effect of Aggregate Angularity
The aggregate angularity could represent the morphological characteristics of aggregates, and significantly affect the skeleton stability, interlocking force, and shear resistance of aggregates. It has been proven that the aggregate angularity is closely related to the strength of asphalt mixtures, especially the permanent deformation resistance and shear resistance [28,29].
Some previous studies about aggregate angularity have been carried out, and some indexes to describe the aggregate angularity have been proposed. These angularity indexes are mainly divided into two categories, one denoting the roundness of aggregate corners and the other one denoting the roundness of the overall outline of the aggregate [34]. However, most of these angularity indexes are in two dimensions. Moreover, the measurement and calculation of the first type of angularity index is very complex, while the second type could be affected by the overall outline of aggregates. Therefore, the angularity indexes mentioned above cannot fully represent the angularity characteristics of aggregates. In this study, the ratio of surface area of irregular aggregate to equivalent ellipsoid was proposed to characterize the aggregate angularity. The equivalent ellipsoid of an aggregate has the same volume as the aggregate. It was considered that the equivalent ellipsoid of aggregates could accurately reflect the overall outline, and this angularity index could significantly reduce the influence of aggregate outline on angularity index, as expressed in Equation (7). The larger AI is, the better the aggregate angularity is.

Effect of Aggregate Angularity
The aggregate angularity could represent the morphological characteristics of aggregates, and significantly affect the skeleton stability, interlocking force, and shear resistance of aggregates. It has been proven that the aggregate angularity is closely related to the strength of asphalt mixtures, especially the permanent deformation resistance and shear resistance [28,29].
Some previous studies about aggregate angularity have been carried out, and some indexes to describe the aggregate angularity have been proposed. These angularity indexes are mainly divided into two categories, one denoting the roundness of aggregate corners and the other one denoting the roundness of the overall outline of the aggregate [34]. However, most of these angularity indexes are in two dimensions. Moreover, the measurement and calculation of the first type of angularity index is very complex, while the second type could be affected by the overall outline of aggregates. Therefore, the angularity indexes mentioned above cannot fully represent the angularity characteristics of aggregates. In this study, the ratio of surface area of irregular aggregate to equivalent ellipsoid was proposed to characterize the aggregate angularity. The equivalent ellipsoid of an aggregate has the same volume as the aggregate. It was considered that the equivalent ellipsoid of aggregates could accurately reflect the overall outline, and this angularity index could significantly reduce the influence of aggregate outline on angularity index, as expressed in Equation (7). The larger AI is, the better the aggregate angularity is.
where AI is the angularity index of aggregates, S is surface area of aggregates, and S ellipse is surface area of equivalent ellipsoid.
As mentioned above, irregular polyhedral aggregates were generated by cutting a cube or a sphere with random planes. Thus, the equivalent ellipsoid of an irregular polyhedral aggregate can be approximated to an equivalent sphere, and then Equation (8) can be further expressed as follows: where S sphere is surface area of equivalent sphere and V is the volume of aggregates.
Coarse aggregates with different angularities were generated by cutting a cube or a sphere with random planes. Some coarse aggregates with different angularities are shown in Figure 10, in which their angularity become better with the increasing number. DE models of asphalt mixtures with different aggregate angularities were then obtained. where Ssphere is surface area of equivalent sphere and V is the volume of aggregates.
Coarse aggregates with different angularities were generated by cutting a cube or a sphere with random planes. Some coarse aggregates with different angularities are shown in Figure 10 In order to calculate the surface area of aggregates, the number of discrete elements on aggregate surface and the number of discrete elements that made up the aggregates were counted. The angularity index of the aggregates was then calculated by Equation (9). The method used to count the number of discrete elements on the aggregate surface and the number of discrete elements that make up aggregates is shown in Figure 11.
where ns is the number of discrete elements on the aggregate surface, nv is the number of discrete elements that make up the aggregate, and R is radius of the discrete element. The angularity of a single aggregate can be obtained by Equation (9). However, there are many aggregates within asphalt mixtures, and the angularities of each aggregate are different. It is necessary to evaluate the overall angularity of the aggregates within asphalt mixtures. In this study, the volumetrically weighted average of the angularities of each aggregate was used to describe the overall angularity of the aggregates within the asphalt mixtures, as expressed in Equation (10).
where is the overall angularity of aggregates, AIi is the angularity of aggregate i, Vi is the volume of aggregate i, and n is the number of aggregates. In order to calculate the surface area of aggregates, the number of discrete elements on aggregate surface and the number of discrete elements that made up the aggregates were counted. The angularity index of the aggregates was then calculated by Equation (9). The method used to count the number of discrete elements on the aggregate surface and the number of discrete elements that make up aggregates is shown in Figure 11. Simulations of uniaxial static creep tests were then conducted on the DE models of asphalt mixtures with different overall aggregate angularities. The results are shown in Figure 12. The permanent deformation performance of asphalt mixtures improved with the increasing aggregate angularity, but the improvement weakened gradually. The axial deformation of the asphalt mixture was the smallest when the aggregates were closest to a cube. This is consistent with the requirement of selecting aggregates with shapes close to a cube. Therefore, aggregates with good angularities should be selected to improve the permanent deformation resistance of asphalt pavement. The angularity of a single aggregate can be obtained by Equation (9). However, there are many aggregates within asphalt mixtures, and the angularities of each aggregate are different. It is necessary to evaluate the overall angularity of the aggregates within asphalt mixtures. In this study, the volumetrically weighted average of the angularities of each aggregate was used to describe the overall angularity of the aggregates within the asphalt mixtures, as expressed in Equation (10).
where AI is the overall angularity of aggregates, AI i is the angularity of aggregate i, V i is the volume of aggregate i, and n is the number of aggregates.
Simulations of uniaxial static creep tests were then conducted on the DE models of asphalt mixtures with different overall aggregate angularities. The results are shown in Figure 12. The permanent deformation performance of asphalt mixtures improved with the increasing aggregate angularity, but the improvement weakened gradually. The axial deformation of the asphalt mixture was the smallest when the aggregates were closest to a cube. This is consistent with the requirement of selecting aggregates with shapes close to a cube. Therefore, aggregates with good angularities should be selected to improve the permanent deformation resistance of asphalt pavement.

Effect of Aggregate Orientation
Previous studies about aggregate orientation have generally defined the long axis of the equivalent ellipse of aggregate in two dimensions as the long axis of aggregate, and used the angle between the long axis and the horizontal direction to represent the aggregate orientation [18][19][20][21]. In this study, the aggregate orientation was defined in three dimensions. The long axis of the equivalent ellipsoid of the aggregate was taken as the long axis, and the angles α, β, and γ between the long axis and X axis, Y axis, and Z axis were used to represent the aggregate orientation in three dimensions, as shown in Figure 13.
In order to investigate the effect of aggregate orientation on the permanent deformation performance of asphalt mixtures, a subroutine was written in PFC3D with the built-in "Fish" language. The long axis of the aggregate was set to be parallel to the YOZ surface, and the aggregate orientation was described by the angle between the long axis of the aggregate and the XOY surface. Some coarse aggregates with different orientations are shown in Figure 14.

Effect of Aggregate Orientation
Previous studies about aggregate orientation have generally defined the long axis of the equivalent ellipse of aggregate in two dimensions as the long axis of aggregate, and used the angle between the long axis and the horizontal direction to represent the aggregate orientation [18][19][20][21]. In this study, the aggregate orientation was defined in three dimensions. The long axis of the equivalent ellipsoid of the aggregate was taken as the long axis, and the angles α, β, and γ between the long axis and X axis, Y axis, and Z axis were used to represent the aggregate orientation in three dimensions, as shown in Figure 13.

Effect of Aggregate Orientation
Previous studies about aggregate orientation have generally defined the long axis of the equivalent ellipse of aggregate in two dimensions as the long axis of aggregate, and used the angle between the long axis and the horizontal direction to represent the aggregate orientation [18][19][20][21]. In this study, the aggregate orientation was defined in three dimensions. The long axis of the equivalent ellipsoid of the aggregate was taken as the long axis, and the angles α, β, and γ between the long axis and X axis, Y axis, and Z axis were used to represent the aggregate orientation in three dimensions, as shown in Figure 13.
In order to investigate the effect of aggregate orientation on the permanent deformation performance of asphalt mixtures, a subroutine was written in PFC3D with the built-in "Fish" language. The long axis of the aggregate was set to be parallel to the YOZ surface, and the aggregate orientation was described by the angle between the long axis of the aggregate and the XOY surface. Some coarse aggregates with different orientations are shown in Figure 14. In order to investigate the effect of aggregate orientation on the permanent deformation performance of asphalt mixtures, a subroutine was written in PFC3D with the built-in "Fish" language. The long axis of the aggregate was set to be parallel to the YOZ surface, and the aggregate orientation was described by the angle between the long axis of the aggregate and the XOY surface. Some coarse aggregates with different orientations are shown in Figure 14.

Effect of Aggregate Orientation
Previous studies about aggregate orientation have generally defined the long axis of the equivalent ellipse of aggregate in two dimensions as the long axis of aggregate, and used the angle between the long axis and the horizontal direction to represent the aggregate orientation [18][19][20][21]. In this study, the aggregate orientation was defined in three dimensions. The long axis of the equivalent ellipsoid of the aggregate was taken as the long axis, and the angles α, β, and γ between the long axis and X axis, Y axis, and Z axis were used to represent the aggregate orientation in three dimensions, as shown in Figure 13.
In order to investigate the effect of aggregate orientation on the permanent deformation performance of asphalt mixtures, a subroutine was written in PFC3D with the built-in "Fish" language. The long axis of the aggregate was set to be parallel to the YOZ surface, and the aggregate orientation was described by the angle between the long axis of the aggregate and the XOY surface. Some coarse aggregates with different orientations are shown in Figure 14. Graded coarse aggregates with different orientations were generated and randomly put into the cylinder space of the asphalt mixture model. The unbalanced forces between the aggregates were eliminated by the "cycle" command, so that the coarse aggregates within the asphalt mixtures could reach equilibrium. The angular velocities of the coarse aggregates were fixed to keep the aggregate orientations constant during modeling.
The effects of the orientations of aggregates of different sizes on the mechanical properties of the asphalt mixtures are different. To address the size effect, 30 cylinder DE models of asphalt mixtures with different orientations of aggregates with different sizes were developed. The aggregate orientations were varied from 0 • , 30 • , 45 • , 60 • , to 90 • . The rest of the aggregates remained spherical. Figure 15 shows the coarse aggregates with different orientations. Graded coarse aggregates with different orientations were generated and randomly put into the cylinder space of the asphalt mixture model. The unbalanced forces between the aggregates were eliminated by the "cycle" command, so that the coarse aggregates within the asphalt mixtures could reach equilibrium. The angular velocities of the coarse aggregates were fixed to keep the aggregate orientations constant during modeling.
The effects of the orientations of aggregates of different sizes on the mechanical properties of the asphalt mixtures are different. To address the size effect, 30 cylinder DE models of asphalt mixtures with different orientations of aggregates with different sizes were developed. The aggregate orientations were varied from 0°, 30°, 45°, 60°, to 90°. The rest of the aggregates remained spherical. Figure 15 shows the coarse aggregates with different orientations. Simulations of uniaxial static creep tests were then conducted on the prepared DE models. The results are shown in Figure 16. It can be seen from Figure 16 that the axial strains of the mixtures increased gradually with the increasing aggregate orientations. The orientations of 19-26.5mm and 16-19 mm aggregates had insignificant effects on the permanent deformation behavior of asphalt mixtures. The orientations of 13.2-16 mm, 9.5-13.2 mm, and 2.36-4.75 mm aggregates had more significant effects than those of 19-26.5 mm and 16-19 mm. The orientation of 4.75-9.5 mm aggregates had the greatest effect on the permanent deformation behavior of the asphalt mixtures. The reason is that the aggregates tended to rotate horizontally under vertical loads when the aggregate orientation was large. Thus, the aggregate skeleton was unstable and led to large axial deformation of the asphalt mixtures. When the aggregate orientation was closer to horizontal, the aggregate skeleton was more stable, and the axial deformation of the mixture specimens was smaller. It has been shown that the percentage pass of 4.75 mm plays an important role in controlling the aggregate skeleton, and the influence is next to that of air void content on the permanent deformation behavior of asphalt mixtures. This may be the reason why the orientation of 4.75-9.5 mm aggregates had such a great influence. Therefore, different kinds of rollers should be used during rolling compaction of asphalt pavement, which could force the aggregate orientation to be more horizontal. Thus, the asphalt pavement would reach a more stable state and its permanent deformation performance could be significantly improved.

Effect of Aggregate Surface Texture
The surface texture of aggregates reflects the roughness of aggregate surface. The rougher the aggregate surface, the better the surface texture. Internal friction between aggregates and adhesion between aggregate and asphalt binder are closely related to the surface texture of aggregates. It has been proven that the surface texture of aggregates has a significant impact on the performance of asphalt mixtures, especially the permanent deformation performance. The randomly created irregular polyhedral aggregates in the DE model could only reflect geometry and angularity, without characterizing the surface texture of aggregates. In this study, the friction coefficient between aggregates was used to indirectly represent the surface texture of aggregates. The higher the friction coefficient, the rougher the surface texture. A series of DE models of asphalt mixture specimens with different aggregate friction coefficients assigned were developed. The internal structures of the developed DE models were the same, and only the aggregate friction coefficients of polyhedron aggregates in the DE models were different. Simulations of uniaxial static creep tests were conducted on the prepared DE models. The results are shown in Figure 17. It can be seen from Figure 16 that the axial strains of the mixtures increased gradually with the increasing aggregate orientations. The orientations of 19-26.5mm and 16-19 mm aggregates had insignificant effects on the permanent deformation behavior of asphalt mixtures. The orientations of 13.2-16 mm, 9.5-13.2 mm, and 2.36-4.75 mm aggregates had more significant effects than those of 19-26.5 mm and 16-19 mm. The orientation of 4.75-9.5 mm aggregates had the greatest effect on the permanent deformation behavior of the asphalt mixtures. The reason is that the aggregates tended to rotate horizontally under vertical loads when the aggregate orientation was large. Thus, the aggregate skeleton was unstable and led to large axial deformation of the asphalt mixtures. When the aggregate orientation was closer to horizontal, the aggregate skeleton was more stable, and the axial deformation of the mixture specimens was smaller. It has been shown that the percentage pass of 4.75 mm plays an important role in controlling the aggregate skeleton, and the influence is next to that of air void content on the permanent deformation behavior of asphalt mixtures. This may be the reason why the orientation of 4.75-9.5 mm aggregates had such a great influence. Therefore, different kinds of rollers should be used during rolling compaction of asphalt pavement, which could force the aggregate orientation to be more horizontal. Thus, the asphalt pavement would reach a more stable state and its permanent deformation performance could be significantly improved.

Effect of Aggregate Surface Texture
The surface texture of aggregates reflects the roughness of aggregate surface. The rougher the aggregate surface, the better the surface texture. Internal friction between aggregates and adhesion between aggregate and asphalt binder are closely related to the surface texture of aggregates. It has been proven that the surface texture of aggregates has a significant impact on the performance of asphalt mixtures, especially the permanent deformation performance. The randomly created irregular polyhedral aggregates in the DE model could only reflect geometry and angularity, without characterizing the surface texture of aggregates. In this study, the friction coefficient between aggregates was used to indirectly represent the surface texture of aggregates. The higher the friction coefficient, the rougher the surface texture. A series of DE models of asphalt mixture specimens with different aggregate friction coefficients assigned were developed. The internal structures of the developed DE models were the same, and only the aggregate friction coefficients of polyhedron aggregates in the DE models were different. Simulations of uniaxial static creep tests were conducted on the prepared DE models. The results are shown in Figure 17. It can be seen from Figure 17 that the permanent deformation performance of the asphalt mixtures was improved with the increasing friction coefficient, but the increment was gradually lowered. The permanent deformation resistance of the asphalt mixtures increased with rougher aggregate surface texture. Therefore, aggregates with good surface texture should be used to improve the permanent deformation performance of asphalt pavement.

Effect of Aggregate Distribution
Aggregates are ideally uniformly distributed within asphalt mixture for optimum performance. However, insufficient mixing force or mixing time may cause non-uniform distribution of aggregates. Aggregate segregation may also occur during paving. Early distresses of asphalt pavement such as cracking, pothole, blooding, and rutting are closely related to the material inhomogeneity caused by variable construction quality. The non-uniform distribution of aggregates greatly influences the internal structure of asphalt mixtures and consequently lowers the permanent deformation performance of asphalt mixtures.
To evaluate the uniformity of aggregate distribution, a cubic space containing asphalt mixtures was constructed by "wall" and the specimen was divided into three parts evenly in both the horizontal and vertical directions, as shown in Figure 18. Coarse aggregates with different sizes were added to each of the three parts, as shown in Figure 19 to create non-uniform distribution. The total volume of aggregates within the three parts remained constant. The distribution percentages of the aggregate volume within each part are shown in Table 4. It can be seen from Figure 17 that the permanent deformation performance of the asphalt mixtures was improved with the increasing friction coefficient, but the increment was gradually lowered. The permanent deformation resistance of the asphalt mixtures increased with rougher aggregate surface texture. Therefore, aggregates with good surface texture should be used to improve the permanent deformation performance of asphalt pavement.

Effect of Aggregate Distribution
Aggregates are ideally uniformly distributed within asphalt mixture for optimum performance. However, insufficient mixing force or mixing time may cause non-uniform distribution of aggregates. Aggregate segregation may also occur during paving. Early distresses of asphalt pavement such as cracking, pothole, blooding, and rutting are closely related to the material inhomogeneity caused by variable construction quality. The non-uniform distribution of aggregates greatly influences the internal structure of asphalt mixtures and consequently lowers the permanent deformation performance of asphalt mixtures.
To evaluate the uniformity of aggregate distribution, a cubic space containing asphalt mixtures was constructed by "wall" and the specimen was divided into three parts evenly in both the horizontal and vertical directions, as shown in Figure 18. Coarse aggregates with different sizes were added to each of the three parts, as shown in Figure 19 to create non-uniform distribution. The total volume of aggregates within the three parts remained constant. The distribution percentages of the aggregate volume within each part are shown in Table 4. It can be seen from Figure 17 that the permanent deformation performance of the asphalt mixtures was improved with the increasing friction coefficient, but the increment was gradually lowered. The permanent deformation resistance of the asphalt mixtures increased with rougher aggregate surface texture. Therefore, aggregates with good surface texture should be used to improve the permanent deformation performance of asphalt pavement.

Effect of Aggregate Distribution
Aggregates are ideally uniformly distributed within asphalt mixture for optimum performance. However, insufficient mixing force or mixing time may cause non-uniform distribution of aggregates. Aggregate segregation may also occur during paving. Early distresses of asphalt pavement such as cracking, pothole, blooding, and rutting are closely related to the material inhomogeneity caused by variable construction quality. The non-uniform distribution of aggregates greatly influences the internal structure of asphalt mixtures and consequently lowers the permanent deformation performance of asphalt mixtures.
To evaluate the uniformity of aggregate distribution, a cubic space containing asphalt mixtures was constructed by "wall" and the specimen was divided into three parts evenly in both the horizontal and vertical directions, as shown in Figure 18. Coarse aggregates with different sizes were added to each of the three parts, as shown in Figure 19 to create non-uniform distribution. The total volume of aggregates within the three parts remained constant. The distribution percentages of the aggregate volume within each part are shown in Table 4.  A quantitative indicator of aggregate distribution was proposed to measure the uniformity of distribution of aggregate. The volume of aggregates within each part was recorded and the coefficient of variation (var) of aggregate distribution was calculated as follows:  A quantitative indicator of aggregate distribution was proposed to measure the uniformity of distribution of aggregate. The volume of aggregates within each part was recorded and the coefficient of variation (var) of aggregate distribution was calculated as follows: var j (14) where V i is the volume of aggregate i within part j, V j is the total volume of aggregates within part j, V is the average aggregate volume of all divided parts, var j is the coefficient of variation of the total volume of aggregates within part j, and var is the coefficient of variation of the total volume of aggregates within specimen. Figure 20 and Table 5 below present the results, in which the axial strain was calculated by measuring the axial deformation and dividing it by the initial axial length.
where Vi is the volume of aggregate i within part j, Vj is the total volume of aggregates within part j, is the average aggregate volume of all divided parts, varj is the coefficient of variation of the total volume of aggregates within part j, and var is the coefficient of variation of the total volume of aggregates within specimen. Figure 20 and Table 5 below present the results, in which the axial strain was calculated by measuring the axial deformation and dividing it by the initial axial length.  As shown in Figure 20 and Table 5, axial deformation of the specimens increased as the coefficient of variation of aggregate distribution increased using both partition methods. Aggregate distribution had a significant impact on the deformation performance of the asphalt mixtures. Nonuniformity in the vertical direction had much more impact on the deformation of the specimen than that in the horizontal direction, as shown in Figure 20.

Conclusions
In this study, three-dimensional discrete element (DE) models of asphalt mixtures composed of coarse aggregates, asphalt mastic, and air voids was developed using PFC3D, and uniaxial static load creep tests were conducted on the prepared models to investigate the effects of aggregate mesostructure on the permanent deformation behavior of the asphalt mixtures. The effects of aggregate angularity, orientation, surface texture, and distribution on the permanent deformation behavior of asphalt mixtures were carefully analyzed. The following conclusions were drawn:  As shown in Figure 20 and Table 5, axial deformation of the specimens increased as the coefficient of variation of aggregate distribution increased using both partition methods. Aggregate distribution had a significant impact on the deformation performance of the asphalt mixtures. Non-uniformity in the vertical direction had much more impact on the deformation of the specimen than that in the horizontal direction, as shown in Figure 20.

Conclusions
In this study, three-dimensional discrete element (DE) models of asphalt mixtures composed of coarse aggregates, asphalt mastic, and air voids was developed using PFC3D, and uniaxial static load creep tests were conducted on the prepared models to investigate the effects of aggregate mesostructure on the permanent deformation behavior of the asphalt mixtures. The effects of aggregate angularity, orientation, surface texture, and distribution on the permanent deformation behavior of asphalt mixtures were carefully analyzed. The following conclusions were drawn: (1) The constructed DE models of asphalt mixtures could accurately capture the mesostructures of aggregates, such as angularity, orientation, surface texture, and distribution. The simulation of the uniaxial creep test conducted using the DEM could precisely estimate the permanent deformation behavior of asphalt mixtures.
(2) The permanent deformation performance of the asphalt mixtures was improved with increasing aggregate angularity. The axial deformation of the asphalt mixtures was smallest when the aggregates were close to cubical. Good aggregate angularity had a positive effect on the permanent deformation performances of asphalt mixtures, especially when the aggregates were nearly cubical.
(3) Aggregate packing was more stable when the orientations tended to be horizontal, which improved the permanent deformation performances of asphalt mixtures. The percentage pass of 4.75 mm played an important role in controlling the aggregate skeleton and had a great impact on the permanent deformation behaviors of the asphalt mixtures.
(4) The permanent deformation performance of the asphalt mixtures was improved with increasing friction coefficient, indicating that aggregates with good surface texture benefitted the permanent deformation performance of the asphalt mixtures.
(5) Axial deformation of the asphalt mixture specimen increased with increasing coefficient of variation of aggregate distribution. Non-uniform distribution of the aggregates had a negative impact on the permanent deformation performance of the asphalt mixtures, especially when the aggregates were distributed non-uniformly in the vertical direction.
In the future research steps, a viscoplastic mesomechanical model will be developed to further describe the viscoplastic behavior of the permanent deformation behavior of asphalt mixtures.

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