STUDY ON CREEP BEHAVIOR OF ASPHALT MIXTURE BASED ON DISCRETE ELEMENT METHOD

A three-dimensional (3D) microstructure-based discrete element (DE) model was developed to study the creep behaviour of high viscoelastic asphalt sand (HVAS) with the uniaxial compression creep tests. The threepoint bending creep tests of asphalt mortar were carried out in order to obtain the parameters of the Burger model, to determine the transformation method of macroscopic parameters and microscopic parameters of the model in theory, to obtain the parameters used in the discrete element model, and then establish the discrete element analysis model for the asphalt mixture. A 3D-DE digital specimen was composed of coarse aggregates, asphalt mortar and air voids, which could also take gradation, irregular shape, random distribution of


Introduction
An asphalt mixture composed of aggregates, asphalt mortar and air voids is characterized by inhomogeneity and discreteness. Continuum mechanics has been used to analyse permanent deformation performance of asphalt mixtures and guide material design for a long time. However, it is inconsistent with the discontinuous characteristics of the asphalt mixture, moreover, there are difficulties in reflecting the actual mechanical properties. A clear relationship between the overall performance and the composed phase is important for understanding the properties of the asphalt mixture. DEM is a numerical simulation method used to study the mechanical behaviour of granular materials, which can be also used to quantify the internal stress of the asphalt mixture during loadings (An, Ling, Geng, & Zhang, 2018;Lin, Wang, & Guo, 2015), considering the anisotropic characteristics of asphalt and aggregate, and simulate the effect of aggregate size, distribution, composition, shape and other factors on the performance of the asphalt mixture. It could provide a new microanalysis solution.
At present, there are two main means to study the asphalt mixtures based on the meso-mechanical model: one is to combine numerical simulation with digital image processing to reconstruct the asphalt mixture model; the other is by means of Discrete Element Method (DEM) or Finite Element Method (FEM). In addition, some scholars began to study viscoelasticity of the asphalt mixture by DEM. The aggregates were modelled as linear elastic, and the binder was modeled as linear viscoelastic or nonlinear viscoelastic in the finite element model (Papagiannakis, Abbas, & Masad, 2002). For each strain level, the viscoelastic behaviour of the binder was described by creating a model and considering the obtained experimental data. By comparing the viscoelastic response prediction and the direct measurement results of the test, the correctness of the asphalt concrete microstructure model in the finite element method and the rationality of the adhesive viscoelastic properties were verified. Through comparisons of experimental results and numerical simulation results, the fracture model of asphalt concrete was established based on the discrete element method (Kim & Buttlar, 2009). The combined effects of aggregate gradation, shape, stiffness, and strength on hot-mix asphalt were analysed by the discrete element method and image processing techniques (Mahmoud, Masad, & Nazarian, 2010). The microstructure of asphalt concretes was captured by using X-Ray Computed Tomography (CT) and using the Discrete Element Method (DEM). The viscoelastic properties of mastics were fitted by the Burger model to Dynamic Shear Rheometer (DSR) data. By comparing the DEM simulation results with the uniaxial creep test results, the DEM model underestimated the creep performance of the mixture in the primary creep stage, but could accurately simulate the second creep stage (Zelelew & Papagiannakis, 2010). The creep responses of asphalt mixture were analysed by the discrete element model considering 3D microstructure; an approach to reduce the computation time with the time-temperature superposition principle was developed (You, Liu, & Dai, 2011).
A meso-scale finite element (FE) model was established to simulate the complex geometry of several types of asphalt mixture with different gradations of aggregates. Different moving wheel loads with different passing velocities were considered and their effects on the mechanical responses of the asphalt mixture were examined (Ziaei-Rad, Nouri, Ziaei-Rad, & Abtahi, 2012). The discrete element program PFC3D was used to fill each particle size aggregate in the simulated asphalt mixture. It was found that the contact force of the larger aggregate was greater than that of the smaller aggregate in the asphalt mixture, and the aggregate size of 4.75 mm was the key particle size of the skeleton asphalt mixture (Chen & Huang, 2012). The two-dimensional discrete element model was used to simulate the permanent deformation characteristics of asphalt concrete, and the actual and virtual simulated asphalt concrete material microscopic results were compared using image processing technology (Fakhri, Kheiry, & Mirghasemi, 2012). At the same time, many dynamic creep tests were carried out. The research results showed that the developed 2D model could better predict the deformation characteristics of asphalt concrete, but the simulation results underestimated the cumulative permanent strain of the material compared with the actual tests. An image-based multi-scale finite element method was developed to predict the mechanical response of asphalt mixtures (Arshadia & Bahia, 2015). Multi-scale effects and four-dimensional levels were considered in the method in order to simulate particle contact at each scale for the presence of mineral particles. The four-dimensional levels were composed of asphalt binder, mastic, mortar, and asphalt mixture. The elastic properties of the binder were obtained through experiments and compiled into ABAQUS user material subroutine for finite element analysis of the asphalt concrete. A uniaxial creep curve and load recovery curve for asphalt mixtures were obtained by 2D image scanning. A comprehensive experimental and theoretical research method proposed by  combined with digital image processing and micromechanical models were used to study low temperature creep stiffness characteristics of asphalt mixtures. The effect of aggregate size on material properties was studied by advanced digital image processing and numerical calculation of two-dimensional images for asphalt mixture . The feasibility of the virtual uniaxial creep test method was verified to accurately predict deformation of the asphalt mixture. The virtual uniaxial creep test could truly reflect the deformation characteristics of different gradations. Preliminary numerical simulation of uniaxial creep test for asphalt mixture was realized (Chen, Zhao & Wan, 2015). The heterogeneous 2D discrete element method was used to simulate the fracture of the asphalt mixture in the semi-circular bending (SCB) test (Renteria, Mahmoud, Yanez, & Burbach, 2017).
A three-dimensional (3D) microstructure-based discrete element (DE) model was developed to study creep behavior of high viscoelastic asphalt sand (HVAS) with the uniaxial compression creep tests. The three-point bending creep tests of asphalt mortar were carried out to obtain the parameters of the Burger model, to determine the transformation method of macroscopic parameters and microscopic parameters of the model in theory, to obtain the parameters used in the discrete element model, and then establish the discrete element analysis model for the asphalt mixture. A 3D-DE digital specimen was composed of coarse aggregates, asphalt mortar and air voids, which could also take gradation, irregular shape, random distribution of aggregate and air voids into consideration, and the boundary conditions of the model were set through the simulation of the uniaxial compression creep tests. An accurate and extensive mapping model of HVAS was built by 3D-PFC (Particle Flow Code), which can provide a simple alternative to the laboratory tests. This method can simulate a series of numerical examples based on different stress levels, coarse aggregate homogenizations, mortar homogenizations and temperatures in a single factor method. Comparison of results of laboratory and numerical tests shows that the 3D-PFC-viscoelastic model can reflect the creep mechanical behaviour of asphalt mixture accurately. It provides the theoretical basis and auxiliary means for analysing the mechanical properties of asphalt mixtures using PFC software. The research on creep behaviour of the asphalt mixture by numerical simulation opens up a new way for research on creep behaviour of the asphalt mixture, it is of considerable theoretical value and has broad application prospects.

