Cross-Functional Test to Explore the Determination Method of Meso-Parameters in the Discrete Element Model of Asphalt Mixtures

In order to obtain more accurate parameters required for the simulation of asphalt mixtures in the discrete element method (DEM), this study carried out a series of cross-functional asphalt mixture experiments to obtain the DEM simulation meso-parameters. By comparing the results of simulation and actual experiments, a method to obtain the meso-parameters of the DEM simulation was proposed. In this method, the numerical aggregate profile was obtained by X-ray CT scanning and the 3D aggregate model was reconstructed in MIMICS. The linear contact parameters of the aggregate and the Burgers model parameters of the asphalt mastic were obtained by nanoindentation technology. The parameters of the parallel bonding model between the aggregate and mastic were determined by the macroscopic tensile adhesion test and shear bond test. The results showed that the meso-parameters obtained by the macroscopic experiment provide a basis for the calibration of DEM parameters to a certain extent. The trends in simulation results are similar to the macro test results. Therefore, the newly proposed method is feasible.


Introduction
With the development of road engineering, various pavement materials and complex mechanical environments increase the difficulty of using traditional macro-mechanics tests to solve engineering problems. Numerical simulation technology brings new prospects for solving complex engineering problems. The discrete element method (DEM) can deal with the mechanical problems of discontinuous media and effectively simulate the discontinuous phenomena such as cracking and separation of materials [1]. The DEM calculation uses the time-stepping algorithm to apply the motion equation on each particle repeatedly. At the beginning of the time step, the force-displacement equations between particles and between particles and the wall are updated to analyze the motion and displacement of each particle in the system, as well as the interaction between them [2]. Particle flow code (PFC) simplifies the DEM to study the mechanical behavior of granular materials from the perspective of the microstructure. In the simulation process, different contact models and corresponding meso-parameters should be appointed for different materials to achieve a relevant macro mechanical response. However, there is no clear mapping between the meso-parameters of the model and the macro mechanical parameters of the object, and the acquisition of accurate meso-parameters with reasonable experiments and calculation methods are the focal point of DEM simulation [3].
Discrete element numerical simulation technology can solve various current difficulties in the actual macro experiment, such as long periods, low cost-effectiveness and dentation technology, X-ray CT technology, the tensile adhesion test and the shear bond test. X-ray CT technology can obtain a realistic profile of the aggregate. Nanoindentation technology also provides technical support for appointing the mechanical parameters of asphalt mastic and aggregate in the in situ nondestructive asphalt concrete. Combined with the macroscopic tensile adhesion test and shear bond test, the mechanical parameters of the asphalt mixture were determined. After that, the simulation and actual test results of the dynamic modulus were compared to verify the accuracy of the proposed method.

Materials
Due to a large amount of fine aggregate particles in continuous-graded asphalt mixture, the interface between the aggregate and asphalt mastic is complex, which brings more difficulties for the nanoindentation test. However, for the gap-graded asphalt mixture, the interface between the aggregate and asphalt mastic is clearer, which is suitable for the nanoindentation test. Therefore, the gap-graded aggregate was selected to prepare the asphalt mixture specimens. The gradation is shown in Table 1. The basalt aggregate and a PG64-22 base asphalt binder were used in this study.

Macro Dynamic Modulus Test
The dynamic modulus test was adopted to measure the macroscopic properties of the mixture according to the test standard AASHTO TP62-03. The specimens were loaded at the frequencies of 25, 10, 5 and 1 Hz with a universal material testing machine UTM-25. The test specimens were cored from standard gyratory compacted mixtures. The gyratory compacted mixtures were prepared according to AASHTO T312. The air void of the mixture was 7%. The temperature of the dynamic modulus test was 25 • C.

Aging of Specimen
Mixtures with four different aging conditions were selected, including non-aging, aging for 2 d, aging for 5 d and aging for 8 d. The aging test was performed on the dynamic modulus specimens based on AASHTO PP2. The specimens aged for 2 d, 5 d and 8 d were obtained by maintaining non-aged specimens at 85 • C for different amounts of time.

