#### Abstract

In order to analyze the torsional shear process of asphalt mixtures in a microscopic view, the numerical simulation of a torsional shear test of an asphalt mixture was carried out by discrete element method. Based on the defects of existing algorithms, the method of random reconstruction of the existing 3D model of the asphalt mixture was improved, and a new reconstruction method was proposed. A 3D numerical model of the asphalt mixture contained irregular-shaped coarse aggregate, mineral gradation, and asphalt mortar; furthermore, the particle algorithm established the air void distribution. Then, the numerical simulation of the asphalt mixture’s torsional shear was completed; in addition, the stress, displacement, and contact of the specimens at each stage were analyzed. The results showed that the stress and displacement in different stages changed greatly with the loading, i.e., the crack generated from a weak point on the surface and then spread to the ends with an oblique angle of about 45°. At the same time, the shear failure process of the asphalt mixture was studied. The virtual test method could accomplish the implementation of the numerical simulation of torsional shear; it also provided a good research method for analysis of the asphalt mixture’s shear failure process.

#### 1. Introduction

With the increasing volume of road traffic and the continuous development of heavy-duty transportation, the ruts and cracking of pavement surface have become the main forms of damage for China’s asphalt pavement. Studies show [1–3] that these damages are associated with the insufficient shear performance of the asphalt mixture. Therefore, it has great significance to analyze the shear failure and shear performance of asphalt pavement.

Researchers have carried out related research, but the existing research has been mainly carried out by means of laboratory tests, which were macroscopic and had the disadvantages of high cost and poor reproducibility [4]. The asphalt mixture was a typical multiphase composite material composed of aggregates, asphalt, and air voids; thus, the internal stress and strain are discontinuous. However, most studies have used continuous homogenization hypothesis, which remains inconsistent with the actual situation. Thus, it proved necessary to develop further research on mechanical properties considering the heterogeneous internal structure of the asphalt mixture in a microscopic view.

Numerical simulation is a method of studying engineering problems through numerical calculation and image displays. In the field of pavement, numerical simulation methods are used to analyze the mesostructure and performance of an asphalt mixture. Numerical analysis methods commonly include the finite element method, the boundary element method, the finite difference method, and the discrete element method. Among them, the finite element method is suitable for a simulation situation when the strain level is small (i.e., small deformation), and the discrete element method has the advantage of simulating large deformation or cracking problems; this method is especially used for the study of discontinuous particulate materials such as asphalt mixtures. For example, the discrete element method was used to simulate the splitting test [5–7], the virtual uniaxial creep test [8], the virtual fatigue test of the trabecular specimen [9], and the prediction of the asphalt mixture dynamic modulus [10], creep stiffness [11], and so on. These studies mainly focused on the splitting, modulus, and fatigue properties of an asphalt mixture, and the content of the shear test simulation of an asphalt mixture was relatively rare.

Based on the related research, the author proposes a new torsional shear test method under normal stress conditions [12]. In order to further analyze the torsional shear development process of the asphalt mixture in the test and analyze the internal variation of the asphalt mixture during the shearing process in a microscopic view, the discrete element method was used to establish the numerical model of the asphalt mixture, and the relevant microscopic parameters were determined; then, the numerical simulation of torsional shear was developed. The flowchart of the research is shown in Figure 1.

#### 2. Construction of 3D Numerical Model of Asphalt Mixture

The internal structure of an asphalt mixture is complex; it includes multiphase composite materials composed of aggregate, asphalt, and air voids. To perform numerical simulations, the primary problem should be solved to establish a numerical model that could reflect the internal characteristics of the asphalt mixture. The traditional method for obtaining the numerical model of an asphalt mixture is to take a photograph with a digital camera of the mixture’s surface or a section of the specimen after cutting; one can also use industrial CT to perform tomographic scanning on the asphalt mixture specimen and obtain a binary image of the sample by digital image processing technologies such as image enhancement technology and image segmentation technology; then, a 2D or 3D numerical model of the mixture specimen can be obtained by digital reconstruction [13]. In this way, a digital model with a high degree of matching with the real specimen could be obtained. However, there are also deficiencies to consider such as high test requirements, high cost, and the inability to implement simulation requirements in some cases. Therefore, some scholars have proposed a new reconstruction generation method using random generation techniques to reconstruct a mixture model [14]; however, this method has defects in the aggregate generation algorithm, e.g., parallel or possible overlap of the cutting faces. Therefore, a new random generation method has been proposed on the base of improvement of the existing random reconstruction algorithm.

