# Numerical Simulation Virtual Test of Torsion Shear for Asphalt Mixture

Authors: Jun Xie; Yougang Yang
Journal: Advances in Materials Science and Engineering
 (2019)

---

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


---

## Body

## 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.Figure 1
The flowchart of the research.


## 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:(1)x−x0−R≥0,x−x0+R≤0,y−y0−R≥0,y−y0+R≤0,z−z0−R≥0,z−z0+R≤0,wherex0,y0, and z0 are the coordinates of cube-shaped centroid O and 2R 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 (nj1,nj2,nj3) is generated for each quadrant cutting plane. The normal vector was generated randomly, where j is equal to 1∼8 andnj1,nj2,andnj3 are calculated by (2)α=urand×π2.0+j−1.0×π2.0,β=urand×π2.0+intj−1.04.0×π2.0,t1=cosα,t2=sinα,t3=cosβ,tt=t12+t22+t32,nj1=t1tt,nj2=t2tt,nj3=t3tt,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 by(3)d=R×q,where q is the cutting control coefficient, ranging between 0 and 1.0.According to the distanced and the normal vector (nj1,nj2,nj3), the coordinates of the intersection Mcxc,yc,zc of the centroid O and the cutting plane were obtained as follows:(4)xc=x0+d×nj1,yc=y0+d×nj2,zc=z0+d×nj3.Then, the plane equation is(5)Pj=x×nj1+y×nj2+z×nj3−xc×nj1+yc×nj2+zc×nj3=0.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:(6)x−x0−R≥0,x−x0+R≤0,y−y0−R≥0,y−y0+R≤0,z−z0−R≥0,z−z0+R≤0,Pj≤0.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 Figure2.Figure 2
Regular arrangement particles.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 Figure3. A single different particle was generated by running once.Figure 3
Irregular aggregate particles model.


### 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:(7)ni=V×JDi4/3πri3,where V=πd2h/4 (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, ri is the average radius of the particles, ri=Di−1+Di/4, and Di−1 and Di are the dimensions of the i−1th and ith 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ρc, the density of asphalt is ρl, the asphalt-aggregate ratio is VV, the design air voids is V=πd2h/4, and the volume of the cylindrical specimen is V=πd2h/4, where d is the diameter and h is the height. The volume of the asphalt is set to h, the volume of the coarse aggregate is Vc, and the volume of the fine aggregate is Vx. 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.Table 1
Percentages of aggregate gradation.
Sieve size (mm)
D
1
D
2
D
3
…
D
n−2
D
n−1
D
n = 2.36 mm
<2.36
Passing percentages (%)
100
P
D2
P
D3
…
P
Dn−2
P
Dn−1
P
2.36
…
Individual percent retain on each sieve (%)
0
100 −PD2
P
D2 − PD3
…
P
Dn−3 − PDn-2
P
Dn−2 − PDn−1
P
Dn−1 − P2.36
P
2.36According to the equation of the asphalt-aggregate ratio calculation, namely,(8)a=mlmc+mx=ρlVlρcVc+ρcVx.The content of asphalt in the numerical specimen is as follows:(9)Vl=aρcρlVc+Vx.According to the volume composition of the whole specimen, there are(10)V1−VV=Vl+Vc+Vx.Substituting equation (9) into equation (10), we obtain(11)Vc+VxV=1−VV1+aρc/ρl.According to the gradation passing percentages table,(12)mxmc=P2.36×mx+mc100−P2.36×mx+mc=P2.36100−P2.36.And, because the density is equal,(13)VxVc=P2.36100−P2.36.Substituting equation (13) into equation (11), we obtain(14)VcV=1−VV100−P2.361+aρc/ρl×100.Combined with the passing percentages table, the formula for the volume fraction of aggregates of each sieve is obtained as follows:(15)JDi=Vc+VxV×PDi−1−PDi100.Therefore, from equations (11) and (15),(16)JDi=PDi−1−PDi100×1−VV1+aρc/ρl.In addition, the density parameterρs of the asphalt mortar could also be derived from equations (8) and (10):(17)ρs=msVs=ρlVl+ρxVxVl+Vx.Substituting equation (9) into equation (17), then(18)ρs=a+1+aVx/Vcaρc/ρl+1+aρc/ρlVx/Vcρc.Taking the asphalt mixture of AC-13 as an example, the gradation is shown in Table2. The density of the asphalt is 1.03 g/cm3, the aggregate density is 2.7 g/cm3, 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.Table 2
Gradation of the asphalt mixture of AC-13.
Sieve size
16
13.2
9.5
4.75
2.36
1.18
0.6
0.3
0.15
0.075
Passing percentage (%)
Upper limit
100
100
85
56
38
28
22
17
12
8
Lower limit
100
95
72
42
28
20
12
8
7
4
Median
100
97.5
78.5
49
33
24
17
12.5
9.5
6Table 3
Aggregate particles of the asphalt mixture of AC-13.
Sieve size
16
13.2
9.5
4.75
2.36
1.18
0.6
0.3
0.15
0.075
Passing percentages (%)
100
97.5
78.5
49
33
24
17
12.5
9.5
6
Individual percent retain on each sieve (%)
0
2.5
19
29.5
16
—
—
—
—
—
Volume fraction
0.000
0.021
0.161
0.250
0.136
—
—
—
—
—
Number
0
10
165
1036
4540
—
—
—
—
—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 Figure4. 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.Figure 4
Arrangement of aggregate particle.