Direct Tensile Test and Uniaxial Compression Test
The dynamic modulus specimens with four different aging conditions were cut into cuboid specimens. The cuboid specimens were fixed on the loading plate with epoxy resin. The test temperature of the direct tensile test was 25 • C and the loading rate was 0.5 mm/min. The process is shown in Figure 1. The direct tensile test was conducted according to AASHTO T314. The uniaxial compression test was conducted according to ASTM D 1074. The loading rate was 50 mm/min. A UTM-25 was used to load the cuboid specimens. Materials 2021, 14, x FOR PEER REVIEW 4 of 16 Figure 1. The process of direct tensile test.

Tensile Adhesion Test
The tensile adhesion test was conducted by an LGZ-2 structural layer material strength tensile instrument based on the test standard EN 12697-48. The instrument is shown in Figure 2. The binder in the aging specimen was extracted and recovered based on the standard AASHTO TP2. Then the binder was spread evenly on a stone slab. After that, the drawing head was pressed against the binder surface and kept warm for a while to make a tensile adhesion test specimen. The specimens are shown in Figure 3a. After the test, the drawing head was pulled up, as shown in Figure 3b. The test temperature was 25 °C. Tensile strength was collected to calibrate the normal bonding strength of the parallel bonding model.

Tensile Adhesion Test
The tensile adhesion test was conducted by an LGZ-2 structural layer material strength tensile instrument based on the test standard EN 12697-48. The instrument is shown in Figure 2. The binder in the aging specimen was extracted and recovered based on the standard AASHTO TP2. Then the binder was spread evenly on a stone slab. After that, the drawing head was pressed against the binder surface and kept warm for a while to make a tensile adhesion test specimen. The specimens are shown in Figure 3a. After the test, the drawing head was pulled up, as shown in Figure 3b. The test temperature was 25 • C. Tensile strength was collected to calibrate the normal bonding strength of the parallel bonding model.

Tensile Adhesion Test
The tensile adhesion test was conducted by an LGZ-2 structural layer material strength tensile instrument based on the test standard EN 12697-48. The instrument is shown in Figure 2. The binder in the aging specimen was extracted and recovered based on the standard AASHTO TP2. Then the binder was spread evenly on a stone slab. After that, the drawing head was pressed against the binder surface and kept warm for a while to make a tensile adhesion test specimen. The specimens are shown in Figure 3a. After the test, the drawing head was pulled up, as shown in Figure 3b. The test temperature was 25 °C. Tensile strength was collected to calibrate the normal bonding strength of the parallel bonding model.

Tensile Adhesion Test
The tensile adhesion test was conducted by an LGZ-2 structural layer material strength tensile instrument based on the test standard EN 12697-48. The instrument is shown in Figure 2. The binder in the aging specimen was extracted and recovered based on the standard AASHTO TP2. Then the binder was spread evenly on a stone slab. After that, the drawing head was pressed against the binder surface and kept warm for a while to make a tensile adhesion test specimen. The specimens are shown in Figure 3a. After the test, the drawing head was pulled up, as shown in Figure 3b. The test temperature was 25 °C. Tensile strength was collected to calibrate the normal bonding strength of the parallel bonding model.

Shear Bond Test
The shear bond test was performed on a UTM-25 based on EN 12697-48. After the extracted binder was spread evenly on a stone slab, another slab was covered on the binder. The test specimen is shown in Figure 4. The test process is shown in Figure 5. The shear strength could provide the basis for the calibration of tangential bonding strength in a parallel bonding model. The test temperature was 25 • C.

Shear Bond Test
The shear bond test was performed on a UTM-25 based on EN 12697-48. After the extracted binder was spread evenly on a stone slab, another slab was covered on the binder. The test specimen is shown in Figure 4. The test process is shown in Figure 5. The shear strength could provide the basis for the calibration of tangential bonding strength in a parallel bonding model. The test temperature was 25 °C.