According to the relevant research [15], when constructing the numerical model, the part with the specified particle size greater than 2.36 mm was regarded as the aggregate part, and the part with the particle size less than 2.36 mm and the mineral powder and asphalt were regarded as the asphalt mortar, plus the air void structure, which together make up the asphalt mixture.

##### 2.1. Algorithm for Generating Irregular Aggregate Particles

In order to generate irregular aggregate particles, according to the existing research, an irregular polyhedron random cutting algorithm was modified to overcome the defects of the existing algorithms. The irregular polyhedron algorithm was based on a cube; then, the cube centroid was designated the point of origin, and a spatial Cartesian coordinate system was established. Due to the random nature of the cutting plane, the coordinate axes of the space Cartesian coordinate system were parallel to the true coordinate system axis, and a random face was generated in each quadrant for cutting. The main steps are as follows.

###### 2.1.1. Fill a Regular Arrangement of Small Particle Size Discrete Units

The program was programmed to fill with smaller sizes of particles in a specified order in a specified cube area and then set the side length of the area to be equal to the aggregate particle size. The mathematical equation of the cube is as follows:where,, and are the coordinates of cube-shaped centroid *O* and 2*R* is the square side length, that is, the particle size of the aggregate.

###### 2.1.2. Generate Eight Random Planes to Cut Cubes Randomly

First, a normal vector () is generated for each quadrant cutting plane. The normal vector was generated randomly, where *j* is equal to 1∼8 and are calculated by where *urand* is a computer internal random function with a range of (0, 1) and *int* is the computer internal rounding function, wherein the distance *d* from the centroid of the cube to the cutting plane is represented bywhere *q* is the cutting control coefficient, ranging between 0 and 1.0.

According to the distance *d* and the normal vector (), the coordinates of the intersection of the centroid *O* and the cutting plane were obtained as follows:

Then, the plane equation is

According to the randomly generated cutting surface, the cube was cut to obtain the irregular polyhedral region. The mathematical equation of the irregular polyhedral region is shown in the following equation:

Judging the relative position, by writing a program, traversing all the small-sized particles arranged regularly in the square, and judging the relative position of the irregular particle with the irregular polyhedron, the particles in the irregular polyhedral area were the aggregate part.

According to the above algorithm, a single aggregate having a particle diameter of 14 mm would be described as an example.

First, the regular arrangement of small ball particles was filled. Use the “BALL” command to fill the spheres in the cube area according to the rule arrangement algorithm, take the radius of the pellets to 0.5 mm, and set the centroid of the cube (0, 0, 0) with a side length of 14.0 mm. The area consists of 2,744 regularly arranged small spherical particles, as shown in Figure 2.

Second, set the cutting control coefficient to generate the cutting surface. Third, judging by the position of the small sphere particles and the polygonal area, the particles not in the area were deleted, and the remaining were the desired aggregate particle model, as shown in Figure 3. A single different particle was generated by running once.

##### 2.2. Establishment of 3D Numerical Model

###### 2.2.1. Determination of the Number of Aggregate Units in Each Sieve

Because the passing percentages of the specification gradation were calculated by the mass characterization, in order to characterize the gradation characteristics of the mixture in the model, assuming that the density of the aggregates was equal, then the gradation could be characterized by the volume of aggregates, and the number of gradation ball units was calculated according to the volume fraction of each sieve of aggregates; that is, the volume of the specimen occupied by certain aggregates was calculated according to the volume fraction and then divided by the average volume of the ball; thus, the number of the aggregate particles was obtained, as shown in equation (7), wherein the radius of the sphere was half of the value of the average value between the upper and lower limit sizes of the aggregate:where (a cylindrical specimen as an example) is the volume of the specimen, *d* is the diameter of the circle on the bottom surface of the cylinder, *h* is the height of the cylinder, is the average radius of the particles, , and and are the dimensions of the ^{th} and ^{th} sieve sizes.

