Numerical Simulation Virtual Test of Torsion Shear for Asphalt Mixture

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 diﬀerent 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.


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][2][3] that these damages are associated with the insufficient shear performance of the asphalt mixture. erefore, 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]. e 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. us, 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][6][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. ese 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. e flowchart of the research is shown in Figure 1.

Construction of 3D Numerical Model of
Asphalt Mixture e 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. e 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. erefore, 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. erefore, 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.

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. e 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. e main steps are as follows.

Fill a Regular Arrangement of Small Particle Size Discrete
Units.
e 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. e mathematical equation of the cube is as follows: where x 0 , y 0 , and z 0 are the coordinates of cube-shaped centroid O and 2R is the square side length, that is, the particle size of the aggregate.

Generate Eight Random Planes to Cut Cubes Randomly.
First, a normal vector (n j1 , n j2 , n j3 ) is generated for each quadrant cutting plane. e normal vector was generated randomly, where j is equal to 1∼8 and n j1 , n j2 , and n j3 are calculated by where urand is a computer internal random function with a range of (0, 1) and int is the computer internal rounding 2 Advances in Materials Science and Engineering function, wherein the distance d from the centroid of the cube to the cutting plane is represented by where q is the cutting control coe cient, ranging between 0 and 1.0. According to the distance d and the normal vector (n j1 , n j2 , n j3 ), the coordinates of the intersection M c (x c , y c , z c ) of the centroid O and the cutting plane were obtained as follows: en, the plane equation is According to the randomly generated cutting surface, the cube was cut to obtain the irregular polyhedral region. e 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 lled. Use the "BALL" command to ll 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. e area consists of 2,744 regularly arranged small spherical particles, as shown in Figure 2.
Second, set the cutting control coe cient to generate the cutting surface. ird, 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 di erent particle was generated by running once.

Determination of the Number of Aggregate Units in
Each Sieve. Because the passing percentages of the specication 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:  erefore, 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 ρ c , the density of asphalt is ρ l , the asphalt-aggregate ratio is VV, the design air voids is V � (πd 2 h/4), and the volume of the cylindrical specimen is V � (πd 2 h/4), where d is the diameter and h is the height. e volume of the asphalt is set to h, the volume of the coarse aggregate is V c , and the volume of the fine aggregate is V x . e passing percentages of the aggregate are shown in Table 1. e 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, e 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: erefore, from equations (11) and (15), In addition, the density parameter ρ s 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.
e 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%. e 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. e 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.

Arrangement of Gradation of Parent Particles.
e 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.

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. en, 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.
ird, the irregular aggregate generation program was loaded, and a gradation check was performed. e 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. erefore, 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. en, the next aggregate was generated, and an irregular aggregate model was obtained, as shown in Figure 6.

Advances in Materials Science and Engineering
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. e 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. us, a numerical model of the asphalt mixture was obtained. e number of asphalt mortar particles to be removed, n_del, is calculated by where VV is the air voids and n total is the total number of small particles in the model. e 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. e final numerical model of the cylindrical specimen is shown in Figure 8.

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.

Parameters of 3D
Model. Discrete elements method mainly included three basic models: stiffness model, sliding model, and bond model [16,17]. e 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. erefore, 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. erefore, 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. rough research, the paper found a new definition method, which was defined by the principle of density discrimination. e 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. e 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   Advances in Materials Science and Engineering 7 parameters of each aggregate particle were sequentially incremented. erefore, 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. en, 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. en, 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. e 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 .
e mesocontact force f n in the model could be expressed as where ε k , ε mk , and ε mc are the strain of corresponding units. e stress σ in the macro Burger's model could be expressed as where E 1 , η 1 , E 2 , and η 2 are the macrofitting parameters of Burger's model. And, because of the equivalence between contact force and stress f n � σA, there were For the same reason, there were us, we can obtain (e) (f)  Advances in Materials Science and Engineering e relationship between the tangential parameters of Burger's model and the macroscopic parameters could be obtained as follows: us, 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.

Test Conditions.
In the numerical test, the conditions were the same as those in the laboratory test, and the controlled strain mode was employed. e 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.

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. e simulation test results are shown in Table 10. e failure state of the model is shown in Figure 12.
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 thirdparty 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°. e 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]. e results are shown in Table 11, and the failure specimens are shown in Figure 14.
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.

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. e 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.       Figure 11: Model of loading. 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.
e 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   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.

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.
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.

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.
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. e 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.

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. 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.
It can be seen in the figure that, as the loading progressed, the displacement of the specimen became larger and larger. e 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.
rough 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.

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. e numerical simulation of the torsional shear test of the asphalt mixture was carried out, and the following conclusions were obtained: (1) e 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) e 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. e 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. e 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. e difference in displacement at the location where the crack occurred was large, which indicated that the internal structure was discontinuous. e 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
e data used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare that there are no conflicts of interest.

Supplementary Materials
e 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. e 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)