Nanoindentation Test
The nanoindentation test was applied to evaluate the mechanical properties of the aggregate phase and mastic phase. The gyratory compacted specimens were cut into small squares and cured with epoxy resin. After polishing, the small square was exposed. Finally, a nanoindentation test specimen was made, which is shown in Figure 6. The specimen was set on the nanoindentation instrument, and the appropriate measurement area was selected. The test point interval of the nanoindentation test was set as 80 μm, as shown in Figure 7. The test temperature was 25 °C.

Shear Bond Test
The shear bond test was performed on a UTM-25 based on EN 12697-48. After the extracted binder was spread evenly on a stone slab, another slab was covered on the binder. The test specimen is shown in Figure 4. The test process is shown in Figure 5. The shear strength could provide the basis for the calibration of tangential bonding strength in a parallel bonding model. The test temperature was 25 °C.

Nanoindentation Test
The nanoindentation test was applied to evaluate the mechanical properties of the aggregate phase and mastic phase. The gyratory compacted specimens were cut into small squares and cured with epoxy resin. After polishing, the small square was exposed. Finally, a nanoindentation test specimen was made, which is shown in Figure 6. The specimen was set on the nanoindentation instrument, and the appropriate measurement area was selected. The test point interval of the nanoindentation test was set as 80 μm, as shown in Figure 7. The test temperature was 25 °C.

Nanoindentation Test
The nanoindentation test was applied to evaluate the mechanical properties of the aggregate phase and mastic phase. The gyratory compacted specimens were cut into small squares and cured with epoxy resin. After polishing, the small square was exposed. Finally, a nanoindentation test specimen was made, which is shown in Figure 6. The specimen was set on the nanoindentation instrument, and the appropriate measurement area was selected. The test point interval of the nanoindentation test was set as 80 µm, as shown in Figure 7. The test temperature was 25 • C.

Shear Bond Test
The shear bond test was performed on a UTM-25 based on EN 12697-48. After the extracted binder was spread evenly on a stone slab, another slab was covered on the binder. The test specimen is shown in Figure 4. The test process is shown in Figure 5. The shear strength could provide the basis for the calibration of tangential bonding strength in a parallel bonding model. The test temperature was 25 °C.

Nanoindentation Test
The nanoindentation test was applied to evaluate the mechanical properties of the aggregate phase and mastic phase. The gyratory compacted specimens were cut into small squares and cured with epoxy resin. After polishing, the small square was exposed. Finally, a nanoindentation test specimen was made, which is shown in Figure 6. The specimen was set on the nanoindentation instrument, and the appropriate measurement area was selected. The test point interval of the nanoindentation test was set as 80 μm, as shown in Figure 7. The test temperature was 25 °C.

Discrete Element Model of Asphalt Mixtures
The coarse aggregate was scanned by a Y. CT PRECISION S scanner provided by School of Materials Science and Engineering, Southeast University. The height interval of CT scanning was 0.1 mm. Due to the great difference between the density of the aggregate and air, the gray values of aggregates in the scanning image were obviously different from that of the background. The aggregate profile was reconstructed by the materialize interactive medical image control system (MIMICS 17.0). The three-dimensional structure model was created by importing the CT section image into the software and interpolating the gray value of the image. The appropriate amount of CT scanning section images was selected, then the aggregate numerical profile was exported, as shown in Figure 8. The DEM software used in this study was PFC 5.0. According to the size of the macro dynamic modulus specimen, the virtual specimen was created in the DEM software. The appropriate amount of aggregate was put in the space according to the corresponding gradation. The aggregate distribution of asphalt mixture in DEM mode is shown in Figure  9.