### 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 Figure5.Figure 5
Regular arrangement particle.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 Figure6.Figure 6
Specimen of irregular aggregate (no air voids).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 then_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 by(19)n_del=intVV∗n_total,where VV is the air voids and n_total is the total number of small particles in the model.The air void distribution is shown in Figure7. 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.Figure 7
Distribution of air voids. (a) Entire model. (b) Vertical section view. (c) Transverse section view.Figure 8
Numerical model of irregular aggregate.
(a)
(b)
(c)


## 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:(1)x−x0−R≥0,x−x0+R≤0,y−y0−R≥0,y−y0+R≤0,z−z0−R≥0,z−z0+R≤0,wherex0,y0, and z0 are the coordinates of cube-shaped centroid O and 2R 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 (nj1,nj2,nj3) is generated for each quadrant cutting plane. The normal vector was generated randomly, where j is equal to 1∼8 andnj1,nj2,andnj3 are calculated by (2)α=urand×π2.0+j−1.0×π2.0,β=urand×π2.0+intj−1.04.0×π2.0,t1=cosα,t2=sinα,t3=cosβ,tt=t12+t22+t32,nj1=t1tt,nj2=t2tt,nj3=t3tt,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 by(3)d=R×q,where q is the cutting control coefficient, ranging between 0 and 1.0.According to the distanced and the normal vector (nj1,nj2,nj3), the coordinates of the intersection Mcxc,yc,zc of the centroid O and the cutting plane were obtained as follows:(4)xc=x0+d×nj1,yc=y0+d×nj2,zc=z0+d×nj3.Then, the plane equation is(5)Pj=x×nj1+y×nj2+z×nj3−xc×nj1+yc×nj2+zc×nj3=0.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:(6)x−x0−R≥0,x−x0+R≤0,y−y0−R≥0,y−y0+R≤0,z−z0−R≥0,z−z0+R≤0,Pj≤0.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 Figure2.Figure 2
Regular arrangement particles.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 Figure3. A single different particle was generated by running once.Figure 3
Irregular aggregate particles model.


