Influence of the Mesoscopic Viscoelastic Contact Model on Characterizing the Rheological Behavior of Asphalt Concrete in the DEM Simulation

)e numerical simulation based on the discrete element method (DEM) is popular to analyze the material behavior of asphalt concrete in recent years because of the advantage of the DEM in characterizing the heterogeneous microstructures. As a type of viscoelastic material, the rheological behavior of asphalt concrete is represented depending on the mesoscopic viscoelastic contact model between two particles in a contact in DEM simulations. However, what is missing in the existing literature studies is analysis of the influence of the mesoscopic viscoelastic contact models. Hence, the existing mesoscopic viscoelastic contact models are employed to build different types of DEM numerical samples of asphalt concrete in this study. Laboratory tests and the corresponding numerical tests at different temperatures and frequencies are implemented to investigate the difference in simulation precision in the case of using different mesoscopic viscoelastic contact models via the rheological index of dynamic modulus and phase angle. )e results show the following: (1) the mesoscopic generalized Maxwell contact model provides the best simulation precision at low temperatures; (2) the mesoscopic generalized Kelvin contact model shows an improved precision at high temperatures; and (3) although the mesoscopic Burgers contact model has the simplest mathematical structure, the simulation precisions are obviously lower than those of the other two contact models, particularly when simulating the phase angle at low temperatures and frequencies. )e results will be beneficial to select the appropriate mesoscopic contact model for the DEM modeling of asphalt concrete according to the loading conditions.


Introduction
Asphalt concrete is one of the most important construction materials in the highway industry [1], which can be treated as a three-level distributed system: the coarse aggregates, the asphalt mastics, and the voids [2,3]. e existence of coarse aggregates leads to the conspicuous granular behavior of asphalt concrete [4,5]. Hence, the discrete element method (DEM) introduced by Cundall and Strack [6,7] for studying granular materials [8] draws engineers' great attention. Since the mid-2000s, the DEM has become a prevailing modeling and simulation method for studying the performance of asphalt concrete, of which the foremost advantage for characterizing some important mechanical and structural performance of asphalt concrete is that the method provides an explicit treatment of the materials at a mesoscopic level.
In the numerical simulation based on the DEM, asphalt concrete is always modeled using three parts: the coarse aggregate elements, the asphalt mastic elements, and the void elements. e performance of the aggregate-mastic interaction is characterized as that of contacts between coarse aggregate elements and asphalt mastic elements. Materials in the DEM are modeled as a system of distinct, interacting, and spherical particles subject to motion and deformation. A force-displacement law updated each time step is acted on the contacts between particles, from which macroscopic characteristics of materials are obtained through the interaction of contacts between particles at the mesoscopic scale. ereupon, complex mechanical characteristics of materials can be elucidated by DEM simulations [9]. e focus of this paper is the force-displacement law (i.e., the mesoscopic contact model) in the DEM.
In the DEM models, the coarse aggregates treated as the densely assembled elastic particles described their interactions by the mesoscopic elastic contact models. e asphalt mastics, containing the asphalt binder, fine aggregates, and mineral powder, are treated as the densely assembled viscoelastic particles with equal radius, and their interactions are described through the mesoscopic viscoelastic contact models. Obviously, the mesoscopic viscoelastic contact models acting on the interactions of asphalt mastic particles play a significant role in characterizing the rheological behavior of asphalt concrete.
In the past decade, the applied mesoscopic viscoelastic contact models are increasingly complex. At first, the simple bond contact model was used to describe the viscoelastic properties of asphalt mastics in the studies of Kim et al. [10][11][12] and Mahmoud et al. [13]. Afterwards, Qian et al. [14] and Yang et al. [15] adopted the cohesive models to reflex the viscosity characteristics of asphalt mastics. Meanwhile, a mesoscopic Burgers viscoelastic contact model (embedded by "Itasca Corporation") instead of the simple bond contact models was applied in analyzing the viscoelastic mechanical response of asphalt concrete [16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34]. Research showed that mesoscopic viscoelastic contact models influenced the asphalt concrete performance significantly, especially the rheological behavior. e Burgers contact models used in these studies are all invoked directly from the software (PFC) produced by Itasca Corporation. Subsequently, Ren and Sun [35] developed a new mesoscopic viscoelastic contact model, namely, the mesoscopic generalized Maxwell contact model, which could describe the low-temperature rheological behavior of asphalt concrete more accurately than the mesoscopic Burgers contact model in the DEM [36][37][38]. Shimizu et al. [39] and Jin et al. [40] put forward the mesoscopic generalized Kelvin contact model in the DEM but had not been applied to asphalt-based materials.
e mesoscopic Burgers contact model (B-model), the mesoscopic generalized Maxwell contact model (M-model), and the mesoscopic generalized Kelvin contact model (Kmodel) are the only three mesoscopic viscoelastic contact models in the DEM simulations at present. However, what is missing in the existing literature studies is analysis of the influence of the three mesoscopic viscoelastic contact models on characterizing the rheological behavior of asphalt concrete. Hence, there is not efficient guidance that selects the appropriate mesoscopic viscoelastic contact model under different loading conditions in the DEM simulations.
In response to the above problem, three types of asphalt concrete DEM models using the B-model, M-model, and K-model are built to investigate the difference in their rheological behaviors compared to the experimental data based on the dynamic modulus and phase angle via the dynamic modulus tests in the case of different temperatures (− 6°C, 4°C, 16°C, 27°C, and 38°C) and loading frequencies (0.1 Hz, 0.5 Hz, 1 Hz, 5 Hz, 10 Hz, and 25 Hz). e results will be beneficial to select the appropriate mesoscopic contact model for the DEM modeling of asphalt concrete according to the loading conditions. Figure 1 shows the numerical sample of asphalt concrete using a 2D DEM modeling technology proposed in our previous study [36] based on the software PFC. e numerical sample contained coarse aggregate elements (red region), asphalt mastic elements (black region), and void elements (white region). e coarse aggregate elements represent the aggregates whose diameters are larger than 1.18 mm. e asphalt mastic elements represent the combination of the asphalt binder and fine aggregates whose diameters are smaller than 1.18 mm. e coarse aggregate elements and asphalt mastic elements are both treated as densely assembled balls with an equal radius (0.5 mm), and there is no ball in the void elements.