Discrete Element Model of Asphalt Mixtures
The coarse aggregate was scanned by a Y. CT PRECISION S scanner provided by School of Materials Science and Engineering, Southeast University. The height interval of CT scanning was 0.1 mm. Due to the great difference between the density of the aggregate and air, the gray values of aggregates in the scanning image were obviously different from that of the background. The aggregate profile was reconstructed by the materialize interactive medical image control system (MIMICS 17.0). The three-dimensional structure model was created by importing the CT section image into the software and interpolating the gray value of the image. The appropriate amount of CT scanning section images was selected, then the aggregate numerical profile was exported, as shown in Figure 8.

Discrete Element Model of Asphalt Mixtures
The coarse aggregate was scanned by a Y. CT PRECISION S scanner provided by School of Materials Science and Engineering, Southeast University. The height interval of CT scanning was 0.1 mm. Due to the great difference between the density of the aggregate and air, the gray values of aggregates in the scanning image were obviously different from that of the background. The aggregate profile was reconstructed by the materialize interactive medical image control system (MIMICS 17.0). The three-dimensional structure model was created by importing the CT section image into the software and interpolating the gray value of the image. The appropriate amount of CT scanning section images was selected, then the aggregate numerical profile was exported, as shown in Figure 8. The DEM software used in this study was PFC 5.0. According to the size of the macro dynamic modulus specimen, the virtual specimen was created in the DEM software. The appropriate amount of aggregate was put in the space according to the corresponding gradation. The aggregate distribution of asphalt mixture in DEM mode is shown in Figure  9. The DEM software used in this study was PFC 5.0. According to the size of the macro dynamic modulus specimen, the virtual specimen was created in the DEM software. The appropriate amount of aggregate was put in the space according to the corresponding gradation. The aggregate distribution of asphalt mixture in DEM mode is shown in Figure 9.

Discrete Element Model of Asphalt Mixtures
The coarse aggregate was scanned by a Y. CT PRECISION S scanner provided by School of Materials Science and Engineering, Southeast University. The height interval of CT scanning was 0.1 mm. Due to the great difference between the density of the aggregate and air, the gray values of aggregates in the scanning image were obviously different from that of the background. The aggregate profile was reconstructed by the materialize interactive medical image control system (MIMICS 17.0). The three-dimensional structure model was created by importing the CT section image into the software and interpolating the gray value of the image. The appropriate amount of CT scanning section images was selected, then the aggregate numerical profile was exported, as shown in Figure 8. The DEM software used in this study was PFC 5.0. According to the size of the macro dynamic modulus specimen, the virtual specimen was created in the DEM software. The appropriate amount of aggregate was put in the space according to the corresponding gradation. The aggregate distribution of asphalt mixture in DEM mode is shown in Figure  9. After the aggregate skeleton with specific gradation was established, the blank space of the virtual specimens was filled with spherical elements to form the mastic phase. The distribution of voids along the height of the specimen was obtained by X-ray scanning, and the numerical model curve of the void phase was established. The specimen was divided into several layers along the height direction, and the height of each layer was 10 mm. The number of spherical elements representing the void phase in each layer was calculated, and the corresponding spherical elements in the mastic phase were deleted. Finally, a three-dimensional specimen was constructed, as shown in Figure 10. After the aggregate skeleton with specific gradation was established, the blank space of the virtual specimens was filled with spherical elements to form the mastic phase. The distribution of voids along the height of the specimen was obtained by X-ray scanning, and the numerical model curve of the void phase was established. The specimen was divided into several layers along the height direction, and the height of each layer was 10 mm. The number of spherical elements representing the void phase in each layer was calculated, and the corresponding spherical elements in the mastic phase were deleted. Finally, a three-dimensional specimen was constructed, as shown in Figure 10.

Meso-Parameters of the Discrete Element Model
In this study, the linear model was selected to characterize the contact between aggregates. The parallel bonding model was used to characterize the contact between the mastic and aggregate, or the internal contact of the asphalt mastic. Since the DEM software (PFC 5.0) can only assign one contact model to a contact, to characterize the internal contact of the asphalt mastic, the Burgers model and the parallel bonding model were distributed among the spherical elements of mastic phase in a certain ratio. The schematic diagram of the contact model in asphalt mixture is shown in Figure 11.