## 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:(1)x−x0−R≥0,x−x0+R≤0,y−y0−R≥0,y−y0+R≤0,z−z0−R≥0,z−z0+R≤0,wherex0,y0, and z0 are the coordinates of cube-shaped centroid O and 2R 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 (nj1,nj2,nj3) is generated for each quadrant cutting plane. The normal vector was generated randomly, where j is equal to 1∼8 andnj1,nj2,andnj3 are calculated by (2)α=urand×π2.0+j−1.0×π2.0,β=urand×π2.0+intj−1.04.0×π2.0,t1=cosα,t2=sinα,t3=cosβ,tt=t12+t22+t32,nj1=t1tt,nj2=t2tt,nj3=t3tt,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 by(3)d=R×q,where q is the cutting control coefficient, ranging between 0 and 1.0.According to the distanced and the normal vector (nj1,nj2,nj3), the coordinates of the intersection Mcxc,yc,zc of the centroid O and the cutting plane were obtained as follows:(4)xc=x0+d×nj1,yc=y0+d×nj2,zc=z0+d×nj3.Then, the plane equation is(5)Pj=x×nj1+y×nj2+z×nj3−xc×nj1+yc×nj2+zc×nj3=0.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:(6)x−x0−R≥0,x−x0+R≤0,y−y0−R≥0,y−y0+R≤0,z−z0−R≥0,z−z0+R≤0,Pj≤0.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 Figure2.Figure 2
Regular arrangement particles.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 Figure3. A single different particle was generated by running once.Figure 3
Irregular aggregate particles model.


## 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:(7)ni=V×JDi4/3πri3,where V=πd2h/4 (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, ri is the average radius of the particles, ri=Di−1+Di/4, and Di−1 and Di are the dimensions of the i−1th and ith 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ρc, the density of asphalt is ρl, the asphalt-aggregate ratio is VV, the design air voids is V=πd2h/4, and the volume of the cylindrical specimen is V=πd2h/4, where d is the diameter and h is the height. The volume of the asphalt is set to h, the volume of the coarse aggregate is Vc, and the volume of the fine aggregate is Vx. 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.Table 1
Percentages of aggregate gradation.
Sieve size (mm)
D
1
D
2
D
3
…
D
n−2
D
n−1
D
n = 2.36 mm
<2.36
Passing percentages (%)
100
P
D2
P
D3
…
P
Dn−2
P
Dn−1
P
2.36
…
Individual percent retain on each sieve (%)
0
100 −PD2
P
D2 − PD3
…
P
Dn−3 − PDn-2
P
Dn−2 − PDn−1
P
Dn−1 − P2.36
P
2.36According to the equation of the asphalt-aggregate ratio calculation, namely,(8)a=mlmc+mx=ρlVlρcVc+ρcVx.The content of asphalt in the numerical specimen is as follows:(9)Vl=aρcρlVc+Vx.According to the volume composition of the whole specimen, there are(10)V1−VV=Vl+Vc+Vx.Substituting equation (9) into equation (10), we obtain(11)Vc+VxV=1−VV1+aρc/ρl.According to the gradation passing percentages table,(12)mxmc=P2.36×mx+mc100−P2.36×mx+mc=P2.36100−P2.36.And, because the density is equal,(13)VxVc=P2.36100−P2.36.Substituting equation (13) into equation (11), we obtain(14)VcV=1−VV100−P2.361+aρc/ρl×100.Combined with the passing percentages table, the formula for the volume fraction of aggregates of each sieve is obtained as follows:(15)JDi=Vc+VxV×PDi−1−PDi100.Therefore, from equations (11) and (15),(16)JDi=PDi−1−PDi100×1−VV1+aρc/ρl.In addition, the density parameterρs of the asphalt mortar could also be derived from equations (8) and (10):(17)ρs=msVs=ρlVl+ρxVxVl+Vx.Substituting equation (9) into equation (17), then(18)ρs=a+1+aVx/Vcaρc/ρl+1+aρc/ρlVx/Vcρc.Taking the asphalt mixture of AC-13 as an example, the gradation is shown in Table2. The density of the asphalt is 1.03 g/cm3, the aggregate density is 2.7 g/cm3, 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.Table 2
Gradation of the asphalt mixture of AC-13.
Sieve size
16
13.2
9.5
4.75
2.36
1.18
0.6
0.3
0.15
0.075
Passing percentage (%)
Upper limit
100
100
85
56
38
28
22
17
12
8
Lower limit
100
95
72
42
28
20
12
8
7
4
Median
100
97.5
78.5
49
33
24
17
12.5
9.5
6Table 3
Aggregate particles of the asphalt mixture of AC-13.
Sieve size
16
13.2
9.5
4.75
2.36
1.18
0.6
0.3
0.15
0.075
Passing percentages (%)
100
97.5
78.5
49
33
24
17
12.5
9.5
6
Individual percent retain on each sieve (%)
0
2.5
19
29.5
16
—
—
—
—
—
Volume fraction
0.000
0.021
0.161
0.250
0.136
—
—
—
—
—
Number
0
10
165
1036
4540
—
—
—
—
—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 Figure4. 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.Figure 4
Arrangement of aggregate particle.