Constitutive Relations of Different Contacts in Numerical
Samples. In the numerical sample, four types of contacts in the numerical sample of asphalt concrete are used to represent four different interactions: contacts within coarse aggregate elements, between adjacent coarse aggregate elements, between coarse aggregate elements and asphalt mastic elements, and within asphalt mastic elements. Obviously, the properties of any contact in the numerical samples can be characterized through getting the corresponding mesoscopic contact models of the two balls (two identical aggregate elements, two different aggregate elements, two identical mastic elements, or one aggregate element and a mastic element) in a contact. As shown in Figures 2-4, the constitutive relations within aggregates and between adjacent aggregates are the same among different types of numerical samples, which are described via two tandem springs acting on the coarse aggregate element. e constitutive relations between aggregates and asphalt mastics are characterized via stringing one spring acting on the coarse aggregate element and one mesoscopic viscoelastic contact model acting on the asphalt mastic element. e constitutive relations within asphalt mastics are represented via two tandem mesoscopic viscoelastic contact models acting on the asphalt mastic element.
As previously mentioned, the B-model is employed as the built-in model of PFC in the existing studies. e M-model and the K-model have been developed by Ren and Sun [35] and Shimizu et al. [39] via the finite difference timedomain method [41], respectively. For a more detailed and rigorous introduction about the mathematic relations of the M-model and K-model, readers can refer to the studies of Ren and Sun [35] and Shimizu et al. [39].

Normal direction
Shear direction a: aggregate b: aggregate Contacts within aggregates and between adjacent aggregates Normal direction Shear direction a: aggregate b: asphalt mastic Contacts between aggregate and asphalt mastic

Normal direction
Shear direction a: aggregate b: asphalt mastic Contacts between aggregate and asphalt mastic Shear direction modulus and phase angle of the asphalt mastic at different frequencies.
e parameter calibration method of the B-model, the M-model, and the K-model has been established by Liu and You [22], Ren and Sun [35], and Jin et al. [40], respectively, which is not detailedly introduced here.
e Young modulus of the coarse aggregates is 59.2 GPa. e gradation of the asphalt mastic used in this study is shown in Table 1.
e dynamic modulus and phase angle of the asphalt mastic are measured at different temperatures (− 6°C, 4°C, 15°C, 27°C, and 38°C) and frequencies (0.1 Hz, 0.5 Hz, 1 Hz, 5 Hz, 10 Hz, and 25 Hz), of which the testing results are shown in Figure 5. e parameters of the numerical samples using the B-model, the M-model, and the K-model are listed in Tables 2-4, respectively. e parameters can be assigned to the discrete elements using the PFC command "prop pa_name pa_value range id id_DE." "prop," "range," and "id" are the built-in commands of PFC. "pa_name" and "pa_value" are the name and the value of the target parameter. "id_DE" is the ID number of the discrete element.