Meso-Parameters of the Discrete Element Model
In this study, the linear model was selected to characterize the contact between aggregates. The parallel bonding model was used to characterize the contact between the mastic and aggregate, or the internal contact of the asphalt mastic. Since the DEM software (PFC 5.0) can only assign one contact model to a contact, to characterize the internal contact of the asphalt mastic, the Burgers model and the parallel bonding model were distributed among the spherical elements of mastic phase in a certain ratio. The schematic diagram of the contact model in asphalt mixture is shown in Figure 11.

Determination of Linear Contact Parameters
There was a linear relationship between the contact force and displacement of the linear contact model, the key parameter of which was the stiffness between the contact entities [20]. The load-depth curve of the aggregate in the mixture is shown in Figure 12.

Determination of Linear Contact Parameters
There was a linear relationship between the contact force and displacement of the linear contact model, the key parameter of which was the stiffness between the contact entities [20]. The load-depth curve of the aggregate in the mixture is shown in Figure 12. The test load was maintained for 200 s at the load level of 2.05 mN. Similar to the typical nanoindentation load-depth curve, the indentation test on the aggregate did not show the negative slope of the unloading section curve. After calculation of the instrument, the elastic modulus of aggregate was determined as 55 GPa, the friction coefficient of the aggregate was determined as 0.5 and the Poisson ratio of the aggregate was determined as 0.25. Figure 11. The schematic diagram of the contact model in the asphalt mixture.

Determination of Linear Contact Parameters
There was a linear relationship between the contact force and displacement of the linear contact model, the key parameter of which was the stiffness between the contact entities [20]. The load-depth curve of the aggregate in the mixture is shown in Figure 12. The test load was maintained for 200 s at the load level of 2.05 mN. Similar to the typical nanoindentation load-depth curve, the indentation test on the aggregate did not show the negative slope of the unloading section curve. After calculation of the instrument, the elastic modulus of aggregate was determined as 55 GPa, the friction coefficient of the aggregate was determined as 0.5 and the Poisson ratio of the aggregate was determined as 0.25.

Determination of Burgers Model Parameters
The Burgers model characterizes the time-varying relationship of the interaction between different particles or between particles and the boundary. The significant influential parameters, which include the elastic coefficient Em and viscosity coefficient ηm of the

Determination of Burgers Model Parameters
The Burgers model characterizes the time-varying relationship of the interaction between different particles or between particles and the boundary. The significant influential parameters, which include the elastic coefficient E m and viscosity coefficient η m of the Maxwell model, as well as the elastic coefficient E k and viscosity coefficient η k of the Kelvin model [21]. In order to obtain the meso-parameters of the Burgers model, the nanoindentation test was carried out on the asphalt mastic phase in the specimens. By studying the creep process of mastic, the constitutive parameters of the Burgers model were calculated, and then the component parameters of the Burgers model were converted into meso-parameters of DEM through equations.
The nanoindentation test was applied to the specimens with four aging conditions at 25 • C, and then the load section of the test curve was analyzed. Finally, the parameters of the Burgers constitutive model under quasi-static uniaxial compression were fitted, as shown in Table 2. The parameters of the Burgers model were converted into the meso-parameters in the DEM. The normal contact parameters were calculated according to Equations (1)-(4). Then, the corresponding tangential contact parameters were obtained according to the relationship between elastic modulus and shear modulus (E = 2(1 + υ)G) [22], as shown in Equations (5)- (8).
where the ν is the Poisson ratio of the material; L is the length of the viscoelastic beam, that is, the center distance of adjacent spherical elements; E 1 and η 1 are the modulus of the spring element and the viscosity of the dashpot element in the Maxwell model, respectively; E 2 and η 2 are the modulus of the spring element and the viscosity of the dashpot element in the Kelvin model, respectively. After calculation, the meso-parameters of the Burgers model required for DEM were obtained. The meso-parameters are shown in Table 3.