Therefore, the key to reflect the gradation is to calculate the volume fraction of each aggregate, which is derived as follows:

Assume that the density of the aggregates is , the density of asphalt is , the asphalt-aggregate ratio is , the design air voids is , and the volume of the cylindrical specimen is , where *d* is the diameter and *h* is the height. The volume of the asphalt is set to , the volume of the coarse aggregate is , and the volume of the fine aggregate is . The passing percentages of the aggregate are shown in Table 1. The volume percentages of the aggregate and the asphalt mortar in the asphalt mixture are then calculated.

According to the equation of the asphalt-aggregate ratio calculation, namely,

The content of asphalt in the numerical specimen is as follows:

According to the volume composition of the whole specimen, there are

Substituting equation (9) into equation (10), we obtain

According to the gradation passing percentages table,

And, because the density is equal,

Substituting equation (13) into equation (11), we obtain

Combined with the passing percentages table, the formula for the volume fraction of aggregates of each sieve is obtained as follows:

Therefore, from equations (11) and (15),

In addition, the density parameter of the asphalt mortar could also be derived from equations (8) and (10):

Substituting equation (9) into equation (17), then

Taking the asphalt mixture of AC-13 as an example, the gradation is shown in Table 2. The density of the asphalt is 1.03 g/cm^{3}, the aggregate density is 2.7 g/cm^{3}, the asphalt-aggregate ratio is 5%, and the air voids is 4%. The diameter of the cylindrical specimen is 100 mm, and the height is 100 mm. According to the above formula, the number of particles of each sieve is calculated, as shown in Table 3.

In the asphalt mixture, the aggregates of different particle sizes were divided into different sieve parts, and the proportion of each sieve of the total mass was different, so the gradation was formed. The number of aggregate particles of each sieve was determined; it is equivalent to determining the proportion of the aggregates, so the number of gradation units was obtained.

###### 2.2.2. Arrangement of Gradation of Parent Particles

The number of aggregates for each sieve was calculated according to the above method; then, we put them into the model within the size of the specimen, as shown in Figure 4. By writing a program, scanning particles of each sieve and calculating the corresponding volume fraction, it was found to be consistent with the actual value, which indicates that the arrangement is correct.

##### 2.3. Generation of 3D Numerical Model

First, a program was written to scan the abovementioned gradation ball unit, and information of each aggregate particle was extracted: coordinates (*x*, *y*, *z*) and radius *r* as the geometric parameters of the generated polygon aggregate.

Then, the original particles were deleted, and the particles with smaller diameter particles were regularly filled in the region in a certain order; that is, all the aggregate particles were arranged neatly in the horizontal and vertical directions, and each of the particles in the middle position was arranged adjacent to one particle in the upper, lower, left, right, front, and back directions, as shown in Figure 5.

Third, the irregular aggregate generation program was loaded, and a gradation check was performed. The method of gradation check was calculated according to the percentage of all the small ball particles in each sieve of the total particles in the model. Because the basic unit of the discrete element model is a single sphere, it is difficult to obtain the same value as the volume fraction of the original aggregates by loading the irregular aggregate algorithm. Therefore, considering that, when the percentage of aggregate particles in the entire cylinder model came within the percentage threshold, the gradation test was successful; otherwise, the aggregate was regenerated until the percentage of the aggregate particles in the total cylinder model was within the percentage threshold. Then, the next aggregate was generated, and an irregular aggregate model was obtained, as shown in Figure 6.

In order to create certain air voids for the above model, a random deletion method was used to remove some particles of the asphalt mortar randomly to characterize the air voids of the mixture. The specific method was to generate a random number in the maximum particle range of the irregular aggregate model and determine whether the particle address with the random number as the serial number was empty, and, if not, determine whether the particle with the random number as the serial number was asphalt mortar; if it was, delete the particles until the *n_*del particles were removed and stopped. Thus, a numerical model of the asphalt mixture was obtained. The number of asphalt mortar particles to be removed, *n_*del, is calculated bywhere *VV* is the air voids and is the total number of small particles in the model.