Laboratory Tests.
e material composition of asphalt concrete used in this study is listed in Table 5.
e dynamic modulus tests of asphalt concrete are carried out at the temperatures of − 6°C, 4°C, 16°C, 27°C, and 38°C and at the frequencies of 0.1 Hz, 0.5 Hz, 1 Hz, 5 Hz, 10 Hz, and 25 Hz, of which the result data are shown in Figure 6. e dynamic modulus-testing apparatus is shown in Figure 7. All the laboratory tests are implemented according to the Chinese testing specification "Standard Test Methods of Bitumen and Bituminous Mixtures for Highway Engineering (JTG E20-2011)," and four samples were tested successfully in each testing condition. e following preparation method (Figure 8) of asphalt concrete samples and asphalt mastic samples is implemented to ensure that the asphalt mastic in asphalt mastic samples is similar to that existing in asphalt concrete samples. Figure 9 plots a diagram of the dynamic modulus numerical test. e same loading condition used in the laboratory tests is applied in the numerical simulations. A certain loading frequency equal to that used in the laboratory test is applied to the upper horizontal wall of the numerical sample in the y-direction to simulate the loading procedure, while the bottom horizontal wall used to support the numerical sample is fixed in all the directions.

Numerical Tests.
To ensure the reliability of simulations, the parallel tests are implemented in numerical simulation just as experiment tests. In experimental parallel tests, the asphalt concrete samples usually have the same raw materials and material composition and the random material distribution. Accordingly, in numerical parallel tests, the numerical models have the same material parameters and material composition and the random microstructure distribution.
Moreover, it should be explained that the same series of parallel numerical samples are employed to simulate the dynamic modulus and phase angle of different types of numerical samples in this study. It is beneficial to compare the difference in the rheological behaviors in the case of using different mesoscopic viscoelastic contact models more clearly without the impact of mesoscopic structures (such as aggregate characteristics and void characteristics) and their distributions. In addition, the influences of the temperatures on the behaviors of asphalt concrete are represented indirectly via the difference in parameters calibrated at different temperatures, as shown in Tables 2-4.      Table 4: Parameters of the numerical samples using the K-model.

Results and Discussion
To evaluate the simulation precision, the relative errors between simulation results and laboratory data are calculated using the following equation: where Error |E| and Error θ are the error of dynamic modulus and phase angle, |E| DEM and |E| Lab are the simulation results and laboratory data of dynamic modulus, and θ DEM and θ Lab are the simulation results and laboratory data of phase angle. Moreover, the curves in the following are given in the linear-log scales to plot the results clearly at low frequencies. Figures 10 and 11 plot the data comparison and error analysis of dynamic modulus.

Dynamic Modulus.
As shown in Figures 10 and 11, the following observations can be made:      Figures 12 and 13 plot the data comparison and error analysis of phase angle.

Phase Angle.
As shown in Figures 12 and 13, the following observations can be made:  at different frequencies are on average 9.25%, 8.69%, 5.05%, 3.51%, and 2.49% at − 6°C, 4°C, 16°C, 27°C, and 38°C, respectively. e maximum errors in the case of using different numerical samples are on average 30.72%, 24.93%, 7.56%, 6.23%, and 5.14% at − 6°C, 4°C, 16°C, 27°C, and 38°C, respectively. In fact, when the temperature exceeds 15°C, the difference in the simulation precision in the case of using the three contact models is not significant. However, the numerical

Concluding Remarks
In this study, the influence of mesoscopic viscoelastic contact models in the DEM on characterizing the rheological behavior of asphalt concrete is investigated based on the comparison between laboratory tests and the corresponding numerical tests via the rheological indexes of dynamic modulus and phase angle at different temperatures (− 6°C, 4°C, 16°C, 27°C, and 38°C) and loading frequencies (0.1 Hz, 0.5 Hz, 1 Hz, 5 Hz, 10 Hz, and 25 Hz). e following observations can be made:    (i) e difference in simulation precision in the case of using the three mesoscopic viscoelastic contact models at different temperatures is more significant than that at different frequencies. (ii) e difference in simulation precision of phase angle in the case of using the three mesoscopic viscoelastic contact models is more obvious than that of the dynamic modulus. (iii) e mesoscopic generalized Maxwell contact model provides the best simulation precision at low temperature. e absolute level of the simulation precision is little changed with the temperature and frequency. (iv) e mesoscopic generalized Kelvin contact model shows an improved simulation precision at high temperatures, and this advantage is more significant with the increased temperature. However, the frequency has less effect on the absolute value of the simulation precision. (v) e simulation precision of the mesoscopic Burgers contact model is obviously lower than that of the other two contact models, particularly when simulating the phase angle at low temperatures and frequencies. However, the simulation precision is acceptable when the temperature exceeds 15°C, particularly for the dynamic modulus. Considering the relatively simple mathematical structure, the mesoscopic Burgers contact model can also be used in the DEM simulations at medium and high temperatures.

Data Availability
e data used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare that they have no conflicts of interest.