Determination of Parallel Bonding Model Parameters
Since the Burgers model in DEM does not have adhesion, the particles seriously scatter when they are loaded under the unconfined test conditions. The parallel bonding model has a certain deformability, which can transfer tensile force, shear force and moment, as well as restrict the rotation of particles to a certain extent and simulate the constitutive behavior of bonding materials. This model could describe the bonding or adhesion strength of the asphalt mastic [23]. The characteristic parameters of the parallel bonding model include the linear contact modulus (emod) and parallel bonding modulus (Pb_emod), normal bonding strength, tangential bonding strength and stiffness ratio [24].
In this study, the parallel bonding parameters were determined according to the following steps: (1) The emod was fixed, and the Pb_emod was calculated under direct tensile conditions. The contact model was only set as the parallel bonding model. The direct tensile simulation test was performed. The emod was set to a small value (1 × 10 5 ), and then the Pb_emod was changed to 1 × 10 9 , 5 × 10 9 , 10 × 10 9 and 20 × 10 9 , respectively. The relationship between the Pb_emod and the simulated tensile elastic modulus E t could be obtained, as shown in Figure 13.
stitutive behavior of bonding materials. This model could describe the bonding or adhesion strength of the asphalt mastic [23]. The characteristic parameters of the parallel bonding model include the linear contact modulus (emod) and parallel bonding modulus (Pb_emod), normal bonding strength, tangential bonding strength and stiffness ratio [24].
In this study, the parallel bonding parameters were determined according to the following steps: (1) The emod was fixed, and the Pb_emod was calculated under direct tensile conditions.
The contact model was only set as the parallel bonding model. The direct tensile simulation test was performed. The emod was set to a small value (1e5), and then the Pb_emod was changed to 1e9, 5e9, 10e9 and 20e9, respectively. The relationship between the Pb_emod and the simulated tensile elastic modulus Et could be obtained, as shown in Figure 13. It can be seen from Figure 13 that there is a positive linear relationship between E t and the Pb_emod. The corresponding relationship is obtained by linear fitting, as shown in Equation (9). E t = 1.3002x + 0.03397 (9) where the E t is the tensile elastic modulus (GPa) and x is the actual value of Pb_emod/1 × 10 9 (GPa).
(2) The Pb_emod was fixed, and the emod was calculated under uniaxial compression conditions.
The uniaxial compression simulation test was performed in the DEM software. The Pb_emod was fixed and the emod was changed. The relationship between the emod and the compressive elastic modulus E c obtained by the uniaxial compression simulation test is shown in Figure 14. It can be seen from Figure 13 that there is a positive linear relationship between Et and the Pb_emod. The corresponding relationship is obtained by linear fitting, as shown in Equation (9). Et = 1.3002x + 0.03397 (9) where the Et is the tensile elastic modulus (GPa) and x is the actual value of Pb_emod/1e9 (GPa).
(2) The Pb_emod was fixed, and the emod was calculated under uniaxial compression conditions.
The uniaxial compression simulation test was performed in the DEM software. The Pb_emod was fixed and the emod was changed. The relationship between the emod and the compressive elastic modulus Ec obtained by the uniaxial compression simulation test is shown in Figure 14.
It can be seen from Figure 14 that the compressive elastic modulus EC is positively correlated with the emod. The linear relationship between the EC and emod is shown in Equation (10).
where the EC is the compressive elastic modulus (GPa) and x is the actual value of emod/1e9 (GPa). (3) The Pb_emod and emod were calibrated by macro tests. The tensile elastic modulus and compression elastic modulus were obtained through the direct tensile test and uniaxial compression test, respectively. The results are shown in Table 4. After that, the Pb_emod and emod were calculated according to Equations (9) It can be seen from Figure 14 that the compressive elastic modulus E C is positively correlated with the emod. The linear relationship between the E C and emod is shown in Equation (10). E C = 1.7825x − 0.0449 (10) where the E C is the compressive elastic modulus (GPa) and x is the actual value of emod/1 × 10 9 (GPa).
(3) The Pb_emod and emod were calibrated by macro tests. The tensile elastic modulus and compression elastic modulus were obtained through the direct tensile test and uniaxial compression test, respectively. The results are shown in Table 4. After that, the Pb_emod and emod were calculated according to Equations (9) and (10), respectively. The Pb_emod and emod results are shown in Table 5. (4) The determination of the stiffness ratio of parallel bonding component and linear contact component.
The stiffness ratios were set to 1, 1.5, 2 and 2.5, and the uniaxial compression test was performed. The strain corresponding to the peak strength was used to calculate the Poisson ratio, respectively. The relationship between the Poisson ratio and the stiffness ratio was obtained by fitting, as shown in Equation (11). µ = 0.1083x + 0.0821 (11) where the µ is the Poisson ratio, and x is the stiffness ratio. The Poisson ratio of asphalt mastic was chosen to be 0.35. The stiffness ratio was calculated by Equation (11), and the result was 2.5.
(5) The determination of normal and tangential bond strength.
The tensile and shear forces between the asphalt mastic and aggregate were obtained by the tensile test and shear bond test, which provided the basis for the calibration of normal and tangential bonding strength in the parallel bonding model. The results of the tensile strength and shear strength are listed in Table 6. The tensile strength decreased with the increase of aging time, while the shear strength increased with the aging time.