Determination of parameters based on DE model
In PFC, the viscoelastic contact parameters of mortar-mortar and mortar-aggregate in HVAS can be characterized by the Burger model. The transformation relationship between macro and micro-parameters is derived in this section.

DE model parameters
70-penetration bitumen was used as asphalt binder for preparation of the specimens. Limestone was used as the aggregate, and the ratio of binder to aggregates was 8.1% by weight. The continuous aggregate gradation, having the nominal maximum size of 2.36 mm, is listed in Table 1. An additional 1% activated rubber crumb and 0.7% TCA (a kind of additive) of the mass of the asphalt mixture were added during the blending process.
Track plate specimens manufactured by wheel rolling in agreement with the standard test method (JTG E20-2011, 2011) were cut into small beams with dimensions of 30 mm × 35 mm × 250 mm. The creep tests were conducted for three hours for each specimen at three temperatures of 25 °C, 40 °C and 60 °C, as shown in Figure 1. The creep compliance was obtained by fitting the test data at different temperatures, and the parameters of the four components of the Burger model are listed in Table 2.

Stiffness model of contact inside asphalt mortar
The normal and tangential models of contacts between asphalt mortar sphere-particles are shown in Figure 2.
Combined with the macroscopic parameters and the microparameters of the asphalt mortar, based on the Burger model, the parameters of the mortar particles can be expressed as where E 1 , η 1 , E 2 and η 2 are the macro Burger's parameters, k mn , c mn , k kn and c kn are the normal micro Burger's parameters of asphalt mortar, k ms , c ms , k ks and c ks are the tangential micro Burger's parameters of asphalt mortar, ν is the Poisson's ratio, and L is the particle thickness, can be understood as the diameter of the particle size, which can be regarded as an elastic beam with two ends at the center of the sphere, the force and bending moment act on both ends of the beam. The beam length is L in PFC and its value is 1, and the unit is mm.
Therefore, the parameters of mortar can be obtained by Eq. (1) based on the Burger model. The micro-parameters for mortar particles include the stiffness parameters, the viscosity parameters and the friction coefficient, which affected the failure modes of asphalt mixture directly, and was set at 0.7 considering actual situations.

