#### Abstract

AC-13 asphalt mixture was taken as the research object to investigate the evolution and distribution laws of force chains. A digital specimen of AC-13 asphalt mixture was reconstructed using the discrete element method (DEM) to simulate the simple performance test (SPT). Next, the force chain information among aggregate particles was extracted to analyze the evolution, probability distribution, and angle distribution of force chains. The results indicate that the AC-13 mesoscopic model reconstructed using the DEM is feasible to simulate the mesoscopic mechanical properties of asphalt mixture by comparing the predicted results and laboratory test results. The spatial distributions of force chains are anisotropic. The probability distributions of normal force chains varying with the loading times are consistent. Furthermore, the probability distribution has the maximum value at the minimum (the ratio of contact force to mean contact force); the peak value appears again at = 1.75 and then gradually decreases and tends to be stable. In addition, the angle distributions of force chains mainly locate near 90° and 270°, and the proportions of strong force chains are slightly greater than 50%, but the maximum proportion is only 51.12%.

#### 1. Introduction

Asphalt mixture is a multiphase composite material, which comprises aggregates, asphalt binder, and air voids. Aggregates that dominate by volume and mass in asphalt mixture are considered as the granular media, which are the main material to form a skeleton structure and resist external loading applied to the asphalt mixture. Therefore, the asphalt mixture exhibits the essential attribute of granular media in a certain sense. And there exists the interaction force among the contacting particles that forms the loading transmission paths in a granular system, so Bouchaud et al. proposed the concept of force chain after observing stress distribution, propagation, and arching within granular media [1, 2]. Currently, the research methods on force chains focus on the following areas: mechanical analysis, laboratory test, and numerical simulation.

Cates et al. argued that fragility was linked to the marginal stability of force chains within granular matter and described that the force chains propagated and fluctuated with scalar model and tensorial model in granular media [3–5]. Tordesillas et al. viewed the force chain buckling in a constrained granular medium from the structural mechanics, investigated the force chain evolution for cohesionless granular systems, and presented a regularized two-dimensional model for the force chain buckling [6–8]. Howell et al. carried out the experiments on a granular system consisting of photoelastic polymer disks subjected to shear and compressive loading and measured and calculated the probability distributions of normal and tangential force chains [9, 10]. Yi et al. studied the force chain transmission in a hexagonal-close-packed array granular system with and without point defect submitted to a concentrated loading using the aluminum-plastic board and carbon paper technique [11, 12]. Tordesillas et al. examined the relationship between force cycles and force chain in a dense granular material under quasistatic biaxial loading using a discrete element simulation and analyzed the spatial patterns of force chain buckling with multiscale characterization [13, 14]. Sun et al. proposed the criteria of strong force chain and angle and analyzed the evolution of force chain structure morphology through simulating the uniaxial compression and biaxial test using the discrete element method [15, 16]. You et al. reconstructed the mesoscopic model of asphalt mixture and provided the force chain distributions within aggregates, within mastic, and between aggregates and mastic [17, 18].

It can be found that the researches on force chains concentrate on the granular material with regular morphology, single particle size, and smooth edges by the literature cited above. But few of them are involved in bond effect. However, the asphalt mixture comprises the aggregates with different morphology, geometric shapes, and broad particle size, and the asphalt binder with bond effect. The interaction between aggregates and asphalt binders changes with time and temperature that makes the anisotropic interaction more complex. Furthermore, its spatial and temporal evolution and distribution laws of force chains have not yet reported.

#### 2. Determination of DEM Viscoelastic Mesoscopic Parameters