The air void distribution is shown in Figure 7. If the total number of particles in the model is 48040 and the air voids is 4%, the number of particles to be deleted *n_*del is 1921. The final numerical model of the cylindrical specimen is shown in Figure 8.

**(a)**

**(b)**

**(c)**

#### 3. Numerical Simulation of Torsional Shear of Asphalt Mixture

##### 3.1. Generation of 3D Numerical Specimen of Asphalt Mixture

According to the above model establishment method, two kinds of asphalt mixture 3D cylindrical specimen with gradation types AC-13 and SMA-13 were generated (the diameter was 100 mm, the height was 100 mm, and the air voids was 4%). As shown in Figure 9, the yellow part is the aggregate phase, the blue part is the asphalt cement phase, and the blank part is the air void phase.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

##### 3.2. Parameters of 3D Model

Discrete elements method mainly included three basic models: stiffness model, sliding model, and bond model [16, 17]. The stiffness model characterized the material elastic characteristics, the sliding model characterized the slip characteristics, and the bond model characterized the strength characteristics of the material. However, the above three models could not characterize the viscoelastic properties of the asphalt mixture, so Burger’s model was introduced. At the same time, the model could reflect the normal action between the units and should also reflect the tangential action between the units. Therefore, the tangential direction between the units was also considered by Burger’s model, as shown in Figure 10.

In addition, there are four types of contact between the asphalt mixture units: contact between internal units of the same aggregates, between different aggregate units, between internal units of asphalt mortar, and between asphalt mortar units and aggregate units. Among them, the contact parameter between the asphalt mortar unit and the aggregate unit was greatly affected by the aggregate characteristics, also it could not be tested through a good measuring method currently. Therefore, the contact parameter between the asphalt mortar unit and the aggregate unit could be considered the same as that of the internal contact of the asphalt mortar by appropriately simplification.

At the same time, there was no suitable command in the previous study to define the contact between adjacent aggregates uniformly, considering that, under actual conditions, the aggregate surface always had a layer of asphalt film covered, so it was mostly regarded as the contact among the aggregates and the asphalt mortar, then with the aggregates. Through research, the paper found a new definition method, which was defined by the principle of density discrimination. The mortar and each aggregate particle were defined at different densities first; then, different locations of the contact were distinguished and parameters were assigned based on this. The operation process was as follows: When defining the parameters, in the initial stage of imparting the particle density parameter, first, the density parameter of the asphalt mortar particles was defined as 1; then, each aggregate particle was sequentially incremented to define different density parameters. By programming the *fish* function, each aggregate was scanned, and the density parameters of each aggregate particle were sequentially incremented. Therefore, the density parameters of the different aggregate particles were inconsistent and also different from those of the asphalt mortar. Subsequently, we scanned all contacts and returned the two physical addresses corresponding to the contact. Then, we determined whether the densities of the two entities A and B were equal; if they were equal, it was judged whether the density of the entity A was 1; if it was, the contact belonged to the internal contact of the asphalt mortar; if it was not 1, it indicated that the contact belonged to the internal contact of the aggregates. Then, if the densities of the two entities were not equal, it was determined whether the density of the entity A or the entity B was 1; if so, the contact was the contact between the asphalt mortar unit and the aggregate unit; otherwise, it was the contact between different aggregates. In this way, by multilevel conditional statement judgment, parameters could be assigned separately. When the contact parameter assignment was completed, new density parameters were assigned to the asphalt particles and aggregate particles to replace the original parameters. In this way, the problem of assigning the parameters could be solved.

Because the aggregate property in the asphalt mixture was relatively stable and was generally considered to have linear elastic properties, a linear elastic stiffness contact model was adopted for the aggregate, and a point bonding contact model was adopted for the bonding of the aggregate. Because the viscoelastic properties of the asphalt mixture were mainly due to the nature of the asphalt mortar, Burger’s viscoelastic contact model needed to be considered for the contact between the asphalt mortar units. The four contact models are shown in Table 4.