Macro Tests of Asphalt Mixture
Under different loading frequencies, the dynamic modulus of the asphalt mixture is shown in Table 7. The higher the frequency, the greater the dynamic modulus of the asphalt mixture. In general, the dynamic modulus of the mixture increased with the deepening of the aging degree, which indicates the aging hardening behavior of the asphalt mixture.

Simulation of the Asphalt Mixture Model in DEM
The size and gradation of the specimens in the simulative dynamic modulus test were the same as those of the macro material experiments. During the simulation, the lower boundary was fixed. In order to simulate the loading process of the macro dynamic modulus test, the moving velocity of the upper boundary was changing in the waveform of a half-sine wave. The schematic diagram of loading in the DEM software is shown in Figure 15. After a trial calculation, the ratio of the parallel bonding model and the Burgers model was taken as 8:2. At this ratio, it was observed that the particles of the spherical elements scattered to a certain extent. boundary was fixed. In order to simulate the loading process of the macro dynamic modulus test, the moving velocity of the upper boundary was changing in the waveform of a half-sine wave. The schematic diagram of loading in the DEM software is shown in Figure  15. After a trial calculation, the ratio of the parallel bonding model and the Burgers model was taken as 8:2. At this ratio, it was observed that the particles of the spherical elements scattered to a certain extent. The stress-strain curve was obtained by monitoring the distance between the upper and lower boundaries and the load stress. Figure 16 shows the simulative stress-strain curves of the non-aged mixture under different loading frequencies. It can be found that the strain decreased with the increase of frequency, and the variation trend was the same as that of the macro test. The dynamic modulus under different frequencies was calculated by taking the ratio of stress and strain amplitude from the last 3-5 cycles, and the results are shown in Table 8. It can be seen that the simulated dynamic modulus increased with the increase of the loading frequency and aging time. The stress-strain curve was obtained by monitoring the distance between the upper and lower boundaries and the load stress. Figure 16 shows the simulative stress-strain curves of the non-aged mixture under different loading frequencies. It can be found that the strain decreased with the increase of frequency, and the variation trend was the same as that of the macro test. The dynamic modulus under different frequencies was calculated by taking the ratio of stress and strain amplitude from the last 3-5 cycles, and the results are shown in Table 8. It can be seen that the simulated dynamic modulus increased with the increase of the loading frequency and aging time. and lower boundaries and the load stress. Figure 16 shows the simulative stress-strain curves of the non-aged mixture under different loading frequencies. It can be found that the strain decreased with the increase of frequency, and the variation trend was the same as that of the macro test. The dynamic modulus under different frequencies was calculated by taking the ratio of stress and strain amplitude from the last 3-5 cycles, and the results are shown in Table 8. It can be seen that the simulated dynamic modulus increased with the increase of the loading frequency and aging time.