The SPT dynamic modulus test was developed to measure the dynamic moduli and phase angles of AC-13 asphalt mixture and mastic specimens. Burger’s model parameters were calculated based on the test results, and the mesoscopic parameter values were determined by verifying the corresponding relationship between Burger’s model parameters and DEM mesoscopic parameters. Aggregates are limestone, asphalt binder is SK70^{#}, and their properties are listed in Tables 1 and 2. The asphalt mastic was obtained from AC-13 asphalt mixture’s aggregate gradation by eliminating all the aggregates bigger than 2.36 mm except the asphalt binder. The gradations of AC-13 asphalt mixture and mastic are listed in Table 3. Optimum asphalt contents of AC-13 asphalt mixture and mastic are 4.49% and 11.92%, respectively.

The dynamic moduli and phase angles of asphalt mixture and mastic were measured at three temperatures of 5°C, 20°C, and 40°C and five frequencies of 1 Hz, 5 Hz, 10 Hz, 20 Hz, and 25 Hz. 5°C and 20 Hz were selected to analyze the mesoscopic responses of asphalt mixture and the characteristics of force chains. And the test results were used to calculate the mesoscopic parameters of asphalt mixture viscoelastic model. The mesoscopic parameters include stiffness within aggregates and adjacent aggregates, stiffness and viscosities within asphalt mastic, and stiffness and viscosities adjacent to asphalt mastic as well as stiffness and viscosities between aggregates and mastic.

#### 3. Reconstruction of Digital Specimen and Verification of Mesoscopic Model

The asphalt mixture specimen prepared by gyratory compaction was drilled core to obtain a new specimen with diameter of 100 mm, which was cut along the diameter direction, and then a cross-sectional image (as shown in Figure 1(a)) was received. The image was processed using the digital image processing technique (Image-Pro Plus 6.0, IPP 6.0), its pixel coordinates were extracted and the DEM was imported to generate a digital specimen of asphalt mixture combining Matlab. The whole digital specimen has 50,500 particles, and so many particles affect the computational efficiency. Consequently, this paper captured part of the digital specimen to conduct a numerical simulation, and the digital specimen with 3060 particles was demonstrated in Figure 1(b). The red parts represent the aggregates, the blue parts represent the asphalt mastic, and the upper and lower parts represent the loading platens. The Haversine loading is applied by the top loading platen and the bottom loading platen is served as a fixed base, which are consistent with the SPT laboratory test.

It is necessary to verify the developed DEM model’s reasonability and the consistency between simulation and laboratory results. The frequency applied by the top loading platen was 20 Hz, and the time-step was 1.285 × 10^{−8} s/step. The stresses and strains were extracted after running for 0.1 s (two periods), and the stress and strain curves were plotted in Figure 2. Dynamic moduli and phase angles from the DEM simulation were calculated by where and are the maximum and minimum values of stresses; and are the corresponding maximum and minimum strains; is the time difference between the adjacent peak stress and strain; and is loading period and 20 Hz corresponds to the loading period of 0.05 s.

The predicted dynamic modulus and phase angle were calculated according to the stresses and strains and (1). It is found that the difference of dynamic modulus between measured and predicted results is 0.784 GPa, and its error is 2.9%. Besides, the difference of phase angle between measured and predicted results is 1.8°. Therefore, there are slight differences between measured and predicted dynamic moduli and phase angles. The reconstructed DEM model has the ability to analyze the mesoscopic responses of AC-13 asphalt mixture.

#### 4. Results and Discussion

##### 4.1. Force Chain Evolution

There were produced mesoscopic responses in aggregates, in mastic, and at the interfaces of aggregates/mastic after applying the Haversine loading at 20 Hz for the AC-13 digital specimen. The anisotropic force chain networks directly reflected the mesoscopic responses within the asphalt mixture and the force chain evolution internal aggregates and internal mastic and at the interfaces of aggregates/mastic. The force chain networks were traced at four loading times (0.125 × 10^{−3} s, 0.25 × 10^{−3} s, 0.375 × 10^{−3} s, and 0.5 × 10^{−3} s, which corresponded to 9728 steps, 19455 steps, 29184 steps, and 38910 steps) and were described in Figure 3. The blue parts represent compressive force chain networks and the red parts represent tensile force chain networks.

**(a) 0.125 × 10−3 s**

**(b) 0.25 × 10−3 s**

**(c) 0.375 × 10−3 s**

**(d) 0.5 × 10−3 s**

Figure 3 shows that the spatial distributions of force chains are anisotropic under the Haversine loading. The thicker lines represent the greater force chains, which are mainly the compressive force chains in vertical direction and the main loading transmission paths. Moreover, the force chains that transfer the smaller loading in the other directions are the tensile force chains and these smaller force chain networks distribute widely. The maximum contact forces are 2.970 × 10^{4} N, 3.332 × 10^{4} N, 3.611 × 10^{4} N, and 3.958 × 10^{4} N at four loading times indicating that their acting forces among particles enhance as the loading time grows. This is attributed to the sensitivity of particles inside the specimen to the Haversine loading while it transfers from top to bottom. The interlocked action among particles strengthens, and the particle system adapts to the loading by means of self-organized behaviors. In addition, the force chains within aggregates are all the compressive force chains (blue parts), and the tensile force chains (red parts) are all within mastic and at the interfaces of aggregate/mastic. The reason is that the unbroken aggregates form the skeleton structure subjected to and transferring the Haversine compressive forces, while the Haversine loading produces the stripping trend between asphalt and aggregates within mastic and at the interfaces of aggregate/mastic.

##### 4.2. Probability Distribution of Force Chains

It is difficult to excavate the quantitative distribution information of force chain networks due to their spatial and temporal heterogeneity. Currently, the macroscopic statistical characteristics of force chains are the main way to describe the mesoscopic information of force chain networks. Probability distribution of force chains is conducted as a macroscopic quantitative indicator, where is the ratio of contact force to mean contact force. The probability distributions of normal force chains is provided in Figure 4 at four loading times.

As can be seen in Figure 4, the probability distributions of normal force chains have no significant differences at four loading times and the maximum probability distribution ( = 0.2495) at 0.5 × 10^{−3} s is greater than the other three loading times. decreases with the increase of when ≤ 0.75; is proportional to when , and the peak occurs at = 1.75; is inversely proportional to when > 1.75 and it remains unchanged basically when ≥ 3.0. The variation laws of illustrate that the probability distributions of smaller normal compressive and tensile force chains are greater, and the probability distribution of larger normal compressive force chains is the largest at = 1.75.

##### 4.3. Angle Distribution of Force Chains

Force chains possess great sensitivity and uncertainty due to their anisotropic spatial extension directions depending on the magnitude and time of external loading. The angles between the force chain direction and horizontal direction are served as the angle distribution quantitative indicators. A statistical analysis is employed to obtain the angle distribution characteristics of force chains in 0°–360° (as shown in Figure 5) by extracting the normal and shear force chains.

**(a) 0.125 × 10−3 s**

**(b) 0.25 × 10−3 s**

**(c) 0.375 × 10−3 s**

**(d) 0.5 × 10−3 s**

Figure 5 presents that these angle distributions of force chains mainly reside nearby 90° and 270°; the angle distribution proportions in 0°–180° are significantly greater than those in 180°–360°, which relate to the Haversine loading directions. The asphalt mixture specimen is applied at the vertical Haversine loading that mainly transmits along the vertical direction. But, the asphalt mixture as a composite material comprises aggregate particles with different particle sizes and geometric shapes, so that the angle distributions mainly distribute along the vertical direction while they extend along the other directions.

Force chains can be distinguished into strong and weak force chains according to the force chain thickness and the magnitude of . They are the strong force chains when ≥ 1, while they are the weak force chains when < 1. It is necessary to analyze the proportions of strong force chains at different loading times (shown in Figure 6) due to the strong force chains bearing most of the loading and as the main loading transmission paths.

There are no significant differences in the proportions of strong force chains in Figure 6. The maximum proportion of strong force chains is 51.12% and the minimum proportion of strong force chains is 50.68%, so the maximum difference is 0.44%. The results mentioned in Figure 6 illustrate that the proportions of strong force chains have no changes with the increase of loading time basically after the loading transmits the whole specimen. The strong force chains support most of external loading and become the main parts used to ensure the stability of the whole granular system. Therefore, the interlocked effect from aggregates, especially coarse aggregates, constitutes the skeleton structure, and the aggregates transfer larger loading to form the strong force chain network structure in AC-13 asphalt mixture. However, the fine aggregates and mineral powder particles play a role of filling effect, which are subjected to smaller or zero loading to form the weak force chain network structure.

#### 5. Conclusions

The differences between the measured and predicted dynamic modulus with the error of 2.9% and phase angle of 1.8° show that the developed viscoelastic model using the discrete element method is feasible. It is suitable to investigate mesoscopic responses of the AC-13 asphalt mixture.

Force chains have a heterogeneous distribution in asphalt mixture. The acting forces among the contacting particles enhance as the loading time extends, the force chains within aggregates are all the compressive force chains, and the force chains within mastic and at the interfaces of aggregate/mastic are mainly the tensile force chains.

The probability distributions of normal force chains are basically consistent for different . The probability distributions decrease firstly and then increase. They decrease after reaching the peak at = 1.75 and then stabilize.

The angle distributions of force chains mainly reside nearby 90° and 270°, and the angles distributions in 0°–180° are significantly greater than those in 180°–360°. The proportions of strong force chains are slightly greater than 50%, but the maximum proportion is only 51.12%.

In this paper, we conduct the quantitative analysis of force chains in asphalt mixture under the conditions of 5°C and 20 Hz. In the further study, the evolution and distribution of force chains with different skeleton types can be analyzed under different temperatures and frequencies.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

#### Acknowledgments

This work was supported by the National Natural Science Foundation of China (nos. 51408047, 51378073, and 51408043) and the Special Fund for Basic Scientific Research of Central College of Chang’an University (nos. 310831161002, 310821153502, and 310821152003). The authors gratefully acknowledge their financial support.