In order to analyze the microscopic parameters of Burger’s model, it can be seen in Figure 9 that this model has a total of eight parameters: *K*_{kn}, *K*_{mn}, *C*_{kn}, and *C*_{mn} and tangential *K*_{ks}, *K*_{ms}, *C*_{ks}, and *C*_{ms}.

The mesocontact force *f*_{n} in the model could be expressed aswhere , , and are the strain of corresponding units.

The stress in the macro Burger’s model could be expressed aswhere , , , and are the macrofitting parameters of Burger’s model.

And, because of the equivalence between contact force and stress , there were

For the same reason, there were

Thus, we can obtain

The relationship between the tangential parameters of Burger’s model and the macroscopic parameters could be obtained as follows:

Thus, as long as the corresponding viscoelastic parameters and linear elastic parameters and the bonding parameters between them were obtained, the microscopic parameters of the numerical model could be obtained. For the four macroscopic parameters of Burger's viscoelastic model of asphalt mortar, the macroscopic parameters could be obtained through experiments. According to the relevant research and the results of the literature [14], the macroscopic parameters of Burger’s model of asphalt mortar of AC-13 and SMA-13 were obtained (Table 5).

Assuming the aggregate elastic modulus was 50 GPa, Poisson’s ratio was 0.25, the tensile strength was 6 MPa, and the shear strength was 15 MPa. According to the above analysis, the normal stiffness and tangential stiffness of the aggregate and the normal strength and tangential strength of the point bond could be obtained. Finally, the microscopic parameters of the asphalt mixture were calculated, as shown in Tables 6∼9 (temperature was 60°C), where *kn* and *ks* are the normal and tangential stiffness of the contact unit, respectively, and *n_*bond and *s_*bond are the normal and tangential stiffness of the contact, respectively.

##### 3.3. Test Conditions

In the numerical test, the conditions were the same as those in the laboratory test, and the controlled strain mode was employed. The specific loading mode uses the two upper and lower loadings “clump” to form two loading plates, fix the bottom of the test piece by fixing the lower loading plate, and give the upper loading plate a constant angular velocity *ω*. A high-strength parallel bond contact model was used between the two loading plates and the specimen to simulate the bonding of the actual specimen with epoxy resin (the loading method is shown in Figure 11). By programming the *fish* function, the unbalanced torque of the loading plate was recorded in every time step during the period of testing; in addition, shear stress also was calculated.

##### 3.4. Torsional Shear Numerical Simulation Test Results

According to the above conditions, the corresponding parameters were set and the torsional numerical simulation tests were carried out on the two types of cylindrical asphalt mixture discrete element models. Four sets of parallel tests were completed. The simulation test results are shown in Table 10. The failure state of the model is shown in Figure 12.

**(a)**

**(b)**

According to the results, the shear strength of AC-13 is 0.360 MPa and the shear strength of SMA-13 is 0.506 MPa; furthermore, it shows that shear performance of SMA-13 is better than that of AC-13 under the same conditions. In order to observe the crack shape more intuitively, the third-party postprocessing of the test was carried out, and the failure cloud images are shown in Figure 13. It can be seen in the figure that the specimens of the two gradation types have similar crack modes, basically occurring along an angle of about 45°.

**(a)**

**(b)**

The torsional shear test of AC-13 and SMA-13 actual specimens under the same conditions was carried out by the method of the literature [12]. The results are shown in Table 11, and the failure specimens are shown in Figure 14.

**(a)**

**(b)**

It could be found that the simulation test results were basically consistent with the real results, and the relative errors were less than 5%. At the same time, from Figures 13 and 14, the failure cracks also show similar laws between the numerical model and the laboratory specimens; in addition, they occurred all along an angle of about 45°, which indicates that the numerical simulation virtual test method is feasible.

#### 4. Analysis of Torsional Shear Numerical Simulation Test Process

To analyze the failure process, the asphalt mixture of AC-13 was used as an example to analyze the stress, displacement, and contact changes at various stages of the test. The tensile stress and the compressive stress are set to red and black, respectively, and the depth of the color indicated the magnitude of the stress; further, the displacement is represented by a black arrow, similar to the vector, the direction of which is the direction of the arrow, and the size is represented by the length of the arrow.