## 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:(7)ni=V×JDi4/3πri3,where V=πd2h/4 (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, ri is the average radius of the particles, ri=Di−1+Di/4, and Di−1 and Di are the dimensions of the i−1th and ith 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ρc, the density of asphalt is ρl, the asphalt-aggregate ratio is VV, the design air voids is V=πd2h/4, and the volume of the cylindrical specimen is V=πd2h/4, where d is the diameter and h is the height. The volume of the asphalt is set to h, the volume of the coarse aggregate is Vc, and the volume of the fine aggregate is Vx. 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.Table 1
Percentages of aggregate gradation.
Sieve size (mm)
D
1
D
2
D
3
…
D
n−2
D
n−1
D
n = 2.36 mm
<2.36
Passing percentages (%)
100
P
D2
P
D3
…
P
Dn−2
P
Dn−1
P
2.36
…
Individual percent retain on each sieve (%)
0
100 −PD2
P
D2 − PD3
…
P
Dn−3 − PDn-2
P
Dn−2 − PDn−1
P
Dn−1 − P2.36
P
2.36According to the equation of the asphalt-aggregate ratio calculation, namely,(8)a=mlmc+mx=ρlVlρcVc+ρcVx.The content of asphalt in the numerical specimen is as follows:(9)Vl=aρcρlVc+Vx.According to the volume composition of the whole specimen, there are(10)V1−VV=Vl+Vc+Vx.Substituting equation (9) into equation (10), we obtain(11)Vc+VxV=1−VV1+aρc/ρl.According to the gradation passing percentages table,(12)mxmc=P2.36×mx+mc100−P2.36×mx+mc=P2.36100−P2.36.And, because the density is equal,(13)VxVc=P2.36100−P2.36.Substituting equation (13) into equation (11), we obtain(14)VcV=1−VV100−P2.361+aρc/ρl×100.Combined with the passing percentages table, the formula for the volume fraction of aggregates of each sieve is obtained as follows:(15)JDi=Vc+VxV×PDi−1−PDi100.Therefore, from equations (11) and (15),(16)JDi=PDi−1−PDi100×1−VV1+aρc/ρl.In addition, the density parameterρs of the asphalt mortar could also be derived from equations (8) and (10):(17)ρs=msVs=ρlVl+ρxVxVl+Vx.Substituting equation (9) into equation (17), then(18)ρs=a+1+aVx/Vcaρc/ρl+1+aρc/ρlVx/Vcρc.Taking the asphalt mixture of AC-13 as an example, the gradation is shown in Table2. The density of the asphalt is 1.03 g/cm3, the aggregate density is 2.7 g/cm3, 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.Table 2
Gradation of the asphalt mixture of AC-13.
Sieve size
16
13.2
9.5
4.75
2.36
1.18
0.6
0.3
0.15
0.075
Passing percentage (%)
Upper limit
100
100
85
56
38
28
22
17
12
8
Lower limit
100
95
72
42
28
20
12
8
7
4
Median
100
97.5
78.5
49
33
24
17
12.5
9.5
6Table 3
Aggregate particles of the asphalt mixture of AC-13.
Sieve size
16
13.2
9.5
4.75
2.36
1.18
0.6
0.3
0.15
0.075
Passing percentages (%)
100
97.5
78.5
49
33
24
17
12.5
9.5
6
Individual percent retain on each sieve (%)
0
2.5
19
29.5
16
—
—
—
—
—
Volume fraction
0.000
0.021
0.161
0.250
0.136
—
—
—
—
—
Number
0
10
165
1036
4540
—
—
—
—
—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 Figure4. 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.Figure 4
Arrangement of aggregate particle.