Stiffness model of contact between aggregate and mortar
The parameters for aggregate-mortar can be expressed by the Burger model, where k mn , c mn , k kn and c kn are the normal micro Burger's parameters of mortar-aggregate, k ms , c ms , k ks and c ks are the tangential micro Burger's parameters of the mortar-aggregate, E is the elastic modulus of the coarse aggregate and G is the shear modulus of the coarse aggregate. The elastic modulus and the Poisson's ratio of the coarse aggregate are set at 55 GPa and 0.25, respectively, based on engineering experience.

Overview of uniaxial static creep test
When a constant stress and temperature are applied to the asphalt mixture, the corresponding range of strain can be divided into three stages: the deceleration creep stage, isokinetic creep stage and accelerated creep stage.

Preparation of specimens
The specimens were proposed to be cylinder or cube prism according to the previous research. In our experiments, each beam was cut into 40 mm × 40 mm × 80 mm prism specimens to reduce the size effect with a reasonable ratio of height to diameter of 2 according to JTG E20-2011 (Standard Test Method of Asphalt and Asphalt Mixtures for Highway Engineering), as shown in Figure 4.

Creep tests and results
According to the test standards (T0728-2000), the uniaxial compression creep tests were conducted under a constant stress 0.6 MPa at temperature 40 °C by SPT (Standard Penetration Test) as shown in Figure 5, and the time-strain curve was obtained as shown in Figure 6.  According to the time-strain curve in Figure 6, it can be seen that the deceleration creep phase and the constant velocity creep phase of the asphalt mixture are better measured in the tests comparing with the typical creep deformation curve. However, due to the limitations of the test conditions, the accelerated creep, i.e. the destruction phase, is not measured. The initial creep rate is close to 3.67·10 -6 per second. The isokinetic creep rate is close to 1.37·10 -7 per second.

Model geometry based on DEM
The digital specimen of the asphalt mixture was composed of coarse aggregates, asphalt mortar and air voids, in which the size of aggregates was determined according to the gradation of HAVS. Using FISH (a computer language built into PFC), a subroutine was developed to generate each file of aggregates into the wall space with the height of 80 mm, width of 40 mm and length of 40 mm, respectively, as shown in Figure7. The distributions of 1.18-2.36 mm particles are irregular clusters.
The subroutine was developed using FISH to determine the target porosity and the specified number of asphalt mortar particles was generated for the digital model. Random generated air voids were used  Figure 7a, the model represents the coarse aggregate, which is characterized by irregular shape. On the one hand, the particle size is large; on the other hand, the influence of irregular shape distribution on the experiment is considered in this paper. In Figure 7b, the model represents the fine aggregate, which is used in large quantity due to its small particle size. Considering the computer simulation efficiency, it is expressed in terms of spherical particles. Figure 7c shows particle distribution with particle size less than 2.36 mm.

Contact model and material properties
There are a number of contact models available. The constitutive model for a particular contact consists of two parts: the stiffness and the Burger model, respectively. The stiffness model can predict the modulus of aggregate in the DEM model under compressive loads and provide a linear elastic relationship between the contact force and relative displacement between aggregate particles, while the Burger model can be used to predict the viscoelastic properties of various asphalt mortars.
Usually, it is difficult to define the micro-parameters for a digital specimen to represent the properties of the material. The input properties can be derived directly from laboratory results, as listed in Table 3.

Trial load and boundary conditions for DE Model
The force and displacement boundaries were generated in the rectangular prism-shaped asphalt mixture specimen. The speed of top wall was adjusted to make the axial stress reach 0.6 MPa using servo control program. A uniaxial compression creep test model is shown in Figure 8.

Numerical simulation verification
The aggregate and mortar domains were modelled separately and then combined. In the previous research, air voids were generated by removing discrete elements of mortar randomly so that the percentage of void volume reached a desired value. Numerical creep tests of the uniaxial compression were conducted, and E was used as input parameters of the asphalt mixture modulus under the constant stress 0.6 MPa at 40 °C. The aggregate modulus was fixed at 55 000 MPa. The creep results were derived in the form of CSV files, and then the curves were drawn as shown in Figure 9.
It should be noted that the failure stage of the creep curve disappeared because of the computational efficiency of threedimensional discrete element software. As shown in Figure 9, the results obtained by the numerical simulation are in good agreement with the creep tests, and the difference between the simulation result and the experimental value can be explained by the fact that the shape and position of aggregates in the digital specimen are different from those in the test specimen.
The results show that the numerical method provides better approximation for the experimental data over almost the whole range of the curve, moreover the PFC-3D subroutine could effectively predict the creep behaviour of the asphalt mixtures.