##### 4.1. Initial Stage

At the start period of the test, the distribution of stress and contact conditions at the beginning of the loading were collected by writing a program, as shown in Figure 15.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

The vertical section graph was obtained by cutting along a plane half the width direction (*x* direction) of the cylindrical specimen, and the transverse section graph was obtained by cutting laterally along the half of the height direction (*z* direction) of the specimen. In fact, the cutting observation of different angles and different positions could be performed by modifying the unit normal vector of the cutting plane, and the cutting observation of different angles and different positions could be performed through one point of the cutting plane. Because the number of pictures was too large, for the sake of convenience of observation, only the pictures of the positions in the width direction and the height direction 1/2 were placed here for analysis. The position described later is the same as this.

As shown in Figures 15(a)–15(c), the compressive stress and the tensile stress are only obvious at the top of the specimen and at the top of the loading plate; they had not yet occurred in other positions, which means that the load had not been transferred, which is consistent with the actual situation.

As shown in Figures 15(d)–15(f), the contact is good at this time, and there is no evidence of damage. At the same time, it can be seen that some positions at the contact appear as blanks, which is the position of air voids.

##### 4.2. Earlier Stage of Loading Period

After running a certain number of steps, that is, after the loading plate rotated a certain displacement (the early loading stage), the corresponding stress, contact, and displacement were also collected according to the programmed procedure, as shown in Figure 16.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

**(i)**

As shown in Figures 16(a)–16(c), after running a certain number of steps, the tensile and compressive stresses had changed greatly at this time; further, they had been transmitted to the bottom of the specimen, and the distribution became uneven. At the same time, it was found that the stress mainly occurred between the aggregates, which indicated that the skeleton of the aggregates in the asphalt mixture was more obvious, which is consistent with the actual situation.

As shown in Figures 16(d)–16(f), although the particle contact changed slightly at this time, the disconnection had not occurred, which indicates that the stress did not reach its failure strength, and the specimen could continue to bear the load.

As shown in Figures 16(g)–16(i), the model particles had a certain displacement and exhibited a rotation form consistent with the loading direction, and the outer edge particles had a larger displacement than the inner edge particles; further, the upper particles are more displaced than the lower particles. It can also be seen from the transverse section graph that the displacement rotation is not a completely consistent standard rotation form, which indicates that the interior of the mixture is not uniform.

##### 4.3. Later Stage of Loading Period (Appearance of Crack)

As the upper loading plate continued to rotate, cracks began to occur at the weak position of the specimen, and various changes of the specimen were collected at this time, as shown in Figure 17.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

**(i)**

As shown in Figures 17(a)–17(c), as the loading plate continues to twist, the stress inside the model increases significantly; at the same time, there are obvious cracks at some positions, and cracks are generated at the position where there is no stress, which indicates that the specimen experienced damage.

As shown in Figures 17(d)–17(f), the contact between the particles at some positions is broken or greatly changes at this time, which indicates that cracks have been generated there.

As shown in Figures 17(g)–17(i), the particles in the specimen have large displacement and still exhibit a rotation change consistent with the loading direction, while the displacement of outer edge particles is more than that of the inner edge. The displacement of the upper particles is larger than that of the lower particles, and there is a large displacement difference at the local position, which indicates that cracks are generated here.

##### 4.4. Ending Stage of Loading Period (Model Failure)

As the upper loading plate continued to rotate, the crack gradually expanded. According to the collected stress, when it dropped rapidly, it could be considered that the specimen had been destroyed, as shown in Figure 18.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

**(i)**

**(j)**

**(k)**

**(l)**

As shown in Figures 18(a)–18(c), as the loading plate continues to rotate, the stress inside the model is redistributed after the crack propagates. At the same time, the angle of the crack is about 45° along the specimen.

As shown in Figures 18(d)–18(f), the change in contact at the crack position before this time is larger and the distance farther, which indicates that the crack had developed to a larger width.