## 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 Figure5.Figure 5
Regular arrangement particle.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 Figure6.Figure 6
Specimen of irregular aggregate (no air voids).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 then_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 by(19)n_del=intVV∗n_total,where VV is the air voids and n_total is the total number of small particles in the model.The air void distribution is shown in Figure7. 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.Figure 7
Distribution of air voids. (a) Entire model. (b) Vertical section view. (c) Transverse section view.Figure 8
Numerical model of irregular aggregate.
(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 Figure9, the yellow part is the aggregate phase, the blue part is the asphalt cement phase, and the blank part is the air void phase.Figure 9
Structure of numerical model. (a) Entire model of AC-13. (b) Vertical section view of AC-13. (c) Transverse section view of AC-13. (d) Entire model of SMA-13. (e) Vertical section view of SMA-13. (f) Transverse section view of SMA-13.
(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.Figure 10
Burger’s model of particle contact. (a) Normal direction and (b) tangent direction.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 thefish 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 Table4.Table 4
Contact models at different locations inside the mixture.
Contact position
Model selection
Aggregate interior
Linear model + point bond model
Between adjacent aggregates
Linear model + point bond model
Asphalt mortar internal
Burger’s model + linear model + point bond model
Between asphalt mortar and aggregate
Burger’s model + linear model + point bond modelIn order to analyze the microscopic parameters of Burger’s model, it can be seen in Figure9 that this model has a total of eight parameters: Kkn, Kmn, Ckn, and Cmn and tangential Kks, Kms, Cks, and Cms.The mesocontact forcefn in the model could be expressed as(20)fn=ε˙mcCmnL=LεkKkn+Lε˙kCkn=εmkKmnL,where εk, εmk, and εmc are the strain of corresponding units.The stressσ in the macro Burger’s model could be expressed as(21)σ=ε˙mcη1=εkE2+ε˙kη2=εmkE1,where E1, η1, E2, and η2 are the macrofitting parameters of Burger’s model.And, because of the equivalence between contact force and stressfn=σA, there were(22)εmkKmnL=AεmkE1,For the same reason, there were(23)ε˙mcCmnL=Aε˙mcη1,Thus, we can obtain(24)Kmn=E1L,Cmn=η1L,Kkn=E2L,Ckn=η2L.The relationship between the tangential parameters of Burger’s model and the macroscopic parameters could be obtained as follows:(25)Kms=E1L21+v,Cms=η1L21+v,Kks=E2L21+v,Cks=η2L21+v.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).Table 5
Burger’s model parameter.
Mixture type
Burger’s model parameters
E
1 (MPa)
E
2 (MPa)
η
1 (MPa·s)
η
2 (MPa·s)
AC-13
3.185
1.756
2784.233
216.724
SMA-13
0.478
0.117
83.514
5.487Assuming 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 Tables6∼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.Table 6
Microscopic parameters of mixtures of AC-13.
Microscopic parameter
Stiffness
Strength
kn
ks
n_bond
s_bond
Aggregate interior
2.5 × 108
1.0 × 108
37.5
93.8
Between adjacent aggregates
2.5 × 108
1.0 × 108
4.1
9.4
Asphalt mortar internal
1.0 × 108
0.4 × 108
4.3
9.4
Between asphalt mortar and aggregate
1.75 × 108
0.7 × 108
4.3
9.4Table 7
Microscopic parameters of AC-13 Burger’s model.
Microscopic parameter
Maxwell
Kelvin
Stiffness
Viscosity
Stiffness
Viscosity
Normal direction
7.9 × 103
6.95 × 106
4.4 × 103
5.43 × 105
Tangent direction
3.16 × 103
2.78 × 106
1.76 × 103
2.17 × 105Table 8
Microscopic parameters of mixtures of SMA-13.
Microscopic parameter
Stiffness
Strength
kn
ks
n_bond
s_bond
Aggregate interior
2.5 × 108
1.0 × 108
37.5
93.8
Between adjacent aggregates
2.5 × 108
1.0 × 108
4.2
9.5
Asphalt mortar internal
5 × 107
2 × 107
4.3
9.6
Between asphalt mortar and aggregate
1.5 × 108
0.6 × 108
4.3
9.6Table 9
Microscopic parameters of SMA-13 Burger’s model.
Microscopic parameter
Maxwell
Kelvin
Stiffness
Viscosity
Stiffness
Viscosity
Normal direction
2.5 × 103
4.18 × 105
2.93 × 103
2.75 × 104
Tangent direction
9.7 × 102
1.67 × 105
1.17 × 103
1.1 × 104


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


### 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 Table10. The failure state of the model is shown in Figure 12.Table 10
Results of torsion shear numerical simulation test.
Numbering
Mix type
Result (MPa)
Average value (MPa)
Coefficient of variation (%)
1
AC-13
0.350
0.360
1.94
2
0.365
3
0.368
4
0.355
5
SMA-13
0.510
0.506
1.23
6
0.502
7
0.495
8
0.509Figure 12
Results of two numerical model failure: (a) AC-13 and (b) SMA-13.
(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 Figure13. 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°.Figure 13
Displacement cloud image of two numerical model failure: (a) AC-13 and (b) SMA-13.
(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.Table 11
Comparison of numerical simulation results and laboratory test results.
Gradation type
Numerical simulation
Laboratory test
Relative error (%)
Shear strength (MPa)
Coefficient of variation
Shear strength (MPa)
Coefficient of variation
AC-13
0.360
1.94
0.345
5.89
4.35
SMA-13
0.506
1.23
0.487
4.32
3.90Figure 14
Failure of laboratory specimens: (a) AC-13 and (b) SMA-13.
(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 Figures13 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.


## 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 Figure9, the yellow part is the aggregate phase, the blue part is the asphalt cement phase, and the blank part is the air void phase.Figure 9
Structure of numerical model. (a) Entire model of AC-13. (b) Vertical section view of AC-13. (c) Transverse section view of AC-13. (d) Entire model of SMA-13. (e) Vertical section view of SMA-13. (f) Transverse section view of SMA-13.
(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.Figure 10
Burger’s model of particle contact. (a) Normal direction and (b) tangent direction.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 thefish 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 Table4.Table 4
Contact models at different locations inside the mixture.
Contact position
Model selection
Aggregate interior
Linear model + point bond model
Between adjacent aggregates
Linear model + point bond model
Asphalt mortar internal
Burger’s model + linear model + point bond model
Between asphalt mortar and aggregate
Burger’s model + linear model + point bond modelIn order to analyze the microscopic parameters of Burger’s model, it can be seen in Figure9 that this model has a total of eight parameters: Kkn, Kmn, Ckn, and Cmn and tangential Kks, Kms, Cks, and Cms.The mesocontact forcefn in the model could be expressed as(20)fn=ε˙mcCmnL=LεkKkn+Lε˙kCkn=εmkKmnL,where εk, εmk, and εmc are the strain of corresponding units.The stressσ in the macro Burger’s model could be expressed as(21)σ=ε˙mcη1=εkE2+ε˙kη2=εmkE1,where E1, η1, E2, and η2 are the macrofitting parameters of Burger’s model.And, because of the equivalence between contact force and stressfn=σA, there were(22)εmkKmnL=AεmkE1,For the same reason, there were(23)ε˙mcCmnL=Aε˙mcη1,Thus, we can obtain(24)Kmn=E1L,Cmn=η1L,Kkn=E2L,Ckn=η2L.The relationship between the tangential parameters of Burger’s model and the macroscopic parameters could be obtained as follows:(25)Kms=E1L21+v,Cms=η1L21+v,Kks=E2L21+v,Cks=η2L21+v.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).Table 5
Burger’s model parameter.
Mixture type
Burger’s model parameters
E
1 (MPa)
E
2 (MPa)
η
1 (MPa·s)
η
2 (MPa·s)
AC-13
3.185
1.756
2784.233
216.724
SMA-13
0.478
0.117
83.514
5.487Assuming 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 Tables6∼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.Table 6
Microscopic parameters of mixtures of AC-13.
Microscopic parameter
Stiffness
Strength
kn
ks
n_bond
s_bond
Aggregate interior
2.5 × 108
1.0 × 108
37.5
93.8
Between adjacent aggregates
2.5 × 108
1.0 × 108
4.1
9.4
Asphalt mortar internal
1.0 × 108
0.4 × 108
4.3
9.4
Between asphalt mortar and aggregate
1.75 × 108
0.7 × 108
4.3
9.4Table 7
Microscopic parameters of AC-13 Burger’s model.
Microscopic parameter
Maxwell
Kelvin
Stiffness
Viscosity
Stiffness
Viscosity
Normal direction
7.9 × 103
6.95 × 106
4.4 × 103
5.43 × 105
Tangent direction
3.16 × 103
2.78 × 106
1.76 × 103
2.17 × 105Table 8
Microscopic parameters of mixtures of SMA-13.
Microscopic parameter
Stiffness
Strength
kn
ks
n_bond
s_bond
Aggregate interior
2.5 × 108
1.0 × 108
37.5
93.8
Between adjacent aggregates
2.5 × 108
1.0 × 108
4.2
9.5
Asphalt mortar internal
5 × 107
2 × 107
4.3
9.6
Between asphalt mortar and aggregate
1.5 × 108
0.6 × 108
4.3
9.6Table 9
Microscopic parameters of SMA-13 Burger’s model.
Microscopic parameter
Maxwell
Kelvin
Stiffness
Viscosity
Stiffness
Viscosity
Normal direction
2.5 × 103
4.18 × 105
2.93 × 103
2.75 × 104
Tangent direction
9.7 × 102
1.67 × 105
1.17 × 103
1.1 × 104


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


## 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 Table10. The failure state of the model is shown in Figure 12.Table 10
Results of torsion shear numerical simulation test.
Numbering
Mix type
Result (MPa)
Average value (MPa)
Coefficient of variation (%)
1
AC-13
0.350
0.360
1.94
2
0.365
3
0.368
4
0.355
5
SMA-13
0.510
0.506
1.23
6
0.502
7
0.495
8
0.509Figure 12
Results of two numerical model failure: (a) AC-13 and (b) SMA-13.
(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 Figure13. 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°.Figure 13
Displacement cloud image of two numerical model failure: (a) AC-13 and (b) SMA-13.
(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.Table 11
Comparison of numerical simulation results and laboratory test results.
Gradation type
Numerical simulation
Laboratory test
Relative error (%)
Shear strength (MPa)
Coefficient of variation
Shear strength (MPa)
Coefficient of variation
AC-13
0.360
1.94
0.345
5.89
4.35
SMA-13
0.506
1.23
0.487
4.32
3.90Figure 14
Failure of laboratory specimens: (a) AC-13 and (b) SMA-13.
(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 Figures13 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 Figure15.Figure 15
Initial stage of the loading period. (a) Stress distribution (entire model). (b) Stress distribution (vertical section). (c) Stress distribution (transverse section). (d) Contact distribution (entire model). (e) Contact distribution (vertical section). (f) Contact distribution (transverse section).
(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 Figures15(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 Figures15(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 Figure16.Figure 16
Earlier stage of the loading period. (a) Stress distribution (entire model). (b) Stress distribution (vertical section). (c) Stress distribution (transverse section). (d) Contact distribution (entire model). (e) Contact distribution (vertical section). (f) Contact distribution (transverse section). (g) Displacement field (entire model). (h) Displacement field (vertical section). (i) Displacement field (transverse section).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)As shown in Figures16(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 Figures16(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 Figures16(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 Figure17.Figure 17
Later stage of the loading period (appearance of crack). (a) Stress distribution (entire model). (b) Stress distribution (vertical section). (c) Stress distribution (transverse section). (d) Contact distribution (entire model). (e) Contact distribution (vertical section). (f) Contact distribution (transverse section). (g) Displacement field (entire model). (h) Displacement field (vertical section). (i) Displacement field (transverse section).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)As shown in Figures17(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 Figures17(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 Figures17(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 Figure18.Figure 18
Ending stage of the loading period (model failure). (a) Stress distribution (entire model). (b) Stress distribution (vertical section). (c) Stress distribution (transverse section). (d) Contact distribution (entire model). (e) Contact distribution (vertical section). (f) Contact distribution (transverse section). (g) Displacement field (entire model). (h) Displacement field (vertical section). (i) Displacement field (transverse section). (j) Crack distribution (whole model). (k) Crack distribution (vertical section). (l) Crack distribution (transverse section).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)As shown in Figures18(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 Figures18(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 Figures18(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 Figures18(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 Figure19.Figure 19
Displacement cloud images of two numerical models. (a) Initial stage. (b) Earlier stage of the loading period. (c) Loading period 1. (d) Loading period 2. (e) Loading period 3. (f) Crack appeared. (g) Crack propagation. (h) End of loading.
(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.


## 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 Figure15.Figure 15
Initial stage of the loading period. (a) Stress distribution (entire model). (b) Stress distribution (vertical section). (c) Stress distribution (transverse section). (d) Contact distribution (entire model). (e) Contact distribution (vertical section). (f) Contact distribution (transverse section).
(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 Figures15(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 Figures15(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 Figure16.Figure 16
Earlier stage of the loading period. (a) Stress distribution (entire model). (b) Stress distribution (vertical section). (c) Stress distribution (transverse section). (d) Contact distribution (entire model). (e) Contact distribution (vertical section). (f) Contact distribution (transverse section). (g) Displacement field (entire model). (h) Displacement field (vertical section). (i) Displacement field (transverse section).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)As shown in Figures16(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 Figures16(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 Figures16(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 Figure17.Figure 17
Later stage of the loading period (appearance of crack). (a) Stress distribution (entire model). (b) Stress distribution (vertical section). (c) Stress distribution (transverse section). (d) Contact distribution (entire model). (e) Contact distribution (vertical section). (f) Contact distribution (transverse section). (g) Displacement field (entire model). (h) Displacement field (vertical section). (i) Displacement field (transverse section).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)As shown in Figures17(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 Figures17(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 Figures17(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 Figure18.Figure 18
Ending stage of the loading period (model failure). (a) Stress distribution (entire model). (b) Stress distribution (vertical section). (c) Stress distribution (transverse section). (d) Contact distribution (entire model). (e) Contact distribution (vertical section). (f) Contact distribution (transverse section). (g) Displacement field (entire model). (h) Displacement field (vertical section). (i) Displacement field (transverse section). (j) Crack distribution (whole model). (k) Crack distribution (vertical section). (l) Crack distribution (transverse section).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)As shown in Figures18(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 Figures18(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 Figures18(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 Figures18(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 Figure19.Figure 19
Displacement cloud images of two numerical models. (a) Initial stage. (b) Earlier stage of the loading period. (c) Loading period 1. (d) Loading period 2. (e) Loading period 3. (f) Crack appeared. (g) Crack propagation. (h) End of loading.
(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.



---

*Source: 8709849-2019-09-04.xml*