Comparison of the Macro Test and Simulation Results
The comparison of simulated and macro dynamic modulus results under different frequencies and aging conditions are shown in Figure 17. In general, the method proposed in this study could simulate the macro properties of the mixture to a certain extent. This indicates that the method has the potentials to simulate the macroscopic mechanical response of asphalt mixtures. However, under the same conditions, the simulated values were lower than the actual values. The reason may be that the size of the spherical elements was too large, which led to the excessive void content of the virtual specimen, and the skeleton was not as tight as the actual specimen. In addition, the spherical particles scattered to a certain extent during the loading process, which may also cause the simulated strain value to be larger than the actual value.

Comparison of the Macro Test and Simulation Results
The comparison of simulated and macro dynamic modulus results under different frequencies and aging conditions are shown in Figure 17. In general, the method proposed in this study could simulate the macro properties of the mixture to a certain extent. This indicates that the method has the potentials to simulate the macroscopic mechanical response of asphalt mixtures. However, under the same conditions, the simulated values were lower than the actual values. The reason may be that the size of the spherical elements was too large, which led to the excessive void content of the virtual specimen, and the skeleton was not as tight as the actual specimen. In addition, the spherical particles scattered to a certain extent during the loading process, which may also cause the simulated strain value to be larger than the actual value.

Conclusions
In this study, the direct tensile test, uniaxial compression test, shear bond test and tensile adhesion test were performed to obtain the parameters of the parallel bonding model. Nanoindentation tests were carried out on the aggregate phase and mastic phase of mixture specimens to obtain the parameters of the linear contact model and the Burgers model. Besides, X-ray CT technology was adopted to obtain a more real profile of the aggregate. Through comprehensive analysis, the following conclusions can be drawn.
The viscoelastic parameters obtained by the nanoindentation test could reflect the viscoelastic properties of the asphalt mastic. The newly proposed method for obtaining the meso-parameters of asphalt mixtures can be applied to the DEM. The meso-parameters obtained by this method can be used to simulate the macroscopic response of asphalt mixtures in DEM software. The simulated dynamic modulus decreased with the increase of frequency and aging time, which was consistent with the actual macro tests. However, it was also found the simulated values were lower than the actual values. In further research, reducing the particle size of the spherical elements will be taken into consideration to further improve the simulation accuracy.

Conclusions
In this study, the direct tensile test, uniaxial compression test, shear bond test and tensile adhesion test were performed to obtain the parameters of the parallel bonding model. Nanoindentation tests were carried out on the aggregate phase and mastic phase of mixture specimens to obtain the parameters of the linear contact model and the Burgers model. Besides, X-ray CT technology was adopted to obtain a more real profile of the aggregate. Through comprehensive analysis, the following conclusions can be drawn.
The viscoelastic parameters obtained by the nanoindentation test could reflect the viscoelastic properties of the asphalt mastic. The newly proposed method for obtaining the meso-parameters of asphalt mixtures can be applied to the DEM. The meso-parameters obtained by this method can be used to simulate the macroscopic response of asphalt mixtures in DEM software. The simulated dynamic modulus decreased with the increase of frequency and aging time, which was consistent with the actual macro tests. However, it was also found the simulated values were lower than the actual values. In further research, reducing the particle size of the spherical elements will be taken into consideration to further improve the simulation accuracy.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

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