Influence of random distribution of the aggregate on simulation results
Three kinds of aggregate position distributions were modelled with the constant stress 0.6 MPa in the creep tests, as shown in Figure 10. The numerical simulations of creep tests were carried out as shown in Figure 11, to compare with the experimental results. Figure 11 shows that the two curves do not coincide with each other causing by the change in aggregate position. Therefore, it is feasible to carry out the macro factor influence analysis for random distribution of the aggregate. The effects of different stress levels (0.4 MPa, 0.6 MPa, 0.8 MPa) on the numerical results were studied by the uniaxial compression creep tests, as shown in Figure 12.
It can be seen from Figure 12 that there is a positive correlation between the creep compliance (the creep strain) and the stress level. The larger the stress level corresponding to the greater creep strain and compliance, the shorter the service life is, which shows that traffic overload has a great influence on the pavement performance and service life of the asphalt mixture. Figure 12 shows that creep and creep compliance curve become steeper with the increase of loading stress, whose variation decreases with time. With the increase of the loading stress level, the time required to reach isokinetic creep stage is reduced. When creep strain is constant, the required loading time reduces with the increase of the stress level.

Influence of homogenization on creep behaviour
This paper attempts to study the effect of homogenization of the aggregate and mortar on the creep behaviour of HVAS combined with Weibull distribution. The probability density function of the Weibull distribution is shown in Eq. (3). , ( where u is the material mechanical parameter, m is the homogeneity of material, and u 0 is a parameter associated with the mean of mechanical parameters. The micro-parameters of the aggregate and asphalt mortar are defined by the Weibull distribution. In the discrete element software, the procedure for the setting of the micro-parameters of the aggregates and asphalt mortar is shown in Figure 13. The u represents the value satisfying the mechanical parameters of the material; m is the shape parameter of the curve of probability density function, the discreteness of reaction mechanics parameters in physics is called the homogeneity of the material, the larger m, the higher the homogeneity of the material is; u 0 is a parameter related to the average value of all mechanical parameters. The difference of m value will affect the shape of the curve. The bigger m value, the higher the probability density function curve is and the narrower it is, which means that the distribution of material parameters is no longer scattered, but rather concentrated.
The digital test was carried out based on uniaxial compression creep test in which the parameters were subjected to the Weibull distribution. The results of creep strain and compliance were obtained with different homogenization, as shown in Figures 14 and 15. Figure 14 shows that the distribution of aggregate stiffness becomes uniform gradually with the increase of the aggregate. There is a negative correlation between the creep strain and the parameter of the aggregate, and the variation of the creep strain decreases with the increase of m. When m approaches infinity, the aggregate stiffness parameters distribute uniformly in the digital specimen, and both the corresponding creep strain and compliance are minimum. In practical engineering, a stable homogeneity of the coarse aggregate is preferable, and the coarse aggregate mixed with large difference should be avoided. Figure 14 shows the influence of the heterogeneity of the coarse aggregate on the creep behaviour of the asphalt mixture, the higher the homogeneity of the coarse aggregate, the more uniform the distribution of the stiffness of coarse aggregate is. The creep strain of the asphalt mixture is negatively correlated with the homogeneity of coarse aggregate, the larger m, the smaller the variation range of creep strain is. When the value of m is close to infinity, the stiffness parameter of coarse aggregate is evenly distributed in the virtual digital specimen, while the creep strain value and creep compliance value are minimum.
To sum up, in practical engineering, the selection of materials should be carefully considered. To avoid mixing coarse aggregate with huge difference in properties, a coarse aggregate with relatively stable homogeneity is the first choice. Figure 15 shows that the creep strain decreases with the increase of m. When m approaches to infinity, the viscoelastic parameters of the asphalt mortar distribute uniformly in the digital specimen, and the corresponding creep strain and the compliance are the minimum. It's necessary to improve the creep resistance of asphalt mixture by considering homogeneity of the asphalt mixture, the mixing time and the method. Figure 15 shows the influence of heterogeneity of asphalt mortar on creep behaviour of asphalt mixture. The creep strain decreases with the increase of homogeneity of asphalt mortar, but the decline is small. When the value of m is close to infinity, the viscoelastic parameters of asphalt mortar are evenly distributed in the virtual digital specimen, and the creep strain and creep compliance of asphalt mixture reach the minimum value. To sum up, the uniformity of the asphalt b) curves of creep compliances with time In the Figure 16, it can be seen that the tails of the curves of the asphalt mixtures become steeper as the temperature rises, and the faster increase corresponds to the steady creep rate. When the creep strain is constant, the loading time decreases with the temperature. Therefore, it is necessary to improve the performance of the asphalt mixture at high temperature to avoid pavement disease and increase the service life in summer. When the loading time is constant, the creep compliance increases with the temperature. When the temperature is constant, the creep compliance increases along with the loading time. Therefore, the creep behaviour of the asphalt mixture is sensitive to the temperature and time under uniaxial loads.