As shown in Figures 18(g)–18(i), the particles inside the specimen have larger displacement, and the position of the crack is already obvious. At the same time, the displacement of the outer edge particles is larger than that of the inner edge particles, and the displacement of the upper particles is larger than that of the lower particles. In addition, it is found that the displacement difference at the position where the crack appeared is particularly large, which indicates that it is discontinuous.

As shown in Figures 18(j)–18(l), after the specimen is loaded, an obvious crack has appeared; the width of the crack is large, and the crack is formed along the oblique direction of the specimen by 45°.

In order to observe the development process of the crack, the test results were exported at regular intervals in the calculation process, and the postcorrelation process was used to obtain the cloud image at different stages, as shown in Figure 19.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

It can be seen in the figure that, as the loading progressed, the displacement of the specimen became larger and larger. The shear failure process indicates that the crack was generated from a weak position on the surface and then spread to both ends at an oblique angle of about 45°; at the same time, the crack gradually spread to the inside of the specimen, resulting in a crack similar to the spiral shape, up to the complete destruction of the specimen. It was observed that the failure surface was not smooth; this was also caused by the unevenness of the mixture.

Through the analysis of the internal conditions of the specimen at different testing stages, the whole process of the asphalt mixture under the torsion action was obtained during the virtual experiment. It could be considered that the numerical simulation virtual test method had good simulation of the torsional shear process of the asphalt mixture; further, it was reasonable and feasible and also provided a good research method for the analysis of asphalt mixture failure process.

#### 5. Conclusions and Recommendations

In order to analyze the torsional shear process of the asphalt mixture in a mesoscopic view, the discrete element method was used to establish the numerical model of an asphalt mixture. The numerical simulation of the torsional shear test of the asphalt mixture was carried out, and the following conclusions were obtained:(1)The existing 3D model random reconstruction method of the asphalt mixture was improved, and a new reconstruction method was proposed, which was an algorithm for obtaining irregular aggregate particles by using a random plane cutting regular hexahedron in eight different quadrants. Based on this and combined with the gradation characteristics of the mixture, a 3D numerical model of an asphalt mixture contained irregular-shaped coarse aggregate, mineral aggregate gradation, and asphalt mortar, and air void distribution was established.(2)The numerical simulation of the torsional shear test of an asphalt mixture was completed, and the stress, displacement, and contact of the specimens at each stage were analyzed. It was found that the stress and displacement of the model showed a large change at different stages. The specific performance was as follows: As the loading progressed, the stress was quickly transmitted to the whole model, getting larger and larger, and there was a phenomenon of stress redistribution after the crack was generated. The displacement of the model particles was generally performed in a rotating mode along the loading direction, the outer edge particles moved larger than the inner edge particles, and the upper particles moved larger than the lower particles. The difference in displacement at the location where the crack occurred was large, which indicated that the internal structure was discontinuous. The shear failure process of the asphalt mixture was that the crack generated from a weak position on the surface, then spread to the both ends at an oblique angle of about 45°, and gradually spread to the inside of the specimen to produce a crack similar to the spiral shape. Finally, as the loading continued, the specimen was completely destroyed.(3)According to the research results, the numerical simulation of the torsional shear test of an asphalt mixture could be realized by the virtual test method, and it provided a good research method for analyzing the failure process of the specimen. In future work, a more refined model could be established, and various load and temperature conditions could be considered to better study the mechanical properties of the asphalt mixture.

#### Data Availability

The data used to support the findings of this study are included within the article.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest.

#### Acknowledgments

The authors would like to acknowledge the financial support of the National Natural Science Foundation of China (Grant no. 51008039) and the financial support of China Scholarship Council (Grant no. 201808430101).

#### Supplementary Materials

The supplementary material is the code of the discrete element model according with the editor’s and reviewer’s comments of 3rd, it cannot be attached as an appendix behind the manuscript because it had too many pages. The section 1 on page 1–5 of the supplementary material is the discrete source program about the “Algorithm for generating irregular aggregate particles,” and section 2 on page 6–36 is the discrete source program about “Generating 3D models,” “Torsional shear numerical simulation test,” and “Loading setting.”* (Supplementary Materials)*